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

 

 

Advertisements