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.
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