make code robust against broken events

SingleGammaPt_samples.txt

@@ -0,0 +1,3 @@
+_SingleGammaPt100Eta1p6_2p8_PhaseIITDRFall17DR-PU200FEVT_93X_upgrade2023_realistic_v2-v1_GEN-SIM-RECO
+_SingleGammaPt35Eta1p6_2p8_PhaseIITDRFall17DR-noPUFEVT_93X_upgrade2023_realistic_v2-v1_GEN-SIM-RECO
+_SingleGammaPt50Eta1p6_2p8_PhaseIITDRFall17DR-PU200FEVT_93X_upgrade2023_realistic_v2-v1_GEN-SIM-RECO

hgcalHelpers.py

@@ -106,8 +106,9 @@ def getHighestEnergyObjectIndex(ref_etaphi, obj_etaphi, obj_energy, deltaR=0.1):
         # Handle the -pi pi transition
         matched_sym = kdtree.query_ball_point([row.eta, row.phi-np.sign(row.phi)*2.*math.pi], deltaR)
         matched = np.unique(np.concatenate((matched, matched_sym))).astype(int)
-        best_match = np.argmax(obj_energy[matched])
-        matched_indices[index] = best_match
+        if len(matched) > 0:
+            best_match = np.argmax(obj_energy[matched])
+            matched_indices[index] = best_match
     return matched_indices
 
 

megaClustering.py

@@ -15,6 +15,7 @@
 maxlayer = 52
 energyWeights = list(repeat(1.02, 28)) + list(repeat(0.86, 12)) + list(repeat(1.12, 12))
 
+
 def getConeRadius(frontRadius, backRadius, z, maxval=9999.):
     depthTerm = backRadius * (abs(z)-320.7)/(407.8-320.7)
     val = frontRadius + depthTerm
@@ -39,13 +40,18 @@ def getMegaClusters(genParticles, multiClusters, layerClusters, recHits, gun_typ
         selectedGen = selectedGen[(selectedGen.energy >= GEN_engpt*.999)]
     # print selectedGen
 
+    if multiClusters.shape[0] == 0:
+        return pd.DataFrame(columns=['pt', 'eta', 'phi', 'energy'])
     # for the mega cluster axis, take highest energy multicluster within dR = 0.1
     bestMultiClusterIndices = hgcalHelpers.getHighestEnergyObjectIndex(selectedGen[['eta', 'phi']], multiClusters[['eta', 'phi']], multiClusters['energy'], 0.1)
     # print bestMultiClusterIndices
 
     megaClusters = []
 
     for idx, genPart in selectedGen.iterrows():
+        if idx not in bestMultiClusterIndices:
+            # this is the case if there's no match between gen and multi clusters
+            continue
         matchedMultiCluster = multiClusters.iloc[[bestMultiClusterIndices[idx]]]
         energySum = 0
         pTSum = 0

resolutionScaleFiller_submit.sh

@@ -2,7 +2,7 @@
 
 SAMPLEFILE=$1
 
-OUTDIR="/eos/user/c/clange/HGCal/ScaleResolution/"
+OUTDIR="/eos/user/c/clange/HGCal/ScaleResolution/_new/"
 PREFIX="root://eoscms.cern.ch//eos/cms/store/cmst3/group/hgcal/CMG_studies/Production/"
 GUNTYPE="pt"
 POSTFIX="NTUP"
@@ -12,13 +12,13 @@ QUEUE="8nh"
 
 for SAMPLE in `cat $SAMPLEFILE`; do
   SAMPLEDIR="${PREFIX}${SAMPLE}/${POSTFIX}/"
-  PTVAL=`echo ${SAMPLE} | gawk 'match($0, /.*SinglePiPt(.*)Eta.*/, arr) { print arr[1] }'`
   PARTICLE=`echo ${SAMPLE} | gawk 'match($0, /.*Single(.*)Pt.*/, arr) { print arr[1] }'`
+  PTVAL=`echo ${SAMPLE} | gawk -v part=".*Single${PARTICLE}Pt(.*)Eta.*" 'match($0, part, arr) { print arr[1] }'`
   TAG=`echo ${SAMPLE} | gawk 'match($0, /.*Fall17DR-(.*)FEVT.*/, arr) { print arr[1] }'`
   PID=0
   if [ "$PARTICLE" == "Pi" ]; then
     PID=211
-  elif [ "$PARTICLE" == "Photon" ]; then
+  elif [ "$PARTICLE" == "Gamma" ]; then
     PID=22
   fi;
   for OBJNAME in $OBJNAMES; do

resolutionScaleFiller_submitIlya.sh

@@ -4,7 +4,7 @@
 
 SAMPLEFILE=$1
 
-OUTDIR="/eos/user/c/clange/HGCal/ScaleResolution/Ilya/"
+OUTDIR="/eos/user/c/clange/HGCal/ScaleResolution/_new/Ilya/"
 PREFIX="root://eoscms.cern.ch//eos/cms/store/cmst3/group/hgcal/CMG_studies/Production/"
 GUNTYPE="pt"
 POSTFIX="NTUP"