Determine URQMD to be K.E. or momomentum

As per redmine #2684, determine wether number in URQMD filename is kinetic energy or momentum of beam.

scripts/README.md

+## check_pEs.sh
+
+Change flag for more or less output.
+
+May be used to determine if the number in a URQMD filename  is the kinetic energy or momentum of the beam. e.g.
+
+For single file
+```
+sh ./scripts/check_pEs.sh urqmd.nini.1.93gev.mbias.00001.root
+```
+
+For all files in directory
+```
+for file in `ls urqmd*root`; do sh scripts/check_pEs.sh $file; done
+```
+

scripts/check_pEs.sh

+#!/bin/bash
+
+# Author: Eoin Clerkin (FAIR GMBH)
+# Gnu Public License v3
+
+# Checks a urqmd root file to determine whether number in the file name is the beam kinetic energy or momentum per nucleon.
+
+# Needs Cern Root installed.
+
+DEBUG=true;   # Verbosity
+
+# Below is the information stored in the header of a typical uqmd input root file
+: << 'HEADER_OF_ROOT_FILE'
+Decayer                       :                                                                                 
+Projectile mass               : 107                                                                             
+Projectile charge             : 47                                                                              
+Projectile momentum (AGeV/c)  : 1.45286                                                                         
+Target mass                   : 107                                                                             
+Target charge                 : 47                                                                              
+Target momentum (AGeV/c)      : -1.45286                                                                        
+Minimal impact parameter (fm) : 0                                                                               
+Maximal impact parameter (fm) : 13.34                                                                           
+Impact parameter weightning   : -1                                                                              
+Minimal azimuthal angle (rad) : 0                                                                               
+Maximal azimuthal angle (rad) : 0                                                                               
+Cross-section (mb)            : 5586.73                                                                         
+Requested number of events    : 10000
+HEADER_OF_ROOT_FILE
+
+# The momenta of the target and projectile are recorded in a center of momentum coordinate frame.
+# It will be the task of this program to convert this into a coordinate frame where the target is
+# stationary and only the projectile moves.
+
+# In order to achieve this, I will convert the momenta (same) into a velocity.
+# I will them subtract the "negative" velocity of target from the projectile to get the velocity of the projectile in the lab frame.
+# I will convert this back to a momentum for the projectile in the lab frame
+# I will determine the kinetic energy from this.
+# I then compare these calculated kinetic energy and momentum of the projectile in the lab frame to the number in the filename.
+
+
+# Extract information from header
+root -b $1 2>&1 1>/tmp/pEs << EOT
+_file0->GetListOfKeys()->Print();
+_file0->Get("run")->Print();
+EOT
+
+
+
+
+
+
+${DEBUG} && cat /tmp/pEs
+
+
+
+
+
+# Some constants for mass of protron and neutron.
+MASS_OF_PROTON=0.938271998 # GeV/c^2
+MASS_OF_NEUTRON=0.93956536 # GeV/c^2
+
+
+# Obtaining information for the projectile from the root file header
+PROJECTILE_MOMENTUM_CMFRAME=`sed -nE 's/Projectile momentum \(AGeV\/c\)[ ]*: //p' /tmp/pEs`
+${DEBUG} && echo "Projectile momentum in CM frame is " ${PROJECTILE_MOMENTUM_CMFRAME} "AGeV/c"
+
+PROJECTILE_Z=`sed -nE 's/Projectile charge [ ]*: //p' /tmp/pEs`
+${DEBUG} && echo "Projectile charge is taken as the number of proton which is " ${PROJECTILE_Z} "."
+
+PROJECTILE_A=`sed -nE 's/Projectile mass [ ]*: //p' /tmp/pEs`
+${DEBUG} && echo "Projectile mass is taken as the  number of proton and neutron " ${PROJECTILE_A} "."
+
+PROJECTILE_MASS_PER_NUCLEON=`echo "((${PROJECTILE_A}-${PROJECTILE_Z})*${MASS_OF_NEUTRON}+${PROJECTILE_Z}*${MASS_OF_PROTON})/${PROJECTILE_A}" | bc -l`
+${DEBUG} && echo "Average mass of a nucleon for the projectile is " ${PROJECTILE_MASS_PER_NUCLEON} " GeV/c^2."
+
+# Here the standard p = gamma mass is inverted. Mass is in A*average_mass*GeV
+PROJECTILE_VGAMMA=`echo "${PROJECTILE_MOMENTUM_CMFRAME}/${PROJECTILE_MASS_PER_NUCLEON}" | bc -l`
+${DEBUG} && echo "Vgamma(V) for projectile is is ${PROJECTILE_VGAMMA}"
+
+# I will use a Newton Rhapson method to obtain the v from this.
+
+# Newton Raphons Method
+# f(x) = x gamma(x) - vGamma = x/sqrt(1-x*x) - vGamma
+# -> f'(x) = 1/sqrt(1-x*x) + x*x/(1-x*x)^(3/2)         NO   = (1/x)*( f(x) + vGamma + (f(x) + vGamma )^3 )
+# x_(n+1) = x_n - f(x_n)/f'(x_n)
+
+PROJECTILE_VELOCITY_CMFRAME=`cat << EOT | bc -ql
+define f(x) {return (x/sqrt(1-x*x) - ${PROJECTILE_VGAMMA});}
+define fp(x) {return (1/sqrt(1-x*x) + (x/sqrt(1-x*x))*(x/(1-x*x)) );}
+x=0.5;
+diff=(f(x)/fp(x));
+while(diff*diff>0.00000000000000001){
+x = x - 0.01*diff;
+diff=(f(x)/fp(x));
+}
+x
+EOT`
+
+${DEBUG} && echo "Ergo velocity of the projectile in the CM frame is ${PROJECTILE_VELOCITY_CMFRAME}c."
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+${DEBUG} && echo "############################################################################################"
+
+
+
+# The simulation will be the target and the projectile moving with the same velocity towards each other.
+# we take the target only for sake of briefity.
+
+TARGET_MOMENTUM_CMFRAME=`sed -nE 's/Target momentum \(AGeV\/c\)[ ]*: //p' /tmp/pEs`
+${DEBUG} && echo "Target momentum in CM frame is " ${TARGET_MOMENTUM_CMFRAME} "AGeV/c"
+
+TARGET_Z=`sed -nE 's/Target charge [ ]*: //p' /tmp/pEs`
+${DEBUG} && echo "Target charge is taken as the number of proton which is " ${TARGET_Z} "."
+
+TARGET_A=`sed -nE 's/Target mass [ ]*: //p' /tmp/pEs`
+${DEBUG} && echo "Target mass is taken as the  number of proton and neutron " ${TARGET_A} "."
+
+TARGET_MASS_PER_NUCLEON=`echo "((${TARGET_A}-${TARGET_Z})*${MASS_OF_NEUTRON}+${TARGET_Z}*${MASS_OF_PROTON})/${TARGET_A}" | bc -l`
+${DEBUG} && echo "Average mass of a nucleon for the target is " ${TARGET_MASS_PER_NUCLEON} " GeV/c^2."
+
+
+# Here the standard p = gamma mass is inverted. Mass is in A*average_mass*GeV
+TARGET_vGamma=`echo "${TARGET_MOMENTUM_CMFRAME}/${TARGET_MASS_PER_NUCLEON}" | bc -l`
+${DEBUG} && echo "Vgamma(V) for Target is is ${TARGET_vGamma}"
+
+# I will use a Newton Rhapson method to obtain the v from this.
+
+# Newton Raphons Method
+# f(x) = x gamma(x) - vGamma = x/sqrt(1-x*x) - vGamma
+# -> f'(x) = 1/sqrt(1-x*x) + x*x/(1-x*x)^(3/2)         NO   = (1/x)*( f(x) + vGamma + (f(x) + vGamma )^3 )
+# x_(n+1) = x_n - f(x_n)/f'(x_n)
+
+TARGET_VELOCITY_CMFRAME=`cat << EOT | bc -ql
+define f(x) {return (x/sqrt(1-x*x) - ${TARGET_vGamma});}
+define fp(x) {return (1/sqrt(1-x*x) + (x/sqrt(1-x*x))*(x/(1-x*x)) );}
+x=0.5;
+diff=(f(x)/fp(x));
+while(diff*diff>0.00000000000000001){
+x = x - 0.01*diff;
+diff=(f(x)/fp(x));
+}
+x
+EOT`
+
+${DEBUG} && echo "Ergo velocity of the Target in the CM frame is ${TARGET_VELOCITY_CMFRAME}c."
+
+${DEBUG} && echo "############################################################################################"
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+# Here summations of velocities in special relativitiy  w = (u + v)/(1+uv)
+# where v is to be found       w*( 1 + u*v ) = u + v
+# w - u = v - u*v*w = v ( 1 - u*w)
+# v = (w-u)/(1-u*w)
+
+PROJECTILE_VELOCITY_LABFRAME=`echo "( ${PROJECTILE_VELOCITY_CMFRAME} - 1*${TARGET_VELOCITY_CMFRAME} )/( 1 - ${TARGET_VELOCITY_CMFRAME}*${PROJECTILE_VELOCITY_CMFRAME} )" | bc -l`
+
+
+${DEBUG} && echo "Velocity of projectile in the CM frame is ${PROJECTILE_VELOCITY_CMFRAME}"
+${DEBUG} && echo "Velocity of projectile in the LAB frame is ${PROJECTILE_VELOCITY_LABFRAME}"
+
+GAMMA=`echo "1/sqrt(1-${PROJECTILE_VELOCITY_LABFRAME}*${PROJECTILE_VELOCITY_LABFRAME})" | bc -l`
+${DEBUG} && echo "GAMMA is $GAMMA"
+
+PROJECTILE_MOMENTUM_LABFRAME=`echo "$GAMMA*${PROJECTILE_MASS_PER_NUCLEON}*${PROJECTILE_VELOCITY_LABFRAME}" | bc -l`
+
+PROJECTILE_MOMENTUM_LABFRAME=`echo "0.005 + ${PROJECTILE_MOMENTUM_LABFRAME}" | bc -l | sed 's_\([0-9]*\.[0-9][0-9]\).*_\1_'`   # To two decimal places.
+
+
+${DEBUG} && echo "***************************"
+${DEBUG} && echo "Projectile in the lab frame"
+
+${DEBUG} && echo "Momentum: ${PROJECTILE_MOMENTUM_LABFRAME} AGeV/c"
+
+PROJECTILE_EK_LABFRAME=`echo "sqrt(${PROJECTILE_MOMENTUM_LABFRAME}*${PROJECTILE_MOMENTUM_LABFRAME} + ${PROJECTILE_MASS_PER_NUCLEON}*${PROJECTILE_MASS_PER_NUCLEON}) - ${PROJECTILE_MASS_PER_NUCLEON} " | bc -l`
+PROJECTILE_EK_LABFRAME=`echo "0.005 + ${PROJECTILE_EK_LABFRAME}" | bc -l | sed 's_\([0-9]*\.[0-9][0-9]\).*_\1_'`   # To two decimal places.
+
+${DEBUG} && echo "Kinetic energy is ${PROJECTILE_EK_LABFRAME} AGeV"
+${DEBUG} && echo "################################################"
+
+
+
+
+NUM_IN_FILENAME=`
+echo $1 \
+| sed 's/.*urqmd\.[a-z]*\.\([0-9\.]*\)gev\..*\.root/\1/'
+`
+
+${DEBUG} && echo "Number extracted from the filename is ${NUM_IN_FILENAME}"
+
+# If NUM_IN_FILE is closer to Kinetic energy then second terms is smaller and thus expression is positive, ergo momentum
+${DEBUG} && echo "(${PROJECTILE_MOMENTUM_LABFRAME}-${NUM_IN_FILENAME})^2 - (${PROJECTILE_EK_LABFRAME}-${NUM_IN_FILENAME})^2 > 0"
+
+
+${DEBUG} && if [ \
+$(echo "(${PROJECTILE_MOMENTUM_LABFRAME}-${NUM_IN_FILENAME})^2 - (${PROJECTILE_EK_LABFRAME}-${NUM_IN_FILENAME})^2 > 0" | bc) -eq 1 ]; then
+echo "Number in the filename best fits the Kinetic energy of the beam.";
+else
+echo "Number in the filename best fits the Momentum of the beam.";
+fi
+
+
+
+
+
+# Final Output.  
+#filename momentum kinetic_energy num_extracted determine
+if ! ${DEBUG}; then
+echo -n $1"    "
+echo -n ${PROJECTILE_MOMENTUM_LABFRAME}"    "
+echo -n ${PROJECTILE_EK_LABFRAME}"    "
+echo -n ${NUM_IN_FILENAME}"    "
+if [ \
+$(echo "(${PROJECTILE_MOMENTUM_LABFRAME}-${NUM_IN_FILENAME})^2 - (${PROJECTILE_EK_LABFRAME}-${NUM_IN_FILENAME})^2 > 0" | bc) -eq 1 ]; then
+echo "kinetic";
+else
+echo "momentum";
+fi
+fi
+