diff --git a/Tracker/TrackerRecAlgs/PairVertex/scripts/CutNtuple.py b/Tracker/TrackerRecAlgs/PairVertex/scripts/CutNtuple.py
new file mode 100644
index 0000000000000000000000000000000000000000..ea69111172c5a38f2c27608b68d0b3a1d11b2c0a
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/PairVertex/scripts/CutNtuple.py
@@ -0,0 +1,55 @@
+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
diff --git a/Tracker/TrackerRecAlgs/PairVertex/scripts/NtupleReader.py b/Tracker/TrackerRecAlgs/PairVertex/scripts/NtupleReader.py
index ba1696276cf4c9f37fb6b1fb5be9bbb657f7f558..08f0888f95902674676d71761872774141bacdaa 100644
--- a/Tracker/TrackerRecAlgs/PairVertex/scripts/NtupleReader.py
+++ b/Tracker/TrackerRecAlgs/PairVertex/scripts/NtupleReader.py
@@ -1,10 +1,12 @@
#!/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=[])
\ No newline at end of file
+ #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")
+
+
+ GeneralPlotter(TestArrDict,dcuts=[{'longTracks':[None,None,2,None],}],lcuts=[],savefig=True)
\ No newline at end of file
diff --git a/Tracker/TrackerRecAlgs/PairVertex/scripts/PairVertexReader.py b/Tracker/TrackerRecAlgs/PairVertex/scripts/PairVertexReader.py
index f585d0e61a539b8015ca4897b51f5a231c9804ac..7bb21a7e3f961033da2b3b2d2c244d7eed784fbe 100644
--- a/Tracker/TrackerRecAlgs/PairVertex/scripts/PairVertexReader.py
+++ b/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]))
diff --git a/Tracker/TrackerRecAlgs/PairVertex/scripts/PlotHelper.py b/Tracker/TrackerRecAlgs/PairVertex/scripts/PlotHelper.py
index 8597d9cab9779c8ff2fd2e04744730a6f2f3efa7..5ebbdf55462c15b5d2f31ce8c326fdadaa41c298 100644
--- a/Tracker/TrackerRecAlgs/PairVertex/scripts/PlotHelper.py
+++ b/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
diff --git a/Tracker/TrackerRecAlgs/PairVertex/scripts/Utilities.py b/Tracker/TrackerRecAlgs/PairVertex/scripts/Utilities.py
index 43eba0369ec626cf87400a803d4ca6d0655dbdce..03d4cba8a24fb251b0dcb9e478860eff90c4fe1f 100644
--- a/Tracker/TrackerRecAlgs/PairVertex/scripts/Utilities.py
+++ b/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=[]