Vertexing

.gitmodules

 [submodule "faser-common"]
 	   path = faser-common
-	   url = https://gitlab.cern.ch/faser/faser-common.git
+	   url = ../../faser/faser-common.git

PhysicsAnalysis/NtupleDumper/scripts/analyzeRun2.py

+#!/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())    

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)

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())    
+    sys.exit(not sc.isSuccess())
\ No newline at end of file

Tracker/TrackerRecAlgs/PairVertex/scripts/PairVertexReader.py

+#!/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()
+   

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
+
+
+

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