Kinematics correlations

benchmarks/dis/analysis/kinematics_correlations.py

+#!/usr/bin/env python
+# coding: utf-8
+
+import os
+import numpy as np
+import uproot as ur
+import awkward as ak
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+import mplhep 
+import argparse
+
+
+
+parser = argparse.ArgumentParser()
+parser.add_argument('--rec_file', type=str, help='Reconstructed track file.')
+parser.add_argument('--ebeam', type=float, help='Electron beam energy.')
+parser.add_argument('--pbeam', type=float, help='Proton (or ion) beam energy.')
+parser.add_argument('--minq2', type=float, help='Minimum four-momentum transfer squared Q2.')
+parser.add_argument('-o', dest='outdir', default='results/dis/', help='Output directory.')
+args = parser.parse_args()
+kwargs = vars(args)
+
+rec_file = args.rec_file
+minq2 = int(args.minq2)
+k = int(args.ebeam)
+p = int(args.pbeam)
+
+
+
+#logarithmically spaced bins in Q2 and 𝑥 
+Q2_bins = [0.40938507,0.64786184,1.025353587,1.626191868,2.581181894,4.091148983,6.478618696,10.25353599,16.26191868,25.81181894,40.91148983,64.78618696,102.5353599,162.6191868,258.1181894,409.1148983,647.8618696,1025.353599,1626.191868,2581.181894,4091.148983,6482.897648]
+x_bins = [4.09385E-05,6.47862E-05,0.000102535,0.000162619,0.000258118,0.000409115,0.000647862,0.001025354,0.001626192,0.002581182,0.004091149,0.006478619,0.010253536,0.016261919,0.025811819,0.04091149,0.064786187,0.10253536,0.162619187,0.25811819,0.409114898,0.647861868,1.025257131]
+
+
+
+#function to construct Q2 correlation plots
+def Q2correlation(minq2,method): #minq2 can be 1,10,100, or 1000; method can be 'E','DA', or 'JB'
+
+    Q2values_Y = method_Q2values_dict['{}'.format(method)]  #Q2 values of the given method, that are mapped onto the y axis 
+    
+    Q2_List_T = Q2values_T['{}'.format(minq2)] #Truth Q2 values for given minq2, mapped along x axis
+    Q2_List_Y = Q2values_Y['{}'.format(minq2)] #method (E/DA/JB) Q2 values for given minq2, mapped along y axis
+
+    T_len = ak.count(Q2_List_T,axis=0) #total number of events in Truth
+    Y_len = ak.count(Q2_List_Y,axis=0) #total number of events in method
+
+    if T_len > Y_len: #if total number of events for Truth is greater
+        Y_boolean = ak.count(Q2_List_Y,axis=-1) >= 1 #boolean to filter ak.Arrays wrt single events in method 
+        Q2_List_T_F = Q2_List_T[Y_boolean] #filtered Truth Q2 values
+        Q2_List_Y_F = Q2_List_Y[Y_boolean] #filtered method Q2 values
+    else: #if total number of events for method is greater
+        T_boolean = ak.count(Q2_List_T,axis=-1) >= 1 #boolean to filter ak.Arrays wrt single events in Truth
+        Q2_List_T_F = Q2_List_T[T_boolean] #filtered Truth Q2 values
+        Q2_List_Y_F = Q2_List_Y[T_boolean] #filtered method Q2 values
+    
+    T_Q2s = np.array(ak.flatten(Q2_List_T_F)) #Truth Q2 values, mapped along x axis
+    Y_Q2s = np.array(ak.flatten(Q2_List_Y_F)) #methos Q2 values, mapped along y axis
+    
+    #2-dimensional histogram, h
+    h, xedges, yedges = np.histogram2d(x=T_Q2s,y=Y_Q2s, bins=[Q2_bins,Q2_bins]) 
+
+    minq2_dict = {'1':2,'10':7,'100':12,'1000':17} #Q2 bin index at which minq2 starts
+    h[0:minq2_dict['{}'.format(minq2)]]=0 #ignore values before minq2
+
+    #normalization of h:
+    col_sum = ak.sum(h,axis=-1) #number of events in each (verticle) column 
+    norm_h = [] #norm_h is the normalized matrix
+    norm_h_text = [] #display labels matrix
+    for i in range(len(col_sum)):
+        if col_sum[i] != 0:
+            norm_c = h[i]/col_sum[i] #normalized column = column values divide by sum of the column
+        else:
+            norm_c = h[i]
+        norm_h.append(norm_c)
+        norm_c_text = [ '%.3f' % elem for elem in norm_c ] #display value to 3 dp
+        norm_h_text.append(norm_c_text)
+    
+    fig = plt.figure()
+    mplhep.hist2dplot(H=norm_h,norm=mpl.colors.LogNorm(),labels=norm_h_text,xbins=Q2_bins,ybins=Q2_bins,cmax=1,cmin=1e-5)
+    plt.yscale('log')
+    plt.xscale('log')
+    fig.set_figwidth(11)
+    fig.set_figheight(11)
+    plt.xlabel('$Q^2$ [$GeV^2$] Truth')
+    plt.ylabel('$Q^2$ [$GeV^2$] {}'.format(method_dict['{}'.format(method)]))
+    plt.title('{}   $Q^2$ correlation   {}x{}   $minQ^2=${}$GeV^2$'.format(method_dict['{}'.format(method)],k,p,minq2))
+    plt.show()
+    plt.savefig(os.path.join(args.outdir, '%gon%g/Q2_correlation_%s_%gx%g_minQ2=%g.png' %(k,p,method,k,p,minq2)))
+
+
+
+#function to construct Bjorken-x correlation plots
+def Xcorrelation(minq2,method): #minq2 can be 1,10,100, or 1000; method can be 'E','DA', or 'JB'
+
+    Xvalues_Y = method_Xvalues_dict['{}'.format(method)] #x values of the given method, that are mapped onto the y axis 
+        
+    X_List_T = Xvalues_T['{}'.format(minq2)] #Truth x values for given minq2, mapped along x axis
+    X_List_Y = Xvalues_Y['{}'.format(minq2)] #method (E/DA/JB) x values for given minq2, mapped along y axis
+
+    T_len = ak.count(X_List_T,axis=0) #total number of events in Truth
+    Y_len = ak.count(X_List_Y,axis=0) #total number of events in method
+
+    if T_len > Y_len: #if total number of events for Truth is greater
+        Y_boolean = ak.count(X_List_Y,axis=-1) >= 1 #boolean to filter ak.Arrays wrt single events in method 
+        X_List_T_F = X_List_T[Y_boolean] #filtered Truth x values
+        X_List_Y_F = X_List_Y[Y_boolean] #filtered method x values
+    else: #if total number of events for method is greater
+        T_boolean = ak.count(X_List_T,axis=-1) >= 1 #boolean to filter ak.Arrays wrt single events in Truth
+        X_List_T_F = X_List_T[T_boolean] #filtered Truth x values
+        X_List_Y_F = X_List_Y[T_boolean] #filtered method x values
+    
+    T_Xs = np.array(ak.flatten(X_List_T_F)) #Truth Bjorken-x values, mapped along x axis
+    Y_Xs = np.array(ak.flatten(X_List_Y_F)) #method Bjorken-x values, mapped along y axis
+    
+    T_x_bool = T_Xs>=minq2/(4*k*p) #boolean to filter x values that satisfy bjorken-x equation for minq2, ebeam and pbeam
+    T_Xs = T_Xs[T_x_bool]
+    Y_Xs = Y_Xs[T_x_bool]
+
+    #2-dimensional histogram, h
+    h, xedges, yedges = np.histogram2d(x=T_Xs,y=Y_Xs, bins=[x_bins,x_bins])
+    
+    #normalization of h:
+    col_sum = ak.sum(h,axis=-1) #number of events in each (verticle) column 
+    norm_h = [] #norm_h is the normalized matrix
+    norm_h_text = [] #display labels matrix
+    for i in range(len(col_sum)):
+        if col_sum[i] != 0:
+            norm_c = h[i]/col_sum[i] #normalized column = column values divide by sum of the column
+        else:
+            norm_c = h[i]
+        norm_h.append(norm_c)
+        norm_c_text = [ '%.2f' % elem for elem in norm_c ] #display value to 2 dp
+        norm_h_text.append(norm_c_text)
+
+    fig = plt.figure()
+    mplhep.hist2dplot(H=norm_h,norm=mpl.colors.LogNorm(),labels=norm_h_text,xbins=x_bins,ybins=x_bins,cmax=1,cmin=1e-5)
+    plt.yscale('log')
+    plt.xscale('log')
+    fig.set_figwidth(11)
+    fig.set_figheight(11)
+    plt.xlabel('x Truth')
+    plt.ylabel('$x$   {}'.format(method_dict['{}'.format(method)]))
+    plt.title('{}   $x$ correlation   {}x{}   $minQ^2=${}$GeV^2$'.format(method_dict['{}'.format(method)],k,p,minq2))
+    plt.show()
+    plt.savefig(os.path.join(args.outdir, '%gon%g/x_correlation_%s_%gx%g_minQ2=%g.png' %(k,p,method,k,p,minq2)))    
+
+
+
+keys = ur.concatenate(rec_file + ':events/' + 'InclusiveKinematicsTruth')
+minq2_1_T =   [keys['InclusiveKinematicsTruth.Q2'],keys['InclusiveKinematicsTruth.x']]
+keys = ur.concatenate(rec_file + ':events/' + 'InclusiveKinematicsElectron')
+minq2_1_E =   [keys['InclusiveKinematicsElectron.Q2'], keys['InclusiveKinematicsElectron.x']]
+keys = ur.concatenate(rec_file + ':events/' + 'InclusiveKinematicsDA')
+minq2_1_DA =  [keys['InclusiveKinematicsDA.Q2'], keys['InclusiveKinematicsDA.x']]
+keys = ur.concatenate(rec_file + ':events/' + 'InclusiveKinematicsJB')
+minq2_1_JB =  [keys['InclusiveKinematicsJB.Q2'], keys['InclusiveKinematicsJB.x']]
+
+Q2values_T = {'1':minq2_1_T[0]}
+Q2values_E = {'1':minq2_1_E[0]}
+Q2values_DA = {'1':minq2_1_DA[0]}
+Q2values_JB = {'1':minq2_1_JB[0]}
+Xvalues_T = {'1':minq2_1_T[1]}
+Xvalues_E = {'1':minq2_1_E[1]}
+Xvalues_DA = {'1':minq2_1_DA[1]}
+Xvalues_JB = {'1':minq2_1_JB[1]}
+
+method_dict = {'e':'Electron','DA':'Double-Angle','JB':'Jacquet-Blondel'}
+method_Q2values_dict = {'e':Q2values_E,'DA':Q2values_DA,'JB':Q2values_JB}
+method_Xvalues_dict = {'e':Xvalues_E,'DA':Xvalues_DA,'JB':Xvalues_JB}
+
+Q2correlation(minq2,'e')
+Xcorrelation(minq2,'e')
+Q2correlation(minq2,'DA')
+Xcorrelation(minq2,'DA')
+Q2correlation(minq2,'JB')
+Xcorrelation(minq2,'JB')

