diff --git a/Tracker/TrackerRecAlgs/PairVertex/scripts/NtupleReader.py b/Tracker/TrackerRecAlgs/PairVertex/scripts/NtupleReader.py
index cb7501eb978684994dacdbe47e52309c65186568..6b38ab88087a105b0cc872b84dba2aa0a8c8b005 100644
--- a/Tracker/TrackerRecAlgs/PairVertex/scripts/NtupleReader.py
+++ b/Tracker/TrackerRecAlgs/PairVertex/scripts/NtupleReader.py
@@ -29,20 +29,19 @@ BASELINE={
'Calo_total_E_EM':[500,None,None,None]
}
-DATAONLY={'distanceToCollidingBCID ':[None,None,0,None],}
-
PHOTONCONV={
'Preshower0_charge':[12,None,None,None],
'Preshower1_charge':[12,None,None,None],
+ 'nClusters0':[None,8,None,None],
+ 'nClusters2':[14,None,None,None],
'nClusters3':[14,None,None,None],
- 'Calo_total_E_EM':[None,5000000,None,None]
-
+ 'Calo_total_E_EM':[None,5000000,None,None]
}
PHOTONCONV_L=[
lambda event: (event.Timing0_charge+event.Timing1_charge>120) or (event.Timing2_charge+event.Timing3_charge>120),
lambda event: (event.Calo_total_E_EM-330)/(sum(event.Track_p0)-max(event.Track_p0))>0.7
-]
+] #Also: BoolOrphanTSegSt3,
#LAMBDAS AND CONDITIONALS
@@ -50,9 +49,6 @@ PHOTONCONV_L=[
REALLAMBDAS=[
lambda event: event.Timing0_charge+event.Timing1_charge+event.Timing2_charge+event.Timing3_charge>0,
lambda event: event.Calo0_charge+event.Calo1_charge+event.Calo2_charge+event.Calo3_charge>0]
-# lambda event: AnyCaloThreshold,
-# lambda event: all([event.Track_nLayers[i]>6 for i in event.longTracks])
-# ]
WAVEFORMLAMBDAS=[
lambda event: event.Timing0_charge+event.Timing1_charge+event.Timing2_charge+event.Timing3_charge>0,
@@ -103,7 +99,7 @@ def BoolGoodTracks(event):
tests+=[sqrt(event.Track_x0[i]**2+event.Track_y0[i]**2)<R and sqrt(event.Track_x1[i]**2+event.Track_y1[i]**2)<R and sqrt(event.Track_X_atVetoNu[i]**2+event.Track_Y_atVetoNu[i]**2)<R]
if all(tests): gt_i+=[i]
- return len(gt_i)>=2
+ return len(gt_i)==event.longTracks
#iterate over all tracks
#count tracks that are good
@@ -515,21 +511,117 @@ def GetepenSep(event):
def testst0(event):
return [i for i in event.TrackSegment_z]
-#
+def FloatPtleInFirst10(event):
+ exclude=[]
+ tp0list=[j for j in event.Track_p0]
+ for i in range(event.longTracks):
+ if tp0list.count(tp0list[i])>1: exclude+=[i]
+
+ t10=[]
+ for j in event.truthParticleBarcode:
+ if len(t10)<=10: t10+=[j]
+
+ doublecheck=[]
+ dci=[]
+ dcount=0
+ dp=[]
+ for i in range(event.longTracks):
+ if i not in exclude:
+ bc = event.t_barcode[i]
+ if bc in doublecheck:
+ #print(event.t_pdg)
+ #print("Double-counted!")
+ dp+=[event.Track_p0[i]]
+ dcount+=1
+ elif i not in doublecheck:
+ doublecheck+=[bc]
+ dci+=[i]
+
+ return dp#dcount/len(event.t_barcode)
+
+def BoolCorrectSignalPhConv(event):
+ return -11 in event.t_pdg and 11 in event.t_pdg and (13 in event.t_pdg or -13 in event.t_pdg)
+
+def BoolCloseSignalPhConv(event):
+ return (-11 in event.t_pdg or 11 in event.t_pdg) and (13 in event.t_pdg or -13 in event.t_pdg)
+
+def FloatGetDiffVertex(event):
+ #are the electrons in the event from the same vertex?
+ #get indicies of electron tracks
+ ei=[i for i in range(event.longTracks) if abs(event.t_pdg[i])==11]
+ dist=sqrt(((event.truth_prod_x[ei[0]]-event.truth_prod_x[ei[1]])**2)+((event.truth_prod_y[ei[0]]-event.truth_prod_y[ei[1]])**2)) #+((event.truth_prod_z[ei[0]]-event.truth_prod_z[ei[1]])**2)
+ return dist
+
+def BoolIsMaxMomMuon(event):
+ plist=[i for i in event.Track_p0]
+ mi=plist.index(max(plist))
+ return event.t_pdg[mi] in [-13,13]
+
+def GetMaxPIndex(event):
+ plist=[i for i in event.Track_p0]
+ return plist.index(max(plist))
+
+# def GetTrackSepEE(event):
+# mi=GetMaxPIndex(event)
+# tri = [i for i in range(event.longTracks) if i!=mi]
+# return
+
+def GetAllNonMuonMom(event):
+ mi=GetMaxPIndex(event)
+ p=[]
+ for i in range(event.longTracks):
+ if i!=mi: p+=[event.Track_p0[i]]
+
+ return p
+
+def GetAllNonMuonMomTruth(event):
+ mi=GetMaxPIndex(event)
+ p=[]
+ if event.t_pdg[mi] not in [-13,13]: return p
+ for i in range(event.longTracks):
+ if i!=mi:
+ p+=[event.truth_P[i]]
+
+ return p
+
+def PolarMomCoord(event):
+ p=[0,0,0]
+
+ for pi in range(event.longTracks):
+ p[0]+=event.Track_px0[pi]
+ p[1]+=event.Track_py0[pi]
+ p[2]+=event.Track_pz0[pi]
+
+ return arctan(sqrt(p[0]**2+p[1]**2)/p[2])
+
+def PolarMomCoordTruth(event):
+ p=[0,0,0]
+
+ for pi in range(event.longTracks):
+ p[0]+=event.t_px[pi]
+ p[1]+=event.t_py[pi]
+ p[2]+=event.t_pz[pi]
+
+ return arctan(sqrt(p[0]**2+p[1]**2)/p[2])
+
+
+
+
+
MainArrDict={
'CutTest':{
'array':[],
'lambda': ApplyAllCuts,
- 'lkwargs':{'dcuts':[],'lcuts':[]}, #{'dcuts':[BASELINE],'lcuts':[BoolGoodTracks,BoolTimingBaseline]},
+ 'lkwargs':{'dcuts':[],'lcuts':[BoolIsMaxMomMuon]}, #{'dcuts':[BASELINE],'lcuts':[BoolGoodTracks,BoolTimingBaseline]},
'plotterfunc': PrintBoolPercent,
'plotkwargs' : {},
'descr':'Check how much a cut cuts.'
},
'BasicHist':{
'array':[],
- 'lambda': testst0,
+ 'lambda': FloatPtleInFirst10,
'plotterfunc': configBasicHist,
- 'plotkwargs' : {},
+ 'plotkwargs' : {'pckwargs':{'xscale':'log'}},
'descr':'Single variable tests.'
},
'HighestMomentumUp':{
@@ -747,6 +839,41 @@ MainArrDict={
'plotkwargs' : {},
'descr':'Single variable tests.'
},
+ 'VertexDiff':{
+ 'array':[],
+ 'lambda': FloatGetDiffVertex,
+ 'plotterfunc': configBasicHist,
+ 'plotkwargs' : {},
+ 'descr':'Single variable tests.'
+ },
+ 'NonMuonMom':{
+ 'array':[],
+ 'lambda': GetAllNonMuonMom,
+ 'plotterfunc': configBasicHist,
+ 'plotkwargs' : {'kwargs':{'range':[10**3,10**8]},'pckwargs':{'xscale':'log','xlabel':"Momentum (MeV)",'ylabel':"#Tracks",'legend':True}},#
+ 'descr':'Look at momentum tracks upstream'
+ },
+ 'AllNonMuonMomTruth':{
+ 'array':[],
+ 'lambda': GetAllNonMuonMomTruth,
+ 'plotterfunc': configBasicHist,
+ 'plotkwargs' : {'kwargs':{'range':[10**3,10**8]},'pckwargs':{'xscale':'log','xlabel':"Momentum (MeV)",'ylabel':"#Tracks",'legend':True}},#
+ 'descr':'Look at momentum tracks upstream'
+ },
+ 'PolarMomCoord':{
+ 'array':[],
+ 'lambda': PolarMomCoord,
+ 'plotterfunc': configBasicHist,
+ 'plotkwargs' : {'kwargs':{'range':[0,0.01]},'pckwargs':{}},
+ 'descr':'Single variable tests.'
+ },
+ 'PolarMomCoordTruth':{
+ 'array':[],
+ 'lambda': PolarMomCoordTruth,
+ 'plotterfunc': configBasicHist,
+ 'plotkwargs' : {'kwargs':{'range':[0,0.01]},'pckwargs':{}},
+ 'descr':'Single variable tests.'
+ },
}
TestArrDict=dict()
@@ -758,8 +885,8 @@ TestArrDict[testname]=MainArrDict[testname]
def getRunFiles(rnum):
#/eos/home-s/sshively/CALYPZONE/run/MC-tuple-110038.root
- dataDir=f"/eos/experiment/faser/sim/mc22/foresee/{rnum}/phy/r0013/" #f"{dataDir}{rnum}/"
- return sorted(glob(f"{dataDir}/Faser*-r0013-PHYS.root"))
+ dataDir=f"/eos/experiment/faser/sim/mc22/foresee/{rnum}/phy/r0014/" #f"{dataDir}{rnum}/"
+ return sorted(glob(f"{dataDir}/Faser*-PHYS.root"))
def getFlukaFile(rnum):
dataDir=f"/eos/experiment/faser/sim/mc22/particle_gun/{rnum}/phy/s0010-r0013/" #f"{dataDir}{rnum}/"
@@ -801,7 +928,7 @@ def GetPDG(t):
#EXECUTE
if __name__ == "__main__":
- #ITERATE OVER SKIMMED DATA
+ #ITERATE OVER SKIMMED DATA (phConv)
# t = TChain("nt")
# nfiles = 0
# fgrab=''
@@ -813,19 +940,19 @@ if __name__ == "__main__":
# fgrab=f
# nfiles += t.Add(f)
- #ITERATE OVER SKIMMED MC
- t = TChain("nt")
- nfiles = 0
- fgrab=''
- for runnum in getRunNums(dataDir='/eos/home-s/sshively/SKIMMED_STUFF/'):
- if runnum=='r0013' or int(runnum) not in [100041,100042]: continue #
- files = sorted(glob(f"/eos/home-s/sshively/SKIMMED_STUFF/{runnum}/test_0.root"))
- for f in files:
- print(f)
- fgrab=f
- nfiles += t.Add(f)
+ #ITERATE OVER SKIMMED MC (phConv)
+ # t = TChain("nt")
+ # nfiles = 0
+ # fgrab=''
+ # for runnum in getRunNums(dataDir='/eos/home-s/sshively/SKIMMED_STUFF/'):
+ # if runnum=='r0013' or int(runnum) not in [100041,100042]: continue #
+ # files = sorted(glob(f"/eos/home-s/sshively/SKIMMED_STUFF/{runnum}/test_0.root"))
+ # for f in files:
+ # print(f)
+ # fgrab=f
+ # nfiles += t.Add(f)
- #ITERATE OVER FULL DATA
+ #ITERATE OVER FULL DATA (phConv)
# t = TChain("nt")
# nfiles = 0
# fgrab=''
@@ -836,7 +963,7 @@ if __name__ == "__main__":
# fgrab=f
# nfiles += t.Add(f)
- #ITERATE OVER FULL MC
+ #ITERATE OVER FULL MC (phConv)
# t = TChain("nt")
# nfiles = 0
# fgrab=''
@@ -846,27 +973,28 @@ if __name__ == "__main__":
# print(f)
# fgrab=f
# nfiles += t.Add(f)
-
- #for runnum in [f'1100{ri}' for ri in [52,60,48,61,62,49,50,41,63,64,42,43,44,45,73,74,75,37,38,39,32]]:#,52,55,35,42,49,36,43,50,37,44,31,45]]:#getRunNums():
- #
+
+ #FULL DP SAMPLE
+ for runnum in [f'1100{ri}' for ri in [52,60,48,61,62,49,50,41,63,64,42,43,44,45,73,74,75,37,38,39,32]]:#,52,55,35,42,49,36,43,50,37,44,31,45]]:#getRunNums():
+ # if int(runnum[-2:]) in [28,33,38,51]: continue
+ if int(runnum[-2:]) != 38: continue
+ # if int(runnum[-2:]) not in [48,49]: continue#
- #if int(runnum[-2:]) in [28,33,38,51]: continue
- #if int(runnum[-2:]) != 38: continue
- #if int(runnum[-2:]) not in [48,49]: continue#
- #if runnum=='r0013' or int(runnum) in [100042,100041]: continue #skip MC files
- # if runnum=='r0013' or int(runnum) not in [100041]: continue
+ files = getRunFiles(runnum) # getAnyFile(runnum) # getRunFiles(runnum) getAnyRunFiles(runnum,"/eos/experiment/faser/phys/2022/r0014/") getFlukaFile('100041')
+ user_input=f'Ntuple_Run{runnum}'
- # files = getRunFiles(runnum) # getAnyFile(runnum) # getRunFiles(runnum) getAnyRunFiles(runnum,"/eos/experiment/faser/phys/2022/r0014/") getFlukaFile('100041')
- # user_input=f'Ntuple_Run{runnum}'
-
- # for f in files:
- # print(f)
- # fgrab=f
- # nfiles += t.Add(f) #t.Add("/eos/experiment/faser/phys/2022/r0011/"+user_input+".root")
+ t = TChain("nt")
+ nfiles = 0
+ fgrab=''
+ for f in files:
+ print(f)
+ fgrab=f
+ nfiles += t.Add(f) #t.Add("/eos/experiment/faser/phys/2022/r0011/"+user_input+".root")
#print(files[0])
- #GeneralPlotter(t,MakePlotDict('CutTest'),dcuts=[],lcuts=[BoolBothDaughtersFiducial])
+ GeneralPlotter(t,MakePlotDict('PolarMomCoord'),dcuts=[BASELINE,{'longTracks':[0,3,None,None],}],lcuts=[BoolTimingBaseline,BoolGoodTracks],label='Reco',showplot=False)
+ GeneralPlotter(t,MakePlotDict('PolarMomCoordTruth'),dcuts=[BASELINE,{'longTracks':[0,3,None,None],}],lcuts=[BoolTimingBaseline,BoolGoodTracks],label='Truth')
#GetProb(t,runnum)
@@ -921,22 +1049,32 @@ if __name__ == "__main__":
#ShowTrackSepVsRecoMom(t)
#3-3-2023
- # x=GetTrackEfficiency(t,dcutN=[PHOTONCONV,{'longTracks':[None,None,3,None]}],dcutD=[PHOTONCONV,{'longTracks':[1,None,None,None]}],lcutN=PHOTONCONV_L+[BoolOrphanTSegSt3,BoolMultiTrack],lcutD=PHOTONCONV_L+[BoolOrphanTSegSt3,BoolMultiTrack])
+ # x=GetTrackEfficiency(t,dcutN=[PHOTONCONV,{'longTracks':[None,None,3,None]}],dcutD=[PHOTONCONV,{'longTracks':[1,None,None,None]}],lcutN=PHOTONCONV_L+[BoolOrphanTSegSt3],lcutD=PHOTONCONV_L+[BoolOrphanTSegSt3])
# print(x)
# GeneralPlotter(t,MakePlotDict('NTSegs'),dcuts=[{'longTracks':[1,None,None,None]},PHOTONCONV],lcuts=PHOTONCONV_L)
# GeneralPlotter(t,MakePlotDict('TSeg_z'),dcuts=[{'longTracks':[1,None,None,None]},PHOTONCONV],lcuts=PHOTONCONV_L)
# GeneralPlotter(t,MakePlotDict('TSeg_z'),dcuts=[{'longTracks':[1,None,None,None]}],lcuts=[])
#3-6-2023
- #CutEfficiencyArrays(t,lambda event:event.Track_p0[1],dcutN=[PHOTONCONV,{'longTracks':[None,None,3,None]}],dcutD=[PHOTONCONV,{'longTracks':[1,None,None,None]}],lcutN=PHOTONCONV_L+[BoolOrphanTSegSt3,BoolMultiTrack],lcutD=PHOTONCONV_L+[BoolOrphanTSegSt3,BoolMultiTrack],title='Max Momentum test eff',bins=15,hmin=0,hmax=500000)
+ # CutEfficiencyArrays(t,lambda event:event.Track_p0[1],dcutN=[PHOTONCONV,{'longTracks':[None,None,3,None]}],dcutD=[PHOTONCONV,{'longTracks':[1,None,None,None]}],lcutN=PHOTONCONV_L+[BoolOrphanTSegSt3],lcutD=PHOTONCONV_L+[BoolOrphanTSegSt3],title='Max Momentum test eff',bins=15,hmin=0,hmax=500000)
#repeat but w/ 2>? Use skimmed events from jamie
#3-7-2023
# GeneralPlotter(t,MakePlotDict('epenSep'),dcuts=[PHOTONCONV,{'longTracks':[None,None,None,3]}],lcuts=PHOTONCONV_L+[BoolOrphanTSegSt3,BoolMultiTrack,BoolHas1TrackIn4St])
# print('Muons: ',COUNTMUON,' NotMu: ',COUNTNOTMUON)
- GeneralPlotter(t,MakePlotDict('BasicHist'),dcuts=[{'longTracks':[1,None,None,None]}],lcuts=[])
+ #GeneralPlotter(t,MakePlotDict('CutTest'),dcuts=[{'longTracks':[1,None,None,None]},PHOTONCONV],lcuts=PHOTONCONV_L+[BoolOrphanTSegSt3])#BoolCorrectSignalPhConv
+ #GeneralPlotter(t,MakePlotDict('VertexDiff'),dcuts=[{'longTracks':[None,None,3,None]},PHOTONCONV],lcuts=PHOTONCONV_L+[BoolOrphanTSegSt3,BoolCorrectSignalPhConv])
#print('Muons: ',COUNTMUON,' NotMu: ',COUNTNOTMUON)
+
+ #3-9-2023
+ # GeneralPlotter(t,MakePlotDict('NonMuonMom'),dcuts=[{'longTracks':[None,None,2,None]},PHOTONCONV],lcuts=PHOTONCONV_L+[BoolOrphanTSegSt3],label='==2 tracks',showplot=False)
+ # GeneralPlotter(t,MakePlotDict('NonMuonMom'),dcuts=[{'longTracks':[None,None,3,None]},PHOTONCONV],lcuts=PHOTONCONV_L+[BoolOrphanTSegSt3],label='==3 tracks')
+
+ # GeneralPlotter(t,MakePlotDict('AllNonMuonMomTruth'),dcuts=[{'longTracks':[None,None,2,None]},PHOTONCONV],lcuts=PHOTONCONV_L+[BoolOrphanTSegSt3],label='==2 tracks',showplot=False)
+ # GeneralPlotter(t,MakePlotDict('AllNonMuonMomTruth'),dcuts=[{'longTracks':[None,None,3,None]},PHOTONCONV],lcuts=PHOTONCONV_L+[BoolOrphanTSegSt3],label='==3 tracks')
+
+
# MODELS IN NEW PARAM SPACE
# 52,60,48,61,62,49,50,41,63,64,42,43,44,45,73,74,75,37,38,39,32
diff --git a/Tracker/TrackerRecAlgs/PairVertex/scripts/PlotHelper.py b/Tracker/TrackerRecAlgs/PairVertex/scripts/PlotHelper.py
index 326e7fff657f3d385843eb38ed6f41f9d52797a2..4a42d281f91589d59fb65b54266ccb81ecdb897b 100644
--- a/Tracker/TrackerRecAlgs/PairVertex/scripts/PlotHelper.py
+++ b/Tracker/TrackerRecAlgs/PairVertex/scripts/PlotHelper.py
@@ -146,10 +146,10 @@ def CorrectBinsForLog(arr,bins,xscale,rrange=None):
if minx==maxx:
maxx+=1
bins=1
- minx=1
- maxx=10.0**4
+ minx=10**3
+ maxx=10.0**8
#bins=25
- bins2 = [10**e for e in range(-5,4,1)]#logspace(start=log10(minx), stop=log10(maxx), num=bins)
+ bins2 = logspace(start=log10(minx), stop=log10(maxx), num=bins)# [10**e for e in range(-5,4,1)]
#print(type(bins2),bins2)
if xscale=='symlog':
bins2 = [0]+list(logspace(start=log10(1), stop=log10(maxx), num=bins-1))