diff --git a/Tracker/TrackerRecAlgs/PairVertex/python/PairVertexAlgConfig.py b/Tracker/TrackerRecAlgs/PairVertex/python/PairVertexAlgConfig.py index 378692f459e389c332c3757bad933102e7ceec2b..46247541be2a987df63ae42c2b2db2a32ed49564 100644 --- a/Tracker/TrackerRecAlgs/PairVertex/python/PairVertexAlgConfig.py +++ b/Tracker/TrackerRecAlgs/PairVertex/python/PairVertexAlgConfig.py @@ -10,11 +10,11 @@ from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg #Command line config from argparse import ArgumentParser -test_file = "/eos/experiment/faser/sim/mc22/foresee/110033/rec/r0013/FaserMC-MC22_FS_Aee_010MeV_1Em5-110033-00000-00009-s0008-r0013-xAOD.root" #'/eos/experiment/faser/sim/mdc/foresee/110001/rec/r0009/FaserMC-MDC_FS_Aee_100MeV_1Em5_shift-110001-00000-00009-s0007-r0009-xAOD.root' +test_file = f"/eos/experiment/faser/sim/mc22/foresee/110033/rec/r0013/" #'/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=1100033) # 110001-4 +parser.add_argument("-r","--runnum",type=int,default=110033) # 110001-4 +parser.add_argument("-i","--inputfile",type=str,default=test_file) parser.add_argument("-v","--version",type=str,default='r0013') # 110001-4 ui=parser.parse_args() @@ -27,12 +27,9 @@ inputIsMC = 'FaserMC' in ui.inputfile #Get/Make input/output file names from ROOT import TFile from glob import glob -dataDir=ui.inputfile #f"/eos/experiment/faser/sim/mdc/foresee/{ui.runnum}/rec/{ui.version}/" +dataDir=f"/eos/experiment/faser/sim/mc22/foresee/{ui.runnum}/rec/{ui.version}/" #ui.inputfile #f"/eos/experiment/faser/sim/mdc/foresee/{ui.runnum}/rec/{ui.version}/" files=sorted(glob(f"{dataDir}/Faser*")) -name=files[0].split(ui.version+"/")[-1].split("_FS_")[-1].split("shift")[0] #Generate output file name from input file name - -#------------------- - +name=files[0].split(dataDir)[-1].split(f'{ui.runnum}')[0].replace('-','_').rstrip('_')#[0]#files[0].split(ui.version+"/")[-1].split("_FS_")[-1].split(f'{ui.runnum}')[0].split() #Generate output file name from input file name def PairVertexAlgCfg(flags, **kwargs): acc = FaserSCT_GeometryCfg(flags) @@ -42,7 +39,7 @@ def PairVertexAlgCfg(flags, **kwargs): #Hist output thistSvc = CompFactory.THistSvc() - thistSvc.Output = [f"HIST DATAFILE='Ntuple_{ui.version}_{name}.root' OPT='RECREATE'"] + thistSvc.Output = [f"HIST DATAFILE='Ntuple_{name}.root' OPT='RECREATE'"] acc.addService(thistSvc) return acc diff --git a/Tracker/TrackerRecAlgs/PairVertex/scripts/PairVertexReader.py b/Tracker/TrackerRecAlgs/PairVertex/scripts/PairVertexReader.py index 8f380a827816cbb782d85531c202e1f90d494642..76947280955f5ea9b0fcb5e52974beeb38a8e26e 100644 --- a/Tracker/TrackerRecAlgs/PairVertex/scripts/PairVertexReader.py +++ b/Tracker/TrackerRecAlgs/PairVertex/scripts/PairVertexReader.py @@ -16,18 +16,6 @@ WaveformCuts2 = [ lambda e: e.Calo0_integral+e.Calo1_integral+e.Calo2_integral+e.Calo3_integral>0, lambda e: e.Timing0_integral+e.Timing1_integral+e.Timing2_integral+e.Timing3_integral>0] -#Q/P METHODS -def GetPBounds(dcuts=[],lcuts=[]): - TotalTracks=0 - PArray=[] - for event in t: - if ApplyAllCuts(event,dcuts,lcuts): - TotalTracks+=event.TrackCount - PArray+=list(event.tracksTruthPmag) - Ilims=[int(i*TotalTracks/10) for i in range(11)] - PArray.sort() - Plims=[0]+[PArray[i] for i in Ilims] - return Plims #Q/P P bounds entered manually because storing them was weird in dictionary generation PH = [191909.7203710638, 248751.19099289246] @@ -157,7 +145,6 @@ def SharedClusterLocation(event): return sharedcountarr[1] - def VertexZCompare(event): val = event.vertexPos[2] countarr=[0,0,0] #[No filter, HasShared, NoShared] @@ -217,17 +204,18 @@ def Chi2DoF(event): if event.tracksDoF[tr]>0: a+=[event.tracksChi2[tr]/event.tracksDoF[tr]] return a - #BOOLEAN FUNCTIONS def BoolDaughtersMatched(event): #both daughter particles are matched to tracks - if event.TrackCount>=2 and len(event.tracksTruthBarcode)>=2: + if len(event.tracksTruthBarcode)>=2: return 2 in event.tracksTruthBarcode and 3 in event.tracksTruthBarcode else: return False def BoolDaughtersFiducial(event): - return event.truthd0_IsFiducial and event.truthd1_IsFiducial + #There should only be one of these but they're stored in an array of doubles like all other truth particle info + #Check that there's only one then call the first element + return event.truthd0_IsFiducial[0] and event.truthd1_IsFiducial[0] def BoolDaughtersUltraFiducial(event): layermapPos=[[0,0,0],[0,0,0],[0,0,0],[0,0,0]] @@ -248,7 +236,7 @@ def BoolDaughtersUltraFiducial(event): def BoolDaughtersMatchedFiducial(event): #passes track cuts and is fiducial - return event.TrackCount>=2 and 2 in event.tracksTruthBarcode and 3 in event.tracksTruthBarcode and event.truthd0_IsFiducial and event.truthd1_IsFiducial + return BoolDaughtersMatched(event) and BoolDaughtersFiducial(event) def BoolDaughtersUnmatchedFiducial(event): #fails track cuts and is not fiducial @@ -262,35 +250,12 @@ def BoolHasSharedCluster(event): return sharedcount>0 def BoolNotHasSharedCluster(event): - #Does the event have shared clusters? - sharedcount=0 - for c in range(len(event.clusterEFracPos)): - if event.clusterETotPos[c]>0 and event.clusterETotNeg[c]>0 and event.clusterEFracPos[c]<0.99 and event.clusterEFracNeg[c]<0.99: sharedcount+=1 - return sharedcount==0 - -def BoolWhatever(event): - boollist=[ - BoolDaughtersMatched(event), - event.TrackCount>=2] - - td2i=None - td3i=None - for i in range(event.TrackCount): - if event.tracksTruthBarcode[i]==2: td2i=i - if event.tracksTruthBarcode[i]==3: td3i=i - if event.tracksDoF[i]==0: boollist+= [False] - else: boollist+= [event.tracksChi2[i]/event.tracksDoF[i]<=10] - - if type(td2i)==int and type(td3i)==int: - boollist+= [event.tracksNLayersHit[td2i]>=7 and event.tracksNLayersHit[td3i]>=7] - else: boollist+= [False] - - - return all(boollist) + #Does the event have no shared clusters? + return not BoolHasSharedCluster(event) def FlexiBool(event,stage=0): + #Sequential bool tests if stage==0: return all([BoolDaughtersMatched(event) and BoolDaughtersFiducial(event)]) - boollist=[BoolDaughtersMatched(event)] #numerator 1 if stage==1: return BoolDaughtersFiducial(event) and BoolTrackCut(event) and BoolHasSharedCluster(event) and BoolDaughtersMatched(event) @@ -305,13 +270,21 @@ def FlexiBool(event,stage=0): if stage==4: return BoolDaughtersFiducial(event) and BoolNotHasSharedCluster(event) #num 3 - if stage==5: return BoolTrackCut(event) + if stage==5: return BoolTrackCut(event) and BoolWaveformCuts(event) #denom 3 if stage==6: return event.vertexPos[2]>-1585 and event.vertexPos[2]<-45 and sqrt(event.vertexPos[0]**2+event.vertexPos[1]**2)<=100 + boollist=[] + if stage==7: return BoolFiducialVolume(event) + if stage>=8: boollist+=[BoolWaveformCuts(event)] + if stage>=9: boollist+=[event.TrackCount==2] + if stage>=10: boollist+=[BoolLayersHit(event)] + if stage>=11: boollist+=[BoolChi2DoF(event)] + return all(boollist) + def BoolChi2DoF(event): - return all([event.tracksDoF[i]>0 and event.tracksChi2[i]/event.tracksDoF[i]<10 for i in range(event.TrackCount)]) + return all([event.tracksDoF[i]>0 and (event.tracksChi2[i]/event.tracksDoF[i])<10 for i in range(event.TrackCount)]) def BoolLayersHit(event): return all([event.tracksNLayersHit[i]>=7 for i in range(event.TrackCount)]) @@ -323,9 +296,16 @@ def BoolRecoCut(event): return event.TrackCount==2 and ApplyAllCuts(event,[WaveformCuts],WaveformCuts2) and BoolChi2DoF(event) and BoolLayersHit(event) def BoolTrackCut(event): - return event.TrackCount==2 and BoolLayersHit(event) #and BoolChi2DoF(event) + return event.TrackCount==2 and BoolLayersHit(event) and BoolChi2DoF(event) + +def BoolWaveformCuts(event): + return ApplyAllCuts(event,[WaveformCuts],WaveformCuts2) + +def BoolFiducialVolume(event): + return event.vertexPos[2]>-1585 and event.vertexPos[2]<-45 and sqrt(event.vertexPos[0]**2+event.vertexPos[1]**2)<=100 -#For storing different plot configurations +#For storing different plot configurations. +#Old stuff gets moved to Storage.py MainArrDict={ 'RecoX':{ 'array':[], @@ -339,13 +319,6 @@ MainArrDict={ 'plotterfunc': configBasicHist, 'plotkwargs' : {'pckwargs':{'xscale':'log'}} }, - 'CountSharedCluster':{ - 'array':[], - 'lambda': CountSharedCluster, - 'plotterfunc': configBasicHist, - 'plotkwargs' : {'bins':25,'kwargs':{'range':[0,25]},'pckwargs':{'xlabel':'# Shared Clusters'}}, - 'descr':'How many clusters have energy deposits from both daughter particles and fail track cuts?' - }, 'CountSharedCluster_Compare':{ 'array':[[],[],[]], 'lambda': CountSharedClusterCompare, @@ -354,160 +327,6 @@ MainArrDict={ 'descr':'How many clusters have energy deposits from both daughter particles?', 'SepHist':'arrays allowed to have different lengths' }, - 'DPDecayVertexZ':{ - 'array':[], - 'lambda': lambda event: event.vertexPos[2], - 'plotterfunc': configBasicHist, - 'plotkwargs' : {}, - 'descr':'Dark photon decay vertex coordinate along beamline' - }, - 'DPDecayVertexZ_Compare':{ - 'array':[[],[],[]], - 'lambda': VertexZCompare, - 'plotterfunc': configMultiHist, - 'plotkwargs' : {'label':['No Cuts','Has Shared Clusters','No Shared Clusters'],'bins':25,'kwargs':{'range':[-1600,0]},'pckwargs':{'xscale':'linear','legend':True,'legloc':'center right'}}, - 'descr':'Compare dp-decay z-coordinate cuts', - 'SepHist':'arrays allowed to have different lengths' - }, - 'NumClusters':{ - 'array':[], - 'lambda': lambda event: event.clusterCount, - 'plotterfunc': configBasicHist, - 'plotkwargs' :{'kwargs':{'range':[0,200]}}, - }, - 'NumClustersBool':{ - 'array':[], - 'lambda': lambda event: event.clusterCount<200, - 'plotterfunc': PrintBoolPercent, - 'plotkwargs' :{}, - 'descr':'How many events have a cluster count in a certain range?' - }, - 'ClusterStation':{ - 'array':[], - 'lambda': lambda event: [i for i in event.clusterLocStation], - 'plotterfunc': configBasicHist, - 'plotkwargs' :{'bins':8,'kwargs':{'range':[0,4]},'pckwargs':{'xlabel':'Station #','ylabel':'#Clusters'}}, - 'descr':'Where are any clusters?' - }, - 'SharedClusterLocation':{ - 'array':[], - 'lambda': SharedClusterLocation, - 'plotterfunc': configBasicHist, - 'plotkwargs' :{'bins':8,'kwargs':{'range':[0,4]},'pckwargs':{'xlabel':'Station #','ylabel':'#Clusters'}}, - 'descr': 'Where are the shared clusters?' - }, - 'FiducialDaughters':{ - 'array':[], - 'lambda': BoolDaughtersFiducial, - 'plotterfunc': PrintBoolPercent, - 'plotkwargs' :{}, - 'descr':'' - }, - 'FiducialMatchedDaughters':{ - 'array':[], - 'lambda': BoolDaughtersMatchedFiducial, - 'plotterfunc': PrintBoolPercent, - 'plotkwargs' :{}, - 'descr':'' - }, - 'BoolDaughtersUltraFiducial':{ - 'array':[], - 'lambda': BoolDaughtersUltraFiducial, - 'plotterfunc': PrintBoolPercent, - 'plotkwargs' :{}, - 'descr':'' - }, - 'BoolDaughtersUltraFiducialMatched':{ - 'array':[], - 'lambda':BoolDaughtersUltraFiducial, - 'plotterfunc': PrintBoolPercent, - 'plotkwargs' :{}, - 'descr':'' - }, - 'BoolWhatever':{ - 'array':[], - 'lambda': BoolWhatever, - 'plotterfunc': PrintBoolPercent, - 'plotkwargs' :{}, - 'descr':'something with efficiency' - }, - 'ClustersPerDaughter':{ - 'array':[], - 'lambda': ClustersPerDaughter, - 'plotterfunc': configBasicHist, - 'plotkwargs' :{'bins':41,'kwargs':{'range':[0,41]},'pckwargs':{'xlabel':'# of Clusters','ylabel':'# Daughter Particle'}}, - 'descr': 'How many clusters does a daughter particle cause?' - }, - 'PercentDaughterClusterShared':{ - 'array':[], - 'lambda': PercentDaughterClusterShared, - 'plotterfunc': configBasicHist, - 'plotkwargs' :{'pckwargs':{'xlabel':'Percent of clusters shared','ylabel':'# Daughter Paricles'}}, - 'descr': 'How many clusters does a daughter particle cause?' - }, - 'Chi2DoF':{ - 'array':[], - 'lambda': Chi2DoF, - 'plotterfunc': configBasicHist, - 'plotkwargs' :{'kwargs':{'range':[0,20]},'pckwargs':{'xlabel':'Chi2/DoF','ylabel':'# Tracks'}}, - 'descr': 'How many clusters does a daughter particle cause?' - }, - # 1-24-23 plots - 'Efficiency0':{ - 'array':[], - 'lambda': FlexiBool, - 'lkwargs':{'stage':0}, - 'plotterfunc': PrintBoolPercent, - 'plotkwargs' :{}, - }, - 'Efficiency1':{ - 'array':[], - 'lambda': FlexiBool, - 'lkwargs':{'stage':1}, - 'plotterfunc': PrintBoolPercent, - 'plotkwargs' :{}, - }, - 'Efficiency2':{ - 'array':[], - 'lambda': FlexiBool, - 'lkwargs':{'stage':2}, - 'plotterfunc': PrintBoolPercent, - 'plotkwargs' :{}, - }, - 'Efficiency3':{ - 'array':[], - 'lambda': FlexiBool, - 'lkwargs':{'stage':3}, - 'plotterfunc': PrintBoolPercent, - 'plotkwargs' :{}, - }, - 'Efficiency4':{ - 'array':[], - 'lambda': FlexiBool, - 'lkwargs':{'stage':4}, - 'plotterfunc': PrintBoolPercent, - 'plotkwargs' :{}, - }, - 'Efficiency5':{ - 'array':[], - 'lambda': FlexiBool, - 'lkwargs':{'stage':5}, - 'plotterfunc': PrintBoolPercent, - 'plotkwargs' :{}, - }, - 'Efficiency6':{ - 'array':[], - 'lambda': FlexiBool, - 'lkwargs':{'stage':6}, - 'plotterfunc': PrintBoolPercent, - 'plotkwargs' :{}, - }, - 'RecoMinMomentum':{ - 'array':[], - 'lambda': lambda event:min(event.tracksPmag), - 'plotterfunc': configBasicHist, - 'plotkwargs' : {'pckwargs':{'xscale':'log'}} - }, } @@ -560,6 +379,7 @@ def GeneralPlotter(arrDict,dcuts=[],lcuts=[],showplot=True,savefig=False,titleex for name in arrDict.keys(): isPrintout= arrDict[name]['plotterfunc'] in [PrintBoolPercent] #Is it text-only output? print(name) + #print(arrDict[name]['descr']) print('LEN:',len(arrDict[name]['array'])) #how many events remain after the cut? if label: arrDict[name]['plotkwargs'].update({'label':label}) @@ -576,41 +396,15 @@ def GeneralPlotter(arrDict,dcuts=[],lcuts=[],showplot=True,savefig=False,titleex #EXECUTE if __name__ == "__main__": - DP_fnames=['Ntuple_Aee_100MeV_1Em5_']# Ntuple_Aee_10MeV_1Em5_ Ntuple_Amm_316MeV_2Em6_ Ntuple_Aee_10MeV_1Em4_ + DP_fnames=['Ntuple_FaserMC_MC22_FS_Aee_100MeV_1Em5']# 'Ntuple_Aee_10MeV_1Em4_' 'Ntuple_Aee_100MeV_1Em5_','Ntuple_Aee_10MeV_1Em5_', ,'Ntuple_Amm_316MeV_2Em6_' .root for f in DP_fnames: user_input=f t = TChain("events") nfiles = 0 - nfiles += t.Add("/eos/home-s/sshively/CALYPZONE/run/BackUp2/"+user_input+".root") + nfiles += t.Add("/eos/home-s/sshively/CALYPZONE/run/"+user_input+".root") print(f) - #1-24-23 Plots - # GeneralPlotter(MakePlotDict('Efficiency0'),dcuts=[{}],lcuts=[]) - # GeneralPlotter(MakePlotDict([f'Efficiency{i}' for i in range(1,7)]),dcuts=[{}],lcuts=[]) - # print() - # GeneralPlotter(MakePlotDict('Efficiency4'),dcuts=[{}],lcuts=[BoolDaughtersFiducial,BoolHasSharedCluster]) - # GeneralPlotter(MakePlotDict('Efficiency4'),dcuts=[{}],lcuts=[BoolDaughtersFiducial,BoolNotHasSharedCluster]) - - # GeneralPlotter(MakePlotDict('Chi2DoF'),dcuts=[{'TrackCount':[0,None,None,None]}],lcuts=[],label='All Events',showplot=False) - # GeneralPlotter(MakePlotDict('Chi2DoF'),dcuts=[{'TrackCount':[0,None,None,None]}],lcuts=[BoolDaughtersMatched,BoolDaughtersFiducial],label='Daughters Matched & Fiducial') - - # GeneralPlotter(MakePlotDict('CountSharedCluster'),dcuts=[],lcuts=[],label='All Events',showplot=False) - # GeneralPlotter(MakePlotDict('CountSharedCluster'),dcuts=[{'TrackCount':[1,None,None,None]}],lcuts=[BoolDaughtersMatched,BoolDaughtersFiducial,BoolChi2DoF, BoolLayersHit],label='Passes Track Cuts',showplot=False) - # GeneralPlotter(MakePlotDict('CountSharedCluster'),dcuts=[{'TrackCount':[None,None,None,2]}],lcuts=[BoolDaughtersMatched,BoolDaughtersFiducial,BoolAntiTrack],label='Fails Track Cuts') - - ####GeneralPlotter(MakePlotDict('ClusterStation'),dcuts=[],lcuts=[],label='All Clusters',showplot=False) #meh - #GeneralPlotter(MakePlotDict('SharedClusterLocation'),dcuts=[],lcuts=[],label='Shared Cllusters') - - #GeneralPlotter(MakePlotDict('DPDecayVertexZ_Compare'),dcuts=[],lcuts=[]) - - #GeneralPlotter(MakePlotDict('RecoMinMomentum'),dcuts=[{'TrackCount':[None,None,2,None]},WaveformCuts],lcuts=[BoolDaughtersMatched,BoolDaughtersFiducial,BoolChi2DoF, BoolLayersHit]+WaveformCuts2) - - #1-25-23 - #GeneralPlotter(MakePlotDict('FiducialDaughters'),dcuts=[{}],lcuts=[]) - # GeneralPlotter(MakePlotDict([f'Efficiency{i}' for i in range(6,7)]),dcuts=[{}],lcuts=[]) - - # GeneralPlotter(MakePlotDict('Chi2DoF'),dcuts=[{}],lcuts=[BoolDaughtersMatched,BoolTrackCut,BoolNotHasSharedCluster],label='No Shared',showplot=False) - # GeneralPlotter(MakePlotDict('Chi2DoF'),dcuts=[{}],lcuts=[BoolDaughtersMatched,BoolTrackCut,BoolHasSharedCluster],label='Shared') + GeneralPlotter(MakePlotDict('RecoMaxMomentum'),dcuts=[{}],lcuts=[]) diff --git a/Tracker/TrackerRecAlgs/PairVertex/scripts/PlotHelper.py b/Tracker/TrackerRecAlgs/PairVertex/scripts/PlotHelper.py index e3da186362c8f90e0f931c1182b1d56070f40fbc..db770deaf55ac8c212fa5fbc848b07971f1632a1 100644 --- a/Tracker/TrackerRecAlgs/PairVertex/scripts/PlotHelper.py +++ b/Tracker/TrackerRecAlgs/PairVertex/scripts/PlotHelper.py @@ -25,9 +25,11 @@ def configHitMap(arrays,bins=(50,50),pckwargs={'xlabel':"X (mm)",'ylabel':"Y (mm if counts: plt.text(50,-185,"#EVENTS: {}".format(len(arrayx))) PlotConfigs(**pckwargs) -def configBasicHist(arr,bins=25,label=None,pckwargs={'xscale':'linear'},kwargs={}): +def configBasicHist(arr,bins=25,label=None,pckwargs={'xscale':'linear'},kwargs={},rrange=None): #1D Histogram - if 'xscale' in pckwargs.keys() and 'log' in pckwargs['xscale']: bins=CorrectBinsForLog(arr,bins,pckwargs['xscale']) + if 'xscale' in pckwargs.keys() and 'log' in pckwargs['xscale']: + if rrange:CorrectBinsForLog(arr,bins,pckwargs['xscale'],rrange=rrange) + else: bins=CorrectBinsForLog(arr,bins,pckwargs['xscale']) n, binsOut, patches = plt.hist(arr,bins=bins,label=label,align="left",histtype=u'step',**kwargs) PlotConfigs(**pckwargs) @@ -88,8 +90,9 @@ def PlotRootRatio(arrlist,bins=15,hmin=0,hmax=10000,title=None): input() #Utilities -def CorrectBinsForLog(arr,bins,xscale): +def CorrectBinsForLog(arr,bins,xscale,rrange=None): minx,maxx = (min(arr),max(arr)) + if rrange!=None: minx,maxx = (rrange[0],rrange[1]) if minx==maxx: maxx+=1 bins=1