benchmarks/dis/config.yml

@@ -12,13 +12,16 @@ dis:generate:
     matrix:
       - EBEAM: 5
         PBEAM: 41
+        MINQ2: 1
       - EBEAM: 10
         PBEAM: 100
+        MINQ2: 1
       - EBEAM: 18
         PBEAM: 275
+        MINQ2: 1
   timeout: 1 hours
   script:
-    - bash benchmarks/dis/get.sh --config dis_${EBEAM}x${PBEAM} --ebeam ${EBEAM} --pbeam ${PBEAM}
+    - bash benchmarks/dis/get.sh --config dis_${EBEAM}x${PBEAM}_minQ2=${MINQ2} --ebeam ${EBEAM} --pbeam ${PBEAM} --minq2 ${MINQ2}
 
 dis:simulate:
   stage: simulate
@@ -28,13 +31,16 @@ dis:simulate:
     matrix:
       - EBEAM: 5
         PBEAM: 41
+        MINQ2: 1
       - EBEAM: 10
         PBEAM: 100
+        MINQ2: 1
       - EBEAM: 18
         PBEAM: 275
+        MINQ2: 1
   timeout: 2 hour
   script:
-    - bash benchmarks/dis/dis.sh --config dis_${EBEAM}x${PBEAM} --ebeam ${EBEAM} --pbeam ${PBEAM}
+    - bash benchmarks/dis/dis.sh --config dis_${EBEAM}x${PBEAM}_minQ2=${MINQ2} --ebeam ${EBEAM} --pbeam ${PBEAM} --minq2 ${MINQ2}
   retry:
     max: 2
     when:

