diff --git a/PhysicsAnalysis/NtupleDumper/scripts/analyzeRun.py b/PhysicsAnalysis/NtupleDumper/scripts/analyzeRun.py
index fcee52fbc227a9746d5b45d285295ccfeae95736..8a47857614ebb70589ed7cb5ab6ca20aef91c878 100755
--- a/PhysicsAnalysis/NtupleDumper/scripts/analyzeRun.py
+++ b/PhysicsAnalysis/NtupleDumper/scripts/analyzeRun.py
@@ -56,8 +56,10 @@ if __name__ == "__main__":
Configurable.configurableRun3Behavior = True
dataDir=f"/eos/experiment/faser/rec/2022/{ptag}/{runno:06d}"
- files=sorted(glob.glob(f"{dataDir}/Faser-Physics*"))
- fileListInitial=files[num*filesPerJob:(num+1)*filesPerJob]
+ if inputIsMC: dataDir=f"/eos/experiment/faser/sim/mdc/foresee/{ui.runnum}/rec/r0009/"
+ files=sorted(glob.glob(f"{dataDir}/Faser*"))
+ if num==-1: num=len(files)
+ fileListInitial=files[0:num]#[num*filesPerJob:(num+1)*filesPerJob]
fileList=[]
for fName in fileListInitial:
try:
@@ -66,7 +68,7 @@ if __name__ == "__main__":
except OSError:
print("Warning bad file: ",fName)
- log.info(f"Analyzing Run {runno} files {num*filesPerJob} to {(num+1)*filesPerJob} (num={num})")
+ log.info(f"Analyzing Run {runno} files {0} to {num} (num={num})")
log.info(f"Got {len(fileList)} files out of {len(fileListInitial)}")
outName=f"Data-tuple-{runno:06d}-{num:05d}-{filesPerJob}.root"
@@ -76,7 +78,7 @@ if __name__ == "__main__":
ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-03" # Always needed; must match FaserVersionS
ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions for now
ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig
- ConfigFlags.Input.isMC = False # Needed to bypass autoconfig
+ ConfigFlags.Input.isMC = inputIsMC # Needed to bypass autoconfig
ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # FASER geometry
ConfigFlags.Common.isOnline = False
ConfigFlags.GeoModel.Align.Dynamic = False
@@ -98,11 +100,12 @@ if __name__ == "__main__":
acc.addService(AthenaEventLoopMgr)
# # Hack to avoid problem with our use of MC databases when isMC = False
- replicaSvc = acc.getService("DBReplicaSvc")
- replicaSvc.COOLSQLiteVetoPattern = ""
- replicaSvc.UseCOOLSQLite = True
- replicaSvc.UseCOOLFrontier = False
- replicaSvc.UseGeomSQLite = True
+ if not inputIsMC:
+ replicaSvc = acc.getService("DBReplicaSvc")
+ replicaSvc.COOLSQLiteVetoPattern = ""
+ replicaSvc.UseCOOLSQLite = True
+ replicaSvc.UseCOOLFrontier = False
+ replicaSvc.UseGeomSQLite = True
# Timing
#acc.merge(MergeRecoTimingObjCfg(ConfigFlags))
@@ -119,7 +122,7 @@ if __name__ == "__main__":
# ConfigFlags.dump()
# Execute and finish
- sc = acc.run(maxEvents=-1)
+ sc = acc.run(maxEvents=ui.events)
# Success should be 0
sys.exit(not sc.isSuccess())
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
index 9614de99c394407659689aea08630ed141cdbbb9..1369684d16f1071c2b5eab6e5853b87459d0e06a 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
@@ -136,6 +136,8 @@ StatusCode NtupleDumperAlg::initialize()
}
m_tree = new TTree("nt", "NtupleDumper tree");
+
+ //META INFORMATION
m_tree->Branch("run", &m_run_number, "run/I");
m_tree->Branch("eventID", &m_event_number, "eventID/I");
m_tree->Branch("eventTime", &m_event_time, "eventTime/I");
@@ -156,6 +158,7 @@ StatusCode NtupleDumperAlg::initialize()
m_tree->Branch("inputBits", &m_inputBits, "inputBits/I");
m_tree->Branch("inputBitsNext", &m_inputBitsNext, "inputBitsNext/I");
+ //WAVEFORMS
addWaveBranches("VetoNu",2,4);
addWaveBranches("VetoSt1",2,6);
addWaveBranches("VetoSt2",1,14);
@@ -172,6 +175,7 @@ StatusCode NtupleDumperAlg::initialize()
addBranch("Preshower_total_nMIP", &m_preshower_total_nMIP);
addBranch("Preshower_total_E_dep", &m_preshower_total_E_dep);
+ //TRACKER
addBranch("nClusters0",&m_station0Clusters);
addBranch("nClusters1",&m_station1Clusters);
addBranch("nClusters2",&m_station2Clusters);
@@ -192,6 +196,7 @@ StatusCode NtupleDumperAlg::initialize()
m_tree->Branch("TrackSegment_py", &m_trackseg_py);
m_tree->Branch("TrackSegment_pz", &m_trackseg_pz);
+ //TrackCollection
m_tree->Branch("longTracks", &m_longTracks, "longTracks/I");
m_tree->Branch("Track_Chi2", &m_Chi2);
m_tree->Branch("Track_nDoF", &m_DoF);
@@ -212,11 +217,13 @@ StatusCode NtupleDumperAlg::initialize()
m_tree->Branch("Track_charge", &m_charge);
m_tree->Branch("Track_nLayers", &m_nLayers);
+ //Which Stations were hit?
m_tree->Branch("Track_InStation0",&m_nHit0);
m_tree->Branch("Track_InStation1",&m_nHit1);
m_tree->Branch("Track_InStation2",&m_nHit2);
m_tree->Branch("Track_InStation3",&m_nHit3);
+ //Extrapolated Track Info
m_tree->Branch("Track_X_atVetoNu", &m_xVetoNu);
m_tree->Branch("Track_Y_atVetoNu", &m_yVetoNu);
m_tree->Branch("Track_ThetaX_atVetoNu", &m_thetaxVetoNu);
@@ -252,15 +259,20 @@ StatusCode NtupleDumperAlg::initialize()
m_tree->Branch("Track_ThetaX_atCalo", &m_thetaxCalo);
m_tree->Branch("Track_ThetaY_atCalo", &m_thetayCalo);
+ //TRUTH
m_tree->Branch("t_pdg", &m_t_pdg);
m_tree->Branch("t_barcode", &m_t_barcode);
m_tree->Branch("t_truthHitRatio", &m_t_truthHitRatio);
+
m_tree->Branch("t_prodVtx_x", &m_t_prodVtx_x);
m_tree->Branch("t_prodVtx_y", &m_t_prodVtx_y);
m_tree->Branch("t_prodVtx_z", &m_t_prodVtx_z);
+
m_tree->Branch("t_decayVtx_x", &m_t_decayVtx_x);
m_tree->Branch("t_decayVtx_y", &m_t_decayVtx_y);
m_tree->Branch("t_decayVtx_z", &m_t_decayVtx_z);
+
+ //Truth track info
m_tree->Branch("t_px", &m_t_px);
m_tree->Branch("t_py", &m_t_py);
m_tree->Branch("t_pz", &m_t_pz);
diff --git a/Tracker/TrackerRecAlgs/PairVertex/CMakeLists.txt b/Tracker/TrackerRecAlgs/PairVertex/CMakeLists.txt
index ff270bb69c9defc70b31e5e50f8dbc93ffa3992f..07e4fc4065838c8e54104fdbb25625c4b0bb2d21 100644
--- a/Tracker/TrackerRecAlgs/PairVertex/CMakeLists.txt
+++ b/Tracker/TrackerRecAlgs/PairVertex/CMakeLists.txt
@@ -9,6 +9,7 @@ atlas_add_component(
FaserActsGeometryLib TrackerSimEvent TrackerIdentifier
TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack
xAODFaserWaveform ScintIdentifier FaserCaloIdentifier
+ FaserActsKalmanFilterLib
)
atlas_install_python_modules(python/*.py)
diff --git a/Tracker/TrackerRecAlgs/PairVertex/scripts/EfficiencyPlotter.py b/Tracker/TrackerRecAlgs/PairVertex/scripts/EfficiencyPlotter.py
new file mode 100644
index 0000000000000000000000000000000000000000..8fc57cc162617b2d5c5b5b1bd6035329d95ea096
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/PairVertex/scripts/EfficiencyPlotter.py
@@ -0,0 +1,53 @@
+from PlotHelper import *
+from Utilities import *
+
+#Efficiency Plotters Seperate from General Plotter system. Maybe put this in its own script
+def CutEfficiencyArrays(lmbda,dcutN=[],dcutD=[],lcutN=[],lcutD=[],title=None,bins=15,hmin=0,hmax=10000):
+ Func = lmbda # (Returns 1 value)
+ narray=[]
+ darray=[]
+ for event in t:
+ if ApplyAllCuts(event,dcutN,lcutN): narray+=[Func(event)]
+ if ApplyAllCuts(event,dcutD,lcutD): darray+=[Func(event)]
+
+ PlotRootRatio([narray,darray],title=f"{user_input} {title}",bins=bins,hmin=hmin,hmax=hmax)
+
+def FiducialEff_vs_Rho2():
+ l= lambda event: (event.vertexPos[0]**2)+(event.vertexPos[1]**2)
+ CutEfficiencyArrays(l,title='BothDFiducial vs TRho^2',bins=20,hmin=0,hmax=10000,lcutN=[lambda event: event.truthd0_IsFiducial[0] and event.truthd1_IsFiducial[0]])
+
+def FiducialEff_vs_VertexZ():
+ l= lambda event: event.vertexPos[2]
+ CutEfficiencyArrays(l,title='BothDFiducial vs TVertexZ',bins=20,hmin=-1600,hmax=0,lcutN=[lambda event: event.truthd0_IsFiducial[0] and event.truthd1_IsFiducial[0]])
+
+EfficiencyDict={
+ 'FiducialEff_vs_Rho2':{
+ 'lambda': lambda event: (event.vertexPos[0]**2)+(event.vertexPos[1]**2),
+ 'kwargs':{'title':'BothDFiducial vs TRho^2',
+ 'bins':20,'hmin':0,'hmax':10000,
+ 'lcutN':[lambda event: event.truthd0_IsFiducial[0] and event.truthd1_IsFiducial[0]]}
+ },
+ 'FiducialEff_vs_Rho2':{
+ 'lambda': lambda event: event.vertexPos[2],
+ 'kwargs':{'title':'BothDFiducial vs TVertexZ',
+ 'bins':20,'hmin':-1600,'hmax':0,
+ 'lcutN':[lambda event: event.truthd0_IsFiducial[0] and event.truthd1_IsFiducial[0]]}
+ }
+}
+
+def DoEfficiencyPlots():
+ for k in EfficiencyDict.keys():
+ print(k)
+ CutEfficiencyArrays(EfficiencyDict[k]['lambda'],**EfficiencyDict[k]['kwargs'])
+
+#EXECUTE
+if __name__ == "__main__":
+ DP_fnames=['Ntuple_Aee_100MeV_1Em5_']# Ntuple_Aee_10MeV_1Em5_ Ntuple_Amm_316MeV_2Em6_ Ntuple_Aee_10MeV_1Em4_
+ 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")
+ print(f)
+
+ DoEfficiencyPlots()
\ No newline at end of file
diff --git a/Tracker/TrackerRecAlgs/PairVertex/scripts/PairVertexReader.py b/Tracker/TrackerRecAlgs/PairVertex/scripts/PairVertexReader.py
index 2468be886114912fa36555c88abaad99e8e77076..f585d0e61a539b8015ca4897b51f5a231c9804ac 100644
--- a/Tracker/TrackerRecAlgs/PairVertex/scripts/PairVertexReader.py
+++ b/Tracker/TrackerRecAlgs/PairVertex/scripts/PairVertexReader.py
@@ -1,587 +1,249 @@
#!/usr/bin/env python
+from PlotHelper import *
+from Utilities import *
-# Set up (Py)ROOT.
-from ROOT import TChain, TH2F, TH1F
-print("Root loaded.")
-import matplotlib.pyplot as plt
-print("MPL loaded.")
-from operator import lt,eq,ne,gt
-from numpy import nan, seterr, inf,arctan,log10,logspace, histogram2d, meshgrid, ones_like,array
-from numpy import sum as npsum
-import numpy as np
-from math import sqrt
-from cmcrameri import cm
-print("Other stuff loaded.")
-seterr(all='ignore')
-missingattr=[]
-
-#attr: [>,<,==,!=] if none --> skip
-CUTDICT={
- 'TrackCount':[None,None,2,None],
- 'VetoNu0_integral':[None,None,0,None],
- 'VetoNu1_integral':[None,None,0,None],
- 'VetoSt10_integral':[None,None,0,None],
- 'VetoSt11_integral':[None,None,0,None],
- 'VetoSt20_integral':[None,None,0,None],
- 'nlayersHit':[6,None,None,None],
- 'playersHit':[6,None,None,None],
- 'niftHit':[None,None,0,None],
- 'piftHit':[None,None,0,None],
- 'pDoF':[0,None,None,None],
- 'nDoF':[0,None,None,None],
- 'Preshower0_integral':[0,None,None,None],
- 'Preshower1_integral':[0,None,None,None]
- }
-
-#CUT DICTIONARIES/LISTS
-LAMBDACUTS=[
- lambda event : len(event.Tracks[0])==2,
- lambda event : event.Tracks[0][0]==-1*event.Tracks[0][1],
- lambda event : event.Timing0_integral+event.Timing1_integral+event.Timing2_integral+event.Timing3_integral>0,
- lambda event : event.Calo0_integral+event.Calo1_integral+event.Calo2_integral+event.Calo3_integral>0,
- lambda event : event.nChi2/event.nDoF <10 and event.pChi2/event.pDoF <10,
- lambda event : sqrt((event.pMom[0]**2)+(event.pMom[1]**2)+(event.pMom[2]**2))>1]
-
-FIDUCIALCUT=[lambda event : event.vertexPos[2]>-1585 and event.vertexPos[2]<-45]
-
-BADCUT=list(LAMBDACUTS)
-#BADCUT[4]=lambda event : event.nChi2/event.nDoF >10 and event.pChi2/event.pDoF >10
-#BADCUT[]
-
-HITCUT= [lambda event : all([abs(i)<200 for i in [event.pPos_S1L0[0],event.pPos_S1L0[1],event.nPos_S1L0[0],event.nPos_S1L0[1]]])]
-
-
-#UTILITY
-def configHitMap(arrayx,arrayy=[],xlabel="X (mm)",ylabel="Y (mm)",xlim=None,ylim=None,title=None,show=False):
- plt.hist2d(arrayx, arrayy, bins=(50, 50), cmap=cm.imola)
- plt.colorbar()
- plt.xlabel(xlabel)
- plt.ylabel(ylabel)
- if xlim: plt.xlim(xlim)
- if ylim: plt.ylim(ylim)#-100,100
- #plt.text(50,-185,"COUNTS: {}".format(len(arrayx)))
- if title: plt.title(title)
- if show: plt.show()
-
-def configHitMap2(arrays,bins=(50,50),xlabel="X (mm)",ylabel="Y (mm)",xlim=None,ylim=None,title=None,show=False):
- arrayx=arrays[0]
- arrayy=arrays[1]
- plt.hist2d(arrayx, arrayy, bins=bins, cmap=cm.imola)
- plt.colorbar()
- plt.xlabel(xlabel)
- plt.ylabel(ylabel)
- if xlim: plt.xlim(xlim)
- if ylim: plt.ylim(ylim)#-100,100
- #plt.text(50,-185,"COUNTS: {}".format(len(arrayx)))
- if title: plt.title(title)
- if show: plt.show()
-
-def configBasicHist(arr,bins=25,xscale='linear',yscale='linear',title=None,xlabel=None,legend=False,kwargs={}):
- if len(arr)==0:
- print("Empty array.")
- return
- if 'log' in xscale:
- minx,maxx = (min(arr),max(arr))
- if minx==maxx:
- maxx+=1
- bins=1
- bins = logspace(start=log10(minx), stop=log10(maxx), num=bins)
- if xscale=='symlog':
- bins = [0]+list(logspace(start=log10(1), stop=log10(maxx), num=24))
- if minx>1: bins = logspace(start=log10(minx), stop=log10(maxx), num=25)
-
- n, bins, patches = plt.hist(arr,bins=bins,align="left",histtype=u'step',**kwargs)
- plt.xscale(xscale)
- #plt.xticks([i for i in range(0,21)])
- if xlabel: plt.xlabel(xlabel)
- plt.yscale(yscale)
- plt.ylabel('Number of Events')
- if legend: plt.legend()
- if title:plt.title(title)
+#Q/P METHODS
+def GetPBounds(dcuts=[],lcuts=[]):
+ TotalTracks=0
+ PArray=[]
+ for event in t:
+ if ApplyAllCuts(event,dcuts,lcuts):
+ TotalTracks+=event.TrackCount
+ PArray+=list(event.tracksTruthPmag)
+ Ilims=[int(i*TotalTracks/10) for i in range(11)]
+ PArray.sort()
+ Plims=[0]+[PArray[i] for i in Ilims]
+ return Plims
+
+#Q/P P bounds entered manually because storing them was weird in dictionary generation
+PH = [191909.7203710638, 248751.19099289246]
+PHr=0.02
+def MakeNewQPDicts(Plims=None):
+ newdict=dict()
+ basename='AllTracksQP'
+ dictbase={
+ 'array':[],
+ 'lambda': GetAllTrackQPReco,
+ 'plotterfunc': configBasicHist,
+ 'plotkwargs' :{'bins':100,'kwargs':{'range':[-1*PHr,PHr]}}} #[2*-10**-3,2*10**-3]
- return n,bins
-
-def configMultiHist(arrlist,bins=25,xscale='linear',yscale='linear',title=None,xlabel=None,label=[],legend=True,kwargs={}):
- megarr=[]
- for a in arrlist: megarr+=a
- if 'log' in xscale:
- minx,maxx = (min(megarr),max(megarr))
- if minx==maxx:
- maxx+=1
- bins=1
- bins = logspace(start=log10(minx), stop=log10(maxx), num=bins)
- if xscale=='symlog':
- bins = [0]+list(logspace(start=log10(1), stop=log10(maxx), num=24))
- if minx>1: bins = logspace(start=log10(minx), stop=log10(maxx), num=25)
-
- for iarr in range(len(arrlist)):
- arr=arrlist[iarr]
- labelused=''
- if len(megarr)==0:
- print(f"Empty array.")
- else:
-
- if len(label)==len(arrlist): labelused=label[iarr]
- else:
- print("Insufficient labels given.")
- plt.hist(arr,bins=bins,align="left",histtype=u'step',label=labelused,**kwargs)
-
- plt.xscale(xscale)
- plt.yscale(yscale)
- if xlabel: plt.xlabel(xlabel)
- plt.ylabel('Number of Events')
- if legend: plt.legend()
- if title:plt.title(title)
-
-def sortList(tlist):
- tlist.sort(key = lambda x: x[4], reverse=True)
- return tlist
-
-def plotLEGO(xarr,yarr):
- hist = TH2F("hist","Ndiff",20,-10,10,20,-10,10)
- print("Fill start.")
- for i in range(len(xarr)):
- hist.Fill(xarr[i],yarr[i])
-
- print("Fill done.")
- hist.Draw("LEGO") #colz
- input()
-
-def getQPNSorted(event,highestrank=1,lowestrank=3): #Get sorted array containing [q,px,py,pz,p-mag], ranked by p-mag
- trackQPN=[] # [(q,px,py,pz,p)]
- for i in range(highestrank-1,lowestrank):
- #Momenum is in MeV by default, so we convert to GeV since 1 MeV = 10^-3 GeV
- pmag = sqrt((event.Tracks[4][i]**2)+(event.Tracks[5][i]**2)+(event.Tracks[6][i]**2))
- trackQPN+=[(event.Tracks[0][i],event.Tracks[4][i]*(10**-3),event.Tracks[5][i]*(10**-3),event.Tracks[6][i]*(10**-3),pmag*(10**-3))]
-
- return sortList(trackQPN)
+ #for i in range(len(Plims)-1):
+ name=f'{basename}_{int(PH[0]*10**-3)}-{int(PH[1]*10**-3)}GeV'
+ newdict[name]=dict(dictbase)
+ newdict[name]['lkwargs']={'tar':[PH[0],PH[1]]}
+ # break #TEMP
-def ApplyCut(event, cut):
- tests=[gt,lt,eq,ne]
- bool1=[]
- global missingattr
- for attr in cut.keys():
- for i in range(len(tests)):
- try:
- if cut[attr][i]!=None:
- value = getattr(event,attr)
- bool1 += [tests[i](value,cut[attr][i]) and value != inf and value !=nan and value !=-1*inf]
- if not all(bool1): break
- except AttributeError:
- if attr not in missingattr:
- print("No attribute: ",attr)
- missingattr+=[attr]
-
- if not all(bool1):
- break
-
- return all(bool1)
+ return newdict
-def ApplyAllCuts(event,dcuts=[CUTDICT],lcuts=LAMBDACUTS): #Requires all True to be True
- #cutList = list of dictionaries
- bool1=True
- boollist=[]
- for cut in dcuts:
- bool1 = ApplyCut(event,cut)
- if not bool1: break
- if bool1:
- for lcut in lcuts:
- boollist+=[lcut(event)]
- if not all(boollist): break
+#DATA LAMBDAS (for retrieving values in complicated ways)
+def RecoVsTruthMomentumMatched(event):
+ return [event.tracksPmag[event.trackD2i-1],event.tracksTruthPmag[event.trackD2i-1]]
- return all(boollist) and bool1
-
-def RootRatio(narray,darray):
- #rarray = np.divide(narray,darray)
- h_numerator=TH1F("num","num",15,0,10000)
- h_denominator=TH1F("den","den",15,0,10000)
- for i in narray: h_numerator.Fill(i)
- for j in darray: h_denominator.Fill(j)
-
- h_ratio = h_numerator.Clone()
- h_ratio.Sumw2()
- h_ratio.Divide(h_denominator)
- h_ratio.SetStats(0)
- h_ratio.SetTitle(f"{user_input} Efficiency vs Rho^2;Rho^2 (mm^2);Efficiency(%)")
- h_ratio.SetMinimum(0.)
- h_ratio.Draw("ep")
- input()
-
-def MakeWFDict():
- WFDict=dict()
- waveforms={'VetoSt1':(0,1), 'VetoSt2':(0,0),'Timing':(0,3),'Preshower':(0,1),'Calo':(0,3)}
- meas = 'integral'
- TODO = list(waveforms.keys())
- for wf in TODO:
- st_r=waveforms[wf]
- for st in range(st_r[0],st_r[1]+1):
- key=f'{wf}{st}_{meas}'
- WFDict[key]={
- 'array':[],
- 'lkwargs':{'k':key},
- 'lambda': lambda event,k: getattr(event, k),
- 'plotterfunc': configBasicHist,
- 'plotkwargs' : {'xscale':'symlog'}
- }
- return WFDict
-
-#SPECIFIC PLOTS
-
-def plotMomentumDistribution(highestrank=1,lowestrank=3):
+def RecoVsTruthQPMatched(event,charge=1):
+ return [RecoVsMCQP(event,charge,'qopr'),RecoVsMCQP(event,charge,'qopm')]
- arr_pmag=[[] for i in range(highestrank,lowestrank+1)]
- for event in t:
- if len(event.Tracks[0])==lowestrank:
- trackQPNsorted=getQPNSorted(event,highestrank=1,lowestrank=lowestrank)
- isSimilarP = True # abs(trackQPNsorted[0][4]-trackQPNsorted[1][4])<100
- if trackQPNsorted[0][0]!=trackQPNsorted[1][0] and isSimilarP:
- for rank in range(highestrank-1,lowestrank):
- arr_pmag[rank]+=[trackQPNsorted[rank][4]]
+def ResidualVSQP_truth(event):
+ return [RecoVsMCQP(event,1,'qopm'),RecoVsMCQP(event,1),RecoVsMCQP(event,-1,'qopm'),RecoVsMCQP(event,-1)]
- for rank in range(highestrank-1,lowestrank):
- rlabel = f'P Rank {rank+1}'
- configBasicHist(arr_pmag[rank],label=rlabel,bins=logspace(start=log10(30), stop=log10(7000), num=25),range=(30,7000), legend=True,kwargs={'histtype':u'step'})
- plt.xscale('log')
- plt.xlabel('Track Momentum Magnitude (GeV) ')
- plt.ylabel('Number of Events')
- plt.title(f"Momentum Distribution of {lowestrank}-track Events ("+user_input+")")
- plt.show()
+def RecoVsMCQP(event,charge=-1,param=None):
+ i=event.trackD2i-1
+ if event.tracksTruthCharge[i]!=charge: i=event.trackD3i-1
-def efficiencyP(pmin=1,pmax=4000,pinc=100):
- arr_eff=[]
- arr_p=[]
- for pcut in range(pmin,pmax,pinc):
- eventcountcut=0
- eventcount=0
- for event in t:
- if ApplyAllCuts(event):
- trackQPNsorted=getQPNSorted(event,highestrank=1,lowestrank=2) #[()]
- if trackQPNsorted[0][0]!=trackQPNsorted[1][0]:
- eventcount+=1
- if trackQPNsorted[0][-1]>pcut and trackQPNsorted[1][-1]>pcut: eventcountcut+=1
- arr_eff+= [100*eventcountcut/eventcount]
- arr_p+= [pcut]
- plt.plot(arr_p,arr_eff)
- plt.ylabel('Cut Efficiency (%)')
- plt.xlabel('P-cut (GeV)')
- plt.title(user_input+'Momentum Cut Efficiency')
- plt.savefig(user_input+'MomCutEff')
- plt.show()
+ q_r=event.tracksCharge[i]
+ p_r=event.tracksPmag[i]
+ q_m=event.tracksTruthCharge[i]
+ p_m=event.tracksTruthPmag[i]
-def dpHitMap():
- arrxp=[]
- arryp=[]
- arrxn=[]
- arryn=[]
- for event in t:
- if ApplyAllCuts(event):
- xp,yp,xn,yn=(event.pPos_S1L0[0],event.pPos_S1L0[1],event.nPos_S1L0[0],event.nPos_S1L0[1])
- if all([abs(i)<200 for i in [xp,yp,xn,yn]]):
- arrxp+=[xp]
- arryp+=[yp]
- arrxn+=[xn]
- arryn+=[yn]
-
- #plt.scatter(arrxp,arryp,color='red',s=0.2)
- configHitMap(arrxp,arryp,title='Positive Track Upstream Hits')
- #plotLEGO(arrxp,arryp)
- plt.savefig(f"{user_input}HitMapP")
- plt.close()
- plt.show()
- #plt.scatter(arrxn,arryn,color='blue',s=0.1)
- configHitMap(arrxn,arryn,title='Negative Track Upstream Hits')
- plt.savefig(f"{user_input}HitMapN")
- #plotLEGO(arrxn,arryn)
- plt.show()
-
-def CutEfficiency(Ndcuts,Nlcuts,Ddcuts=[],Dlcuts=[]):
- ecount=0
- ccount=0
- for event in t:
- if ApplyAllCuts(event,Ddcuts,Dlcuts): ecount+=1
- if ApplyAllCuts(event,Ndcuts,Nlcuts): ccount+=1
- return ccount,ecount,100*ccount/ecount
-
-def CutEfficiencyArrays(xlab='Z (mm)'):
- Func = lambda event: (event.vertexPos[0]**2)+(event.vertexPos[1]**2) #event.vertexPos[2]#
- narray=[]
- darray=[]
- for event in t:
- if ApplyAllCuts(event,[CUTDICT],LAMBDACUTS): narray+=[Func(event)]
- if ApplyAllCuts(event,[],FIDUCIALCUT): darray+=[Func(event)]
-
- # if 'Z' not in xlab:
- # dn,dbin = configBasicHist(darray,bins=20,kwargs={'range':(0,10000),'label':'Fiducial Volume'},legend=True)
- # nn,nbin = configBasicHist(narray,bins=20,xlabel='Rho^2 (mm)',kwargs={'range':(0,10000),'label':'Selection Cuts'},legend=True)
- # else:
- # dn,dbin = configBasicHist(darray,bins=20,kwargs={'range':(-1600,0),'label':'Fiducial Volume'},legend=True)
- # nn,nbin = configBasicHist(narray,bins=20,xlabel=xlab,kwargs={'range':(-1600,0),'label':'Selection Cuts'},legend=True)
-
- # plt.title(user_input+"Efficiency vs "+xlab)
- # plt.show()
-
- RootRatio(narray,darray)
-
- # x=[]
- # y=[]
- # err=[]
- # for i in range(len(nn)):
- # x+=[nbin[i]]
- # y+=[100*nn[i]/dn[i]]
- # err+=[100*(1/dn[i])*sqrt((nn[i]*(dn[i]-nn[i]))/dn[i])]
- # #plt.plot(x,y)
- # plt.errorbar(x,y,yerr=err)
- # plt.ylim(bottom=0)
- # plt.xlabel(xlab)#Rho^2 (mm)
- # plt.ylabel("Efficiency")
- # plt.title(user_input+"Efficiency vs"+xlab)
- # plt.show()
-
-def EfficiencyPlotter(minx,maxx,stepx=1):
- arr=[]
- arrx=[]
- for i in range(minx+stepx,maxx,stepx):
- lcut=lambda event: ((event.vertexPos[0]**2)+(event.vertexPos[1]**2))<i
- arr+=[CutEfficiency([CUTDICT],LAMBDACUTS,dlcuts=[lcut])]
- arrx+=[i]
- plt.plot(arrx,arr)
- plt.ylabel("Fiducial Event Count/Cut-Applied Event Count")
- plt.show()
-
-#FOR PLOT DICTIONARY (complex lambda functions)
-def TwoMomentum(event,sorte=True):
- ppmag=sqrt((event.pMom[0]**2)+(event.pMom[1]**2)+(event.pMom[2]**2))*(10**-3)
- npmag=sqrt((event.nMom[0]**2)+(event.nMom[1]**2)+(event.nMom[2]**2))*(10**-3)
- final = [ppmag,npmag]
- if sorte: final.sort()
- return final
-
-def AngularDistribution(event):
- px1,py1,pz1 = (event.pMom[0],event.pMom[1],event.pMom[2])
- px2,py2,pz2 = (event.nMom[0],event.nMom[1],event.nMom[2])
- return [(arctan(py1/pz1)-arctan(py2/pz2))*1000,(arctan(px1/pz1)-arctan(px2/pz2))*1000]
-
-def PAngleX(event):
- px1,py1,pz1 = (event.pMom[0],event.pMom[1],event.pMom[2])
- px2,py2,pz2 = (event.nMom[0],event.nMom[1],event.nMom[2])
- return [arctan(px1/pz1)*1000,arctan(px2/pz2)*1000]
-
-def PAngleY(event):
- px1,py1,pz1 = (event.pMom[0],event.pMom[1],event.pMom[2])
- px2,py2,pz2 = (event.nMom[0],event.nMom[1],event.nMom[2])
- return [arctan(py1/pz1)*1000,arctan(py2/pz2)*1000]
- #plotLEGO(ndiff_arrx,ndiff_arry)
-
-def UpstreamHitMap(event,charge=1):
- xp,yp,xn,yn=(event.pPos_S1L0[0],event.pPos_S1L0[1],event.nPos_S1L0[0],event.nPos_S1L0[1])
- if charge==1: return [xp,yp]
- elif charge==-1: return [xn,yn]
-
-
-def TrackSeparation(event):
- xp,yp,xn,yn=(event.pPos_S1L0[0],event.pPos_S1L0[1],event.nPos_S1L0[0],event.nPos_S1L0[1])
- if all([abs(i)<200 for i in [xp,yp,xn,yn]]):
- return sqrt((xp-xn)**2+(yp-yn)**2)
- else: return nan
-
-def SeparationDistribution(event):
- xp,yp,zp,xn,yn,zn=(event.pPos_S1L0[0],event.pPos_S1L0[1],event.pPos_S1L0[2],event.nPos_S1L0[0],event.nPos_S1L0[1],event.nPos_S1L0[2])
- if all([abs(i)<200 for i in [xp,yp,xn,yn]]):
- return [(arctan(yp/zp)-arctan(yn/zn))*1000,(arctan(xp/zp)-arctan(xn/zn))*1000]
- else: return [-101,-101] #outside range, therefore not plotted
-
-def Chi2DoF(event):
- try: return [event.pChi2/event.pDoF,event.nChi2/event.nDoF]
- except ZeroDivisionError: return [nan,nan]
-
-def SumCalo(event):
- Etot=0
- for st in range(0,4):
- aname=f'Calo{st}_integral'
- value = getattr(event, aname)
- bool1= value != inf and value !=nan and value !=-1*inf
- if bool1: Etot+=value
- else: return nan
- return Etot
-
-#PLOT DICTIONARY FORMAT
-# arrDict={
-# 'ExampleName':{
-# 'array':[],
-# 'lambda': lambda event: event.pChi2/event.pDoF,
-# 'plotterfunc': configBasicHist,
-# 'plotkwargs' : {'edgecolor':'black','linewidth':1.2, 'rwidth':0.9}
-# }
-# }
-# Veto St1-2, Ch0-1 - VetoStXY_
-# Timing 0-3 - TimingX_
-# Preshower 0-1 - PreshowerX_
-# Calo 0-3 - CaloX_
-# Sum of integrals in calo wave channels --> Calo Edep sort of (missing some callibration constants)
-
-
-#templambda = lambda ti : lambda event : event.vertexPos[2]<ti
-#EfficiencyPlotter(templambda,0,1000)
-
+ if param=='qr': return q_r
+ if param=='pr': return p_r
+ if param=='qm': return q_m
+ if param=='pm': return p_m
+ if param=='qopr': return q_r/p_r
+ if param=='qopm': return q_m/p_m
+ return ((q_r/p_r)-(q_m/p_m))/event.tracksSigmaQovP[i]
+
+def NormalizedResidual2(event):
+ return [RecoVsMCQP(event,1),RecoVsMCQP(event,-1)]
+
+def GetAllTrackQPReco(event,tar=1000):
+ measarr = []
+ for ti in range(event.TrackCount):
+ q_r=event.tracksCharge[ti]
+ p_r=event.tracksPmag[ti]
+ p_m=event.tracksTruthPmag[ti]
+ if type(tar)==list:
+ tar=PH
+ if tar[0]<p_m and p_m<tar[1]: measarr+=[(q_r/p_r)*10**3]
+ #if tar[0]<p_m and p_m<tar[1]: measarr+=[(q_r/p_r)*10**3]
+ #elif p_m>tar*0.9 and p_m<tar*1.2: measarr+=[(q_r/p_r)*10**3]
+
+ return measarr
+
+def GetAllTrackQPTruth(event):
+ measarr = []
+ for ti in range(len(event.truthd0_charge)):
+ q_0=event.truthd0_charge[ti]
+ p_0=event.truthd0_P[ti]
+ measarr+=[(q_0/p_0)*10**3]
+
+ for ti in range(len(event.truthd1_charge)):
+ q_1=event.truthd1_charge[ti]
+ p_1=event.truthd1_P[ti]
+ measarr+=[(q_1/p_1)*10**3]
+
+ return measarr
+
+def GetAllTrackMomentum(event):
+ measarr = []
+
+ for ti in range(event.TrackCount):
+ p_r=event.tracksPmag[ti]
+ measarr+=[p_r*10**-3]
+
+ return finalval
+
+
+#For storing different plot configurations
MainArrDict={
- 'CaloSum':{
+ 'RecoX':{
'array':[],
- 'lambda': SumCalo,
+ 'lambda': lambda event: event.tracksX[0],
'plotterfunc': configBasicHist,
- 'plotkwargs' : {'xscale':'log'}
+ 'plotkwargs' : {}
},
- 'TracksPerEvent':{
+ 'RecoMaxMomentum':{
'array':[],
- 'lambda': lambda event: len(event.Tracks[0]),
+ 'lambda': lambda event:max(event.tracksPmag),
'plotterfunc': configBasicHist,
- 'plotkwargs' : {}
+ 'plotkwargs' : {'pckwargs':{'xscale':'log'}}
},
- 'TrackSeparation':{
+ 'TruthMaxMomentum':{
'array':[],
- 'lambda': TrackSeparation,
+ 'lambda': lambda event:max(event.tracksTruthPmag),
'plotterfunc': configBasicHist,
- 'plotkwargs' : {'bins':20,'xlabel':'Distance (mm)','kwargs':{'range':(0,10)}}
- },
- 'TwoMomentum':{
- 'array':[[],[]],
- 'lambda': TwoMomentum,
- 'plotterfunc': configMultiHist,
- 'plotkwargs' : {'label':['Positive','Negative'],'xscale':'log','xlabel':'Track Momentum Magnitude (GeV) '}
- },
- 'AngularDistribution':{
- 'array':[[],[]],
- 'lambda': AngularDistribution,
- 'plotterfunc': configMultiHist,
- 'plotkwargs' : {'label':['Momentum-Y','Momentum-X'],'xlabel':'N-diff (mrad)','kwargs':{'range':(-35,35)}}
- },
- 'PAngleX':{
- 'array':[[],[]],
- 'lambda': PAngleX,
- 'plotterfunc': configMultiHist,
- 'plotkwargs' : {'label':['Positive','Negative'],'xlabel':'Angle (mrad)','kwargs':{'range':(-35,35)}}
+ 'plotkwargs' : {'pckwargs':{'xscale':'log'}}
},
- 'PAngleY':{
+ 'RecoVsTruthMaxMomentum':{
'array':[[],[]],
- 'lambda': PAngleY,
+ 'lambda': lambda event:[max(event.tracksPmag),max(event.tracksTruthPmag)],
'plotterfunc': configMultiHist,
- 'plotkwargs' : {'label':['Positive','Negative'],'xlabel':'Angle (mrad)','kwargs':{'range':(-35,35)}}
+ 'plotkwargs' : {'label':['Reco','Truth'],'pckwargs':{'xscale':'log','legend':True}}
},
- 'SeparationDistribution':{
+ 'RecoVsTruthMaxMomentumScatter':{
'array':[[],[]],
- 'lambda': SeparationDistribution,
- 'plotterfunc': configMultiHist,
- 'plotkwargs' : {'label':['Y','X'],'kwargs':{'range':(-100,100)}}
+ 'lambda': lambda event:[max(event.tracksPmag),max(event.tracksTruthPmag)],
+ 'plotterfunc': config2DScatter,
+ 'plotkwargs' : {}
},
- 'Chi2DoF':{
- 'array':[[],[]],
- 'lambda': Chi2DoF,
+ 'QP_reco Vs QP_truth':{
+ 'array':[[],[],[],[]],
+ 'lambda': RecoVsTruthQPMatched,
'plotterfunc': configMultiHist,
- 'plotkwargs' : {'bins':20,'label':['Positive','Negative'],'kwargs':{'range':(0,11)}}
+ 'plotkwargs' : {'label':['Reco +','Truth +','Reco -','Truth -'],'kwargs':{'range':[-10**-5,10**-5]},'pckwargs':{'xscale':'linear','xlabel':"Q/P",'legend':True}}
},
- 'DoF':{
+ 'QP_truth_vs_Residual':{
'array':[[],[]],
- 'lambda': lambda event:[event.pDoF,event.nDoF],
- 'plotterfunc': configMultiHist,
- 'plotkwargs' : {'bins':30,'label':['Positive','Negative'],'kwargs':{'range':(0,31)}}
+ 'lambda': ResidualVSQP_truth,
+ 'plotterfunc': config2DScatter,
+ 'plotkwargs' : {'kwargs':{},'pckwargs':{'xlabel':"Q/P_truth",'ylabel':"Residual"}}
},
- 'LayersHit':{
+ 'NormalizedResidual':{
'array':[[],[]],
- 'lambda': lambda event: [event.playersHit,event.nlayersHit],
+ 'lambda': NormalizedResidual2,
'plotterfunc': configMultiHist,
- 'plotkwargs' : {'bins':5,'label':['Positive','Negative'],'kwargs':{'range':(5,10)}}
+ 'plotkwargs' : {'bins':20,'kwargs':{'range':[-10,11]},'label':['Positive','Negative'],'pckwargs':{'xscale':'linear','xlabel':"Residual",'legend':True}}
},
- 'VertexZ':{
+ 'AllTracksQP_reco':{
'array':[],
- 'lambda': lambda event: event.vertexPos[2],
+ 'lambda': GetAllTrackQPReco,
'plotterfunc': configBasicHist,
- 'plotkwargs' : {'xlabel':'Z (mm)'}
+ 'plotkwargs' :{'bins':100,'kwargs':{'range':[2*-10**-3,2*10**-3]}},
},
- 'TrueVertexRho':{
+ 'AllTracksQP_truth':{
'array':[],
- 'lambda': lambda event: (event.vertexPos[0]**2)+(event.vertexPos[1]**2),
+ 'lambda': GetAllTrackQPTruth,
'plotterfunc': configBasicHist,
- 'plotkwargs' : {'xlabel':'Rho^2 (mm^2)'}
+ 'plotkwargs' :{'kwargs':{'range':[2*-10**-3,2*10**-3]}},
},
- 'Rho2vsVertexZ':{
- 'array':[[],[]],
- 'lambda': lambda event:[(event.vertexPos[0]**2)+(event.vertexPos[1]**2),event.vertexPos[2]],
- 'plotterfunc': configHitMap2,
- 'plotkwargs' : {'bins':(10,10),'xlabel':'Rho^2 (mm^2)','ylabel':'Vertex-Z (mm)'}
+ 'AllMomentum_Reco':{
+ 'array':[],
+ 'lambda': GetAllTrackMomentum,
+ 'plotterfunc': configBasicHist,
+ 'plotkwargs' : {'kwargs':{'range':[0,4000]}}
},
- 'UpstreamHitsPositive':{
- 'array':[[],[]],
- 'lambda': UpstreamHitMap,
- 'plotterfunc': configHitMap2,
- 'plotkwargs' : {'xlim':(-100,100),'ylim':(-100,100)},
- 'lkwargs':{'charge':1}
+ 'AllTracksQP_reco':{
+ 'array':[],
+ 'lambda': GetAllTrackQPReco,
+ 'plotterfunc': configBasicHist,
+ 'plotkwargs' :{'bins':50,'kwargs':{'range':[2*-10**-3,2*10**-3]}},
+ 'lkwargs' : {'tar':2000*10**3}
},
- 'UpstreamHitsNegative':{
- 'array':[[],[]],
- 'lambda': UpstreamHitMap,
- 'plotterfunc': configHitMap2,
- 'plotkwargs' : {'xlim':(-100,100),'ylim':(-100,100)},
- 'lkwargs':{'charge':-1}
+ 'AllTracksTruthPmag':{
+ 'array':[],
+ 'lambda': lambda event: list(event.tracksTruthPmag),
+ 'plotterfunc': configBasicHist,
+ 'plotkwargs' :{'bins':50,'pckwargs':{'xscale':'log'}},
},
- 'UpstreamXYPositive':{
- 'array':[[],[]],
- 'lambda': lambda event: (event.pPos_S1L0[0],event.pPos_S1L0[1]),
- 'plotterfunc': configMultiHist,
- 'plotkwargs' : {'label':['Positive','Negative']},
+ 'NormalizedResidualSum':{
+ 'array':[],
+ 'lambda': NormalizedResidual2,
+ 'plotterfunc': configBasicHist,
+ 'plotkwargs' : {'bins':20,'kwargs':{'range':[-10,11]},'pckwargs':{'xscale':'linear','xlabel':"Residual"}}
},
-}
-MainArrDict.update(MakeWFDict()) #Fill with wf plots
+ }
+
+
-#For testing/plotting purposes, plotting specific members of MainArrDict
TestArrDict=dict()
-testname='TrueVertexRho'
+testname='NormalizedResidualSum' #Any key from MainArrDict
TestArrDict[testname]=MainArrDict[testname]
# for k in MainArrDict.keys():
-# if 'Calo' in k:
+# if 'Check' in k:
# TestArrDict[k]=MainArrDict[k]
-def GeneralPlotter(arrDict,dcuts=[CUTDICT],lcuts=LAMBDACUTS,showplot=True,savefig=False):
+def GeneralPlotter(arrDict,dcuts=[],lcuts=[],showplot=True,savefig=False):
nans=0
+
for event in t:
if ApplyAllCuts(event,dcuts,lcuts):
for name in arrDict.keys():
value=nan
- if arrDict[name]['plotterfunc'] in [configMultiHist,configHitMap2]:
- if 'lkwargs' in list(arrDict[name].keys()): value = arrDict[name]['lambda'](event,**arrDict[name]['lkwargs'])
- else: value = arrDict[name]['lambda'](event)
- checkReal = all([v != inf and v !=nan and v !=-1*inf for v in value])
- if checkReal:
+ if 'lkwargs' in list(arrDict[name].keys()):
+ value = arrDict[name]['lambda'](event,**arrDict[name]['lkwargs'])
+ else: value = arrDict[name]['lambda'](event)
+ 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 checkReal and type(value)==list:
+ if arrDict[name]['plotterfunc'] in [configBasicHist]: arrDict[name]['array']+=value
+ else:
for iv in range(len(value)):
arrDict[name]['array'][iv]+=[value[iv]]
- else: nans+=1
- else:
- if 'lkwargs' in list(arrDict[name].keys()): value = arrDict[name]['lambda'](event,**arrDict[name]['lkwargs'])
- else: value = arrDict[name]['lambda'](event)
- checkReal = value != inf and value !=nan and value !=-1*inf
- if checkReal: arrDict[name]['array']+=[value]
- else: nans+=1
-
+ elif checkReal: arrDict[name]['array']+=[value]
+ else: nans+=1
+
print(f"# of NaNs:{nans}")
for name in arrDict.keys():
+ isPrintout= arrDict[name]['plotterfunc'] in [PrintBoolPercent]
+ print(name)
+ print('LEN:',len(arrDict[name]['array']))
arrDict[name]['plotterfunc']( arrDict[name]['array'], **arrDict[name]['plotkwargs'])
- plt.title(user_input+name)
- if savefig: plt.savefig(user_input+name)
- if showplot: plt.show()
-
-
-
-#Iterate over all models and do a thing
-DP_fnames=['Ntuple_Amm_316MeV_2Em6_','Ntuple_Aee_10MeV_1Em5_','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/"+user_input+".root")
- print(f)
-
- #GeneralPlotter(TestArrDict,dcuts=[CUTDICT],lcuts=LAMBDACUTS)
- #EfficiencyPlotter(0,200,stepx=10)
- #print(CutEfficiency([CUTDICT],LAMBDACUTS,Ddcuts=[],Dlcuts=FIDUCIALCUT))
- CutEfficiencyArrays(xlab='Rho^2 (mm^2)')
- #EfficiencyTest()
-
+ if not isPrintout:
+ plt.title(user_input+name)
+ if savefig: plt.savefig(user_input+name)
+ if showplot: plt.show()
+
+#EXECUTE
+if __name__ == "__main__":
+ DP_fnames=['Ntuple_Aee_100MeV_1Em5_']# Ntuple_Aee_10MeV_1Em5_ Ntuple_Amm_316MeV_2Em6_ Ntuple_Aee_10MeV_1Em4_
+ 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")
+ print(f)
+
+ GeneralPlotter(TestArrDict,dcuts=[{}],lcuts=[BoolDaughterCut])
+ #print(GetPBounds(lcuts=[BoolDaughterCut]))
diff --git a/Tracker/TrackerRecAlgs/PairVertex/scripts/PlotHelper.py b/Tracker/TrackerRecAlgs/PairVertex/scripts/PlotHelper.py
new file mode 100644
index 0000000000000000000000000000000000000000..8597d9cab9779c8ff2fd2e04744730a6f2f3efa7
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/PairVertex/scripts/PlotHelper.py
@@ -0,0 +1,102 @@
+from numpy import log10,logspace
+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'):
+ plt.xscale(xscale)
+ plt.yscale(yscale)
+ if xticks: plt.xticks(xticks)
+ if xlabel: plt.xlabel(xlabel)
+ if ylabel: plt.ylabel(ylabel)
+ if legend: plt.legend(loc=legloc)
+ if title:plt.title(title)
+ if show: plt.show()
+
+def configHitMap(arrays,bins=(50,50),pckwargs={'xlabel':"X (mm)",'ylabel':"Y (mm)"},counts=False,kwargs={}):
+ #2D Histogram with colors
+ arrayx=arrays[0]
+ arrayy=arrays[1]
+ plt.hist2d(arrayx, arrayy, bins=bins, cmap=cm.imola,**kwargs)
+ plt.colorbar()
+ if counts: plt.text(50,-185,"#EVENTS: {}".format(len(arrayx)))
+ PlotConfigs(**pckwargs)
+
+def configBasicHist(arr,bins=25,pckwargs={'xscale':'linear'},kwargs={}):
+ #1D Histogram
+ if 'log' in pckwargs['xscale']: bins=CorrectBinsForLog(arr,bins,pckwargs['xscale'])
+ n, binsOut, patches = plt.hist(arr,bins=bins,align="left",histtype=u'step',**kwargs)
+ PlotConfigs(**pckwargs)
+
+ return n,binsOut
+
+def configMultiHist(arrlist,bins=25,label=[],pckwargs={'xscale':'linear','legloc':'upper left'},kwargs={}):
+ #1D Histogram for comparison
+ megarr=[]
+ for a in arrlist: megarr+=a
+
+ if 'log' in pckwargs['xscale']: bins=CorrectBinsForLog(megarr,bins,pckwargs['xscale'])
+ for iarr in range(len(arrlist)):
+ arr=arrlist[iarr]
+ labelused=''
+ if len(arr)==0:
+ print(f"Empty array.")
+ else:
+ if len(label)==len(arrlist): labelused=label[iarr]
+ else:
+ print("Insufficient labels given.")
+ plt.hist(arr,bins=bins,align="left",histtype=u'step',label=labelused,**kwargs)
+
+ PlotConfigs(**pckwargs)
+
+def config2DScatter(arrlist,pckwargs=dict(),kwargs={}):
+ x=arrlist[0]
+ y=arrlist[1]
+ # plt.ylim(-25,25)
+ # plt.xlim(-10**-5,10**-5)
+ plt.scatter(x,y,**kwargs)
+ PlotConfigs(**pckwargs)
+
+def PrintBoolPercent(arr,kwargs={}):
+ #For an array of boolians, print out the % of elements that are "true".
+ n=0.0
+ d=len(arr)
+ for x in arr:
+ if x: n+=1.0
+
+ finalval= (n/d)*100
+ print("Calculating...")
+ print(f"{n}/{d}")
+ print('Percent:'+'{:.2f}'.format(finalval)+'%')
+
+def PlotRootRatio(arrlist,bins=15,hmin=0,hmax=10000,title=None):
+
+ narray,darray=arrlist[0],arrlist[1]
+ h_numerator=TH1F("num","num",bins,hmin,hmax)
+ h_denominator=TH1F("den","den",bins,hmin,hmax)
+ for i in narray: h_numerator.Fill(i)
+ for j in darray: h_denominator.Fill(j)
+
+ h_ratio = TEfficiency(h_numerator,h_denominator)
+ c = h_ratio.CreateGraph()
+ c.GetYaxis().SetLimits(0.0, 1.05)
+ if title: h_ratio.SetTitle(title)
+ h_ratio.Draw("ep")
+ input()
+
+#Utilities
+def CorrectBinsForLog(arr,bins,xscale):
+ minx,maxx = (min(arr),max(arr))
+ if minx==maxx:
+ maxx+=1
+ bins=1
+ bins = logspace(start=log10(minx), stop=log10(maxx), num=bins)
+ if xscale=='symlog':
+ bins = [0]+list(logspace(start=log10(1), stop=log10(maxx), num=24))
+ if minx>1: bins = logspace(start=log10(minx), stop=log10(maxx), num=25)
+
+ return bins
+
+
+
+
diff --git a/Tracker/TrackerRecAlgs/PairVertex/scripts/Utilities.py b/Tracker/TrackerRecAlgs/PairVertex/scripts/Utilities.py
new file mode 100644
index 0000000000000000000000000000000000000000..ff711357eed1ed411749dcf9ed2e800546f358c3
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/PairVertex/scripts/Utilities.py
@@ -0,0 +1,65 @@
+#UTILITY
+from ROOT import TChain, TH2F, TH1F
+print("Root loaded.")
+import matplotlib.pyplot as plt
+print("MPL loaded.")
+from operator import lt,eq,ne,gt
+from numpy import nan, seterr, inf,arctan, array
+from numpy import sum as npsum
+import numpy as np
+from math import sqrt
+print("Other stuff loaded.")
+
+user_input=''
+
+seterr(all='ignore')
+missingattr=[]
+
+#CUTS
+#attr: [>,<,==,!=] if none --> skip
+BASICCUT={
+ 'TrackCount':[0,None,None,None]
+ }
+
+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=[]
+ global missingattr
+ for attr in cut.keys():
+ for i in range(len(tests)):
+ try:
+ if cut[attr][i]!=None:
+ value = getattr(event,attr)
+ bool1 += [tests[i](value,cut[attr][i]) and value != inf and value !=nan and value !=-1*inf]
+ if not all(bool1): break
+ except AttributeError:
+ if attr not in missingattr:
+ print("No attribute: ",attr)
+ missingattr+=[attr]
+
+ if not all(bool1):
+ break
+
+ return all(bool1)
+
+def ApplyAllCuts(event,dcuts=[],lcuts=[]): #Requires all True to be True
+ #cutList = list of dictionaries
+ bool1=True
+ boollist=[]
+ for cut in dcuts:
+ bool1 = ApplyCut(event,cut)
+ if not bool1: break
+ if bool1:
+ for lcut in lcuts:
+ boollist+=[lcut(event)]
+ if not all(boollist): break
+
+ return all(boollist) and bool1
+
diff --git a/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.cxx b/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.cxx
index 8009ce9a72e0efb119b0e5cf52c52d99a2add749..a44ba24c44f2b516d4f9a7769d8bad8b7405810e 100644
--- a/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.cxx
+++ b/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.cxx
@@ -14,13 +14,14 @@
#include "FaserCaloIdentifier/EcalID.h"
#include "TrackerRIO_OnTrack/FaserSCT_ClusterOnTrack.h"
#include <cmath>
+#include <numeric>
#include <TTree.h>
#include <TBranch.h>
#include <typeinfo>
-
+constexpr float NaN = std::numeric_limits<double>::quiet_NaN();
namespace Tracker {
@@ -89,6 +90,10 @@ StatusCode PairVertexAlg::initialize()
ATH_CHECK(m_preshowerContainer.initialize());
ATH_CHECK(m_ecalContainer.initialize());
+ ATH_CHECK(m_fiducialParticleTool.retrieve());
+ ATH_CHECK(m_trackTruthMatchingTool.retrieve());
+ ATH_CHECK(m_truthParticleContainer.initialize());
+
ATH_CHECK(detStore()->retrieve(m_idHelper, "FaserSCT_ID"));
ATH_CHECK(detStore()->retrieve(m_detMgr, "SCT"));
ATH_CHECK(detStore()->retrieve(m_vetoNuHelper, "VetoNuID"));
@@ -98,31 +103,39 @@ StatusCode PairVertexAlg::initialize()
ATH_CHECK(detStore()->retrieve(m_ecalHelper, "EcalID"));
m_tree = new TTree("events", "Events");
- m_tree->Branch("Tracks", &m_tracks);
+
m_tree->Branch("TrackCount", &m_trackCount, "trackCount/I");
- //POSITIVE
- m_tree->Branch("pPos_S1L0", &m_ppos_s1l0); // [x,y,z]
- m_tree->Branch("pMom", &m_pmom); // [px,py,pz]
- m_tree->Branch("pChi2", &m_pChi2); // float
- m_tree->Branch("pDoF", &m_pDoF); // double
- m_tree->Branch("pstationMap", &m_pstationMap); // set (int)
- m_tree->Branch("playerMap", &m_playerMap); // vec(vec(int))
- m_tree->Branch("playersHit", &m_playersHit); // int
- m_tree->Branch("piftHit", &m_piftHit); // int
-
- //NEGATIVE
- m_tree->Branch("nPos_S1L0", &m_npos_s1l0); // [x,y,z]
- m_tree->Branch("nMom", &m_nmom); // [px,py,pz]
- m_tree->Branch("nChi2", &m_nChi2); // float
- m_tree->Branch("nDoF", &m_nDoF); // double
- m_tree->Branch("nstationMap", &m_nstationMap); // set (int)
- m_tree->Branch("nlayerMap", &m_nlayerMap); // vec(vec(int))
- m_tree->Branch("nlayersHit", &m_nlayersHit); // int
- m_tree->Branch("niftHit", &m_niftHit); // int
-
- //CALO
- m_tree->Branch("Calo_total_Edep", &m_Calo_Total_Edep);
+ m_tree->Branch("tracksX", &m_tracksX);
+ m_tree->Branch("tracksY", &m_tracksY);
+ m_tree->Branch("tracksZ", &m_tracksZ);
+ m_tree->Branch("tracksPx", &m_tracksPx);
+ m_tree->Branch("tracksPy", &m_tracksPy);
+ m_tree->Branch("tracksPz", &m_tracksPz);
+ m_tree->Branch("tracksPmag", &m_tracksPmag);
+
+ m_tree->Branch("tracksCharge", &m_tracksCharge);
+ m_tree->Branch("tracksSigmaQovP", &m_sigmaQoverP);
+ m_tree->Branch("tracksChi2", &m_tracksChi2);
+ m_tree->Branch("tracksDoF", &m_tracksDoF);
+ m_tree->Branch("trackPositiveIndicies", &m_trackPositiveIndicies);
+ m_tree->Branch("trackNegativeIndicies", &m_trackNegativeIndicies);
+
+ m_tree->Branch("tracksTruthBarcode", &m_tracksTruthBarcode);
+ m_tree->Branch("tracksTruthPdg", &m_tracksTruthPdg);
+ m_tree->Branch("tracksTruthCharge", &m_tracksTruthCharge);
+ m_tree->Branch("tracksTruthPmag", &m_tracksTruthPmag);
+ m_tree->Branch("tracksIsFiducial", &m_tracksIsFiducial);
+
+ m_tree->Branch("trackD2i", &m_trackD2i); //Index of track with daughter particle Barcode 2
+ m_tree->Branch("trackD3i", &m_trackD3i); //Index of track with daughter particle Barcode 3
+
+ m_tree->Branch("tracksNLayersHit", & m_tracksnLayersHit);
+ m_tree->Branch("tracksLayersIFT", & m_tracksLayersIFT);
+ m_tree->Branch("tracksIFTHit", & m_tracksIFThit);
+
+ m_tree->Branch("tracksLayerMap", &m_layerMap);
+
//WAVEFORM
addWaveBranches("VetoSt1",2,6);
@@ -133,7 +146,44 @@ StatusCode PairVertexAlg::initialize()
addWaveBranches("Calo",4,0);
//TRUTH
- m_tree->Branch("vertexPos", &m_vertexPos); // [x,y,z]
+ m_tree->Branch("tTrkBarcode", &m_truthBarcode);
+ m_tree->Branch("tTrkPdg", &m_truthPdg);
+
+ m_tree->Branch("vertexPos", &m_vertexPos); // [x,y,z]
+ m_tree->Branch("isFiducial", &m_isFiducial);
+
+ m_tree->Branch("truthM_P", &m_truthM_P);
+ m_tree->Branch("truthM_px", &m_truthM_px);
+ m_tree->Branch("truthM_py", &m_truthM_py);
+ m_tree->Branch("truthM_pz", &m_truthM_pz);
+ m_tree->Branch("truthM_x", &m_truthM_x);
+ m_tree->Branch("truthM_y", &m_truthM_y);
+ m_tree->Branch("truthM_z", &m_truthM_z);
+
+ //Electron
+ m_tree->Branch("truthd0_P", &m_truthd0_P);
+ m_tree->Branch("truthd0_px", &m_truthd0_px);
+ m_tree->Branch("truthd0_py", &m_truthd0_py);
+ m_tree->Branch("truthd0_pz", &m_truthd0_pz);
+ m_tree->Branch("truthd0_x", &m_truthd0_x);
+ m_tree->Branch("truthd0_y", &m_truthd0_y);
+ m_tree->Branch("truthd0_z", &m_truthd0_z);
+ m_tree->Branch("truthd0_IsFiducial", &m_truthd0_IsFiducial);
+ m_tree->Branch("truthd0_charge", &m_truthd0_charge);
+
+ //Positron
+ m_tree->Branch("truthd1_P", &m_truthd1_P);
+ m_tree->Branch("truthd1_px", &m_truthd1_px);
+ m_tree->Branch("truthd1_py", &m_truthd1_py);
+ m_tree->Branch("truthd1_pz", &m_truthd1_pz);
+ m_tree->Branch("truthd1_x", &m_truthd1_x);
+ m_tree->Branch("truthd1_y", &m_truthd1_y);
+ m_tree->Branch("truthd1_z", &m_truthd1_z);
+ m_tree->Branch("truthd1_IsFiducial", &m_truthd1_IsFiducial);
+ m_tree->Branch("truthd1_charge", &m_truthd1_charge);
+
+
+
ATH_CHECK(histSvc()->regTree("/HIST/events", m_tree));
return StatusCode::SUCCESS;
@@ -159,9 +209,75 @@ StatusCode PairVertexAlg::execute(const EventContext &ctx) const
FillWaveBranches(*triggerContainer);
FillWaveBranches(*preshowerContainer);
FillWaveBranches(*ecalContainer);
- FillWaveBranches(*vetoNuContainer);
+ FillWaveBranches(*vetoNuContainer);
+
+ //FILL WITH TRUTHPARTICLECONTAINER
+ SG::ReadHandle<xAOD::TruthParticleContainer> truthParticleContainer { m_truthParticleContainer, ctx };
+ if (truthParticleContainer.isValid() && truthParticleContainer->size() > 0)
+ {
+ for (auto particle : *truthParticleContainer)
+ {
+ if ( particle->barcode() <= 3)
+ {
- // TRACK COLLECTION
+ if ( particle->barcode() == 1) // mother particle (A')
+ {
+
+ m_truthM_P.push_back(particle->p4().P());
+ m_truthM_px.push_back(particle->p4().X());
+ m_truthM_py.push_back(particle->p4().Y());
+ m_truthM_pz.push_back(particle->p4().Z());
+ if (m_fiducialParticleTool->isFiducial(particle->barcode())) {ATH_MSG_ERROR("DP has hits! Something is wrong");}
+ if ( particle->hasProdVtx()) {
+ m_truthM_x.push_back(particle->prodVtx()->x());
+ m_truthM_y.push_back(particle->prodVtx()->y());
+ m_truthM_z.push_back(particle->prodVtx()->z());
+ } else {
+ m_truthM_x.push_back(NaN);
+ m_truthM_y.push_back(NaN);
+ m_truthM_z.push_back(NaN);
+ }}
+
+ if ( particle->charge() < 0) // negative daughter particle
+ {
+ m_truthd0_P.push_back(particle->p4().P());
+ m_truthd0_px.push_back(particle->p4().X());
+ m_truthd0_py.push_back(particle->p4().Y());
+ m_truthd0_pz.push_back(particle->p4().Z());
+ m_truthd0_charge.push_back(particle->charge());
+ m_truthd0_IsFiducial.push_back(m_fiducialParticleTool->isFiducial(particle->barcode()));
+
+ if ( particle->hasProdVtx()) {
+ m_truthd0_x.push_back(particle->prodVtx()->x());
+ m_truthd0_y.push_back(particle->prodVtx()->y());
+ m_truthd0_z.push_back(particle->prodVtx()->z());
+ } else {
+ m_truthd0_x.push_back(NaN);
+ m_truthd0_y.push_back(NaN);
+ m_truthd0_z.push_back(NaN);
+ }}
+
+ if ( particle->charge() > 0) // positive daughter particle
+ {
+ m_truthd1_P.push_back(particle->p4().P());
+ m_truthd1_px.push_back(particle->p4().X());
+ m_truthd1_py.push_back(particle->p4().Y());
+ m_truthd1_pz.push_back(particle->p4().Z());
+ m_truthd1_charge.push_back(particle->charge());
+ m_truthd1_IsFiducial.push_back(m_fiducialParticleTool->isFiducial(particle->barcode()));
+
+ if ( particle->hasProdVtx()) {
+ m_truthd1_x.push_back(particle->prodVtx()->x());
+ m_truthd1_y.push_back(particle->prodVtx()->y());
+ m_truthd1_z.push_back(particle->prodVtx()->z());
+ } else {
+ m_truthd1_x.push_back(NaN);
+ m_truthd1_y.push_back(NaN);
+ m_truthd1_z.push_back(NaN);
+ }}
+ }}}
+
+ //FILL WITH TRACKCOLLECTION
SG::ReadHandle<TrackCollection> tracks { m_trackCollectionKey, ctx };
ATH_CHECK(tracks.isValid());
for (auto trk : *tracks)
@@ -204,14 +320,31 @@ StatusCode PairVertexAlg::execute(const EventContext &ctx) const
m_numberOfNoUpstreamTracks++;
ATH_MSG_WARNING("Unable to find track parameters on any surface");
continue;
- //End error catching; assign state to upstream if everything else is ok.
+ }//End error catching; assign state to upstream if everything else is ok.
+
+
+ //Truth info
+ auto [truthParticle, hitCount] = m_trackTruthMatchingTool->getTruthParticle(trk);
+ if (truthParticle != nullptr) {
+ m_tracksTruthPdg.push_back(truthParticle->pdgId());
+ m_tracksTruthBarcode.push_back(truthParticle->barcode());
+ m_tracksIsFiducial.push_back(m_fiducialParticleTool->isFiducial(truthParticle->barcode()));
+ m_tracksTruthCharge.push_back(truthParticle->charge());
+ m_tracksTruthPmag.push_back(truthParticle->p4().P());
+ if (truthParticle->barcode()==2) m_trackD2i=m_trackCount;
+ if (truthParticle->barcode()==3) m_trackD3i=m_trackCount;
+ }
+ else {
+ m_tracksTruthPdg.push_back(NaN);
+ m_tracksTruthBarcode.push_back(NaN);
+ m_tracksIsFiducial.push_back(NaN);
+ m_tracksTruthCharge.push_back(NaN);
+ m_tracksTruthPmag.push_back(NaN);
}
-
- // Check for hit in the three downstream stations
+ //LAYERMAPPING
std::vector<std::vector<int>> layerMap={{0,0,0},{0,0,0},{0,0,0},{0,0,0}};
std::set<int> stationMap;
- std::vector<double> PH_pos_S1L0={0,0,0}; //x,y,z
for (auto measurement : *(trk->measurementsOnTrack())) {
const Tracker::FaserSCT_ClusterOnTrack* cluster = dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*>(measurement);
@@ -221,64 +354,45 @@ StatusCode PairVertexAlg::execute(const EventContext &ctx) const
int layer = m_idHelper->layer(id);
stationMap.emplace(station);
layerMap[station][layer]+=1;
- }
- }
-
+ }}
int layersHit=0;
int iftHit=0;
for (int st=0;st<4;st++) {
for (int lr=0;lr<3;lr++){
if (layerMap[st][lr]>0) layersHit++;
if (st==0 && layerMap[st][lr]>0) iftHit++;
- }
- }
+ }}
if (stationMap.count(1) == 0 || stationMap.count(2) == 0 || stationMap.count(3) == 0) continue;
+
+ m_tracksnLayersHit.push_back(layersHit);
+ m_tracksIFThit.push_back(iftHit);
+
+ //TRACKPARAMETERS
const Trk::TrackParameters* upstreamParameters = trk->trackParameters()->front();
const Trk::TrackParameters* downstreamParameters = trk->trackParameters()->back();
+ auto covariance = *upstreamParameters->covariance();
+ double sigmaQoverP = sqrt(covariance(4,4));
+ m_sigmaQoverP.push_back(sigmaQoverP);
auto momentum = upstreamParameters->momentum();
double charge = upstreamParameters->charge();
auto position = upstreamParameters->position();
-
- m_tracks[0].push_back(charge);
- m_tracks[1].push_back(position.x());
- m_tracks[2].push_back(position.y());
- m_tracks[3].push_back(position.z());
- m_tracks[4].push_back(momentum.x());
- m_tracks[5].push_back(momentum.y());
- m_tracks[6].push_back(momentum.z());
-
- if (charge>0){
- m_ppos_s1l0[0]=position.x();
- m_ppos_s1l0[1]=position.y();
- m_ppos_s1l0[2]=position.z();
- m_pmom[0]=momentum.x();
- m_pmom[1]=momentum.y();
- m_pmom[2]=momentum.z();
- m_pChi2=trk->fitQuality()->chiSquared();
- m_pDoF=trk->fitQuality()->numberDoF();
- m_pstationMap=stationMap;
- m_playerMap=layerMap;
- m_playersHit=layersHit;
- m_piftHit=iftHit;
- }
- else{
- m_npos_s1l0[0]=position.x();
- m_npos_s1l0[1]=position.y();
- m_npos_s1l0[2]=position.z();
- m_nmom[0]=momentum.x();
- m_nmom[1]=momentum.y();
- m_nmom[2]=momentum.z();
- m_nChi2=trk->fitQuality()->chiSquared();
- m_nDoF=trk->fitQuality()->numberDoF();
- m_nstationMap=stationMap;
- m_nlayerMap=layerMap;
- m_nlayersHit=layersHit;
- m_niftHit=iftHit;
-
- }
+ m_tracksCharge.push_back(charge);
+ if (charge>0){ m_trackPositiveIndicies.push_back(m_trackCount);}
+ if (charge>0) {m_trackNegativeIndicies.push_back(m_trackCount);}
+ m_tracksChi2.push_back(trk->fitQuality()->chiSquared());
+ m_tracksDoF.push_back(trk->fitQuality()->numberDoF());
+
+ m_tracksX.push_back(position.x());
+ m_tracksY.push_back(position.y());
+ m_tracksZ.push_back(position.z());
+
+ m_tracksPx.push_back(momentum.x());
+ m_tracksPy.push_back(momentum.y());
+ m_tracksPz.push_back(momentum.z());
+ m_tracksPmag.push_back(momentum.mag());
++m_trackCount;
@@ -305,7 +419,6 @@ StatusCode PairVertexAlg::execute(const EventContext &ctx) const
m_vertexPos[2]=v->position().z();
break;
}
-
m_tree->Fill();
@@ -325,26 +438,39 @@ StatusCode PairVertexAlg::finalize()
void PairVertexAlg::clearTree() const
{
- for (int i = 0; i < 7; i++) {m_tracks[i].clear();}
m_trackCount=0;
+ m_tracksPmag.clear();
+ m_sigmaQoverP.clear();
+
+ m_tracksX.clear();
+ m_tracksY.clear();
+ m_tracksZ.clear();
+ m_tracksPx.clear();
+ m_tracksPy.clear();
+ m_tracksPz.clear();
+ m_tracksPmag.clear();
+
+ m_tracksCharge.clear();
+ m_sigmaQoverP.clear();
+ m_tracksChi2.clear();
+ m_tracksDoF.clear();
+
+ m_trackPositiveIndicies.clear();
+ m_trackNegativeIndicies.clear();
+ m_trackD2i=-1;
+ m_trackD3i=-1;
+
+ m_tracksTruthBarcode.clear();
+ m_tracksTruthPdg.clear();
+ m_tracksTruthCharge.clear();
+ m_tracksTruthPmag.clear();
+ m_tracksIsFiducial.clear();
+
+ m_tracksnLayersHit.clear();
+ m_tracksLayersIFT.clear();
+ m_tracksIFThit.clear();
+
- m_pstationMap.clear();
- m_playerMap={{},{},{},{}};
- m_ppos_s1l0= {-1,-1,-1};
- m_pmom= {-1,-1,-1};
- m_pChi2=0;
- m_pDoF=0;
- m_playersHit=0;
- m_piftHit=0;
-
- m_nstationMap.clear();
- m_nlayerMap={{},{},{},{}};
- m_npos_s1l0= {-1,-1,-1};
- m_nmom= {-1,-1,-1};
- m_nChi2=0;
- m_nDoF=0;
- m_nlayersHit=0;
- m_niftHit=0;
for(int ii=0;ii<15;ii++) {
m_wave_localtime[ii]=0;
@@ -359,6 +485,23 @@ void PairVertexAlg::clearTree() const
m_wave_status[ii]=0;
}
m_vertexPos= {-1,-1,-1};
+ m_isFiducial.clear();
+ m_truthPdg.clear();
+ m_truthBarcode.clear();
+
+ m_truthd0_P.clear();
+ m_truthd0_px.clear();
+ m_truthd0_py.clear();
+ m_truthd0_pz.clear();
+ m_truthd0_IsFiducial.clear();
+ m_truthd0_charge.clear();
+
+ m_truthd1_P.clear();
+ m_truthd1_px.clear();
+ m_truthd1_py.clear();
+ m_truthd1_pz.clear();
+ m_truthd1_IsFiducial.clear();
+ m_truthd1_charge.clear();
}
diff --git a/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.h b/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.h
index d0700d3a5bf4b4e891a55a7b04b102de07117f79..3495c0ed4f9b1b96aee1e364c214336ef3eb328e 100644
--- a/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.h
+++ b/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.h
@@ -12,6 +12,11 @@ Copyright (C) 2022 CERN for the benefit of the FASER collaboration
#include "xAODFaserWaveform/WaveformHit.h"
#include "xAODFaserWaveform/WaveformHitContainer.h"
+#include "FaserActsKalmanFilter/IFiducialParticleTool.h"
+#include "FaserActsKalmanFilter/ITrackTruthMatchingTool.h"
+#include "xAODTruth/TruthParticleContainer.h"
+#include "xAODTruth/TruthEventContainer.h"
+
#include <vector>
#include "TTree.h"
#include <TH1.h>
@@ -75,6 +80,10 @@ namespace Tracker
SG::ReadHandleKey<xAOD::WaveformHitContainer> m_preshowerContainer { this, "PreshowerContainer", "PreshowerWaveformHits", "Preshower hit container name" };
SG::ReadHandleKey<xAOD::WaveformHitContainer> m_ecalContainer { this, "EcalContainer", "CaloWaveformHits", "Ecal hit container name" };
+ ToolHandle<IFiducialParticleTool> m_fiducialParticleTool {this, "FiducialParticleTool", "FiducialParticleTool"};
+ ToolHandle<ITrackTruthMatchingTool> m_trackTruthMatchingTool {this, "TrackTruthMatchingTool", "TrackTruthMatchingTool"};
+ SG::ReadHandleKey<xAOD::TruthParticleContainer> m_truthParticleContainer { this, "ParticleContainer", "TruthParticles", "Truth particle container name." };
+
mutable std::atomic<size_t> m_numberOfEvents{0};
mutable std::atomic<size_t> m_numberOfPositiveTracks{0};
mutable std::atomic<size_t> m_numberOfNegativeTracks{0};
@@ -86,41 +95,39 @@ namespace Tracker
mutable std::atomic<size_t> m_numberOfNoUpstreamTracks{0};
mutable TTree* m_tree;
- mutable std::vector<std::vector<double>> m_tracks={{},{},{},{},{},{},{}}; //charge, x, y, z, px, py, pz, bS0,bS1,bS2,bS3
- mutable int m_trackCount=0;
-
- //POSITIVE
- mutable std::vector<std::vector<int>> m_playerMap={{},{},{},{}};
- mutable std::set<int> m_pstationMap;
- mutable std::vector<double> m_ppos_s1l0 = {-1,-1,-1}; // {x,y,z}
- mutable std::vector<double> m_pmom;
- mutable float m_pChi2;
- mutable double m_pDoF;
- mutable int m_pnLayers;
- mutable int m_pnLayersIFT;
- mutable int m_playersHit;
- mutable int m_piftHit;
-
- //NEGATIVE
- mutable std::vector<std::vector<int>> m_nlayerMap={{},{},{},{}};
- mutable std::set<int> m_nstationMap;
- mutable std::vector<double> m_npos_s1l0 = {-1,-1,-1}; // {x,y,z}
- mutable std::vector<double> m_nmom;
- mutable float m_nChi2;
- mutable double m_nDoF;
- mutable int m_nlayersHit;
- mutable int m_niftHit;
-
- //CALORIMETER
- //StringProperty m_CaloConfig { this, "CaloConfig", "Low_gain", "Configuration found at http://aagaard.web.cern.ch/aagaard/FASERruns.html (spaces replaced with '_')" };
- mutable float m_Calo0_Edep;
- mutable float m_Calo1_Edep;
- mutable float m_Calo2_Edep;
- mutable float m_Calo3_Edep;
- mutable float m_Calo_Total_Edep;
- mutable float m_MIP_sim_Edep_calo;
- mutable float m_MIP_sim_Edep_preshower;
+ //TRACK COLLECTION
+ mutable int m_trackCount=0;
+
+ mutable std::vector<double> m_tracksX;
+ mutable std::vector<double> m_tracksY;
+ mutable std::vector<double> m_tracksZ;
+ mutable std::vector<double> m_tracksPx;
+ mutable std::vector<double> m_tracksPy;
+ mutable std::vector<double> m_tracksPz;
+ mutable std::vector<double> m_tracksPmag;
+
+ mutable std::vector<double> m_tracksCharge;
+ mutable std::vector<double> m_sigmaQoverP;
+ mutable std::vector<double> m_tracksChi2;
+ mutable std::vector<double> m_tracksDoF;
+
+ mutable std::vector<double> m_tracksTruthBarcode;
+ mutable std::vector<double> m_tracksTruthPdg;
+ mutable std::vector<double> m_tracksTruthCharge;
+ mutable std::vector<double> m_tracksTruthPmag;
+ mutable std::vector<bool> m_tracksIsFiducial;
+
+ mutable std::vector<int> m_tracksnLayersHit;
+ mutable std::vector<int> m_tracksLayersIFT;
+ mutable std::vector<int> m_tracksIFThit;
+
+ mutable std::vector<int> m_trackPositiveIndicies;
+ mutable std::vector<int> m_trackNegativeIndicies;
+ mutable int m_trackD2i; //which track is daughter; -1 if not matched
+ mutable int m_trackD3i;
+
+ mutable std::vector<std::vector<int>> m_layerMap={{},{},{},{}};
//WAVEFORM
mutable float m_wave_localtime[15];
@@ -133,8 +140,58 @@ namespace Tracker
mutable float m_wave_baseline_rms[15];
mutable unsigned int m_wave_status[15];
- //EXTRAPOLATION
+ //TRUTH
+ DoubleProperty m_minMomentum { this, "MinMomentum", 50000.0, "Write out all truth particles with a momentum larger MinMomentum"};
+
mutable std::vector<double> m_vertexPos = {-1,-1,-1};
+ mutable std::vector<bool> m_isFiducial; // track is fiducial if there are simulated hits for stations 1 - 3 and the distance from the center is smaller than 100 mm
+
+ mutable std::vector<int> m_truthParticleBarcode; // vector of barcodes of all truth particles with a momentum larger 50 GeV
+ mutable std::vector<int> m_truthParticleMatchedTracks; // vector of number of tracks to which a truth particle is matched to
+ mutable std::vector<bool> m_truthParticleIsFiducial; // vector of boolean showing whether a truth particle is fiducial
+
+ mutable std::vector<int> m_truthBarcode;
+ mutable std::vector<int> m_truthPdg;
+ mutable std::vector<double> m_truthPmag;
+ mutable std::vector<double> m_truthCharge;
+
+ mutable std::vector<double> m_truthM_P;
+
+ mutable std::vector<double> m_truthM_px;
+ mutable std::vector<double> m_truthM_py;
+ mutable std::vector<double> m_truthM_pz;
+
+ mutable std::vector<double> m_truthM_x;
+ mutable std::vector<double> m_truthM_y;
+ mutable std::vector<double> m_truthM_z;
+
+
+ mutable std::vector<double> m_truthd0_P;
+
+ mutable std::vector<double> m_truthd0_px;
+ mutable std::vector<double> m_truthd0_py;
+ mutable std::vector<double> m_truthd0_pz;
+
+ mutable std::vector<double> m_truthd0_x;
+ mutable std::vector<double> m_truthd0_y;
+ mutable std::vector<double> m_truthd0_z;
+
+ mutable std::vector<double> m_truthd0_IsFiducial;
+ mutable std::vector<double> m_truthd0_charge;
+
+ mutable std::vector<double> m_truthd1_P;
+
+ mutable std::vector<double> m_truthd1_px;
+ mutable std::vector<double> m_truthd1_py;
+ mutable std::vector<double> m_truthd1_pz;
+
+ mutable std::vector<double> m_truthd1_x;
+ mutable std::vector<double> m_truthd1_y;
+ mutable std::vector<double> m_truthd1_z;
+
+ mutable std::vector<double> m_truthd1_IsFiducial;
+ mutable std::vector<double> m_truthd1_charge;
+
};