it works now

Tracker/TrackerRecAlgs/PairVertex/scripts/CutNtuple.py

+import ROOT
+from glob import glob
+
+dataDir="/eos/experiment/faser/phys/2022/r0011/"
+RunDirs = sorted(glob(dataDir+'*'))
+
+outFile = ROOT.TFile.Open('ReNtuple.root',"RECREATE")
+outtree = 0
+
+fcount=0
+
+for runnum in RunDirs:
+    files = sorted(glob(f"{dataDir}{runnum}/Faser*"))
+    for f in files:
+        inFile = ROOT.TFile.Open(f,"READ")
+        intree = inFile.Get('nt')
+        if fcount==0: outtree=intree.CloneTree(0)
+        outtree.CopyEntries(intree)
+        inFile.Close()
+        break
+    break
+
+outFile.Write()
+
+# t = TChain("nt") 
+# nfiles = 0
+# for f in files:
+#     nfiles += t.Add(f) #t.Add("/eos/experiment/faser/phys/2022/r0011/"+user_input+".root") 
+#     print(f)
+
+
+# #----------------------------
+# # Get old file, old tree and set top branch address
+# ofilename="/eos/experiment/faser/phys/2022/r0011/" #....
+# oldfile= ROOT.TFile.Open(ofilename,"READ")
+# oldtree = oldfile.Get("nt")
+
+# nentries = oldtree.GetEntries()
+
+# event = 0
+# oldtree.SetBranchAddress("event", &event)
+
+# # Create a new file + a clone of old tree in new file
+# TFile newfile("small.root", "recreate");
+# auto newtree = oldtree->CloneTree(0);
+
+# for (auto i : ROOT::TSeqI(nentries)) {
+#     oldtree->GetEntry(i);
+#     if (event->GetNtrack() > 605)
+#         newtree->Fill();
+#     event->Clear();
+# }
+
+# #newtree.Print();
+# newfile.Write();
\ No newline at end of file

Tracker/TrackerRecAlgs/PairVertex/scripts/NtupleReader.py

 #!/usr/bin/env python
 from PlotHelper import *
 from Utilities import *
+from glob import glob
 
 #CUTS
+#attr: [>,<,==,!=] if none --> skip
 REALCUT={
-    'longTracks':[2,None,None,None],}
+    'longTracks':[1,None,None,None],}
 #     'Track_nLayers':[6,None,None,None],
 #     'Track_nDoF':[0,None,None,None],
 #     'Track_Chi2':[0,None,None,None],
@@ -12,16 +14,20 @@ REALCUT={
 #     'Preshower1_charge':[0,None,None,None]   
 # }
 
+#LAMBDAS AND CONDITIONALS
+
 REALLAMBDAS=[
-    lambda event: event.Timing0_charge+Timing1_charge+Timing2_charge+Timing3_charge>0,
-    lambda event: event.Calo0_charge+Calo1_charge+Calo2_charge+Calo3_charge>0,
-    lambda event: all([event.Track_nLayers[i]>6 for i in event.longTracks])
-]
+    lambda event: event.Timing0_charge+event.Timing1_charge+event.Timing2_charge+event.Timing3_charge>0,
+    lambda event: event.Calo0_charge+event.Calo1_charge+event.Calo2_charge+event.Calo3_charge>20000]
+#     lambda event: AnyCaloThreshold,
+#     lambda event: all([event.Track_nLayers[i]>6 for i in event.longTracks])
+# ]
 
 def GeneralPlotter(arrDict,dcuts=[],lcuts=[],showplot=True,savefig=False):
     nans=0
-
+    ecount=0
     for event in t:
+        
         if ApplyAllCuts(event,dcuts,lcuts):
             for name in arrDict.keys():
                 value=nan
@@ -39,6 +45,8 @@ def GeneralPlotter(arrDict,dcuts=[],lcuts=[],showplot=True,savefig=False):
                             arrDict[name]['array'][iv]+=[value[iv]]
                 elif checkReal: arrDict[name]['array']+=[value]
                 else: nans+=1