benchmarks/dis/dis.sh

@@ -23,6 +23,7 @@ echo "Running the DIS benchmarks"
 ## - CONFIG:   The specific generator configuration
 ## - EBEAM:    The electron beam energy
 ## - PBEAM:    The ion beam energy
+export REQUIRE_MINQ2=true
 source ${LOCAL_PREFIX}/bin/parse_cmd.sh $@
 
 ## To run the reconstruction, we need the following global variables:
@@ -110,6 +111,7 @@ cat << EOF > ${CONFIG}
   "output_prefix": "${RESULTS_PATH}/${PLOT_TAG}",
   "ebeam": ${EBEAM},
   "pbeam": ${PBEAM},
+  "minq2": ${MINQ2},
   "test_tag": "${BEAM_TAG}"
 }
 EOF
@@ -120,6 +122,12 @@ if [[ "$?" -ne "0" ]] ; then
   exit 1
 fi
 
+python benchmarks/dis/analysis/kinematics_correlations.py --rec_file ${REC_FILE} --ebeam ${EBEAM} --pbeam ${PBEAM} --minq2 ${MINQ2}
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running kinematics_correlations script"
+  exit 1
+fi
+
 CONFIG="${TMP_PATH}/${PLOT_TAG}.raw.json"
 cat << EOF > ${CONFIG}
 {
@@ -128,6 +136,7 @@ cat << EOF > ${CONFIG}
   "output_prefix": "${RESULTS_PATH}/${PLOT_TAG}",
   "ebeam": ${EBEAM},
   "pbeam": ${PBEAM},
+  "minq2": ${MINQ2},
   "test_tag": "${BEAM_TAG}"
 }
 EOF
