diff --git a/.gitmodules b/.gitmodules
index 82b8e1319b74172b9299e014eceb3803ae142cde..5f4e81b19df741d282252ef76cb5841311e2477a 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -1,3 +1,3 @@
[submodule "faser-common"]
path = faser-common
- url = https://gitlab.cern.ch/faser/faser-common.git
+ url = ../../faser/faser-common.git
diff --git a/PhysicsAnalysis/NtupleDumper/scripts/analyzeRun2.py b/PhysicsAnalysis/NtupleDumper/scripts/analyzeRun2.py
new file mode 100644
index 0000000000000000000000000000000000000000..011cdbe9f8ae24e3503e34e314e14801ed7092fc
--- /dev/null
+++ b/PhysicsAnalysis/NtupleDumper/scripts/analyzeRun2.py
@@ -0,0 +1,137 @@
+#!/usr/bin/env python3
+
+"""
+ Copyright (C) 2002-2022 CERN for the benefit of the FASER collaboration
+"""
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg
+
+
+def NtupleDumperAlgCfg(flags, OutName, **kwargs):
+ # Initialize GeoModel
+ from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
+ acc = FaserGeometryCfg(flags)
+
+ acc.merge(MagneticFieldSvcCfg(flags))
+ # acc.merge(FaserActsTrackingGeometrySvcCfg(flags))
+ # acc.merge(FaserActsAlignmentCondAlgCfg(flags))
+
+ actsExtrapolationTool = CompFactory.FaserActsExtrapolationTool("FaserActsExtrapolationTool")
+ actsExtrapolationTool.MaxSteps = 10000
+ actsExtrapolationTool.TrackingGeometryTool = CompFactory.FaserActsTrackingGeometryTool("TrackingGeometryTool")
+
+ NtupleDumperAlg = CompFactory.NtupleDumperAlg("NtupleDumperAlg",**kwargs)
+ NtupleDumperAlg.ExtrapolationTool = actsExtrapolationTool
+ acc.addEventAlgo(NtupleDumperAlg)
+
+ thistSvc = CompFactory.THistSvc()
+ thistSvc.Output += [f"HIST2 DATAFILE='{OutName}' OPT='RECREATE'"]
+ acc.addService(thistSvc)
+
+ return acc
+
+if __name__ == "__main__":
+
+ import glob
+ import sys
+ import ROOT
+
+ runno=int(sys.argv[1])
+ num=int(sys.argv[2])
+ filesPerJob=int(sys.argv[3])
+ run_config=str(sys.argv[4])
+
+ ptag="p0008"
+
+ from AthenaCommon.Logging import log, logging
+ from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+ from AthenaCommon.Configurable import Configurable
+ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+ from AthenaConfiguration.TestDefaults import defaultTestFiles
+ from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+ from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+ # from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+ # Set up logging and new style config
+ log.setLevel(INFO)
+ Configurable.configurableRun3Behavior = True
+
+ #fName = dataDir+"FaserMC-MDC_FS_Aee_100MeV_1Em5_shift-110001-00000-00009-s0007-r0009-xAOD.root"
+
+
+ dataDir=f"/eos/experiment/faser/sim/mdc/foresee/{ptag}/rec/r0009/"
+ files=sorted(glob.glob(f"{dataDir}/FaserMC-MDC*"))
+ fileListInitial=files[num*filesPerJob:(num+1)*filesPerJob]
+ fileList=[]
+ for fName in fileListInitial:
+ try:
+ fh=ROOT.TFile(fName)
+ fileList.append(fName)
+ 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"Got {len(fileList)} files out of {len(fileListInitial)}")
+
+ outName=f"Data-tuple-{run_config}-{runno:06d}-{num:05d}-{filesPerJob}.root"
+
+ log.info(f"Analyzing Run {runno} files {num*filesPerJob} to {(num+1)*filesPerJob} (num={num})")
+ log.info(f"Got {len(fileList)} files out of {len(fileListInitial)}")
+
+ outName=f"Data-tuple-{run_config}-{runno:06d}-{num:05d}-{filesPerJob}.root"
+
+ # outName = "MC_MDC_darkphoton_tuple.root"
+
+ # Configure
+ ConfigFlags.Input.Files = fileList
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # 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 = True # Needed to bypass autoconfig
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # FASER geometry
+ ConfigFlags.Common.isOnline = False
+ ConfigFlags.GeoModel.Align.Dynamic = False
+ ConfigFlags.Beam.NumberOfCollisions = 0.
+
+ ConfigFlags.Detector.GeometryFaserSCT = True
+
+ ConfigFlags.lock()
+
+ # Core components
+ acc = MainServicesCfg(ConfigFlags)
+ acc.merge(PoolReadCfg(ConfigFlags))
+
+ # algorithm
+ acc.merge(NtupleDumperAlgCfg(ConfigFlags, outName, UseFlukaWeights=True, CaloConfig=run_config))
+
+ AthenaEventLoopMgr = CompFactory.AthenaEventLoopMgr()
+ AthenaEventLoopMgr.EventPrintoutInterval=1000
+ 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
+
+ # Timing
+ #acc.merge(MergeRecoTimingObjCfg(ConfigFlags))
+
+ # Dump config
+ # logging.getLogger('forcomps').setLevel(VERBOSE)
+ # acc.foreach_component("*").OutputLevel = VERBOSE
+ # acc.foreach_component("*ClassID*").OutputLevel = INFO
+ # acc.getCondAlgo("FaserSCT_AlignCondAlg").OutputLevel = VERBOSE
+ # acc.getCondAlgo("FaserSCT_DetectorElementCondAlg").OutputLevel = VERBOSE
+ # acc.getService("StoreGateSvc").Dump = True
+ # acc.getService("ConditionStore").Dump = True
+ # acc.printConfig(withDetails=True)
+ # ConfigFlags.dump()
+
+ # Execute and finish
+ sc = acc.run(maxEvents=-1)
+
+ # Success should be 0
+ sys.exit(not sc.isSuccess())
diff --git a/Tracker/TrackerRecAlgs/PairVertex/CMakeLists.txt b/Tracker/TrackerRecAlgs/PairVertex/CMakeLists.txt
index 85053a1f49bd551ea136c5cd61f6cee40543ed26..ff270bb69c9defc70b31e5e50f8dbc93ffa3992f 100644
--- a/Tracker/TrackerRecAlgs/PairVertex/CMakeLists.txt
+++ b/Tracker/TrackerRecAlgs/PairVertex/CMakeLists.txt
@@ -5,7 +5,10 @@ atlas_add_component(
src/PairVertexAlg.h
src/PairVertexAlg.cxx
src/component/PairVertex_entries.cxx
- LINK_LIBRARIES AthenaBaseComps StoreGateLib GeneratorObjects FaserActsGeometryLib TrackerSimEvent TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives
+ LINK_LIBRARIES AthenaBaseComps StoreGateLib GeneratorObjects
+ FaserActsGeometryLib TrackerSimEvent TrackerIdentifier
+ TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack
+ xAODFaserWaveform ScintIdentifier FaserCaloIdentifier
)
atlas_install_python_modules(python/*.py)
diff --git a/Tracker/TrackerRecAlgs/PairVertex/python/PairVertexAlgConfig.py b/Tracker/TrackerRecAlgs/PairVertex/python/PairVertexAlgConfig.py
index f98af1629a39bfc590ef3ab0929da4b849ad5ce6..4243fe99559a9ff4f16753f0e53c51c9d8ea1d2a 100644
--- a/Tracker/TrackerRecAlgs/PairVertex/python/PairVertexAlgConfig.py
+++ b/Tracker/TrackerRecAlgs/PairVertex/python/PairVertexAlgConfig.py
@@ -7,19 +7,38 @@ from AthenaConfiguration.ComponentFactory import CompFactory
from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg
from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg
+#Command line config
+from argparse import ArgumentParser
+
+test_file = '/eos/experiment/faser/sim/mdc/foresee/110001/rec/r0009/FaserMC-MDC_FS_Aee_100MeV_1Em5_shift-110001-00000-00009-s0007-r0009-xAOD.root'
+parser=ArgumentParser(description="EventLooper")
+parser.add_argument("-i","--inputfile",type=str,default=test_file)
+parser.add_argument("-e", "--events", type=int, default=1000) # Max events processed in alg. Set to -1 to process all events.
+parser.add_argument("-r","--runnum",type=int,default=110001) # 110001-4
+ui=parser.parse_args()
+
+inputIsMC = 'FaserMC' in ui.inputfile
+
+#Get/Make input/output file names
+from ROOT import TFile
+from glob import glob
+dataDir=f"/eos/experiment/faser/sim/mdc/foresee/{ui.runnum}/rec/r0009/"
+files=sorted(glob(f"{dataDir}/Faser*"))
+name=files[0].split("r0009/")[-1].split("MDC_FS")[-1].split("shift")[0] #Generate output file name from input file name
+
+#-------------------
+
+
def PairVertexAlgCfg(flags, **kwargs):
acc = FaserSCT_GeometryCfg(flags)
acc.merge(MagneticFieldSvcCfg(flags))
- # acc.merge(FaserActsTrackingGeometrySvcCfg(flags))
- # acc.merge(FaserActsAlignmentCondAlgCfg(flags))
-
- actsExtrapolationTool = CompFactory.FaserActsExtrapolationTool("FaserActsExtrapolationTool")
- actsExtrapolationTool.MaxSteps = 1000
- actsExtrapolationTool.TrackingGeometryTool = CompFactory.FaserActsTrackingGeometryTool("TrackingGeometryTool")
-
PairVertexAlg = CompFactory.Tracker.PairVertexAlg("PairVertexAlg",**kwargs)
- PairVertexAlg.ExtrapolationTool = actsExtrapolationTool
acc.addEventAlgo(PairVertexAlg)
+
+ #Hist output
+ thistSvc = CompFactory.THistSvc()
+ thistSvc.Output = [f"HIST DATAFILE='Ntuple{name}.root' OPT='RECREATE'"]
+ acc.addService(thistSvc)
return acc
if __name__ == "__main__":
@@ -35,17 +54,15 @@ if __name__ == "__main__":
from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
# Set up logging and new style config
- log.setLevel(DEBUG)
+ log.setLevel(INFO)
Configurable.configurableRun3Behavior = True
# Configure
- ConfigFlags.Input.Files = [
- '/eos/experiment/faser/sim/mdc/foresee/110004/rec/r0007/FaserMC-MDC_FS_Amm_316MeV_2Em6-110004-00000-s0003-r0007-xAOD.root'
- ]
+ ConfigFlags.Input.Files = files # [ui.inputfile]
ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # 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 = True # 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
@@ -62,29 +79,16 @@ if __name__ == "__main__":
# algorithm
acc.merge(PairVertexAlgCfg(ConfigFlags))
- # # 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
-
- # Timing
- #acc.merge(MergeRecoTimingObjCfg(ConfigFlags))
-
- # Dump config
- # logging.getLogger('forcomps').setLevel(VERBOSE)
- # acc.foreach_component("*").OutputLevel = VERBOSE
- # acc.foreach_component("*ClassID*").OutputLevel = INFO
- # acc.getCondAlgo("FaserSCT_AlignCondAlg").OutputLevel = VERBOSE
- # acc.getCondAlgo("FaserSCT_DetectorElementCondAlg").OutputLevel = VERBOSE
- # acc.getService("StoreGateSvc").Dump = True
- # acc.getService("ConditionStore").Dump = True
- # acc.printConfig(withDetails=True)
- # ConfigFlags.dump()
+ # Hack to avoid problem with our use of MC databases when isMC = False
+ if not inputIsMC:
+ replicaSvc = acc.getService("DBReplicaSvc")
+ replicaSvc.COOLSQLiteVetoPattern = ""
+ replicaSvc.UseCOOLSQLite = True
+ replicaSvc.UseCOOLFrontier = False
+ replicaSvc.UseGeomSQLite = True
# Execute and finish
- sc = acc.run(maxEvents=-1)
+ sc = acc.run(maxEvents=ui.events)
# Success should be 0
- sys.exit(not sc.isSuccess())
\ No newline at end of file
+ sys.exit(not sc.isSuccess())
\ No newline at end of file
diff --git a/Tracker/TrackerRecAlgs/PairVertex/scripts/PairVertexReader.py b/Tracker/TrackerRecAlgs/PairVertex/scripts/PairVertexReader.py
new file mode 100644
index 0000000000000000000000000000000000000000..2468be886114912fa36555c88abaad99e8e77076
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/PairVertex/scripts/PairVertexReader.py
@@ -0,0 +1,587 @@
+#!/usr/bin/env python
+
+# 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)
+
+ 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)
+
+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=[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
+
+ 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):
+
+ 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]]
+
+ 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 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()
+
+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)
+
+MainArrDict={
+ 'CaloSum':{
+ 'array':[],
+ 'lambda': SumCalo,
+ 'plotterfunc': configBasicHist,
+ 'plotkwargs' : {'xscale':'log'}
+ },
+ 'TracksPerEvent':{
+ 'array':[],
+ 'lambda': lambda event: len(event.Tracks[0]),
+ 'plotterfunc': configBasicHist,
+ 'plotkwargs' : {}
+ },
+ 'TrackSeparation':{
+ 'array':[],
+ 'lambda': TrackSeparation,
+ '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)}}
+ },
+ 'PAngleY':{
+ 'array':[[],[]],
+ 'lambda': PAngleY,
+ 'plotterfunc': configMultiHist,
+ 'plotkwargs' : {'label':['Positive','Negative'],'xlabel':'Angle (mrad)','kwargs':{'range':(-35,35)}}
+ },
+ 'SeparationDistribution':{
+ 'array':[[],[]],
+ 'lambda': SeparationDistribution,
+ 'plotterfunc': configMultiHist,
+ 'plotkwargs' : {'label':['Y','X'],'kwargs':{'range':(-100,100)}}
+ },
+ 'Chi2DoF':{
+ 'array':[[],[]],
+ 'lambda': Chi2DoF,
+ 'plotterfunc': configMultiHist,
+ 'plotkwargs' : {'bins':20,'label':['Positive','Negative'],'kwargs':{'range':(0,11)}}
+ },
+ 'DoF':{
+ 'array':[[],[]],
+ 'lambda': lambda event:[event.pDoF,event.nDoF],
+ 'plotterfunc': configMultiHist,
+ 'plotkwargs' : {'bins':30,'label':['Positive','Negative'],'kwargs':{'range':(0,31)}}
+ },
+ 'LayersHit':{
+ 'array':[[],[]],
+ 'lambda': lambda event: [event.playersHit,event.nlayersHit],
+ 'plotterfunc': configMultiHist,
+ 'plotkwargs' : {'bins':5,'label':['Positive','Negative'],'kwargs':{'range':(5,10)}}
+ },
+ 'VertexZ':{
+ 'array':[],
+ 'lambda': lambda event: event.vertexPos[2],
+ 'plotterfunc': configBasicHist,
+ 'plotkwargs' : {'xlabel':'Z (mm)'}
+ },
+ 'TrueVertexRho':{
+ 'array':[],
+ 'lambda': lambda event: (event.vertexPos[0]**2)+(event.vertexPos[1]**2),
+ 'plotterfunc': configBasicHist,
+ 'plotkwargs' : {'xlabel':'Rho^2 (mm^2)'}
+ },
+ '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)'}
+ },
+ 'UpstreamHitsPositive':{
+ 'array':[[],[]],
+ 'lambda': UpstreamHitMap,
+ 'plotterfunc': configHitMap2,
+ 'plotkwargs' : {'xlim':(-100,100),'ylim':(-100,100)},
+ 'lkwargs':{'charge':1}
+ },
+ 'UpstreamHitsNegative':{
+ 'array':[[],[]],
+ 'lambda': UpstreamHitMap,
+ 'plotterfunc': configHitMap2,
+ 'plotkwargs' : {'xlim':(-100,100),'ylim':(-100,100)},
+ 'lkwargs':{'charge':-1}
+ },
+ 'UpstreamXYPositive':{
+ 'array':[[],[]],
+ 'lambda': lambda event: (event.pPos_S1L0[0],event.pPos_S1L0[1]),
+ 'plotterfunc': configMultiHist,
+ 'plotkwargs' : {'label':['Positive','Negative']},
+ },
+}
+MainArrDict.update(MakeWFDict()) #Fill with wf plots
+
+#For testing/plotting purposes, plotting specific members of MainArrDict
+TestArrDict=dict()
+testname='TrueVertexRho'
+TestArrDict[testname]=MainArrDict[testname]
+# for k in MainArrDict.keys():
+# if 'Calo' in k:
+# TestArrDict[k]=MainArrDict[k]
+
+def GeneralPlotter(arrDict,dcuts=[CUTDICT],lcuts=LAMBDACUTS,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:
+ 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
+
+ print(f"# of NaNs:{nans}")
+ for name in arrDict.keys():
+ 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()
+
diff --git a/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.cxx b/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.cxx
index 99058166c7252ddbff852b0a09ec8d49a15e54fb..8009ce9a72e0efb119b0e5cf52c52d99a2add749 100644
--- a/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.cxx
+++ b/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.cxx
@@ -3,49 +3,167 @@
#include "FaserActsGeometryInterfaces/IFaserActsTrackingGeometryTool.h"
#include "GeoPrimitives/CLHEPtoEigenConverter.h"
#include "Acts/Surfaces/PerigeeSurface.hpp"
-#include "TrackerIdentifier/FaserSCT_ID.h"
#include "TrackerReadoutGeometry/SCT_DetectorManager.h"
#include "TrackerReadoutGeometry/SiDetectorElement.h"
+#include "Identifier/Identifier.h"
+#include "TrackerIdentifier/FaserSCT_ID.h"
+#include "ScintIdentifier/VetoNuID.h"
+#include "ScintIdentifier/VetoID.h"
+#include "ScintIdentifier/TriggerID.h"
+#include "ScintIdentifier/PreshowerID.h"
+#include "FaserCaloIdentifier/EcalID.h"
+#include "TrackerRIO_OnTrack/FaserSCT_ClusterOnTrack.h"
#include <cmath>
+#include <TTree.h>
+#include <TBranch.h>
+
+#include <typeinfo>
+
+
namespace Tracker {
+
PairVertexAlg::PairVertexAlg(const std::string &name, ISvcLocator *pSvcLocator)
- : AthReentrantAlgorithm(name, pSvcLocator) {}
+ : AthReentrantAlgorithm(name, pSvcLocator), AthHistogramming(name), m_idHelper(nullptr), m_histSvc("THistSvc/THistSvc", name) {}
+
+void PairVertexAlg::addBranch(const std::string &name,
+ float* var) {
+ m_tree->Branch(name.c_str(),var,(name+"/F").c_str());
+}
+void PairVertexAlg::addBranch(const std::string &name,
+ unsigned int* var) {
+ m_tree->Branch(name.c_str(),var,(name+"/I").c_str());
+}
+void PairVertexAlg::addWaveBranches(const std::string &name,
+ int nchannels,
+ int first) {
+ for(int ch=0;ch<nchannels;ch++) {
+ std::string base=name+std::to_string(ch)+"_";
+ addBranch(base+"time",&m_wave_localtime[first]);
+ addBranch(base+"peak",&m_wave_peak[first]);
+ addBranch(base+"width",&m_wave_width[first]);
+ addBranch(base+"integral",&m_wave_integral[first]);
+ addBranch(base+"raw_peak",&m_wave_raw_peak[first]);
+ addBranch(base+"raw_charge",&m_wave_raw_integral[first]);
+ addBranch(base+"baseline",&m_wave_baseline_mean[first]);
+ addBranch(base+"baseline_rms",&m_wave_baseline_rms[first]);
+ addBranch(base+"status",&m_wave_status[first]);
+ first++;
+ }
+}
+
+void PairVertexAlg::FillWaveBranches(const xAOD::WaveformHitContainer &wave) const {
+ for (auto hit : wave) {
+ if (waveformHitOK(hit) && (hit->hit_status()&2)==0) { // dont store secoondary hits as they can overwrite the primary hit
+ int ch=hit->channel();
+ //m_wave_localtime[ch]=hit->localtime()+m_clock_phase;
+ m_wave_peak[ch]=hit->peak();
+ m_wave_width[ch]=hit->width();
+ m_wave_integral[ch]=hit->integral();
+
+ m_wave_raw_peak[ch]=hit->raw_peak();
+ m_wave_raw_integral[ch]=hit->raw_integral();
+ m_wave_baseline_mean[ch]=hit->baseline_mean();
+ m_wave_baseline_rms[ch]=hit->baseline_rms();
+ m_wave_status[ch]=hit->hit_status();
+ }
+ }
+}
+
+bool PairVertexAlg::waveformHitOK(const xAOD::WaveformHit* hit) const
+{
+ if (hit->status_bit(xAOD::WaveformStatus::THRESHOLD_FAILED) || hit->status_bit(xAOD::WaveformStatus::SECONDARY)) return false;
+ return true;
+}
+
StatusCode PairVertexAlg::initialize()
{
ATH_CHECK(m_mcEventCollectionKey.initialize());
ATH_CHECK(m_trackCollectionKey.initialize());
- ATH_CHECK(m_extrapolationTool.retrieve());
+ ATH_CHECK(m_vetoNuContainer.initialize());
+ ATH_CHECK(m_vetoContainer.initialize());
+ ATH_CHECK(m_triggerContainer.initialize());
+ ATH_CHECK(m_preshowerContainer.initialize());
+ ATH_CHECK(m_ecalContainer.initialize());
+
ATH_CHECK(detStore()->retrieve(m_idHelper, "FaserSCT_ID"));
ATH_CHECK(detStore()->retrieve(m_detMgr, "SCT"));
+ ATH_CHECK(detStore()->retrieve(m_vetoNuHelper, "VetoNuID"));
+ ATH_CHECK(detStore()->retrieve(m_vetoHelper, "VetoID"));
+ ATH_CHECK(detStore()->retrieve(m_triggerHelper, "TriggerID"));
+ ATH_CHECK(detStore()->retrieve(m_preshowerHelper, "PreshowerID"));
+ 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);
+
+ //WAVEFORM
+ addWaveBranches("VetoSt1",2,6);
+ addWaveBranches("VetoSt2",1,14);
+ addWaveBranches("VetoNu",2,4);
+ addWaveBranches("Timing",4,8);
+ addWaveBranches("Preshower",2,12);
+ addWaveBranches("Calo",4,0);
+
+ //TRUTH
+ m_tree->Branch("vertexPos", &m_vertexPos); // [x,y,z]
+
+ ATH_CHECK(histSvc()->regTree("/HIST/events", m_tree));
return StatusCode::SUCCESS;
}
StatusCode PairVertexAlg::execute(const EventContext &ctx) const
{
+ clearTree();
m_numberOfEvents++;
- const Acts::GeometryContext gctx =
- m_extrapolationTool->trackingGeometryTool()->getNominalGeometryContext().context();
+ //WAVEFORM & CALO
+ SG::ReadHandle<xAOD::WaveformHitContainer> vetoNuContainer { m_vetoNuContainer, ctx };
+ ATH_CHECK(vetoNuContainer.isValid());
+ SG::ReadHandle<xAOD::WaveformHitContainer> vetoContainer { m_vetoContainer, ctx };
+ ATH_CHECK(vetoContainer.isValid());
+ SG::ReadHandle<xAOD::WaveformHitContainer> triggerContainer { m_triggerContainer, ctx };
+ ATH_CHECK(triggerContainer.isValid());
+ SG::ReadHandle<xAOD::WaveformHitContainer> preshowerContainer { m_preshowerContainer, ctx };
+ ATH_CHECK(preshowerContainer.isValid());
+ SG::ReadHandle<xAOD::WaveformHitContainer> ecalContainer { m_ecalContainer, ctx };
+ ATH_CHECK(ecalContainer.isValid());
- std::vector<double> z_positions {};
-
- SG::ReadHandle<McEventCollection> mcEvents {m_mcEventCollectionKey, ctx};
- ATH_CHECK(mcEvents.isValid());
- if (mcEvents->size() != 1)
- {
- ATH_MSG_ERROR("There should be exactly one event in the McEventCollection.");
- return StatusCode::FAILURE;
- }
+ FillWaveBranches(*vetoContainer);
+ FillWaveBranches(*triggerContainer);
+ FillWaveBranches(*preshowerContainer);
+ FillWaveBranches(*ecalContainer);
+ FillWaveBranches(*vetoNuContainer);
+ // TRACK COLLECTION
SG::ReadHandle<TrackCollection> tracks { m_trackCollectionKey, ctx };
ATH_CHECK(tracks.isValid());
-
- const Trk::TrackParameters* positive {nullptr};
- const Trk::TrackParameters* negative {nullptr};
-
for (auto trk : *tracks)
{
const Trk::TrackParameters* upstream {nullptr};
@@ -56,6 +174,8 @@ StatusCode PairVertexAlg::execute(const EventContext &ctx) const
continue;
}
auto parameters = trk->trackParameters();
+
+ //Error Catching
if (parameters == nullptr || parameters->empty())
{
m_numberOfNoParameterTracks++;
@@ -84,69 +204,112 @@ 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.
}
- auto momentum = upstream->momentum();
- double charge = upstream->charge();
- if (charge > 0 || momentum.mag() > positive->momentum().mag())
- {
- positive = upstream;
+
+
+ // Check for hit in the three downstream stations
+ 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);
+ if (cluster != nullptr) {
+ Identifier id = cluster->identify();
+ int station = m_idHelper->station(id);
+ int layer = m_idHelper->layer(id);
+ stationMap.emplace(station);
+ layerMap[station][layer]+=1;
+ }
}
- else if (charge < 0 || momentum.mag() > negative->momentum().mag())
- {
- negative = upstream;
+
+ 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;
+ const Trk::TrackParameters* upstreamParameters = trk->trackParameters()->front();
+ const Trk::TrackParameters* downstreamParameters = trk->trackParameters()->back();
+
+
+ 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_trackCount;
+
}
- if (positive != nullptr) m_numberOfPositiveTracks++;
- if (negative != nullptr) m_numberOfNegativeTracks++;
-
- if (positive != nullptr && negative != nullptr) m_numberOfOppositeSignPairs++;
-
- // for (const HepMC::GenParticle* particle : mcEvents->front()->particle_range())
- // {
- // if ((std::abs(particle->pdg_id()) != 13)) continue;
- // const HepMC::FourVector& vertex = particle->production_vertex()->position();
- // if (vertex.z() > 0) continue;
- // const HepMC::FourVector& momentum = particle->momentum();
- // double phi = momentum.phi();
- // double theta = momentum.theta();
- // double charge = particle->pdg_id() > 0 ? -1 : 1;
- // double abs_momentum = momentum.rho() * m_MeV2GeV;
- // double qop = charge / abs_momentum;
- // // The coordinate system of the Acts::PlaneSurface is defined as
- // // T = Z = normal, U = X x T = -Y, V = T x U = x
- // // Acts::BoundVector pars;
- // // pars << -vertex.y(), vertex.x(), phi, theta, qop, vertex.t();
- // Acts::BoundVector pars = Acts::BoundVector::Zero();
- // pars[Acts::eBoundLoc0] = -vertex.y();
- // pars[Acts::eBoundLoc1] = vertex.x();
- // pars[Acts::eBoundPhi] = phi;
- // pars[Acts::eBoundTheta] = theta;
- // pars[Acts::eBoundQOverP] = qop;
- // pars[Acts::eBoundTime] = vertex.t();
-
- // auto startSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
- // Acts::Vector3(0, 0, vertex.z()), Acts::Vector3(0, 0, 1));
- // auto targetSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
- // Acts::Vector3(0, 0, z_mean), Acts::Vector3(0, 0, 1));
- // Acts::BoundTrackParameters startParameters(
- // std::move(startSurface), pars, charge);
- // ATH_MSG_DEBUG("vertex: " << vertex.x() << ", " << vertex.y() << ", " << vertex.z());
- // ATH_MSG_DEBUG("vertex momentum: " << momentum.x() * m_MeV2GeV << ", " << momentum.y() * m_MeV2GeV << ", " << momentum.z() * m_MeV2GeV);
- // std::unique_ptr<const Acts::BoundTrackParameters> targetParameters =
- // m_extrapolationTool->propagate(ctx, startParameters, *targetSurface);
- // if (targetParameters)
- // {
- // Acts::Vector3 targetPosition = targetParameters->position(gctx);
- // Acts::Vector3 targetMomentum = targetParameters->momentum();
- // ATH_MSG_DEBUG("vertex: " << vertex.x() << ", " << vertex.y() << ", " << vertex.z());
- // ATH_MSG_DEBUG("origin: " << targetPosition.x() << ", " << targetPosition.y() << ", " << targetPosition.z());
- // ATH_MSG_DEBUG("vertex momentum: " << momentum.x() * m_MeV2GeV << ", " << momentum.y() * m_MeV2GeV << ", " << momentum.z() * m_MeV2GeV);
- // ATH_MSG_DEBUG("origin momentum: " << targetMomentum.x() << ", " << targetMomentum.y() << ", " << targetMomentum.z());
- // }
- // }
- return StatusCode::SUCCESS;
+ // MC TRUTH
+ SG::ReadHandle<McEventCollection> mcEvents {m_mcEventCollectionKey, ctx};
+ ATH_CHECK(mcEvents.isValid());
+ if (mcEvents->size() != 1)
+ {
+ ATH_MSG_ERROR("There should be exactly one event in the McEventCollection.");
+ return StatusCode::FAILURE;
+ }
+ const HepMC::GenEvent* theEvent = mcEvents->at(0);
+ auto pVertices = theEvent->vertices_begin();
+ auto pVerticesEnd = theEvent->vertices_end();
+
+ for (; pVertices != pVerticesEnd; ++pVertices)
+ {
+ const HepMC::GenVertex* v = *pVertices;
+ m_vertexPos[0]=v->position().x();
+ m_vertexPos[1]=v->position().y();
+ m_vertexPos[2]=v->position().z();
+ break;
+ }
+
+
+ m_tree->Fill();
+
+ return StatusCode::SUCCESS;
}
StatusCode PairVertexAlg::finalize()
@@ -160,4 +323,46 @@ StatusCode PairVertexAlg::finalize()
return StatusCode::SUCCESS;
}
-} // Tracker
+void PairVertexAlg::clearTree() const
+{
+ for (int i = 0; i < 7; i++) {m_tracks[i].clear();}
+ m_trackCount=0;
+
+ 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;
+ m_wave_peak[ii]=0;
+ m_wave_width[ii]=0;
+ m_wave_integral[ii]=0;
+
+ m_wave_raw_peak[ii]=0;
+ m_wave_raw_integral[ii]=0;
+ m_wave_baseline_mean[ii]=0;
+ m_wave_baseline_rms[ii]=0;
+ m_wave_status[ii]=0;
+ }
+ m_vertexPos= {-1,-1,-1};
+
+}
+
+} // end Tracker
+
+
+
diff --git a/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.h b/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.h
index d15fddf6e288828d143218fa6d2f686a50076ba0..d0700d3a5bf4b4e891a55a7b04b102de07117f79 100644
--- a/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.h
+++ b/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.h
@@ -2,15 +2,34 @@
Copyright (C) 2022 CERN for the benefit of the FASER collaboration
*/
+
#pragma once
#include "GeoPrimitives/GeoPrimitives.h"
#include "AthenaBaseComps/AthReentrantAlgorithm.h"
-#include "FaserActsGeometryInterfaces/IFaserActsExtrapolationTool.h"
#include "GeneratorObjects/McEventCollection.h"
#include "TrkTrack/TrackCollection.h"
+#include "xAODFaserWaveform/WaveformHit.h"
+#include "xAODFaserWaveform/WaveformHitContainer.h"
+
+#include <vector>
+#include "TTree.h"
+#include <TH1.h>
+#include <TProfile.h>
+#include "AthenaBaseComps/AthHistogramming.h"
+#include "StoreGate/ReadHandleKey.h"
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/ToolHandle.h"
+
class FaserSCT_ID;
+class VetoNuID;
+class TTree;
+class VetoID;
+class TriggerID;
+class PreshowerID;
+class EcalID;
+
namespace TrackerDD
{
class SCT_DetectorManager;
@@ -18,7 +37,7 @@ namespace TrackerDD
namespace Tracker
{
- class PairVertexAlg : public AthReentrantAlgorithm
+ class PairVertexAlg : public AthReentrantAlgorithm,AthHistogramming
{
public:
PairVertexAlg(const std::string& name, ISvcLocator* pSvcLocator);
@@ -26,16 +45,35 @@ namespace Tracker
virtual StatusCode initialize() override;
virtual StatusCode execute(const EventContext& ctx) const override;
virtual StatusCode finalize() override;
-
+ const ServiceHandle <ITHistSvc> &histSvc() const;
+
private:
+
+ void clearTree() const;
+ void addBranch(const std::string &name,float* var);
+ void addBranch(const std::string &name,unsigned int* var);
+ void addWaveBranches(const std::string &name, int nchannels, int first);
+ void FillWaveBranches(const xAOD::WaveformHitContainer &wave) const;
+ bool waveformHitOK(const xAOD::WaveformHit* hit) const;
+
+ ServiceHandle<ITHistSvc> m_histSvc;
+
double m_MeV2GeV = 1e-3;
const FaserSCT_ID* m_idHelper {nullptr};
+ const VetoNuID* m_vetoNuHelper;
+ const VetoID* m_vetoHelper;
+ const TriggerID* m_triggerHelper;
+ const PreshowerID* m_preshowerHelper;
+ const EcalID* m_ecalHelper;
+
const TrackerDD::SCT_DetectorManager* m_detMgr {nullptr};
SG::ReadHandleKey<McEventCollection> m_mcEventCollectionKey { this, "McEventCollection", "TruthEvent" };
-
SG::ReadHandleKey<TrackCollection> m_trackCollectionKey { this, "TrackCollection", "CKFTrackCollection" };
-
- ToolHandle<IFaserActsExtrapolationTool> m_extrapolationTool { this, "ExtrapolationTool", "FaserActsExtrapolationTool" };
+ SG::ReadHandleKey<xAOD::WaveformHitContainer> m_vetoNuContainer { this, "VetoNuContainer", "VetoNuWaveformHits", "VetoNu hit container name" };
+ SG::ReadHandleKey<xAOD::WaveformHitContainer> m_vetoContainer { this, "VetoContainer", "VetoWaveformHits", "Veto hit container name" };
+ SG::ReadHandleKey<xAOD::WaveformHitContainer> m_triggerContainer { this, "TriggerContainer", "TriggerWaveformHits", "Trigger hit container name" };
+ 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" };
mutable std::atomic<size_t> m_numberOfEvents{0};
mutable std::atomic<size_t> m_numberOfPositiveTracks{0};
@@ -46,6 +84,62 @@ namespace Tracker
mutable std::atomic<size_t> m_numberOfNullParameters{0};
mutable std::atomic<size_t> m_numberOfParametersWithoutSurface{0};
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;
+
+
+ //WAVEFORM
+ mutable float m_wave_localtime[15];
+ mutable float m_wave_peak[15];
+ mutable float m_wave_width[15];
+ mutable float m_wave_integral[15];
+ mutable float m_wave_raw_peak[15];
+ mutable float m_wave_raw_integral[15];
+ mutable float m_wave_baseline_mean[15];
+ mutable float m_wave_baseline_rms[15];
+ mutable unsigned int m_wave_status[15];
+
+ //EXTRAPOLATION
+ mutable std::vector<double> m_vertexPos = {-1,-1,-1};
+
};
+
+ inline const ServiceHandle<ITHistSvc>& PairVertexAlg::histSvc() const {return m_histSvc;}
+
+ }
+
-}
\ No newline at end of file