When I some years ago began working with nuclear cross section calculations, via the reaction codes EMPIRE and TALYS, my initial thought was; “No, I can’t start plotting all these cross section plots individually. There has to be a better way”. I think this was when I really started to code in Bash.
data-plot script: terminal input
data-plot script: gnuplot output
data-plot script: gnuplot output
There was an easy solution. It took some scripting, but the solution was there. It didn’t take long either before I have my first working version of the cross section plotting program. It was a simple Bash script tied in with the plotting capabilities of gnuplot. Working in Linux gives some perks, especially through Bash with its easy syntax, wide availability to different sets of tools (for example find, grep, and sed), and the ability to run program scripts inherently. The requirements for my program grew, and I continued developing. Revision after revision the kinks and problems have ebbed away and a solid script has emerged over the years.
Why this solution?
As usual in coding the existing program really didn’t cut it, and I needed functions that wasn’t supplied by the existing solutions. One such solution is EMPIRE’s native ZVD plotting program. The advantages with this script are the following:
- The ability to plot both EMPIRE, TALYS, EXFOR, TENDL, and your own experimental and theoretical cross section data side by side.
- Plotting from a database based on a calculation search strings, which enables multiple plotting in one go. Note that plotting a massive numer of data sets might result in the plot legend (key in gnuplot) overflowing your gnuplot window.
- The option to plot product cross sections for metastable and/or ground states.
- The ability to plot data on natural targets from standard EMPIRE and TALYS calculations. I.e. the script will automatically correct the isotope cross section based on the natural occurrence of said element.
- The option to plot yields as well as all cross sections from a certain reaction.
- An inbuilt ability to plot a certain energy range.
- Excellent options for data output, and the possibility to continue analysing the data via gnuplot or other scripts.
Working notes on the code
There are still some issues with the script, mainly concerning the speed of the program. This can partly be traced down to a incomparability between the EXFOR data format (the EXFOR database is the EMPIRE copy of the latest X4toc4 master version) and gnuplot’s input formats. The problem stems from the different notations of scientific number formats, and my workaround isn’t extremly fast. Another speed-related problem of the program will occur when/if you have amassed large libraries of nuclear data. It boils down to the simple fact that it takes time to search through massive data libraries. The solution to this is optimising the current processes as well as developing parallell processes. Sadly this isn’t really the strenght of Bash. The development in these areas will however continue.
There are also some database related issues and functionalities in development. One is the implementation of an automatic TENDL database. Right now you will need to locate the necessary data files from the TENDL database and save them to your own file system in order to search and plot TENDL cross section data.
License / tips?
Feel free to use the program and contact me if you have any questions regarding it. If you find it helpful, please let me know and I will contine to update the code also here.
Requirements
What are needed in order to run the program?
- Bash. If you want to run it in some other shell it should be possible after some minor modifications to the code.
- gnuplot (preferably 5.0+). I highly recommend using the script with gnuplot v5.0+ in order to have an responsive based legend, whereby you can display and hide plots by clicking the corresponding legend key. The script also works on older versions of gnuplot, but you might need to change the output terminal from wxt to qt in code line 1,221.
- data-plot-bash. The placement of the script on your hard drive doesn’t matter as long as you have set up your bash aliases correctly in your bashrc file.
- data-plot-ref. This should be saved in the same directory as the script.
The code
#!/bin/bash
#
# This script runs Gnuplot for several EMPIRE and/or TALYS data files
# Gnuplot plots data for total cross section comparison with EXFOR, TENDL
# and experimental data files.
#
# Script name: data-plot
# Author: Alex Björkholm 2013-2016
# Date: 2016.08.16
################################ WORKING NOTES ################################
# NOTE: Always check with the EXFOR database for missing data. https://www-nds.iaea.org/exfor/exfor.htm
# NOTE: TENDL database not implemented
################################# DECLARATIONS ################################
echo -e "\n
---------------------------\n
CROSS SECTION PLOT \n
---------------------------\n"
IFS=$'\n'
declare -a data_plot
declare -a input_data_bases
declare -a data_bases
declare -a input_search
declare -a dir_empire
declare -a dir_talys
declare -a dir_tendl
declare -a dir_exper
declare -a dir_exper
declare -a dir_gnuplot
declare -a file_paths_empire
declare -a file_paths_talys
declare -a file_paths_tendl
declare -a file_paths_exfor
declare -a file_paths_exper
declare -a input_reaction
declare -a input_products
declare -a natural_isotopes
declare -a targets
declare -a targetsAA
declare -a targetsZZ
declare -a targetsNN
declare -a targetsXX
declare -a targetsPP
declare -a targetsMM
declare -a temp_targets
declare -a temp_targets_MM
declare -a temp_targets_PP
declare -a beams
declare -a beamsAA
declare -a beamsZZ
declare -a beamsNN
declare -a beamsXX
declare -a temp_beams
declare -a products
declare -a productsAA
declare -a productsZZ
declare -a productsNN
declare -a productsXX
declare -a productsMM
declare -a input_products
declare -a temp_products_MM
declare -a temp_products
declare -a reactions
declare -a temp_reactions
declare -a emissions
declare -a delete_array
cd $HOME
################################## VARIABLES ##################################
# Paths $RESEARCHDIR, $RESEARCHPROGRAMDIR, $EMPIREDIR, $TALYSDIR must be defined in bash
# Directory for this program and its reference data
programdir="$RESEARCHPROGRAMDIR/Cross section - Data plot scripts - Bash Gnuplot"
reference="$programdir/data-plot-ref.dat"
# Directories where EMPIRE or TALYS calculations are stored. Mandatory format for subdirectories is directories starting with "target,beam" with subsequent tags separated with "-". For example "75Br,h", "208Pb,p-run", "209Bi,a-omp-rot", and "202Hg,10B-omp".
dir_empire+=($(find $RESEARCHDIR -maxdepth 2 -type d -iwholename $RESEARCHDIR\*Calculations\ Empire\*))
dir_talys+=($(find $RESEARCHDIR -maxdepth 2 -type d -iwholename $RESEARCHDIR\*Calculations\ Talys\*))
# Directories where EXFOR and TENDL data are stored.
# Exfor data (download from http://www-nds.iaea.org/x4toc4-master/) should be sorted with the parseC4.py script in Empire to achieve the correct structure.
# Tendl data (download from ftp://ftp.nrg.eu/pub/www/talys/tendl2014/tendl2014.html) should be stored in the structure provided from the Tendle database.
dir_tendl+=($(find $RESEARCHDIR -maxdepth 3 -type d -iwholename $RESEARCHDIR\*Experimental\ and\ theoretical\ data\*tendl\*))
dir_exfor="$EMPIREDIR/EXFOR"
#dir_tendl="$TALYSDIR/TENDL"
# Directories where additional data in tab separated format is stored. Mandatory format for subdirectories is directories starting with "target,beam" with subsequent tags separated with "-". For example "75Br,h", "209Bi,a-omp-rot", and "202Hg,10B-omp".
dir_exper+=($(find $RESEARCHDIR -maxdepth 3 -type d -iwholename $RESEARCHDIR\*Experimental\ and\ theoretical\ data\*exper\* -not -iwholename \*chisq\*))
dir_gnuplot+=($(find $RESEARCHDIR -maxdepth 3 -type d -iwholename $RESEARCHDIR\*Experimental\ and\ theoretical\ data\*gnuplot\*))
# Directories where output from plots are saved.
mkdir -p $HOME/gnuplot-plots
gnuplot_plot="$HOME/gnuplot-plots"
# Directory for temp files used by GNUPLOT.
mkdir -p $HOME/gnuplot-temp
gnuplot_temp="$HOME/gnuplot-temp"
################################## FUNCTIONS ##################################
manual () {
echo -e "This script runs Gnuplot for several EMPIRE and/or TALYS data files.
Program data-plot input string:
1. Cross section data to plot
product options: xs, direct, preequilibrium, compound, all-xs
production options: activity, \#isotopes, yield, production
emission options: xs, yield, all-xs
2. System(s) to plot
options: exfor, exper, gnuplot, tendl, talys, empire
multiple options can be chosen by writing option1_option2_option3
3. Reaction input
options: '75As(h,x)76Br-75Br', 'natCr(p,x)52Mn', '52Cr(p,x)52mMn', '52Cr(p,x)52gMn' etc.
multiple options can be chosen by writing option1_option2_option3
product note: set outgoing particles to 'x' (calculated automatically)
4. Search string for EMPIRE and TALYS calculations
multiple options can be chosen by writing option1_option2_option3
wildcard notation * can be used to get the optimal selection
5. Graphical terminal type
options: edit/no-edit/eps/epslatex/svg/table
6. Optional x-axis range in the form '5-33'
if no range is defined the x-axis will autoscale to data
Example: data-plot xs talys_empire 'natCr(p,x)52Mn' 'optimal' no-edit 5-40
data-plot xs exfor_exper '75As(h,x)76Br-75Br' 'calculation' edit 5-50
data-plot xs exfor_exper '75As(h,x)all' 'omp*HFB' edit 5-50\n"
exit 2 ; }
natural_conversion () {
local input=$1
unset AAA; unset NAA; unset NAP
if [[ $input != *[0-9]* ]] && [[ $input == nat[A-Z]* ]] ; then
NZZ=$(echo $input | sed 's/nat//g')
zZZ=$(awk -v src=$NZZ '$2 == src { print $1 }' $reference)
ZZZ=($(echo ${zZZ} | sed 's/^0*//'))
while IFS=' ' read -r -a line; do
# Extracting isotope atomic numbers for natural target
AAA+=($(echo ${line[1]} | sed 's/^0*//'))
# Extracting isotope percentage for natural target
PP+=($( bc -l <<< "scale=6;${line[2]} / 100" | sed 's/^\./0./'))
done < $TALYSDIR/structure/abundance/z${zZZ}
fi
if [[ $AAA != *[0-9]* ]] || [[ $ZZZ != *[0-9]* ]]; then
echo -e "ERROR: natural_conversion function for input '$input'\n\n"; exit 2;
fi ; }
nuclide_conversion () {
local input=$1
# If format of the beam is n, p, h, a etc., resolving ZZZ and AAA
if [[ $input != *[0-9]* ]]; then
AAA=$(awk -v src=$input '$1 == src { print $5 }' $reference)
NZZ=${input}
zZZ=$(awk -v src=$input '$1 == src { print $3 }' $reference)
if [[ $input != n ]]; then ZZZ=($(echo ${zZZ} | sed 's/^0*//')); fi
if [[ $input == n ]]; then ZZZ=($(echo ${zZZ})); fi
fi
# If format of target is ZZZAAA (083209), resolving ZZZ and AAA
if [[ $input != *[a-z]* ]]; then
AAA=$(echo ${input:3:7} | sed 's/^0*//')
ZZZ=$(echo ${input:0:3} | sed 's/^0*//')
zZZ=$(echo ${input:0:3})
NZZ=$(awk -v src=$zZZ '$1 == src { print $2 }' $reference)
fi
# If format of the nuclide is AAAXX (209Bi), resolving ZZZ and AAA
if [[ $input == *[0-9][A-Z]* ]]; then
AAA=$(echo $input | sed 's/[^0-9]*//g')
NZZ=$(echo $input | sed 's/^[0-9]*//g')
zZZ=$(awk -v src=$(echo $NZZ | sed 's/^m//g') '$2 == src { print $1 }' $reference)
ZZZ=($(echo ${zZZ} | sed 's/^0*//'))
fi
if [[ $AAA != *[0-9]* ]] || [[ $ZZZ != *[0-9]* ]]; then
echo -e "ERROR: nuclide_conversion function for input '$input'\n\n"; exit 2
fi
NNN=$(($AAA - $ZZZ)); }
p-n-emission () {
local particle_emission=$1
emission_number=0 ; NOUT=0 ; POUT=0
if [ $particle_emission == gamma ]; then particle_emission=g; else
particle_emission=$(echo $particle_emission | grep -o '[0-9]*[a-z]')
for each_emission in ${particle_emission[@]}; do
emission_particle=$(echo $each_emission | grep -o '[a-z]')
emission_number=$(echo $each_emission | grep -o '[0-9]*')
if [ -z $emission_number ]; then emission_number=1; fi
n=$(awk -v src=$emission_particle '$1 == src { print $4 }' $reference)
p=$(awk -v src=$emission_particle '$1 == src { print $3 }' $reference)
NOUT=$(($NOUT + $emission_number * $n))
POUT=$(($POUT + $emission_number * $p))
done
fi; }
particle-emission () {
NOUT=$(($TNN + $BNN - $PNN))
POUT=$(($TZZ + $BZZ - $PZZ))
particle_emission=${NOUT}n${POUT}p
if [ $particle_emission == 1n1p ]; then particle_emission=d; fi
if [ $particle_emission == 2n1p ]; then particle_emission=t; fi
if [ $particle_emission == 1n0p ]; then particle_emission=n; fi
if [ $particle_emission == 0n0p ]; then particle_emission=gamma; fi
if [[ $particle_emission == *n0p ]]; then particle_emission=$(echo $particle_emission | sed 's/0p//g'); fi
if [[ $particle_emission == *-* ]]; then particle_emission=none; fi; }
MT_number () {
local MT=()
declare -a MT
# See endf-6 manual p. 306-316. Here simplified since we're only interested in the produced nuclei
if [ $POUT == 0 -a $NOUT == 0 ];then MT+=(2); fi
if [ $POUT == 0 -a $NOUT == 1 ];then MT+=(4); fi
if [ $POUT == 1 -a $NOUT == 3 ];then MT+=(11); fi
if [ $POUT == 0 -a $NOUT == 2 ];then MT+=(16); fi
if [ $POUT == 0 -a $NOUT == 3 ];then MT+=(17); fi
if [ $POUT == 2 -a $NOUT == 3 ];then MT+=(22); fi
if [ $POUT == 6 -a $NOUT == 7 ];then MT+=(23); fi
if [ $POUT == 2 -a $NOUT == 4 ];then MT+=(24); fi
if [ $POUT == 2 -a $NOUT == 5 ];then MT+=(25); fi
if [ $POUT == 1 -a $NOUT == 1 ];then MT+=(28); fi
if [ $POUT == 4 -a $NOUT == 5 ];then MT+=(29); fi
if [ $POUT == 4 -a $NOUT == 6 ];then MT+=(30); fi
if [ $POUT == 1 -a $NOUT == 2 ];then MT+=(32); fi
if [ $POUT == 1 -a $NOUT == 3 ];then MT+=(33); fi
if [ $POUT == 2 -a $NOUT == 2 ];then MT+=(34); fi
if [ $POUT == 3 -a $NOUT == 4 ];then MT+=(35); fi
if [ $POUT == 5 -a $NOUT == 7 ];then MT+=(36); fi
if [ $POUT == 0 -a $NOUT == 4 ];then MT+=(37); fi
if [ $POUT == 1 -a $NOUT == 2 ];then MT+=(41); fi
if [ $POUT == 1 -a $NOUT == 3 ];then MT+=(42); fi
if [ $POUT == 2 -a $NOUT == 1 ];then MT+=(44); fi
if [ $POUT == 3 -a $NOUT == 3 ];then MT+=(45); fi
if [ $POUT == 0 -a $NOUT == 0 ];then MT+=(103); fi
if [ $POUT == 1 -a $NOUT == 0 ];then MT+=(103); fi
if [ $POUT == 1 -a $NOUT == 1 ];then MT+=(104); fi
if [ $POUT == 1 -a $NOUT == 2 ];then MT+=(105); fi
if [ $POUT == 2 -a $NOUT == 1 ];then MT+=(106); fi
if [ $POUT == 2 -a $NOUT == 2 ];then MT+=(107); fi
if [ $POUT == 4 -a $NOUT == 4 ];then MT+=(108); fi
if [ $POUT == 6 -a $NOUT == 6 ];then MT+=(109); fi
if [ $POUT == 2 -a $NOUT == 0 ];then MT+=(111); fi
if [ $POUT == 3 -a $NOUT == 2 ];then MT+=(112); fi
if [ $POUT == 5 -a $NOUT == 6 ];then MT+=(113); fi
if [ $POUT == 5 -a $NOUT == 5 ];then MT+=(114); fi
if [ $POUT == 2 -a $NOUT == 1 ];then MT+=(115); fi
if [ $POUT == 2 -a $NOUT == 2 ];then MT+=(116); fi
if [ $POUT == 3 -a $NOUT == 3 ];then MT+=(117); fi
if [[ -z ${MT} ]] && [[ $POUT -ge 0 ]] && [[ $NOUT -ge 0 ]];then MT+=(5); fi
MTnum=(${MT[@]}); }
sort_function_underscore () {
if [ -z $1 ]; then echo -e "ERROR: missing input in 'sort_function_underscore'"; else
for each in $(echo $1 | sed 's/_/ /g' | tr ' ' '\n' | sort -u); do echo $each; done; fi ; }
sort_function_bar () {
if [ -z $1 ]; then echo -e "ERROR: missing input in 'sort_function_bar'"; else
for each in $(echo $1 | sed 's/-/ /g' | tr ' ' '\n' | sort -u); do echo $each; done; fi ; }
############################# DEFINITION OF INPUT #############################
# CONTROLLING AND RESOLVING PARAMETERS
# Setting output and checking input for product/emission
data_plot=($(sort_function_underscore $1))
input_data_bases=($(sort_function_underscore $2))
input_reaction=($(sort_function_underscore $3))
input_search=($(sort_function_underscore $4))
# Setting output and checking input for plot axis
for each in ${data_plot[@]}; do
if [[ ${input_reaction} == *x\)* ]]; then
case $each in
man ) manual; ;;
xs | direct | preequilibrium | compound | all-xs ) echo -e "Plotting product cross section $each\n"; ;;
activity | \#isotopes | yield | production) echo -e "Plotting medical nuclide production $each\n"; ;;
* ) echo -e "Input cross section '$each' incorrect for product plot\n"; exit 2; ;;
esac; fi
if [[ ${input_reaction} != *x\)* ]]; then
case $each in
man ) manual; ;;
xs | yield | all-xs ) echo -e "Plotting emission cross section $each\n"; ;;
* ) echo -e "Input cross section '$each' incorrect for emission plot\n"; exit 2; ;;
esac; fi
done
# Setting output and checking input for systems
for each_system in ${input_data_bases[@]}; do
case $each_system in
empire | exfor | exper | talys | tendl | gnuplot ) data_bases+=(${each_system}); echo -e "Plotting ${each_system^^} data"; ;;
* ) echo -e "Input '${each_system^^}' for data system to plot incorrect\n"; exit 2; ;;
esac
done
case $5 in
edit | no-edit | eps | epslatex | svg | table ) gnuplot_output=$5; ;;
man ) manual; ;;
* ) echo -e "Input '$5' incorrect. Choose edit/no-edit/eps/epslatex/svg or 'man' for manual.\n"; exit 2; ;;
esac
case $6 in
*-* ) gnuplot_x_range=$6; echo; ;;
*) echo -e "Input '$6' for x-axis range missing. Setting auto scale\n"; ;;
esac
######################## RESOLVING REACTION PARAMETERS ########################
for each_input_reaction in ${input_reaction[@]}; do
# Setting 'reactions' array for input target and beam
input_target=($(echo $each_input_reaction | awk -F"(" '{print $1}'))
input_beam=($(echo $each_input_reaction | awk -F"(" '{print $2}' | awk -F"," '{print $1}'))
# Setting natural target isotopes
if [[ ${input_target} == *nat* ]]; then
natural_conversion ${input_target}
for (( N = 0; N < ${#PP[*]}; N++ )); do
temp_targets_PP+=($(echo ${PP[N]}))
temp_targets_MM+=($(echo nat))
temp_targets+=($(echo ${AAA[N]}${NZZ}))
done
else
# Non natural target isotopes and
if [[ ${input_target} == *[0-9][mg][A-Z][a-z] ]]; then
if [[ ${input_target} == *[0-9][g][A-Z]* ]]; then temp_targets_MM+=($(echo g)); fi
if [[ ${input_target} == *[0-9][m][A-Z]* ]]; then temp_targets_MM+=($(echo m)); fi
input_target=$(echo ${input_target} | sed 's/^\(.\{2\}\)\(.\{1\}\)\(.*\)/\1\3/')
else
temp_targets_MM+=($(echo n))
fi
temp_targets+=(${input_target})
temp_targets_PP+=($(echo 1.000000))
fi
for (( N = 0; N < ${#temp_targets[*]}; N++ )); do
temp_beams+=(${input_beam})
done
# Setting 'reactions' array for input emission and products
input_emissions=($(sort_function_bar $(echo $each_input_reaction | awk -F"," '{print $2}' | awk -F")" '{print $1}')))
input_products=($(sort_function_bar $(echo $each_input_reaction | awk -F"," '{print $2}' | awk -F")" '{print $2}')))
# Product dependent sequence
if [[ ${input_emissions[@]} == *x* ]] && [[ ${input_products[@]} != x ]]; then
if [[ ${input_products[@]} == *all* ]]; then
# For all products input
for (( N = 0; N < ${#temp_targets[*]}; N++ )); do
for (( n = 0; n < 10; n++ )) ; do
for (( p = 0; p < 5; p++ )) ; do
nuclide_conversion ${temp_beams[N]}
BAA=${AAA}; BZZ=${ZZZ}
nuclide_conversion ${temp_targets[N]}
TAA=${AAA}; TZZ=${ZZZ}
AAA=$(($TAA + $BAA - $p - $n)); if [ $AAA -le 1 ]; then AAA=1; fi
ZZZ=$(($TZZ + $BZZ - $p)); if [ $ZZZ -le 1 ]; then ZZZ=1; fi
if [ $ZZZ -lt 100 ]; then zZZ="0$ZZZ"; else zZZ="$ZZZ"; fi
if [ $ZZZ -lt 10 ]; then zZZ="00$ZZZ"; fi
NZZ=$(awk -v src="$zZZ" '$1 ~ src { print $2 }' $reference)
temp_products+=(${AAA}${NZZ})
temp_products_MM+=($(echo NaN))
done
done
done
else
# For specific products input
for each_product in ${input_products[@]}; do
if [[ ${each_product} == *[0-9][mg][A-Z][a-z] ]]; then
if [[ ${each_product} == *[0-9][g][A-Z]* ]]; then temp_products_MM+=($(echo g)); fi
if [[ ${each_product} == *[0-9][m][A-Z]* ]]; then temp_products_MM+=($(echo m)); fi
each_product=$(echo ${each_product} | sed 's/^\(.\{2\}\)\(.\{1\}\)\(.*\)/\1\3/')
else
temp_products_MM+=($(echo NaN))
fi
temp_products+=(${each_product})
done
fi
# Resolving reaction arrays
for (( N = 0; N < ${#temp_targets[*]}; N++ )); do
for (( M = 0; M < ${#temp_products[*]}; M++ )); do
# Resolving targets
nuclide_conversion ${temp_targets[N]}
targetsAA+=(${AAA}); targetsZZ+=(${ZZZ})
targetsNN+=(${NNN}); targetsXX+=(${NZZ})
targetsPP+=(${temp_targets_PP[N]})
targetsMM+=(${temp_targets_MM[N]})
targets+=(${temp_targets[N]})
# Resolving beams
nuclide_conversion ${temp_beams[N]}
beamsAA+=(${AAA}); beamsZZ+=(${ZZZ})
beamsNN+=(${NNN})
beams+=(${temp_beams[N]})
# Resolving products
nuclide_conversion ${temp_products[M]}
productsAA+=(${AAA}); productsZZ+=(${ZZZ})
productsNN+=(${NNN}); productsXX+=(${NZZ})
productsMM+=(${temp_products_MM[M]})
products+=(${temp_products[M]})
# Resolving reactions
reactions+=($(echo ${temp_targets[N]},${temp_beams[N]}))
done
done
else
# Emission dependent sequence
# For all emission input
if [[ $input_emissions != *all* ]] && [[ ${input_products} == x ]]; then
echo -e "This option is not yet implemented \n"; exit 2
# For specific emission input
else
echo -e "This option is not yet implemented \n"; exit 2
fi
fi
done
echo -e "Reactions to be plotted"
for (( N = 0; N < ${#reactions[*]}; N++ )); do
printf " ${targets[N]}(${beams[N]},x)${products[N]} \t\t target ocurrence: %10s \t\t target state: ${targetsMM[N]} \t\t product state: ${productsMM[N]}\n" ${targetsPP[N]}
done
##################### EMPIRE DATA AND PLOT SEQUENCE ###########################
echo -e "\nSorting EMPIRE data"
# Excluding irrelevant EMPIRE top directories based in reaction product and target
for (( N = 0; N < ${#reactions[*]}; N++ )); do
for (( M = 0; M < ${#dir_empire[*]}; M++ )); do
if [[ ${dir_empire[M]} == *${productsXX[N]}* ]] || [[ ${dir_empire[M]} == *${targetsXX[N]}* ]]; then
delete_array+=(${dir_empire[M]}); fi
done
done
dir_empire=($(for each in ${delete_array[@]}; do echo $each; done | sort -u)); unset delete_array
#for (( N = 0; N < ${#dir_empire[*]}; N++ )); do echo ${dir_empire[N]}; done
# Excluding file duplicities of EMPIRE
if [[ ${input_emissions[@]} == *x* ]] && [[ ${data_bases[@]} == *empire* ]]; then
temp_reactions=($(for each in ${reactions[@]}; do echo $each; done | sort -u))
# Finding files EMPIRE sequence
for each_dir in ${dir_empire[@]}; do
for search in ${input_search[@]}; do
for (( N = 0; N < ${#temp_reactions[*]}; N++ )); do
# Finding files for EMPIRE (product specification inside EMPIRE files)
if [[ $data_plot == *xs* ]]; then
file_paths_empire+=($(find ${each_dir} -type f -iwholename ${each_dir}\*${temp_reactions[N]}\*${search}\*.xsc -not -iwholename \*outputdir\*))
fi
# Finding data files for medical isotope production
if [[ $data_plot == production ]] || [[ $data_plot == activity ]] || [[ $data_plot == yield ]] || [[ $data_plot == \#isotopes ]]; then
echo -e "EMPIRE: no support for '$data_plot' data"
fi
if [[ ${#file_paths_empire[*]} == $no_files ]] || [[ ${#file_paths_empire[*]} == 0 ]]; then
echo -e "EMPIRE: invalid reaction '${temp_reactions[N]}' and '${search}' in '$(basename ${each_dir})'"
else no_files=${#file_paths_empire[*]}; fi
done
done
done
fi
file_paths_empire=($(for file_path in ${file_paths_empire[@]}; do echo $file_path; done | sort -u))
#for file_path in ${file_paths_empire[@]}; do echo "'$file_path'"; done
# Checking if need to run sequense and correct cross section
if [ ${#file_paths_empire[*]} != 0 ]; then
if [ $data_plot == xs -o $data_plot == all-xs ]; then
for (( N = 0; N < ${#file_paths_empire[*]}; N++ )) ; do
# For product specific reactions. All reactions; resolved in reaction sequence
# Sequence based on extraction of all reaction (emission) channels from empire data file
if [[ ${input_emissions[@]} == *x* ]] && [[ ${input_products[@]} != *all* ]]; then
column_no=`cat ${file_paths_empire[$N]} | head -2 | tail -1 | wc -w`
for (( C = 8; C < $column_no; C++ )) ; do
for (( M = 0; M < ${#reactions[*]}; M++ )) ; do
dir_name=$(basename $(dirname ${file_paths_empire[$N]}))
# Finding the correct products to plot
reaction=$(echo ${dir_name} | awk -F"-" '{print $1}')
if [[ ${reactions[M]} == ${reaction} ]]; then
particle_emission=$(cat ${file_paths_empire[$N]} | head -2 | tail -1 | tr -s ' ' | cut -d ' ' -f$((C+1)) | sed 's/)//g; s/(z,//g')
p-n-emission $particle_emission
BAA=${beamsAA[M]}; BZZ=${beamsZZ[M]}
TAA=${targetsAA[M]}; TZZ=${targetsZZ[M]}
AAA=$(($TAA + $BAA - $POUT - $NOUT))
ZZZ=$(($TZZ + $BZZ - $POUT))
if [ $ZZZ -le 100 ]; then zZZ="0$ZZZ"; else zZZ="$ZZZ"; fi
if [ $ZZZ -lt 10 ]; then zZZ="00$ZZZ"; fi
NZZ=$(awk -v src="$zZZ" '$1 ~ src { print $2 }' $reference)
# Excluding irrelevant products
if [[ ${products[M]} == ${AAA}${NZZ} ]]; then
# Setting target with correct reaction abundance
if [[ ${input_target[@]} == *nat* ]]; then
columns="(\$1):(\$${C}*${targetsPP[M]})"
else
columns="(\$1):(\$${C})"
fi
[ -n "$cmd_emp" ] && cmd_emp=${cmd_emp}, || cmd_emp=''
cmd_emp="${cmd_emp} '${file_paths_empire[$N]}' using ${columns} w lp pi 2 lc palette frac 0.$((3 * $N / 5 + 2 * $C)) title 'EMPIRE-$dir_name-${AAA}${NZZ}'"
fi
fi
done
done
fi
if [[ ${input_emissions[@]} != *x* ]]; then
for (( N = 2; N <= 5; N++ )) ; do
[ -n "$cmd_emp" ] && cmd_emp=${cmd_emp}, || cmd_emp=''
dir_name=$(basename $(dirname $(dirname ${file_paths_empire[$N]})))
xs_name=$(cat ${file_paths_empire[$N]} | head -2 | tail -1 | tr -s ' ' | cut -d ' ' -f$((N+1)))
cmd_emp="${cmd_emp} '${file_paths_empire[$N]}' using 1:$N w lp pi 2 lc palette frac 0.$((3 * $N / 5)) title 'EMPIRE-${xs_name,,}-$dir_name'"
done
fi
done
fi
fi
echo EMPIRE: sorting done
##################### TALYS DATA AND PLOT SEQUENCE ############################
echo -e "\nSorting TALYS data"
# Excluding irrelevant TALYS top directories based in reaction product and target
for (( N = 0; N < ${#reactions[*]}; N++ )); do
for (( M = 0; M < ${#dir_talys[*]}; M++ )); do
if [[ ${dir_talys[M]} == *${productsXX[N]}* ]] || [[ ${dir_talys[M]} == *${targetsXX[N]}* ]]; then
delete_array+=(${dir_talys[M]}); fi
done
done
dir_talys=($(for each in ${delete_array[@]}; do echo $each; done | sort -u)); unset delete_array
#for (( N = 0; N < ${#dir_talys[*]}; N++ )); do echo ${dir_talys[N]}; done
# Finding files TALYS sequence
if [[ ${input_emissions[@]} == *x* ]] && [[ ${data_bases[@]} == *talys* ]]; then
for each_dir in ${dir_talys[@]}; do
for search in ${input_search[@]}; do
for (( N = 0; N < ${#reactions[*]}; N++ )); do
# Finding product specific, natural target, files
if [[ ${input_target[@]} == *nat* ]]; then
search_target=$(echo nat${targetsXX[N]}); search_beam=$(echo ${beams[N]})
else
search_target=$(echo ${targets[N]}); search_beam=$(echo ${beams[N]})
fi
# Finding data files for medical isotope production
if [[ $data_plot == production ]] || [[ $data_plot == activity ]] || [[ $data_plot == yield ]] || [[ $data_plot == \#isotopes ]]; then
file_paths_talys+=($(find ${each_dir} -type f -iwholename \*${search_target},${search_beam}\*${search}\*/Y\*${productsZZ[N]}\*${productsAA[N]}\*.\* -not -name \*prev\*))
fi
# Finding product specific files
if [[ $data_plot == xs ]]; then
file_paths_talys+=($(find ${each_dir} -type f -iwholename \*${search_target},${search_beam}\*${search}\*/rp\*${productsZZ[N]}\*${productsAA[N]}.\* -not -name \*prev\*))
fi
if [[ ${#file_paths_talys[*]} == $no_files ]] || [[ ${#file_paths_talys[*]} == 0 ]]; then
echo -e "TALYS: invalid $data_plot product '${products[N]}' for reaction '${reactions[N]}' and '${search}' in '$(basename ${each_dir})'"
else no_files=${#file_paths_talys[*]}; fi
done
done
done
fi
file_paths_talys=($(for file_path in ${file_paths_talys[@]}; do echo $file_path; done | sort -u))
#for file_path in ${file_paths_talys[@]}; do echo "'$file_path'"; done
# Checking if need to run sequense
if [ ${#file_paths_talys[*]} != 0 ]; then
# Defining plotting function for TALYS
cmd_talys () {
# Resolving reaction parameters from file to get correct plot titles
dir_name=$(basename $(dirname $(dirname ${file_paths_talys[$N]})))
base_name=$(basename ${file_paths_talys[$N]} | sed 's/.tot//g')
product=$(echo $base_name | sed 's/^[A-Z]*//g;s/^[a-z]*//g' | awk -F"." '{print $1}')
nuclide_conversion ${product}
# Setting target with correct title/plot for natural targets
target=$(echo $dir_name | awk -F"-" '{print $1}' | awk -F"," '{print $1}')
base_target=$(echo $base_name | sed 's/.L[0-9][0-9]//' | awk -F"." '{print $2}' | sed 's/^0*//')
if [[ ! -z $base_target ]]; then
target="${target/nat/nat(${base_target})}" # Replaces nat with $base_target
dir_name="${dir_name/nat/nat(${base_target})}" # Replaces nat with $base_target
fi
# Setting product with correct title/plot for ground/meta state
if [[ ${base_name} == *.L[0-9][0-9]* ]] && [[ ${productsMM[@]} != *[g-m]* ]]; then
continue
else
if [[ $base_name == *[0-9].L00* ]]; then state="g"; fi
if [[ $base_name == *[0-9].L[0-9][1-9] ]]; then state="m"; fi
fi
# Setting target with correct reaction abundance
for (( M = 0; M < ${#reactions[*]}; M++ )) ; do
if [[ "nat(${targetsAA[$M]})${targetsXX[$M]}" == ${target} ]]; then targetPP=${targetsPP[M]}; fi
if [[ "nat${targetsXX[$M]}" == ${target} ]]; then targetPP=($(echo 1.000000)); fi
if [[ ${targets[$M]} == ${target} ]]; then targetPP=($(echo 1.000000)); fi
done
if [ $data_plot != all-xs -a $data_plot != production ]; then
if [ $data_plot != xs ]; then data_plot=${data_plot^}; fi
column_no=$(awk -v src=${data_plot} '{ for(i=1;i<=NF;i++){if ($i ~ src) {print i } } }' ${file_paths_talys[$N]})
if [[ ${file_paths_talys[$N]} == *Y*.tot ]]; then column_no=$(( column_no / 2)); else column_no=$(( column_no - 1)); fi
fi
# Setting correct column to plot in files
if [[ ${file_paths_talys[$N]} == *Y*.tot ]]; then
column_name=$(cat ${file_paths_talys[$N]} | head -6 | tail -1 | tr -s ' ' | cut -d ' ' -f$(( column_no * 2 )))
column_names=$(cat ${file_paths_talys[$N]} | head -6 | tail -1 | tr -s ' ' | cut -d ' ' -f$(( column_no * 2 + 1)))
ylabel="$column_name $column_names"
else
column_name=$(cat ${file_paths_talys[$N]} | head -5 | tail -1 | tr -s ' ' | cut -d ' ' -f$(( column_no + 1)))
fi
# Setting final variables
columns="(\$1):(\$${column_no}*$targetPP)"
title=$(echo TALYS-${column_name,,}-$dir_name-${AAA}${state}${NZZ})
[ -n "$cmd_tal" ] && cmd_tal=${cmd_tal}, || cmd_tal=''
cmd_tal="${cmd_tal} '${file_paths_talys[$N]}' using ${columns} w lp pi 2 lc palette frac 0.$((3 + 3 * $N / 5)) title '${title}'"; }
for (( N = 0; N < ${#file_paths_talys[*]}; N++ )) ; do
# Plotting all cross sections/medical isotope production choises
if [ $data_plot == all-xs -o $data_plot == production ]; then
if [ $data_plot == production ]; then MAXC=5; fi
if [ $data_plot == all-xs ]; then MAXC=6; fi
for (( column_no = 2; column_no < $MAXC ; column_no++ )) ; do cmd_talys; done
# Plotting specific cross section/medical isotope production choise
else
cmd_talys
fi
done
fi
echo TALYS: sorting done
##################### TENDL DATA AND PLOT SEQUENCE ############################
echo -e "\nSorting TENDL data"
# Excluding irrelevant TALYS top directories based in reaction product and target
for (( N = 0; N < ${#reactions[*]}; N++ )); do
for (( M = 0; M < ${#dir_tendl[*]}; M++ )); do
if [[ ${dir_tendl[M]} == *${productsXX[N]}* ]] || [[ ${dir_tendl[M]} == *${targetsXX[N]}* ]]; then
delete_array+=(${dir_tendl[M]}); break; fi
done
done
dir_tendl=($(for each in ${delete_array[@]}; do echo $each; done | sort -u)); unset delete_array
#for (( N = 0; N < ${#dir_tendl[*]}; N++ )); do echo ${dir_tendl[N]}; done
# Finding files TENDL sequence
if [[ ${input_emissions[@]} == *x* ]] && [[ ${data_bases[@]} == *tendl* ]]; then
for each_dir in ${dir_tendl[@]}; do
for (( N = 0; N < ${#reactions[*]}; N++ )); do
# Finding product specific files
if [[ $data_plot == xs ]]; then
file_paths_tendl+=($(find ${each_dir} -type f -iwholename ${each_dir}/\*${reactions[N]}/\*${products[N]}\*.tot))
fi
# Finding data files for medical isotope production
if [[ $data_plot == production ]] || [[ $data_plot == activity ]] || [[ $data_plot == nuclides ]]; then
file_paths_tendl+=($(find ${each_dir} -type f -iwholename ${each_dir}/\*${reactions[N]}/Y\*${products[N]}\*.tot))
fi
if [[ ${#file_paths_tendl[*]} == $no_files ]] || [[ ${#file_paths_tendl[*]} == 0 ]]; then
echo -e "TENDL: invalid $data_plot product '${products[N]}' for reaction '${reactions[N]}'"
else no_files=${#file_paths_tendl[*]}; fi
done
done
fi
file_paths_tendl=($(for file_path in ${file_paths_tendl[@]}; do echo $file_path; done | sort -u))
#for file_path in ${file_paths_tendl[@]}; do echo "'$file_path'"; done
# Checking if need to run sequense and correct cross section xs or all-xs
if [[ $data_plot == *xs* ]] && [[ ${data_bases[@]} == *tendl* ]]; then
if [ ${#file_paths_tendl[*]} != 0 ]; then
for (( N = 0; N < ${#file_paths_tendl[*]}; N++ )) ; do
# Resolving reaction parameters from file to get correct plot titles
dir_name=$(basename $(dirname $(dirname ${file_paths_tendl[$N]})))
base_name=$(basename ${file_paths_tendl[$N]} | sed 's/.tot//g')
# Setting target with correct reaction abundance
target=$(echo $base_name | awk -F"," '{print $2}' | awk -F"(" '{print $1}' | sed 's/ //g')
for (( M = 0; M < ${#reactions[*]}; M++ )) ; do
if [[ ${targets[$M]} == ${target} ]]; then
targetPP=${targetsPP[M]}
break
fi
done
# Setting product with correct title/plot for ground/meta state
if [[ ${base_name} == *[0-9[0-9][g-m][A-Z][a-z]* ]] && [[ ${productsMM[@]} != *[g-m]* ]]; then
continue
fi
[ -n "$cmd_ten" ] && cmd_ten=${cmd_ten}, || cmd_ten=''
columns="(\$1):(\$2*$targetPP) w lp pi 1"
cmd_ten="${cmd_ten} '${file_paths_tendl[$N]}' using ${columns} lc palette frac 0.$((6 + 3 * $N / 5)) title '$base_name'"
done
fi
fi
echo TENDL: sorting done
##################### EXFOR DATA AND PLOT SEQUENCE ############################
echo -e "\nSorting EXFOR data"
# Finding files EXFOR sequence
if [[ ${input_emissions[@]} == *x* ]] && [[ $data_plot == xs ]] && [[ ${data_bases[@]} == *exfor* ]]; then
for each_dir in ${dir_exfor}; do
for (( N = 0; N < ${#reactions[*]}; N++ )); do
# Finding product specific files
if [[ ${beams[N]} != *[0-9]* ]]; then
file_paths_exfor+=($(find $dir_exfor -type f -iwholename $dir_exfor/${beams[N]}\*/\*${targetsZZ[N]}\*${targetsXX[N]}\*${targetsAA[N]}.c4))
else
file_paths_exfor+=($(find $dir_exfor -type f -iwholename $dir_exfor/other/\*${targetsZZ[N]}\*${targetsXX[N]}\*${targetsAA[N]}.c4))
fi
if [[ ${#file_paths_exfor[*]} == $no_files ]] || [[ ${#file_paths_exfor[*]} == 0 ]]; then
echo -e "EXFOR: invalid $data_plot product '${products[N]}' for reaction '${reactions[N]}'"
else no_files=${#file_paths_exfor[*]}; fi
done
done
fi
file_paths_exfor=($(for file_path in ${file_paths_exfor[@]}; do echo $file_path; done | sort -u))
#for file_path in ${file_paths_exfor[@]}; do echo "'$file_path'"; done
# EXFOR SEQUENCE
# Checking if need to run sequense and correct cross section
if [[ $data_plot == *xs* ]] && [[ ${data_bases[@]} == *exfor* ]]; then
if [[ ${#file_paths_exfor[*]} != 0 ]]; then
for (( N = 0; N < ${#file_paths_exfor[*]}; N++ )) ; do
# Copying file so that the official EXFOR file keeps intact
cp ${file_paths_exfor[$N]} $gnuplot_temp
file_paths_exfor[$N]=$gnuplot_temp/$(basename ${file_paths_exfor[$N]})
# Producing plots for all products
for (( M = 0; M < ${#reactions[*]}; M++ )) ; do
NOUT=$((${targetsNN[M]} + ${beamsNN[M]} - ${productsNN[M]}))
POUT=$((${targetsZZ[M]} + ${beamsZZ[M]} - ${productsZZ[M]}))
# Number of individual data sets in EXFOR file
# each_MT_sets is the first line for each data set
MT_set=(); MT_number
for each_MT in ${MTnum[@]}; do
MT_set+=($(grep -n " "$each_MT"[[:alpha:]]*[[:blank:]]*]" ${file_paths_exfor[$N]} | grep "C4" | awk -F":" '{print $1}'))
done
for each_MT_set in ${MT_set[@]};do
each_MT_linei=$(($each_MT_set + 4)) # First line of data to plot in each data set
# Last line of data to plot in each data set. Taken from listen number of data sets in EXFOR file
each_MT_linef=$(($(sed -n $(($each_MT_linei - 3))p ${file_paths_exfor[$N]}| awk '{print $2}') + $each_MT_linei))
# Lines omitted by gnuplot and corrected counting from zero
unused_lines=$(($(head -n ${each_MT_linei} ${file_paths_exfor[$N]} | grep "\#" | wc -l) + 1))
# Omitting irrelevant data sets
not_xs_data_set=$(sed -n "$each_MT_linei"p ${file_paths_exfor[$N]} | sed 's/^\(.\{59\}\)\(.\{25\}\)\(.*\)/\2/')
if [[ -z "${not_xs_data_set// }" ]]; then
for (( P = $each_MT_linei; P < ${each_MT_linef}; P++ )); do
dataline=$(sed -n ${P}p ${file_paths_exfor[$N]})
# Setting correct scientific format for data column 8
column08pm=$(echo ${dataline} | sed 's/^\(.\{56\}\)\(.\{1\}\)\(.*\)/\2/')
if [[ ${column08pm// } == + ]] || [[ ${column08pm// } == - ]]; then
sed "${P} s/^\(.\{56\}\)\(.\{1\}\)\(.*\)/\1E\2\3/" -i ${file_paths_exfor[$N]}
else
sed "${P} s/^\(.\{58\}\)\(.\{1\}\)\(.*\)/\1 \3/" -i ${file_paths_exfor[$N]}
fi
# Setting empty value in column 7 to zero
column07=$(echo ${dataline} | sed 's/^\(.\{50\}\)\(.\{8\}\)\(.*\)/\2/')
if [[ -z "${column07// }" ]]; then
sed "${P} s/^\(.\{50\}\)\(.\{1\}\)\(.*\)/\10\3/" -i ${file_paths_exfor[$N]}; fi
# setting correct scientific format for data column 7
column07pm=$(echo ${dataline} | sed 's/^\(.\{47\}\)\(.\{1\}\)\(.*\)/\2/')
if [[ ${column07pm// } == + ]] || [[ ${column07pm// } == - ]]; then
sed "${P} s/^\(.\{47\}\)\(.\{1\}\)\(.*\)/\1E\2\3/" -i ${file_paths_exfor[$N]}
else
sed "${P} s/^\(.\{49\}\)\(.\{1\}\)\(.*\)/\1\2\2\3/" -i ${file_paths_exfor[$N]}
fi
# Setting correct scientific format for data column 6
column06pm=$(echo ${dataline} | sed 's/^\(.\{38\}\)\(.\{1\}\)\(.*\)/\2/')
if [[ ${column06pm// } == + ]] || [[ ${column06pm// } == - ]]; then
sed "${P} s/^\(.\{38\}\)\(.\{1\}\)\(.*\)/\1E\2\3/" -i ${file_paths_exfor[$N]}
else
sed "${P} s/^\(.\{40\}\)\(.\{1\}\)\(.*\)/\1\2\2\3/" -i ${file_paths_exfor[$N]}
fi
# Setting empty value in column 5 to zero
column05=$(echo ${dataline} | sed 's/^\(.\{32\}\)\(.\{8\}\)\(.*\)/\2/')
if [[ -z "${column05// }" ]]; then
sed "${P} s/^\(.\{32\}\)\(.\{1\}\)\(.*\)/\10\3/" -i ${file_paths_exfor[$N]}; fi
# setting correct scientific format for data column 5
column05pm=$(echo ${dataline} | sed 's/^\(.\{29\}\)\(.\{1\}\)\(.*\)/\2/')
if [[ ${column05pm// } == + ]] || [[ ${column05pm// } == - ]]; then
sed "${P} s/^\(.\{29\}\)\(.\{1\}\)\(.*\)/\1E\2\3/" -i ${file_paths_exfor[$N]}
else
sed "${P} s/^\(.\{31\}\)\(.\{1\}\)\(.*\)/\1\2\2\3/" -i ${file_paths_exfor[$N]}
fi
done
# Energy plotted to Cross section; with energy error lower and upper limit ; with cross section error lower and upper limit
# Setting target with correct reaction abundance
columns="(\$5*1E-6):(\$7*1E+3*${targetsPP[M]}):(\$6*1E-6):(\$8*1E+3*${targetsPP[M]}) w xyerrorbars"
every="every ::$(($each_MT_linei - $unused_lines))::$(($each_MT_linef - $unused_lines))"
AUTHOR=$(sed -n "$each_MT_linei"p ${file_paths_exfor[$N]} | sed 's/^\(.\{98\}\)\(.*\)/\2/' | awk -F"(" '{print $1}')
AUTHOR=$(echo ${AUTHOR[@],,} | sed -r 's/\<./\U&/g' | sed 's/([0-9]*)//g;s/[[:space:]]//g;s/\./\. /g;s/,Et. Al./ et al./g' )
YEAR=$(sed -n "$each_MT_linei"p ${file_paths_exfor[$N]} | awk -F")" '{print $1}' | awk -F"(" '{print $2}')
if [[ $YEAR > 20 ]]; then YEAR=$(echo 19${YEAR}); else YEAR=$(echo 20${YEAR}); fi
TITLE=$(echo ${AUTHOR} \(${YEAR}\), ${targets[M]}\(${beams[M]},${NOUT}n${POUT}p\)${products[M]})
TITLE=$(echo $TITLE | sed 's/ / /g;s/0p//g;s/1p/p/g;s/0n//g;s/1n/n/g;s/protons/p/g;s/alphas/a/g;s/deutrons/d/g;s/helions/h/g;')
# Plotting
if [ -n "${columns}" ]; then
[ -n "$cmd_exf" ] && cmd_exf=${cmd_exf}, || cmd_exf=''
cmd_exf="${cmd_exf} '${file_paths_exfor[$N]}' ${every} using ${columns} lc palette frac 0.$((0 + 4 * $P / 2)) title '$TITLE'"
fi
fi
done
done
done
fi
fi
echo EXFOR: sorting done
######################## EXPER DATA AND PLOT SEQUENCE #########################
echo -e "\nSorting EXPER data"
# Finding files EXPER/EXFOR sequence
if [[ ${data_bases[@]} == *exper* ]]; then
for each_dir in ${dir_exper}; do
for (( N = 0; N < ${#reactions[*]}; N++ )); do
# Finding product specific files
if [[ $data_plot == xs ]]; then
if [[ ${targetsMM[N]} == nat ]]; then
file_paths_exper+=($(find ${each_dir} -type f -iwholename ${each_dir}/\*nat${targetsXX[N]}\*${beams[N]}/\*${products[N]}\*))
else
file_paths_exper+=($(find ${each_dir} -type f -iwholename ${each_dir}/\*${reactions[N]}/\*${products[N]}\*))
fi
fi
# Finding data files for medical isotope production
if [[ $data_plot == production ]] || [[ $data_plot == activity ]] || [[ $data_plot == nuclides ]]; then
file_paths_exper+=($(find ${each_dir} -type f -iwholename ${each_dir}/\*${reactions[N]}/Y\*${products[N]}\*))
fi
if [[ ${#file_paths_exper[*]} == $no_files ]] || [[ ${#file_paths_exper[*]} == 0 ]]; then
echo -e "EXPER: invalid $data_plot product '${products[N]}' for reaction '${reactions[N]}'"
else no_files=${#file_paths_exper[*]}; fi
done
done
fi
file_paths_exper=($(for file_path in ${file_paths_exper[@]}; do echo $file_path; done | sort -u))
#for file_path in ${file_paths_exper[@]}; do echo "'$file_path'"; done
# EXPER/EXFOR SEQUENCE
# Checking if need to run sequense and correct cross section xs or all-xs
if [[ $data_plot == *xs* ]] && [[ ${data_bases[@]} == *exper* ]]; then
if [ ${#file_paths_exper[*]} != 0 ]; then
for (( N = 0; N < ${#file_paths_exper[*]}; N++ )) ; do
# Resolving reaction parameters from file to get correct plot titles
dir_name=$(basename $(dirname $(dirname ${file_paths_exper[$N]})))
base_name=$(basename ${file_paths_exper[$N]} | sed 's/.tot//g;s/.x4//g')
# Setting target with correct reaction abundance
target=$(echo $base_name | awk -F"," '{print $2}' | awk -F"(" '{print $1}' | sed 's/ //g')
for (( M = 0; M < ${#reactions[*]}; M++ )) ; do
if [[ ${targets[$M]} == ${target} ]]; then
targetPP=${targetsPP[M]}; break
else
targetPP=$((1))
fi
done
[ -n "$cmd_exp" ] && cmd_exp=${cmd_exp}, || cmd_exp=''
columns="(\$1*1E-6):(\$4*1E+3*$targetPP):((\$1*1E-6)-(\$2*1E-6)):((\$1*1E-6)+(\$3*1E-6)):(((\$4*1E+3)-(\$5*1E+3))*$targetPP):(((\$4*1E+3)+(\$6*1E+3))*$targetPP) w xyerrorbars"
cmd_exp="${cmd_exp} '${file_paths_exper[$N]}' using ${columns} lc palette frac 0.$((8 + 3 * $N / 5)) title '$base_name'"
done
fi
fi
echo EXPER: sorting done
####################### GNUPLOT DATA AND PLOT SEQUENCE ########################
echo -e "\nSorting GNUPLOT data"
# Finding files GNUPLOT sequence
if [[ ${data_bases[@]} == *gnuplot* ]]; then
for each_dir in ${dir_gnuplot[@]}; do
for (( N = 0; N < ${#reactions[*]}; N++ )); do
# Finding product specific files
if [[ $data_plot == xs ]] && [[ ${input_emissions[@]} == *x* ]]; then
file_paths_gnuplot+=($(find ${each_dir} -type f -iwholename ${each_dir}/\*${reactions[N]}/\*${products[N]}\*))
fi
if [[ ${#file_paths_gnuplot[*]} == $no_files ]] || [[ ${#file_paths_gnuplot[*]} == 0 ]]; then
echo -e "GNUPLOT: invalid $data_plot product '${products[N]}' for reaction '${reactions[N]}'"
else no_files=${#file_paths_gnuplot[*]}; fi
done
done
fi
file_paths_gnuplot=($(for file_path in ${file_paths_gnuplot[@]}; do echo $file_path; done | sort -u))
#for file_path in ${file_paths_gnuplot[@]}; do echo "'$file_path'"; done
# Checking if need to run sequense and correct cross section xs or all-xs
if [[ $data_plot == *xs* ]] && [[ ${data_bases[@]} == *gnuplot* ]]; then
if [ ${#file_paths_gnuplot[*]} != 0 ]; then
for (( N = 0; N < ${#file_paths_gnuplot[*]}; N++ )) ; do
# Resolving reaction parameters from file to get correct plot titles
dir_name=$(basename $(dirname $(dirname ${file_paths_gnuplot[$N]})))
base_name=$(basename ${file_paths_gnuplot[$N]} | sed 's/.tot//g')
# Setting target with correct reaction abundance
target=$(echo $base_name | awk -F"," '{print $2}' | awk -F"(" '{print $1}' | sed 's/ //g')
for (( M = 0; M < ${#reactions[*]}; M++ )) ; do
if [[ ${targets[$M]} == ${target} ]]; then
targetPP=${targetsPP[M]}; break
else
targetPP=$((1))
fi
done
[ -n "$cmd_gnu" ] && cmd_gnu=${cmd_gnu}, || cmd_gnu=''
columns="(\$1):(\$2*$targetPP) w lp pi 1"
cmd_gnu="${cmd_gnu} '${file_paths_gnuplot[$N]}' using ${columns} lc palette frac 0.$((9 + 3 * $N / 5)) title '$base_name'"
done
fi
fi
echo GNUPLOT: sorting done
######################### GNUPLOT PLOTLINE PARAMETERS #########################
# Checking plottable files, printing statistics
echo
echo -e "Found $(echo $cmd_emp | grep -o "title" | wc -l) EMPIRE data series to plot"
echo -e "Found $(echo $cmd_tal | grep -o "title" | wc -l) TALYS data series to plot"
echo -e "Found $(echo $cmd_ten | grep -o "title" | wc -l) TENDL data series to plot"
echo -e "Found $(echo $cmd_exf | grep -o "title" | wc -l) EXFOR data series to plot"
echo -e "Found $(echo $cmd_exp | grep -o "title" | wc -l) EXPER/EXFOR data series to plot"
echo -e "Found $(echo $cmd_gnu | grep -o "title" | wc -l) GNUPLOT data series to plot\n"
# MANAGING TITLES
# Managing plot titles
TIME=`date '+%Y.%m.%d-%H.%M'`
plot="Cross section data | ${data_plot[@]} | ${input_reaction[@]} | ${input_data_bases[@]^^} | $TIME"
plot_title="Cross section data: ${input_reaction[@]}| ${input_data_bases[@]^^} $data_plot cross section(s)"
# Managing plot labels and titels
if [[ ${input_reaction} == *x* ]]; then
data_plot=${data_plot,}
if [ $data_plot == production -o $data_plot == activity -o $data_plot == yield -o $data_plot == \#isotopes ]; then
xlabel="Time (h)"
plot=$(echo $plot | sed 's/Cross section data/Nuclide production data/g')
plot_title=$(echo $plot_title | sed 's/Cross section data/Nuclide production data/g')
plot_title=$(echo $plot_title | sed 's/cross section(s)//g')
else
xlabel="Energy (MeV)"; ylabel="Cross section (mb)"
fi
else
xlabel="Emission energy (MeV)"; ylabel="Cross section (mb)"
fi
# Cleaning plot titles
plot=$(echo $plot | sed 's/_/ /g;s/-/,/g;s/|/-/g' | sed -r 's/,+/,/g')
plot_title=$(echo $plot_title | sed 's/_/ /g;s/-/,/g;s/|/,/g' | sed -r 's/,+/,/g')
# MERGING CMD COMMANDS
# Merging commands for theoretical and experimental plots
if [ ! -z $cmd_emp ]; then cmd=($cmd_emp); fi
if [ ! -z $cmd_tal ]; then
if [ -z $cmd ]; then cmd=($cmd_tal)
else cmd=$(echo $cmd, $cmd_tal); fi; fi
if [ ! -z $cmd_ten ]; then
if [ -z $cmd ]; then cmd=($cmd_ten)
else cmd=$(echo $cmd, $cmd_ten); fi; fi
if [ ! -z $cmd_exf ]; then
if [ -z $cmd ]; then cmd=($cmd_exf)
else cmd=$(echo $cmd, $cmd_exf); fi; fi
if [ ! -z $cmd_exp ]; then
if [ -z $cmd ]; then cmd=($cmd_exp)
else cmd=$(echo $cmd, $cmd_exp); fi; fi
if [ ! -z $cmd_gnu ]; then
if [ -z $cmd ]; then cmd=($cmd_gnu)
else cmd=$(echo $cmd, $cmd_gnu); fi; fi
# Checking for errors in cmd
if [ -z $cmd ]; then
echo -e "Found no data to plot. Check reaction parameters. \n"; exit 2
else echo -e "Plotting..."
fi
# GENERATING GNUPLOT SCRIPT
# Setting Gnuplot script parameter for printing
if [ $gnuplot_output == eps -o $gnuplot_output == epslatex -o $gnuplot_output == svg -o $gnuplot_output == table ]; then
if [ $gnuplot_output == eps ]; then
command_print+=("set terminal postscript enhanced color font 'lmodern,10'
set output '$gnuplot_plot/$plot.eps'
replot
");fi
if [ $gnuplot_output == epslatex ]; then
command_print+=("set terminal epslatex size 8.0,4.0 color colortext
unset title
set output '$gnuplot_plot/$plot.tex'
replot
");fi
if [ $gnuplot_output == svg ]; then
command_print+=("set terminal svg size 800,400 fname 'lmodern' fsize 10
unset title
set output '$gnuplot_plot/$plot.svg'
replot
");fi
if [ $gnuplot_output == table ]; then
command_print+=("set table
set output '$gnuplot_plot/$plot.table'
replot
");fi
else
command_print="#set terminal svg size 800,400 fname 'lmodern' fsize 10
#set output '$gnuplot_plot/$plot.svg'
#replot"
fi
if [ -z $gnuplot_x_range ]; then
xrange="#set xrange [ 0 : 40 ]"
else
x_min=`echo $gnuplot_x_range | awk -F"-" '{print $1}'`
x_max=`echo $gnuplot_x_range | awk -F"-" '{print $2}'`
xrange="set xrange [ $x_min : $x_max ]"
fi
# Writing Gnuplot script and parameters
cat <<EOF > $gnuplot_temp/plotfile.gnu
set terminal wxt size 800,400 persist
set autoscale # scale axes automatically
$xrange # scale x-axis
unset log # remove any log-scaling
unset label # remove any previous labels
set xtic auto # set xtics automatically
set ytic auto # set ytics automatically
set title "$plot_title"
set xlabel "$xlabel"
set ylabel "$ylabel"
set key outside left bottom horizontal center Left reverse
set key noinvert samplen 0.2 spacing 1 width 0 height 2
unset colorbox
# MATLAB jet color palette
# line styles
set style line 1 lt 1 lc rgb '#000080' #
set style line 2 lt 1 lc rgb '#0000ff' #
set style line 3 lt 1 lc rgb '#0080ff' #
set style line 4 lt 1 lc rgb '#00ffff' #
set style line 5 lt 1 lc rgb '#80ff80' #
set style line 6 lt 1 lc rgb '#ffff00' #
set style line 7 lt 1 lc rgb '#ff8000' #
set style line 8 lt 1 lc rgb '#ff0000' #
set style line 9 lt 1 lc rgb '#800000' #
# palette
set palette defined (0 0.0 0.0 0.5, 1 0.0 0.0 1.0, 2 0.0 0.5 1.0, 3 0.0 1.0 1.0, \
4 0.5 1.0 0.5, 5 1.0 1.0 0.0, 6 1.0 0.5 0.0, 7 1.0 0.0 0.0, 8 0.5 0.0 0.0 )
plot $cmd
set table
set output '$gnuplot_temp/data-table.table'
replot
unset table
replot
#In case for building an output file
${command_print[@]}
EOF
# Writing Gnuplot plotline
cat <<EOF > $gnuplot_temp/plotline.gnu
plot $cmd
EOF
# Run plot command
case $gnuplot_output in
edit ) gnuplot $gnuplot_temp/plotfile.gnu - ; ;;
no-edit | eps | epslatex | svg | table ) gnuplot $gnuplot_temp/plotfile.gnu ; ;;
esac
echo -e "\nScript finished\n\n"
unset IFS
The reference
To view the code, you will need to expand the box below.
# File name: data-plot-ref.dat # Reference file. Comparison between nuclide names and Z-number # #X Name Z N A #Z X Name g gamma 0 0 0 n neutron 0 1 1 p proton 1 0 1 d deutron 1 1 2 t triton 1 2 3 h helion 2 1 3 a alpha 2 2 4 001 H Hydrogen 002 He Helium 003 Li Lithium 004 Be Beryllium 005 B Boron 006 C Carbon 007 N Nitrogen 008 O Oxygen 009 F Fluorine 010 Ne Neon 011 Na Sodium 012 Mg Magnesium 013 Al Aluminium 014 Si Silicon 015 P Phosphorus 016 S Sulfur 017 Cl Chlorine 018 Ar Argon 019 K Potassium 020 Ca Calcium 021 Sc Scandium 022 Ti Titanium 023 V Vanadium 024 Cr Chromium 025 Mn Manganese 026 Fe Iron 027 Co Cobalt 028 Ni Nickel 029 Cu Copper 030 Zn Zinc 031 Ga Gallium 032 Ge Germanium 033 As Arsenic 034 Se Selenium 035 Br Bromine 036 Kr Krypton 037 Rb Rubidium 038 Sr Strontium 039 Y Yttrium 040 Zr Zirconium 041 Nb Niobium 042 Mo Molybdenum 043 Tc Technetium 044 Ru Ruthenium 045 Rh Rhodium 046 Pd Palladium 047 Ag Silver 048 Cd Cadmium 049 In Indium 050 Sn Tin 051 Sb Antimony 052 Te Tellurium 053 I Iodine 054 Xe Xenon 055 Cs Caesium 056 Ba Barium 057 La Lanthanum 058 Ce Cerium 059 Pr Praseodymium 060 Nd Neodymium 061 Pm Promethium 062 Sm Samarium 063 Eu Europium 064 Gd Gadolinium 065 Tb Terbium 066 Dy Dysprosium 067 Ho Holmium 068 Er Erbium 069 Tm Thulium 070 Yb Ytterbium 071 Lu Lutetium 072 Hf Hafnium 073 Ta Tantalum 074 W Tungsten 075 Re Rhenium 076 Os Osmium 077 Ir Iridium 078 Pt Platinum 079 Au Gold 080 Hg Mercury 081 Tl Thallium 082 Pb Lead 083 Bi Bismuth 084 Po Polonium 085 At Astatine 086 Rn Radon 087 Fr Francium 088 Ra Radium 089 Ac Actinium 090 Th Thorium 091 Pa Protactinium 092 U Uranium 093 Np Neptunium 094 Pu Plutonium 095 Am Americium 096 Cm Curium 097 Bk Berkelium 098 Cf Californium 099 Es Einsteinium 100 Fm Fermium 101 Md Mendelevium 102 No Nobelium 103 Lr Lawrencium 104 Rf Rutherfordium 105 Db Dubnium 106 Sg Seaborgium 107 Bh Bohrium 108 Hs Hassium 109 Mt Meitnerium 110 Ds Darmstadtium 111 Rg Roentgenium 112 Cn Copernicium 113 Uut Ununtrium 114 Uuq Ununquadium 115 Uup Ununpentium 116 Uuh Ununhexium 118 Uuo Ununoctium
Advertisements