@@ -145,6 +154,7 @@ cat << EOF > ${CONFIG}
   "output_prefix": "${RESULTS_PATH}/${PLOT_TAG}",
   "ebeam": ${EBEAM},
   "pbeam": ${PBEAM},
+  "minq2": ${MINQ2},
   "test_tag": "${BEAM_TAG}"
 }
 EOF
@@ -162,6 +172,7 @@ cat << EOF > ${CONFIG}
   "output_prefix": "${RESULTS_PATH}/${PLOT_TAG}",
   "ebeam": ${EBEAM},
   "pbeam": ${PBEAM},
+  "minq2": ${MINQ2},
   "test_tag": "${BEAM_TAG}"
 }
 EOF

benchmarks/dis/get.sh

@@ -21,6 +21,7 @@ pushd ${PROJECT_ROOT}
 ## - CONFIG:   The specific generator configuration --> not currenlty used FIXME
 ## - EBEAM:    The electron beam energy --> not currently used FIXME
 ## - PBEAM:    The ion beam energy --> not currently used FIXME
+export REQUIRE_MINQ2=true
 source parse_cmd.sh $@
 
 ## To run the generator, we need the following global variables:
@@ -56,7 +57,7 @@ fi
 ## =============================================================================
 ## Step 3: Copy the file (about 180 lines per event in DIS NC files)
 nlines=$((180*${JUGGLER_N_EVENTS}))
-DATA_URL=S3/eictest/ATHENA/EVGEN/DIS/NC/${EBEAM}x${PBEAM}/minQ2=1/pythia8NCDIS_${EBEAM}x${PBEAM}_minQ2=1_beamEffects_xAngle=-0.025_hiDiv_vtxfix_1.hepmc
+DATA_URL=S3/eictest/ATHENA/EVGEN/DIS/NC/${EBEAM}x${PBEAM}/minQ2=${MINQ2}/pythia8NCDIS_${EBEAM}x${PBEAM}_minQ2=${MINQ2}_beamEffects_xAngle=-0.025_hiDiv_vtxfix_1.hepmc
 mc config host add S3 https://dtn01.sdcc.bnl.gov:9000 ${S3_ACCESS_KEY} ${S3_SECRET_KEY}
 mc head -n ${nlines} ${DATA_URL} | sanitize_hepmc3 > ${TMP_PATH}/${GEN_TAG}.hepmc
 if [[ "$?" -ne "0" ]] ; then
@@ -67,7 +68,7 @@ fi
 ## =============================================================================
 ## Step 4: Finally, move relevant output into the artifacts directory and clean up
 ## =============================================================================
-echo "Moving generator output into ${INPUT_PATH}"
+echo "Moving generator output to ${INPUT_PATH}/${GEN_TAG}.hepmc"
 mv ${TMP_PATH}/${GEN_TAG}.hepmc ${INPUT_PATH}/${GEN_TAG}.hepmc
 ## this step only matters for local execution
 echo "Cleaning up"