+        if ecount%100000==0:print(ecount,' events processed')
+        ecount+=1
 
     print(f"# of NaNs:{nans}")
     for name in arrDict.keys():
@@ -65,7 +73,7 @@ MainArrDict={
     'CutTest':{
         'array':[],
         'lambda': ApplyAllCuts,
-        'lkwargs':{'dcuts':[REALCUT],'lcuts':[]}, 
+        'lkwargs':{'dcuts':[REALCUT],'lcuts':REALLAMBDAS}, 
         'plotterfunc': PrintBoolPercent,
         'plotkwargs' : {},
         'descr':'Check how much a cut cuts.'
@@ -76,22 +84,44 @@ MainArrDict={
         'plotterfunc': configBasicHist,
         'plotkwargs' : {},
         'descr':'Single variable tests.'
+    },
+    'HighestMomentumUp':{
+        'array':[[],[]],
+        'lambda': lambda event: sorted([p for p in event.Track_p0]),
+        'plotterfunc': configMultiHist,
+        'plotkwargs' : {'label':['Lowest','Highest'],'pckwargs':{'xscale':'symlog','xlabel':"Momentum (MeV)",'legend':True}},
+        'descr':'Look at 2 highest momentum tracks upstream'
+    },
+    'HighestMomentumDown':{
+        'array':[[],[]],
+        'lambda': lambda event: sorted([p for p in event.Track_p1]),
+        'plotterfunc': configMultiHist,
+        'plotkwargs' : {'label':['Lowest','Highest'],'pckwargs':{'xscale':'symlog','xlabel':"Momentum (MeV)",'legend':True}},
+        'descr':'Look at 2 highest momentum tracks downstream'
     },}
 
 TestArrDict=dict()
-testname='BasicHist' #Any key from MainArrDict
+testname='HighestMomentumDown' #Any key from MainArrDict
 TestArrDict[testname]=MainArrDict[testname]
 
+def getRunFiles(rnum,dataDir="/eos/experiment/faser/phys/2022/r0011/"):
+    dataDir=f"{dataDir}{rnum}/"
+    return sorted(glob(f"{dataDir}/Faser*"))
+
+def getRunNums(dataDir="/eos/experiment/faser/phys/2022/r0011/"):
+    return [rn.split(dataDir)[1] for rn in sorted(glob(f"{dataDir}*"))]
 
 
 #EXECUTE
 if __name__ == "__main__":
-    DP_fnames=['008301/Faser-Physics-008301-00000-00049-r0011-PHYS']
-    for f in DP_fnames:
-        user_input=f
-        t = TChain("nt") 
-        nfiles = 0
-        nfiles += t.Add("/eos/experiment/faser/phys/2022/r0011/"+user_input+".root") 
-        print(f)
-        
-        GeneralPlotter(TestArrDict,dcuts=[{}],lcuts=[])
+    #runnum = '008301'
+    t = TChain("nt") 
+    nfiles = 0
+    for runnum in getRunNums():
+        files = getRunFiles(runnum)
+        for f in files:
+            print(f)
+            nfiles += t.Add(f) #t.Add("/eos/experiment/faser/phys/2022/r0011/"+user_input+".root") 
+            
\ No newline at end of file
+    
+    GeneralPlotter(TestArrDict,dcuts=[{'longTracks':[None,None,2,None],}],lcuts=[],savefig=True)
\ No newline at end of file

Tracker/TrackerRecAlgs/PairVertex/scripts/PairVertexReader.py

@@ -103,6 +103,82 @@ def GetAllTrackMomentum(event):
 
     return finalval
 
+def CountSharedCluster(event,filter=None):
+    #How many clusters have energy contributions from both daughter particles?
+    sharedcount=0
+    for c in range(len(event.clusterEFracPos)):
+        if event.clusterETotPos[c]>0 and event.clusterETotNeg[c]>0 and event.clusterEFracPos[c]<0.99 and event.clusterEFracNeg[c]<0.99: sharedcount+=1
+    if (filter and filter(event)) or not filter:
+        return sharedcount
+    else: return nan
+
+def CountSharedClusterCompare(event):
+    #How many clusters have energy contributions from both daughter particles?
+    sharedcount=0 
+    sharedcountarr=[0,0,0] #[No filter, pass, fail]
+    for c in range(len(event.clusterEFracPos)):
+        if event.clusterETotPos[c]>0 and event.clusterETotNeg[c]>0 and event.clusterEFracPos[c]<0.99 and event.clusterEFracNeg[c]<0.99: 
+            sharedcount+=1
+    
+    filters = [lambda e:True,BoolDaughterCut2,BoolDaughterCut3]
+    for i in range(len(filters)):
+        if filters[i](event): sharedcountarr[i]=sharedcount
+        else: sharedcountarr[i]=nan
+
+    return sharedcountarr
+
+def VertexZCompare(event):
+    val = event.vertexPos[2]
+    countarr=[0,0,0] #[No filter, HasShared, NoShared]    
+    filters = [lambda e:True,BoolHasSharedCluster,BoolNotHasSharedCluster]
+    for i in range(len(filters)):
+        if filters[i](event): countarr[i]=val
+        else: countarr[i]=nan
+
+    return countarr
+
+def ClusterLocation(event):
+    starr=[]
+    for loc in event.clusterLocation:
+        starr+=[int(str(loc).split()[0])]
+    return starr
+
+def MakeCompareArray(event,getvaluelambda,filters):
+    val = getvaluelambda(event)
+    countarr=[0 for i in range(len(filters))]
+    for fi in range(len(filters)):
+        if filters[i](event): countarr[i]=val
+        else: countarr[i]=nan
+
+    return countarr
+
+def BoolDaughterCut(event):
+    #both daughter particles are matched to tracks
+    if event.TrackCount>=2 and len(event.tracksTruthBarcode)>=2:
+        return 2 in event.tracksTruthBarcode and 3 in event.tracksTruthBarcode
+    else: return False
+
+def BoolDaughterCut2(event):
+    #passes track cuts and is fiducial
+    return event.TrackCount>=2 and  2 in event.tracksTruthBarcode and 3 in event.tracksTruthBarcode and event.truthd0_IsFiducial and event.truthd1_IsFiducial
+
+def BoolDaughterCut3(event):
+    #fails track cuts and is not fiducial
+    return any([event.TrackCount<2, 2 not in event.tracksTruthBarcode, 3 not in event.tracksTruthBarcode]) and event.truthd0_IsFiducial and event.truthd1_IsFiducial
+
+def BoolHasSharedCluster(event):
+    #Does the event have shared clusters?
+    sharedcount=0
+    for c in range(len(event.clusterEFracPos)):
+        if event.clusterETotPos[c]>0 and event.clusterETotNeg[c]>0 and event.clusterEFracPos[c]<0.99 and event.clusterEFracNeg[c]<0.99: sharedcount+=1
+    return sharedcount>0
+
+def BoolNotHasSharedCluster(event):
+    #Does the event have shared clusters?
+    sharedcount=0
+    for c in range(len(event.clusterEFracPos)):
+        if event.clusterETotPos[c]>0 and event.clusterETotNeg[c]>0 and event.clusterEFracPos[c]<0.99 and event.clusterEFracNeg[c]<0.99: sharedcount+=1
+    return sharedcount==0
 
 #For storing different plot configurations
 MainArrDict={
@@ -191,12 +267,68 @@ MainArrDict={
         'plotterfunc': configBasicHist,
         'plotkwargs' : {'bins':20,'kwargs':{'range':[-10,11]},'pckwargs':{'xscale':'linear','xlabel':"Residual"}} 
     },
+    'CountSharedCluster':{
+        'array':[],
+        'lambda': CountSharedCluster,
+        'plotterfunc': configBasicHist,
+        'plotkwargs' : {'bins':25,'kwargs':{'range':[0,26]}},
+        'descr':'How many clusters have energy deposits from both daughter particles and fail track cuts?'
+    },
+    'CountSharedCluster_PassTrackCuts':{
+        'array':[],
+        'lambda': CountSharedCluster,
+        'lkwargs':{'filter':BoolDaughterCut2},
+        'plotterfunc': configBasicHist,
+        'plotkwargs' : {'bins':25,'kwargs':{'range':[0,26]}},
+        'descr':'How many clusters have energy deposits from both daughter particles and pass track cuts?'
+    },
+    'CountSharedCluster_FailTrackCuts':{
+        'array':[],
+        'lambda': CountSharedCluster,
+        'lkwargs':{'filter':BoolDaughterCut3},
+        'plotterfunc': configBasicHist,
+        'plotkwargs' : {'bins':25,'kwargs':{'range':[0,26]}},
+        'descr':'How many clusters have energy deposits from both daughter particles and fail track cuts?'
+    },
+    'CountSharedCluster_Compare':{
+        'array':[[],[],[]],
+        'lambda': CountSharedClusterCompare,
+        'plotterfunc': configMultiHist,
+        'plotkwargs' : {'label':['No Cuts','Passes Track Cuts','Fails Track Cuts'],'bins':25,'kwargs':{'range':[0,26]},'pckwargs':{'xscale':'linear','legend':True}},
+        'descr':'How many clusters have energy deposits from both daughter particles?',
+        'SepHist':'arrays allowed to have different lengths'
+    },
+    'DPDecayVertexZ':{
+        'array':[],
+        'lambda': lambda event: event.vertexPos[2],
+        'plotterfunc': configBasicHist,
+        'plotkwargs' : {},
+        'descr':'Dark photon decay vertex coordinate along beamline'
+    },
+    'DPDecayVertexZ_Compare':{
+        'array':[[],[],[]],
+        'lambda': VertexZCompare,
+        'plotterfunc': configMultiHist,
+        'plotkwargs' : {'label':['No Cuts','Has Shared Clusters','No Shared Clusters'],'bins':25,'kwargs':{'range':[-1600,0]},'pckwargs':{'xscale':'linear','legend':True,'legloc':'center right'}},
+        'descr':'Compare dp-decay z-coordinate cuts',
+        'SepHist':'arrays allowed to have different lengths'
+    },
+    'AllClusterLocations':{
+        'array':[],
+        'lambda': ClusterLocation,
+        'plotterfunc': configBasicHist,
+        'plotkwargs' :{'bins':5,'kwargs':{'range':[0,5]}},
+    },
     }
+
+    #energy deposited hist
+    #frac of energy deposited hist
+    #
     
 
 
 TestArrDict=dict()
-testname='NormalizedResidualSum' #Any key from MainArrDict
+testname='DPDecayVertexZ_Compare' #Any key from MainArrDict
 TestArrDict[testname]=MainArrDict[testname]
 # for k in MainArrDict.keys():
 #     if 'Check' in k:
@@ -215,6 +347,8 @@ def GeneralPlotter(arrDict,dcuts=[],lcuts=[],showplot=True,savefig=False):
                 checkReal=True
                 if type(value)==list: checkReal = all([v != inf and v !=nan and v !=-1*inf for v in value])
                 else: checkReal = value != inf and value !=nan and value !=-1*inf
+                #if 'SepHist' in list(arrDict[name].keys()): checkReal=False
+                
 
                 if checkReal and type(value)==list:
                     if arrDict[name]['plotterfunc'] in [configBasicHist]: arrDict[name]['array']+=value
@@ -222,6 +356,16 @@ def GeneralPlotter(arrDict,dcuts=[],lcuts=[],showplot=True,savefig=False):
                         for iv in range(len(value)):
                             arrDict[name]['array'][iv]+=[value[iv]]
                 elif checkReal: arrDict[name]['array']+=[value]
+                elif 'SepHist' in list(arrDict[name].keys()):
+                    # if type(value[0])==list:
+                    #     for vi in range(len(value)):
+                    #         checkList=[v != inf and v !=nan and v !=-1*inf for v in value[vi]]
+                    #         for ivch in len(checkList):
+                    #             if checkList[ivch]: arrDict[name]['array'][vi]+=[value[vi][ivch]]
+                    checkList=[v != inf and v !=nan and v !=-1*inf for v in value]
+                    for ich in len(checkList):
+                        if checkList[ich]: arrDict[name]['array'][ich]+=[value[ich]]
+
                 else: nans+=1
 
     print(f"# of NaNs:{nans}")
@@ -237,13 +381,13 @@ def GeneralPlotter(arrDict,dcuts=[],lcuts=[],showplot=True,savefig=False):
 
 #EXECUTE
 if __name__ == "__main__":
-    DP_fnames=['Ntuple_Aee_100MeV_1Em5_']#  Ntuple_Aee_10MeV_1Em5_ Ntuple_Amm_316MeV_2Em6_ Ntuple_Aee_10MeV_1Em4_
+    DP_fnames=['Ntuple_Amm_316MeV_2Em6_']#  Ntuple_Aee_10MeV_1Em5_ Ntuple_Amm_316MeV_2Em6_ Ntuple_Aee_10MeV_1Em4_ Ntuple_Aee_100MeV_1Em5_
     for f in DP_fnames:
         user_input=f
         t = TChain("events") 
         nfiles = 0
-        nfiles += t.Add("/eos/home-s/sshively/CALYPZONE/run/BackUp2/"+user_input+".root") 
+        nfiles += t.Add("/eos/home-s/sshively/CALYPZONE/run/"+user_input+".root") 
         print(f)
         
-        GeneralPlotter(TestArrDict,dcuts=[{}],lcuts=[BoolDaughterCut])
+        GeneralPlotter(TestArrDict,dcuts=[{}],lcuts=[])
         #print(GetPBounds(lcuts=[BoolDaughterCut]))

Tracker/TrackerRecAlgs/PairVertex/scripts/PlotHelper.py

@@ -3,7 +3,7 @@ from cmcrameri import cm
 import matplotlib.pyplot as plt
 from ROOT import TEfficiency, TH1F
 
-def PlotConfigs(xscale='linear',yscale='linear',xlim=None,ylim=None,title=None,show=False,xticks=None,xlabel=None,ylabel='#Events',legend=False,legloc='best'):
+def PlotConfigs(xscale='linear',yscale='linear',xlim=None,ylim=None,title=None,show=False,xticks=None,xlabel=None,ylabel='#Events',legend=False,legloc='best',save=False):
     plt.xscale(xscale)
     plt.yscale(yscale)
     if xticks: plt.xticks(xticks)
@@ -12,6 +12,9 @@ def PlotConfigs(xscale='linear',yscale='linear',xlim=None,ylim=None,title=None,s
     if legend: plt.legend(loc=legloc)
     if title:plt.title(title)
     if show: plt.show()
+    #filename = title.replace(" ", "_").replace("(", "").replace(")", "")
+    #if save: plt.savefig(filename)
+
 
 def configHitMap(arrays,bins=(50,50),pckwargs={'xlabel':"X (mm)",'ylabel':"Y (mm)"},counts=False,kwargs={}):
     #2D Histogram with colors

Tracker/TrackerRecAlgs/PairVertex/scripts/Utilities.py

@@ -22,15 +22,6 @@ BASICCUT={
     }
 
 
-
-
-def BoolDaughterCut(event):
-    if event.TrackCount>=2 and len(event.tracksTruthBarcode)>=2:
-        return 2 in event.tracksTruthBarcode and 3 in event.tracksTruthBarcode
-    else: return False
-
-
-
 def ApplyCut(event, cut):
     tests=[gt,lt,eq,ne]
     bool1=[]