diff --git a/DrawHits/drawHits.py b/DrawHits/drawHits.py index 3748b9f..a140072 100755 --- a/DrawHits/drawHits.py +++ b/DrawHits/drawHits.py @@ -254,6 +254,8 @@ def fillHistoByDef(tree,hDef,extraCuts): if hDef.vetoMType(mType): continue is1D = hDef.getParameter('yNbins',mType)==None + is2D = hDef.getParameter('yNbins',mType)!=None and hDef.getParameter('zNbins',mType)==None + is3D = hDef.getParameter('yNbins',mType)!=None and hDef.getParameter('zNbins',mType)!=None isProfile = hDef.getParameter('profile',mType)!=None and hDef.getParameter('profile',mType) effCuts = hDef.getParameter('effCuts',mType) variable = hDef.getParameter('variable',mType) @@ -278,6 +280,7 @@ def fillHistoByDef(tree,hDef,extraCuts): tree.Project(hName+"_2",variable, \ cutString(extraCuts,hDef.getParameter('baseCuts',mType),"moduleType=="+str(mType),effCuts)) histos[mType][3] = ROOT.TEfficiency(histos[mType][1],histos[mType][0]) + histos[mType][3].SetMarkerStyle(20) else: # always keep final histogram in 4th position histos[mType][3] = histos[mType][0] @@ -290,7 +293,7 @@ def fillHistoByDef(tree,hDef,extraCuts): cutString(extraCuts,hDef.getParameter('baseCuts',mType),"moduleType=="+str(mType))) # always keep final histogram in 4th position histos[mType][3] = histos[mType][0] - else: + elif is2D: nby = hDef.getParameter('yNbins',mType) ymin = hDef.getParameter('yMin',mType) ymax = hDef.getParameter('yMax',mType) @@ -310,13 +313,27 @@ def fillHistoByDef(tree,hDef,extraCuts): else: # always keep final histogram in 4th position histos[mType][3] = histos[mType][0] + elif is3D: + nby = hDef.getParameter('yNbins',mType) + ymin = hDef.getParameter('yMin',mType) + ymax = hDef.getParameter('yMax',mType) + nbz = hDef.getParameter('zNbins',mType) + zmin = hDef.getParameter('zMin',mType) + zmax = hDef.getParameter('zMax',mType) + histos[mType] = [ ROOT.TH3F(hName+"_1",hName+"_1",nbx,xmin,xmax,nby,ymin,ymax,nbz,zmin,zmax), \ + None, None, None ] + tree.Project(hName+"_1",variable, \ + cutString(extraCuts,hDef.getParameter('baseCuts',mType),"moduleType=="+str(mType))) + assert effCuts==None + # always keep final histogram in 4th position + histos[mType][3] = histos[mType][0] #print("Ending for ",hDef.name,hName,hTitle) savedDir.cd() return histos -def drawHistoByDef(histos,hDef,logY=False,same=False): +def drawHistoByDef(histos,hDef,logY=False,logZ=False,same=False): result = { 'cnv' : None, 'histos' : histos, 'pave' : None } savedDir = ROOT.gDirectory @@ -344,6 +361,8 @@ def drawHistoByDef(histos,hDef,logY=False,same=False): if hDef.vetoMType(mType): continue is1D = hDef.getParameter('yNbins',mType)==None + is2D = hDef.getParameter('yNbins',mType)!=None and hDef.getParameter('zNbins',mType)==None + is3D = hDef.getParameter('yNbins',mType)!=None and hDef.getParameter('zNbins',mType)!=None isProfile = hDef.getParameter('profile',mType)!=None and hDef.getParameter('profile',mType) effCuts = hDef.getParameter('effCuts',mType) variable = hDef.getParameter('variable',mType) @@ -362,6 +381,7 @@ def drawHistoByDef(histos,hDef,logY=False,same=False): xmax = hDef.getParameter('xMax',mType) ytitle = hDef.getParameter('yTitle',mType) if hDef.getParameter('yTitle',mType) else "" + ztitle = hDef.getParameter('zTitle',mType) if hDef.getParameter('zTitle',mType) else "" if is1D and ( not isProfile ): ymin = hDef.getParameter('yMin',mType) if hDef.getParameter('yMin',mType)!=None else 0. ymax = hDef.getParameter('yMax',mType) if hDef.getParameter('yMax',mType)!=None else 1.05 @@ -372,13 +392,16 @@ def drawHistoByDef(histos,hDef,logY=False,same=False): histos[mType][2].GetXaxis().SetTitle(xtitle) histos[mType][2].GetYaxis().SetTitle(ytitle) histos[mType][3] = ROOT.TEfficiency(histos[mType][1],histos[mType][0]) - histos[mType][3].SetMarkerSize(0.3) + histos[mType][3].SetMarkerStyle(20) histos[mType][3].Draw("same Z") else: histos[mType][0].SetTitle(hTitle) histos[mType][0].GetXaxis().SetTitle(xtitle) histos[mType][0].GetYaxis().SetTitle(ytitle) histos[mType][0].Draw("same" if same else "") + fitFunc = hDef.getParameter('fit') + if fitFunc!=None: + histos[mType][0].Fit(fitFunc,"Q","same") elif isProfile: histos[mType][0].SetTitle(hTitle) histos[mType][0].GetXaxis().SetTitle(xtitle) @@ -387,7 +410,7 @@ def drawHistoByDef(histos,hDef,logY=False,same=False): #histos[mType][0].SetFillColor(ROOT.TColor.GetColorBright(ROOT.kGray)) #histos[mType][0].SetFillColor(ROOT.kGray) histos[mType][0].Draw("same" if same else "") - else: + elif is2D: assert not same zmin = hDef.getParameter('zMin',mType) if hDef.getParameter('zMin',mType)!=None else 0. zmax = hDef.getParameter('zMax',mType) if hDef.getParameter('zMax',mType)!=None else 1.05 @@ -403,8 +426,20 @@ def drawHistoByDef(histos,hDef,logY=False,same=False): histos[mType][0].GetXaxis().SetTitle(xtitle) histos[mType][0].GetYaxis().SetTitle(ytitle) histos[mType][0].Draw("ZCOL") + elif is3D: + assert not same + zmin = hDef.getParameter('zMin',mType) if hDef.getParameter('zMin',mType)!=None else 0. + zmax = hDef.getParameter('zMax',mType) if hDef.getParameter('zMax',mType)!=None else 1.05 + assert effCuts==None + histos[mType][0].SetTitle(hTitle) + histos[mType][0].GetXaxis().SetTitle(xtitle) + histos[mType][0].GetYaxis().SetTitle(ytitle) + histos[mType][0].GetZaxis().SetTitle(ztitle) + histos[mType][0].Draw() if logY or hDef.getParameter('logY',mType): ROOT.gPad.SetLogy(1) + if is2D and ( logZ or hDef.getParameter('logZ',mType) ): + ROOT.gPad.SetLogz(1) ROOT.gPad.Update() #cnv.cd() @@ -484,10 +519,13 @@ def addHistogram(varString,cuts,effCuts=None,name='userHist'): type=str, default='*') parser.add_argument('--vetoedHistograms', help='comma-separated names of histogram definitions not to be used', type=str, default='') -parser.add_argument('--logY', help='use log scale', action='store_true', default=False) +parser.add_argument('--logY', help='use log scale for y axis', action='store_true', default=False) +parser.add_argument('--logZ', help='use log scale for z axis', action='store_true', default=False) parser.add_argument('--printTree', '-p', help='print TTree contents', action='store_true', default=False) parser.add_argument('--listHistograms', '-l', help='list predefined and selected histograms', \ action='store_true', default=False) +parser.add_argument('--zone', '-z', help='restrict to zone in OT (barrel, tilted, endcap)', type=str, \ + choices=['barrel','tilted','endcap'], default=None) parser.add_argument('file', help='input file', type=str, nargs='+', default=None) args = parser.parse_args() outputFormats = [ ] @@ -498,6 +536,19 @@ def addHistogram(varString,cuts,effCuts=None,name='userHist'): fitResiduals = args.fitResiduals.split(",") if args.fitResiduals else [ ] selectedHistoNames = args.selectedHistograms.split(",") vetoedHistoNames = args.vetoedHistograms.split(",") +# +# add cut for zone definition? +# +if args.zone=="barrel": + zoneCuts = "detNormal.Rho()>0.99" +elif args.zone=="endcap": + zoneCuts = "detNormal.Rho()<0.01" +elif args.zone=="tilted": + zoneCuts = "detNormal.Rho()>0.05&&detNormal.Rho()<0.095" +else: + zoneCuts = "" +args.cuts = cutString(args.cuts,zoneCuts) + # # load histogram definitions # @@ -591,7 +642,8 @@ def addHistogram(varString,cuts,effCuts=None,name='userHist'): for hName in allHDefs.byCanvas[cName]: print("Processing histogram",hName,"in canvas",cName) cHistos[hName] = fillHistoByDef(simHitTree,allHDefs.byCanvas[cName][hName],extraCuts) - allObjects.append(drawHistoByDef(cHistos[hName],allHDefs.byCanvas[cName][hName],logY=args.logY,same=same)) + allObjects.append(drawHistoByDef(cHistos[hName],allHDefs.byCanvas[cName][hName], \ + logY=args.logY,logZ=args.logZ,same=same)) same = True if args.output!=None: c = allObjects[-1]['cnv'] @@ -599,6 +651,8 @@ def addHistogram(varString,cuts,effCuts=None,name='userHist'): basename = os.path.join(args.output,c.GetName()) if args.sampleName!=None: basename += "_" + args.sampleName + if args.zone!=None: + basename += "_" + args.zone print(basename) for fmt in outputFormats: c.SaveAs(basename+"."+fmt) diff --git a/HitAnalyzer/test/histogramDefinition.py b/DrawHits/drawHitsConfiguration.py similarity index 51% rename from HitAnalyzer/test/histogramDefinition.py rename to DrawHits/drawHitsConfiguration.py index 8bfc321..c8b8095 100644 --- a/HitAnalyzer/test/histogramDefinition.py +++ b/DrawHits/drawHitsConfiguration.py @@ -1,20 +1,30 @@ # # Definitions for histograms from a configuration file # +import yaml from fnmatch import fnmatch class HistogramDefinition: ''' A single histogram definition. ''' - reqGenFields = [ 'canvasName', 'histogramName', 'histogramTitle', 'variable', 'baseCuts' ] + # fields that have to be present in the general section + reqGenFields = [ 'canvasName' ] + # fields that have to be present in a histogram section reqHistFields = [ ] - requiredFields = reqGenFields + reqHistFields - optGenFields = [ 'effCuts', 'logY' ] - optHistFields = [ 'xNbins', 'xMin', 'xMax', 'xTitle', 'yTitle', 'yNbins', 'yMin', 'yMax', \ - 'zMin', 'zMax', 'display', 'profile' ] - optionalFields = optGenFields + optHistFields - allFields = requiredFields + optionalFields - allHistFields = reqHistFields + optHistFields + requiredFields = list(set(reqGenFields + reqHistFields)) + # optional fields in the general section + optGenFields = [ 'variable', 'baseCuts', 'effCuts', 'logY', 'logZ' ] + # optional fields in the histogram section + optHistFields = [ 'variable', 'histogramName', 'histogramTitle', 'fit', \ + 'xNbins', 'xMin', 'xMax', 'xTitle', 'yTitle', 'zTitle', \ + 'yNbins', 'yMin', 'yMax', \ + 'zNbins', 'zMin', 'zMax', 'display', 'profile' ] + optionalFields = list(set(optGenFields + optHistFields)) + allFields = list(set(requiredFields + optionalFields)) + allHistFields = list(set(reqHistFields + optHistFields)) + # fields that cannot be present for a single module type + vetoMtypeFields = [ 'variable', 'baseCuts', 'effCuts', 'profile' ] + vetoMtypeFields = [ 'profile' ] def __init__(self,name,inputDict): ''' Define histogram and drawing parameters from dictionary. @@ -98,6 +108,9 @@ def getParameter(self,name,mType=None): if ( mTName in self.parameters ) and ( name in self.parameters[mTName] ): result = self.parameters[mTName][name] if result!=None: + # check for parameters that can only be general + if name in HistogramDefinition.vetoMtypeFields: + raise Exception('Parameter '+name+' cannot be present in the section for an individual module type') return self.parameters[mTName][name] # # not found: use general parameter @@ -106,6 +119,11 @@ def getParameter(self,name,mType=None): return self.parameters[name] return None + def __call__(self,name,mType=None): + ''' Make access to parameters easier - use function call + ''' + return self.getParameter(name,mType) + def vetoMType(self,mType): ''' Check for an mType entry with display = False ''' @@ -151,29 +169,124 @@ def __getitem__(self,name): return self.allDefinitions[name] return None -def loadHistogramDefinitions(configName,selectedNames=[],vetoedNames=[]): - ''' Load histogram definitions from a configuration file. The configuration file +class VarMaskCombinations: + ''' Combinations of a variable and a mask for use with RDF.Histo*. + ''' + def __init__(self): + # + # variable+mask name indexed by variable+selection string + # + self.varMaskDefinitions = { } + + def __call__(self,rdf,varName,selection): + ''' Return the name of a variable+mask combination. Define if necessary. + Arguments: + rdf ........ RDataFrame to be used for definition + varName .... name of the variable + selection .. selection string to be used as mask + ''' + # + # use variable name and normalized selection string (remove spaces) + # + selNorm = selection.replace(" ","") + varMask = "(" + varName + ")["+selNorm+"]" + if varMask in self.varMaskDefinitions: + # combination exists (assume that it's also defined in the RDataFrame) + varMaskName = self.varMaskDefinitions[varMask] + assert varMaskName in rdf.GetDefinedColumnNames() + else: + # create name for combination and define it in the RDataFrame + n = len(self.varMaskDefinitions) + varMaskName = "varMask{:03d}".format(n) + self.varMaskDefinitions[varMask] = varMaskName + rdf = rdf.Define(varMaskName,varMask) + #print("+++ defined new var + mask combinations:",varMaskName,varMask) + #print(rdf.GetDefinedColumnNames()) + #print("+++") + return rdf,varMaskName + +class RDFWrapper: + ''' Wrapper of RDataFrame with possibility to define variable+mask combinations + ''' + def __init__(self,rdf): + # + # variable+mask name indexed by variable+selection string + # + self.rdf = rdf + self.varMaskDefinitions = { } + + def __call__(self): + ''' Easy access to RDataFrame + ''' + return self.rdf + + def defineVarMask(self,varName,selection): + ''' Return the name of a variable+mask combination. Define if necessary. + Arguments: + rdf ........ RDataFrame to be used for definition + varName .... name of the variable + selection .. selection string to be used as mask + ''' + # + # use variable name and normalized selection string (remove spaces) + # + selNorm = selection.replace(" ","") + varMask = "(" + varName + ")["+selNorm+"]" + if varMask in self.varMaskDefinitions: + # combination exists (assume that it's also defined in the RDataFrame) + varMaskName = self.varMaskDefinitions[varMask] + assert varMaskName in self.rdf.GetDefinedColumnNames() + else: + # create name for combination and define it in the RDataFrame + n = len(self.varMaskDefinitions) + varMaskName = "varMask{:03d}".format(n) + self.varMaskDefinitions[varMask] = varMaskName + self.rdf = self.rdf.Define(varMaskName,varMask) + #print("+++ defined new var + mask combinations:",varMaskName,varMask) + #print(rdf.GetDefinedColumnNames()) + #print("+++") + return varMaskName + + +def loadConfiguration(configName,selectedNames=[],vetoedNames=[]): + ''' Load variable and histogram definitions from a configuration file. The configuration file + is expected to contain a "variables" and a "histograms" section. The "histogram" section defines dictionaries with the definitions as variable. The name of the variable serves as name of the HistogramDefinition. Arguments: configName ..... name of the module to import from / python file name selectedNames .. explicit list of histogram names to be imported (filename wildcard syntax) vetoedNames .... explicit list of histogram names to be skipped (filename wildcard syntax) + Result: + Tuple with dictionary of variable definitions and HistogramDefinitions object ''' # load histogram definitions # - result = HistogramDefinitions() + vDefs = { } + hDefs = HistogramDefinitions() if configName==None: - return result + return ( vDefs, hDefs ) - moduleName = configName[:-3] if configName.endswith(".py") else configName - module = __import__(moduleName) - for n in dir(module): - if n.startswith('__'): - continue - hDict = getattr(module,n) - #print(n,type(hDict)) - #sys.exit() + with open(configName,"rt") as yamlFile: + allDicts = yaml.load(yamlFile,Loader=yaml.Loader) + yamlFile.close() + # + # store variable definitions as read from configuration + # + if 'variables' in allDicts: + vDefs = allDicts['variables'] + # + # backward compatibility for histograms: treat full input as "histograms" section + # in case neither variables nor histograms sections are defined + # + if 'histograms' in allDicts: + histDicts = allDicts['histograms'] + else: + assert not ( 'variables' in allDicts ) + histDicts = allDicts + # + for n,hDict in histDicts.items(): + # assert type(hDict)==dict # # check if in list of histograms to be displayed @@ -199,8 +312,8 @@ def loadHistogramDefinitions(configName,selectedNames=[],vetoedNames=[]): # add histogram # hDef = HistogramDefinition(n,hDict) - result.add(hDef) + hDefs.add(hDef) print("Added",hDef.getParameter('canvasName')) - - return result + + return ( vDefs, hDefs ) diff --git a/DrawHits/drawHitsRDF.py b/DrawHits/drawHitsRDF.py new file mode 100755 index 0000000..5852301 --- /dev/null +++ b/DrawHits/drawHitsRDF.py @@ -0,0 +1,692 @@ +#! /bin/env python3 +import sys,os +from drawHitsConfiguration import * +import ROOT +import argparse +from fnmatch import fnmatch + +def divideRatios(cnv,nx,ny,widthRatios,heightRatios,canvasTopMargin=0.,canvasLeftMargin=0.): + result = [ ] + cnv.cd() + + wrs = len(widthRatios) + hrs = len(heightRatios) + nxl = min(nx,wrs) + nyl = min(ny,hrs) + + if wrs==0: nxl = nx + if hrs==0: nyl = ny + + pn = 1 + xr = 0. + yr = 0. + x = 0. + y = 1. + # Check the validity of the margins + if canvasTopMargin<0 or canvasTopMargin>1: + raise Exception("DivideRatios: The canvas top margin must be >= 0 and <= 1") + else: + y = 1.- canvasTopMargin + if canvasLeftMargin <0 or canvasLeftMargin >1: + raise Exception("DivideRatios", "The canvas left margin must be >= 0 and <= 1") + + # Check the validity of the ratios + sumOfHeightRatios = canvasTopMargin + if hrs: + for i in range(nyl): + yr = heightRatios[i] + sumOfHeightRatios = sumOfHeightRatios + yr + if yr<0 or yr>1: + raise Exception("DivideRatios", "Y ratios plus the top margin must be >= 0 and <= 1") + + if sumOfHeightRatios>1.: + raise Exception("DivideRatios", "The sum of Y ratios plus the top margin must be <= 1 %g",sumOfHeightRatios) + sumOfWidthRatios = canvasLeftMargin + if wrs: + for j in range(nxl): + xr = widthRatios[j] + sumOfWidthRatios = sumOfWidthRatios + xr + if xr <0 or xr >1: + raise Exception("DivideRatios", "X ratios must be >= 0 and <= 1") + if sumOfWidthRatios>1.: + raise Exception("DivideRatios", "The sum of X ratios must be <= 1 %g ",sumOfWidthRatios) + + # Create the pads according to the ratios + for i in range(nyl): + x = canvasLeftMargin + if hrs: + yr = heightRatios[i] + else: + yr = 1./nyl + for j in range(nxl): + if wrs: + xr = widthRatios[j] + else: + xr = 1./nxl + x1 = max(0., x) + y1 = max(0., y - yr) + x2 = min(1., x + xr) + y2 = min(1., y) + pad = ROOT.TPad(cnv.GetName()+str(pn),cnv.GetName()+str(pn),x1, y1, x2 ,y2) + result.append(pad) + pad.SetNumber(pn) + pad.Draw() + x = x + xr + pn += 1 + y = y - yr + + return result + +def fitHistogram(mType,h): + f1name = "f1"+str(mType) + f1 = ROOT.TF1(f1name,"gaus(0)") + f1.SetParameter(0,h.GetMaximum()) + f1.SetParLimits(0,0.,2*h.GetMaximum()) + f1.SetParameter(1,0.) + f1.SetParameter(2,h.GetRMS()/10.) + h.Fit(f1name) + f2name = "f2"+str(mType) + f2 = ROOT.TF1(f2name,"gaus(0)+gaus(3)") + f2.SetParameter(0,f1.GetParameter(0)) + f2.SetParameter(1,f1.GetParameter(1)) + f2.SetParameter(2,f1.GetParameter(2)) + f2.SetParameter(3,f1.GetParameter(0)/100.) + f2.SetParameter(4,f1.GetParameter(1)) + f2.SetParameter(5,5*f1.GetParameter(2)) + f2.SetParLimits(5,f1.GetParameter(2),10*f1.GetParameter(2)) + h.Fit(f2name) + ROOT.gPad.SetLogy(1) + ROOT.gPad.Update() + + +def cutString(*cuts): + if len(cuts)>0 and cuts[0]!=None: + return "&&".join([ c for c in cuts if ( c!=None and c.strip()!="" ) ]) + return None + +def cutLines(cuts,maxChar=40): + result = [ ] + line = "" + for ic,c in enumerate(cuts): + line += c + if ic<(len(cuts)-1): + line += "&&" + if len(line)>maxChar: + result.append(line) + line = "" + if line!="": + result.append(line) + return result + +# def drawString(pave,title,s=""): +# t = pave.AddText(title) +# t.SetTextFont(43) +# t.SetTextSize(0.06) +# t.SetTextSize(16) +# t.SetTextAlign(13) +# #yt += 14 +# if s!="": +# #t = pave.AddText("") +# t = pave.AddText(" "+s) +# t.SetTextFont(43) +# t.SetTextSize(14) +# t.SetTextAlign(13) +# #t = pave.AddText("") + +# def drawCuts(pave,title,cuts): +# #t = pave.AddText(title) +# #t = pave.AddText("") +# lines = [ ] +# line = "" +# for ic,c in enumerate(cuts): +# line += c +# if ic<(len(cuts)-1): +# line += "&&" +# if len(line)>40: +# lines.append(line) +# line = "" +# if line!="": +# lines.append(line) +# for l in lines: +# if l!="": +# t = pave.AddText(" "+l) +# t.SetTextFont(43) +# t.SetTextSize(0.04) +# t.SetTextSize(14) +# t.SetTextAlign(13) +# #l = " " + c +# #if ic<(len(cuts)-1): +# # l += " &&" +# ## if len(l +# #if +# #t = pave.AddText(l) + + +def drawCutPave(cnv,ic,variable,cuts,effcuts=None): + indBaseCuts = cuts.split("&&") + indEffCuts = None if effcuts==None else effcuts.split("&&") + cnv.cd(ic+3) + #hpave = 3*0.06+(len(indBaseCuts)+1)*0.04 + #if indEffCuts!=None: + # hpave += (len(indEffCuts)+2)*0.04 + + allLines = [ ] + allLines.append(('hdr','Variable(s)')) + allLines.append(('txt',variable)) + allLines.append(('hdr','Basic selection')) + for l in cutLines(indBaseCuts): + if l!="": + allLines.append(('txt',l)) + if indEffCuts!=None: + allLines.append(('hdr','Efficiency selection')) + for l in cutLines(indEffCuts): + if l!="": + allLines.append(('txt',l)) + + hdrPixels = 16 + txtPixels = 14 + #pave = ROOT.TPaveText(0.05,max(0,1.0-hpave),0.95,1.0) + hpave = ((hdrPixels+txtPixels+2)*len([ x for x in allLines if x[0]=='hdr' ]) + \ + (txtPixels+2)*len([ x for x in allLines if x[0]=='txt' ])) / (ROOT.gPad.VtoPixel(0.)-ROOT.gPad.VtoPixel(1.)) + if hpave>1.: + scale = 1./hpave + hpave = 1 + else: + scale = 1 + + pave = ROOT.TPaveText(0.,1.-hpave,1.,1.) + pave.SetBorderSize(0) + pave.SetFillStyle(0) + pave.SetTextFont(43) + #pave.SetMargin(0) + #pave.SetTextSize(0.04) + pave.SetTextSize(int(scale*txtPixels)) + pave.SetTextAlign(13) + nhdr = 0 + for x,y in allLines: + if x=='hdr': + if nhdr>0: + t = pave.AddText("") + t.SetTextSize(int(scale*txtPixels)) + t = pave.AddText(y) + t.SetTextFont(63) + t.SetTextSize(int(scale*hdrPixels)) + nhdr += 1 + else: + t = pave.AddText(" "+y) + t.SetTextSize(int(scale*txtPixels)) + #drawString(pave,"Variable(s)",variable) + #drawString(pave,"Basic selection") + #drawCuts(pave,"Basic selection",indBaseCuts) + #if indEffCuts!=None: + # #t = pave.AddText("") + # drawString(pave,"Efficiency selection") + # drawCuts(pave,"Efficiency selection",indEffCuts) + + pave.Draw() + ROOT.gPad.Update() + return pave + +#def createHistoByDef(rdf,hDef,extraCuts,varMaskCombs): +def createHistoByDef(rdfWrapper,hDef,extraCuts): + histos = { } + + for mType in range(23,26): + #print("Checking mType",mType,"for",hDef.name) + # + # draw histogram? + # + if hDef.vetoMType(mType): + continue + is1D = hDef.getParameter('yNbins',mType)==None + is2D = hDef.getParameter('yNbins',mType)!=None and hDef.getParameter('zNbins',mType)==None + is3D = hDef.getParameter('yNbins',mType)!=None and hDef.getParameter('zNbins',mType)!=None + isProfile = hDef.getParameter('profile',mType)!=None and hDef.getParameter('profile',mType) + effCuts = hDef.getParameter('effCuts',mType) + variable = hDef.getParameter('variable',mType) + ## + #cnv.cd(ic) + #ROOT.gPad.SetGridx(1) + #ROOT.gPad.SetGridy(1) + #if not is1D: + # ROOT.gPad.SetRightMargin(0.125) + hName = hDef.getParameter('histogramName',mType) + str(mType) + hTitle = hDef.getParameter('histogramTitle',mType) + " module type " +str(mType) + nbx = hDef.getParameter('xNbins',mType) + xmin = hDef.getParameter('xMin',mType) + xmax = hDef.getParameter('xMax',mType) + #print("Starting for ",hDef.name,hName,hTitle) + if is1D and ( not isProfile ): + cuts = cutString(extraCuts,hDef.getParameter('baseCuts',mType),"moduleType=="+str(mType)) +#!# rdf,varMask = varMaskCombs(rdf,variable,cuts) + varMask = rdfWrapper.defineVarMask(variable,cuts) + model = (hName+"_1",hName+"_1",nbx,xmin,xmax) + histos[mType] = [ rdfWrapper().Histo1D(model,varMask), None, None, None ] + if effCuts!=None: + cuts = cutString(extraCuts,hDef.getParameter('baseCuts',mType),"moduleType=="+str(mType),effCuts) +#!# rdf,varMask = varMaskCombs(rdf,variable,cuts) + varMask = rdfWrapper.defineVarMask(variable,cuts) + model = (hName+"_2",hName+"_2",nbx,xmin,xmax) + histos[mType][1] = rdfWrapper().Histo1D(model,varMask) + #sys.exit() + elif isProfile: + ymin = hDef.getParameter('yMin',mType) + ymax = hDef.getParameter('yMax',mType) + v2,v1 = variable.split(":") + cuts = cutString(extraCuts,hDef.getParameter('baseCuts',mType),"moduleType=="+str(mType)) +#!# rdf,varMask1 = varMaskCombs(rdf,v1,cuts) + varMask1 = rdfWrapper.defineVarMask(v1,cuts) +#!# rdf,varMask2 = varMaskCombs(rdf,v2,cuts) + varMask2 = rdfWrapper.defineVarMask(v2,cuts) + model = (hName+"_1",hName+"_1",nbx,xmin,xmax,ymin,ymax) + histos[mType] = [ rdfWrapper().Profile1D(model,varMask1,varMask2), None, None, None ] + #histos[mType] = [ rdf.Profile1D((hName+"_1",hName+"_1",nbx,xmin,xmax,ymin,ymax,'S'), \ + # v1+"["+cutString(extraCuts,hDef.getParameter('baseCuts',mType),"moduleType=="+str(mType))+"]",v2), \ + # None, None, None ] + elif is2D: + nby = hDef.getParameter('yNbins',mType) + ymin = hDef.getParameter('yMin',mType) + ymax = hDef.getParameter('yMax',mType) + v2,v1 = variable.split(":") + cuts = cutString(extraCuts,hDef.getParameter('baseCuts',mType),"moduleType=="+str(mType)) +#!# rdf,varMask1 = varMaskCombs(rdf,v1,cuts) + varMask1 = rdfWrapper.defineVarMask(v1,cuts) +#!# rdf,varMask2 = varMaskCombs(rdf,v2,cuts) + varMask2 = rdfWrapper.defineVarMask(v2,cuts) + model = (hName+"_1",hName+"_1",nbx,xmin,xmax,nby,ymin,ymax) + histos[mType] = [ rdfWrapper().Histo2D(model,varMask1,varMask2), None, None, None ] + if effCuts!=None: + cuts = cutString(extraCuts,hDef.getParameter('baseCuts',mType),"moduleType=="+str(mType),effCuts) +#!# rdf,varMask1 = varMaskCombs(rdf,v1,cuts) + varMask1 = rdfWrapper.defineVarMask(v1,cuts) +#!# rdf,varMask2 = varMaskCombs(rdf,v2,cuts) + varMask2 = rdfWrapper.defineVarMask(v2,cuts) + model = (hName+"_2",hName+"_2",nbx,xmin,xmax,nby,ymin,ymax) + histos[mType][1] = rdfWrapper().Histo2D(model,varMask1,varMask2) + elif is3D: + nby = hDef.getParameter('yNbins',mType) + ymin = hDef.getParameter('yMin',mType) + ymax = hDef.getParameter('yMax',mType) + nbz = hDef.getParameter('zNbins',mType) + zmin = hDef.getParameter('zMin',mType) + zmax = hDef.getParameter('zMax',mType) + v3,v2,v1 = variable.split(":") + cuts = cutString(extraCuts,hDef.getParameter('baseCuts',mType),"moduleType=="+str(mType)) +#!# rdf,varMask1 = varMaskCombs(rdf,v1,cuts) + varMask1 = rdfWrapper.defineVarMask(v1,cuts) +#!# rdf,varMask2 = varMaskCombs(rdf,v2,cuts) + varMask2 = rdfWrapper.defineVarMask(v2,cuts) +#!# rdf,varMask3 = varMaskCombs(rdf,v3,cuts) + varMask3 = rdfWrapper.defineVarMask(v3,cuts) + model = (hName+"_1",hName+"_1",nbx,xmin,xmax,nby,ymin,ymax,nbz,zmin,zmax) + histos[mType] = [ rdfWrapper().Histo3D(model,varMask1,varMask2,varMask3), None, None, None ] + assert effCuts==None + #print("Ending for ",hDef.name,hName,hTitle) + + return histos + + + +def fillHistoByDef(tree,hDef,extraCuts,histos): + #histos = { } + + savedDir = ROOT.gDirectory + ROOT.gROOT.cd() + + ic = 0 + for mType in range(23,26): + #print("Checking mType",mType,"for",hDef.name) + ic += 1 + # + # draw histogram? + # + if hDef.vetoMType(mType): + continue + is1D = hDef.getParameter('yNbins',mType)==None + is2D = hDef.getParameter('yNbins',mType)!=None and hDef.getParameter('zNbins',mType)==None + is3D = hDef.getParameter('yNbins',mType)!=None and hDef.getParameter('zNbins',mType)!=None + isProfile = hDef.getParameter('profile',mType)!=None and hDef.getParameter('profile',mType) + effCuts = hDef.getParameter('effCuts',mType) + variable = hDef.getParameter('variable',mType) + ## + #cnv.cd(ic) + hName = hDef.getParameter('histogramName',mType) + str(mType) + hTitle = hDef.getParameter('histogramTitle',mType) + " module type " +str(mType) + nbx = hDef.getParameter('xNbins',mType) + xmin = hDef.getParameter('xMin',mType) + xmax = hDef.getParameter('xMax',mType) + #print("Starting for ",hDef.name,hName,hTitle) + if is1D and ( not isProfile ): + if effCuts!=None: + histos[mType][3] = ROOT.TEfficiency(histos[mType][1].GetValue(),histos[mType][0].GetValue()) + histos[mType][3].SetMarkerStyle(20) + else: + # always keep final histogram in 4th position + histos[mType][3] = histos[mType][0] + elif isProfile: + # always keep final histogram in 4th position + histos[mType][3] = histos[mType][0] + elif is2D: + if effCuts!=None: + # always keep final histogram in 4th position + histos[mType][3] = histos[mType][1].Divide(histos[mType][0].GetValue()) + else: + # always keep final histogram in 4th position + histos[mType][3] = histos[mType][0] + elif is3D: + # always keep final histogram in 4th position + histos[mType][3] = histos[mType][0] + #print("Ending for ",hDef.name,hName,hTitle) + + savedDir.cd() + return histos + + +def drawHistoByDef(histos,hDef,logY=False,logZ=False,same=False): + result = { 'cnv' : None, 'histos' : histos, 'pave' : None } + + savedDir = ROOT.gDirectory + ROOT.gROOT.cd() + cnvName = hDef.getParameter('canvasName') + cnv = ROOT.TCanvas(cnvName,cnvName,1500,800) + result['cnv'] = cnv + #cnv.Divide(3,2) + result['pads'] = divideRatios(cnv,3,2,[1./3.,1./3.,1/3.],[2/3.,1/3.]) + + #print(hDef['canvasName'],'is',is1D) + + ic = 0 + #hEffVs = { } + for mType in range(23,26): + #print("Checking mType",mType,"for",hDef.name) + ic += 1 + cnv.cd(ic) + #padNameUp = cnvName+str(mType)+'up' + #result[padNameUp] = ROOT.TPad(padNameUp,padNameUp,(ic-1)/3.,1/3.,ic/3.,1.) + #result[padNameUp].Draw() + # + # draw histogram? + # + if hDef.vetoMType(mType): + continue + is1D = hDef.getParameter('yNbins',mType)==None + is2D = hDef.getParameter('yNbins',mType)!=None and hDef.getParameter('zNbins',mType)==None + is3D = hDef.getParameter('yNbins',mType)!=None and hDef.getParameter('zNbins',mType)!=None + isProfile = hDef.getParameter('profile',mType)!=None and hDef.getParameter('profile',mType) + effCuts = hDef.getParameter('effCuts',mType) + variable = hDef.getParameter('variable',mType) + # + #result[padNameUp].cd() + ROOT.gPad.SetGridx(1) + ROOT.gPad.SetGridy(1) + if not is1D: + ROOT.gPad.SetRightMargin(0.125) + hName = hDef.getParameter('histogramName',mType) + str(mType) + hTitle = hDef.getParameter('histogramTitle',mType) + " module type " +str(mType) + + xtitle = hDef.getParameter('xTitle',mType) if hDef.getParameter('xTitle',mType) else variable + nbx = hDef.getParameter('xNbins',mType) + xmin = hDef.getParameter('xMin',mType) + xmax = hDef.getParameter('xMax',mType) + + ytitle = hDef.getParameter('yTitle',mType) if hDef.getParameter('yTitle',mType) else "" + ztitle = hDef.getParameter('zTitle',mType) if hDef.getParameter('zTitle',mType) else "" + if is1D and ( not isProfile ): + ymin = hDef.getParameter('yMin',mType) if hDef.getParameter('yMin',mType)!=None else 0. + ymax = hDef.getParameter('yMax',mType) if hDef.getParameter('yMax',mType)!=None else 1.05 + if effCuts!=None: + if not same: + histos[mType][2] = ROOT.gPad.DrawFrame(xmin,ymin,xmax,ymax) + histos[mType][2].SetTitle(hTitle) + histos[mType][2].GetXaxis().SetTitle(xtitle) + histos[mType][2].GetYaxis().SetTitle(ytitle) + histos[mType][3] = ROOT.TEfficiency(histos[mType][1].GetValue(),histos[mType][0].GetValue()) + histos[mType][3].SetMarkerStyle(20) + histos[mType][3].Draw("same Z") + else: + histos[mType][0].SetTitle(hTitle) + histos[mType][0].GetXaxis().SetTitle(xtitle) + histos[mType][0].GetYaxis().SetTitle(ytitle) + histos[mType][0].Draw("same" if same else "") + fitFunc = hDef.getParameter('fit') + if fitFunc!=None: + histos[mType][0].Fit(fitFunc,"Q","same") + elif isProfile: + histos[mType][0].SetTitle(hTitle) + histos[mType][0].GetXaxis().SetTitle(xtitle) + histos[mType][0].GetYaxis().SetTitle(ytitle) + histos[mType][0].SetMarkerSize(0.5) + #histos[mType][0].SetFillColor(ROOT.TColor.GetColorBright(ROOT.kGray)) + #histos[mType][0].SetFillColor(ROOT.kGray) + histos[mType][0].Draw("same" if same else "") + elif is2D: + assert not same + zmin = hDef.getParameter('zMin',mType) if hDef.getParameter('zMin',mType)!=None else 0. + zmax = hDef.getParameter('zMax',mType) if hDef.getParameter('zMax',mType)!=None else 1.05 + if effCuts!=None: + histos[mType][1].SetTitle(hTitle) + histos[mType][1].GetXaxis().SetTitle(xtitle) + histos[mType][1].GetYaxis().SetTitle(ytitle) + histos[mType][1].SetMinimum(zmin) + histos[mType][1].SetMaximum(zmax) + histos[mType][1].Draw("ZCOL") + else: + histos[mType][0].SetTitle(hTitle) + histos[mType][0].GetXaxis().SetTitle(xtitle) + histos[mType][0].GetYaxis().SetTitle(ytitle) + histos[mType][0].Draw("ZCOL") + elif is3D: + assert not same + zmin = hDef.getParameter('zMin',mType) if hDef.getParameter('zMin',mType)!=None else 0. + zmax = hDef.getParameter('zMax',mType) if hDef.getParameter('zMax',mType)!=None else 1.05 + assert effCuts==None + histos[mType][0].SetTitle(hTitle) + histos[mType][0].GetXaxis().SetTitle(xtitle) + histos[mType][0].GetYaxis().SetTitle(ytitle) + histos[mType][0].GetZaxis().SetTitle(ztitle) + histos[mType][0].Draw() + if logY or hDef.getParameter('logY',mType): + ROOT.gPad.SetLogy(1) + if is2D and ( logZ or hDef.getParameter('logZ',mType) ): + ROOT.gPad.SetLogz(1) + ROOT.gPad.Update() + + #cnv.cd() + #padNameDown = cnvName+str(mType)+'down' + #result[padNameDown] = ROOT.TPad(padNameDown,padNameDown,(ic-1)/3.,0.,ic/3.,1/3.) + #result[padNameDown].Draw() + result['pave'+str(mType)] = drawCutPave(cnv,ic,hDef.getParameter('variable',mType), \ + cutString(extraCuts,hDef.getParameter('baseCuts',mType)), \ + cutString(hDef.getParameter('effCuts',mType))) + + #for k,v in result.items(): + # print(k,v) + savedDir.cd() + return result + +def addHistogram(varString,cuts,effCuts=None,name='userHist'): + extraHDict = { } + # split into string defining the variable(s) and (1 or 2) axis definition(s) + fields1 = varString.split(";") + assert len(fields1)<=3 + extraHDict['variable'] = fields1[0] + #extraHDict['canvasName'] = "cEffArg" + #extraHDict['histogramName'] = "hEffArg" + extraHDict['histogramTitle'] = name + # x-axis + fields2 = fields1[1].split(",") + assert len(fields2)==3 + extraHDict['xNbins'] = int(fields2[0]) + extraHDict['xMin'] = float(fields2[1]) + extraHDict['xMax'] = float(fields2[2]) + # check for info on y axis (== presence of 2nd variable) + if len(fields1)==3: + assert ":" in extraHDict['variable'] + fields3 = fields1[2].split(",") + extraHDict['yNbins'] = int(fields3[0]) + extraHDict['yMin'] = float(fields3[1]) + extraHDict['yMax'] = float(fields3[2]) + extraHDict['xTitle'] = extraHDict['variable'].split(":")[1] + extraHDict['yTitle'] = extraHDict['variable'].split(":")[0] + else: + extraHDict['yMin'] = 0. + extraHDict['yMax'] = 1.05 + extraHDict['xTitle'] = extraHDict['variable'] + extraHDict['yTitle'] = 'efficiency' if effCuts!=None else 'events/bin' + extraHDict['baseCuts'] = cuts + if effCuts!=None: + extraHDict['effCuts'] = effCuts + return HistogramDefinition(name,extraHDict) + + + +parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter) +parser.add_argument('--definitions', '-d', help='python module with dictionaries defining efficiency histograms', \ + type=str, default=None) +parser.add_argument('--histogram', \ + help='definition of extra histogram (format ;,,[;,,)', \ + action='append', type=str, default=[]) +parser.add_argument('--batch', '-b', help='ROOT batch mode', action='store_true', default=False) +parser.add_argument('--dxMax', help='max. local dx for efficiency plots', type=float, default=0.0075) +parser.add_argument('--cuts', '-c', help="basic cut string", type=str, default="") +parser.add_argument('--effCuts', '-e', help="basic cut string", type=str, default=None) +parser.add_argument('--output', '-o', help='output directory for graphic output', type=str, default=None) +parser.add_argument('--formats', help='comma-separated list of extensions for output files', \ + type=str, default="pdf,png") +parser.add_argument('--sampleName', help='sample label for output', type=str, default=None) +parser.add_argument('--fitResiduals', '-f', \ + help='comma-separated list of names of histogram sets with residuals to be fit', \ + type=str, default=None) +parser.add_argument('--selectedHistograms', help='comma-separated names of histogram definitions to be used', \ + type=str, default='*') +parser.add_argument('--vetoedHistograms', help='comma-separated names of histogram definitions not to be used', + type=str, default='') +parser.add_argument('--logY', help='use log scale for y axis', action='store_true', default=False) +parser.add_argument('--logZ', help='use log scale for z axis', action='store_true', default=False) +parser.add_argument('--printTree', '-p', help='print TTree contents', action='store_true', default=False) +parser.add_argument('--listHistograms', '-l', help='list predefined and selected histograms', \ + action='store_true', default=False) +parser.add_argument('--zone', '-z', help='restrict to zone in OT (barrel, tilted, endcap)', type=str, \ + choices=['barrel','tilted','endcap'], default=None) +parser.add_argument('file', help='input file', type=str, nargs='+', default=None) +args = parser.parse_args() +outputFormats = [ ] +if args.output!=None: + assert os.path.isdir(args.output) + print("***",args.output) + outputFormats = [ x.strip() for x in args.formats.strip().split(",") ] +fitResiduals = args.fitResiduals.split(",") if args.fitResiduals else [ ] +selectedHistoNames = args.selectedHistograms.split(",") +vetoedHistoNames = args.vetoedHistograms.split(",") +# +# add cut for zone definition? +# +if args.zone=="barrel": + zoneCuts = "detNormalT>0.99" +elif args.zone=="endcap": + zoneCuts = "detNormalT<0.01" +elif args.zone=="tilted": + zoneCuts = "detNormalT>0.05&&detNormalT<0.095" +else: + zoneCuts = "" +args.cuts = cutString(args.cuts,zoneCuts) + +# +# load histogram definitions +# +allVDefs,allHDefs = loadConfiguration(args.definitions,selectedHistoNames,vetoedHistoNames) +if args.listHistograms: + hnames = sorted(allHDefs.allHistoNames) + for hn in hnames: + print(hn) + sys.exit() + +for ih,h in enumerate(args.histogram): + allHDefs.add(addHistogram(h,args.cuts,args.effCuts,name="userH"+str(ih+1))) + + +#extraCuts = "abs(particleType)==13" +#extraCuts = "tof<12.5" +extraCuts = args.cuts + +if args.batch: + ROOT.gROOT.SetBatch(1) +ROOT.gROOT.ProcessLine(".L setTDRStyle.C") +ROOT.gROOT.ProcessLine(".L floatMod.C+") +ROOT.setTDRStyle() +ROOT.gStyle.SetOptStat(0) +ROOT.gStyle.SetOptFit(1) +ROOT.gStyle.SetTitleFont(42) +ROOT.gStyle.SetTitleFontSize(0.03) +ROOT.gStyle.SetTitleX(0.12) +ROOT.gStyle.SetTitleY(1.00) +ROOT.gStyle.SetTitleAlign(13) +ROOT.gStyle.SetTitleBorderSize(0) +#ROOT.gStyle.SetTitleFillColor(0) +ROOT.gStyle.SetOptTitle(1) +#tf = ROOT.TFile(args.file[0]) +#simHitTree = simHitTree = tf.Get("analysis").Get("SimHitTree") +simHitTree = ROOT.TChain("analysis/SimHitTree") +for fn in args.file: + simHitTree.Add(fn) +#print(type(tchain)) +#tchain.Print() +if args.printTree: + simHitTree.Print() + sys.exit() +# +# create RDataFrame and define additional columns +# +#simHitRDF = ROOT.RDataFrame(simHitTree) +simHitRDF = ROOT.RDataFrame("analysis/SimHitTree",args.file) +for k,v in allVDefs.items(): + simHitRDF = simHitRDF.Define(k,v) + +canvases = [ ] +histos = { } +paves = [ ] +#varMaskCombs = VarMaskCombinations() +simHitRDFW = RDFWrapper(simHitRDF) + +allHistos = { } +for cName in allHDefs.canvasNames(): + allHistos[cName] = {} + for hName in allHDefs.byCanvas[cName]: + print("Processing histogram",hName,"in canvas",cName) + try: + allHistos[cName][hName] = \ + createHistoByDef(simHitRDFW,allHDefs.byCanvas[cName][hName],extraCuts) +# simHitRDF,allHistos[cName][hName] = \ +# createHistoByDef(simHitRDF,allHDefs.byCanvas[cName][hName],extraCuts,varMaskCombs) + except: + print("Exception for",cName,hName) + raise +# keep reference to RDF with all definitions +simHitRDF = simHitRDFW() + +allObjects = [ ] +for cName in allHDefs.canvasNames(): + same = False + cHistos = { } + for hName in allHDefs.byCanvas[cName]: + print("Processing histogram",hName,"in canvas",cName) + #print(allHistos[cName][hName]) + fillHistoByDef(simHitTree,allHDefs.byCanvas[cName][hName],extraCuts,allHistos[cName][hName]) + allObjects.append(drawHistoByDef(allHistos[cName][hName],allHDefs.byCanvas[cName][hName], \ + logY=args.logY,logZ=args.logZ,same=same)) + same = True + if args.output!=None: + c = allObjects[-1]['cnv'] + #for c in [ x['cnv'] for x in allObjects ]: + basename = os.path.join(args.output,c.GetName()) + if args.sampleName!=None: + basename += "_" + args.sampleName + if args.zone!=None: + basename += "_" + args.zone + print(basename) + for fmt in outputFormats: + c.SaveAs(basename+"."+fmt) + +# yMin = min([ x.GetMinimum() for x in cHistos +#sys.exit() + diff --git a/DrawHits/drawHitsResRDF.yaml b/DrawHits/drawHitsResRDF.yaml new file mode 100644 index 0000000..b403ed7 --- /dev/null +++ b/DrawHits/drawHitsResRDF.yaml @@ -0,0 +1,150 @@ +# +# definitions of new variables +# +variables: + pathX: 'path.fCoordinates.fX' + pathY: 'path.fCoordinates.fY' + pathZ: 'path.fCoordinates.fZ' + dxdz: 'pathX/pathZ' + alpha: 'atan2(abs(pathX),abs(pathZ))/3.1415*180' + localPosX: 'localPos.fCoordinates.fX' + localPosY: 'localPos.fCoordinates.fY' + rhLocalPosX: 'rhLocalPos.fCoordinates.fX' + rhLocalPosY: 'rhLocalPos.fCoordinates.fY' + rhLocalErrX: 'rhLocalErr.fCoordinates.fX' + rhLocalErrY: 'rhLocalErr.fCoordinates.fY' + globRho2: 'globalDir.fCoordinates.fX*globalDir.fCoordinates.fX+globalDir.fCoordinates.fY*globalDir.fCoordinates.fY' + globR2: 'globRho2+globalDir.fCoordinates.fZ*globalDir.fCoordinates.fZ' + globPt: 'pabs*sqrt(sqrt(globRho2/globR2))' + detNormalT: 'detNormal.fCoordinates.fX*detNormal.fCoordinates.fX+detNormal.fCoordinates.fY*detNormal.fCoordinates.fY' +# +# +# +histograms: + + pullX: + histogramTitle: 'pulls (x)' + variable: '(localPosX-rhLocalPosX)/rhLocalErrX' + baseCuts: 'tof<12.5&&hasRecHit>0' + xTitle: 'pull local x' + xNbins: 200 + xMax: 5 + xMin: -5 + fit: 'gaus' + logY: true + pullXW1: + histogramTitle: 'pulls (x)' + variable: '(localPosX-rhLocalPosX)/rhLocalErrX' + baseCuts: 'tof<12.5&&hasRecHit>0&&clusterSize==1' + xTitle: 'pull local x' + xNbins: 200 + xMax: 5 + xMin: -5 + fit: 'gaus' + logY: true + pullXW2: + histogramTitle: 'pulls (x)' + variable: '(localPosX-rhLocalPosX)/rhLocalErrX' + baseCuts: 'tof<12.5&&hasRecHit>0&&clusterSize==2' + xTitle: 'pull local x' + xNbins: 200 + xMax: 5 + xMin: -5 + fit: 'gaus' + logY: true +# +# +# + resX: + histogramTitle: 'residuals (x)' + variable: 'localPosX-rhLocalPosX' + baseCuts: 'tof<12.5&&hasRecHit>0' + xTitle: '#Delta x [cm]' + xNbins: 300 + xMax: 0.075 + xMin: -0.075 + fit: 'gaus' + logY: true +# +# +# + resX1: + histogramTitle: 'residuals (x) - cluster size 1' + variable: 'localPosX-rhLocalPosX' + baseCuts: 'tof<12.5&&hasRecHit>0&&clusterSize==1' + xTitle: '#Delta x [cm]' + xNbins: 300 + xMax: 0.075 + xMin: -0.075 + fit: 'gaus' + logY: true +# +# +# + resX2: + histogramTitle: 'residuals (x) - cluster size 2' + variable: 'localPosX-rhLocalPosX' + baseCuts: 'tof<12.5&&hasRecHit>0&&clusterSize==2' + xTitle: '#Delta x [cm]' + xNbins: 800 + xMax: 0.040 + xMin: -0.040 + fit: 'gaus' + logY: true +# +# +# + resX3P: + histogramTitle: 'residuals (x) - cluster size >2' + variable: 'localPosX-rhLocalPosX' + baseCuts: 'tof<12.5&&hasRecHit>0&&clusterSize>2' + xTitle: '#Delta x [cm]' + xNbins: 300 + xMax: 0.075 + xMin: -0.075 + logY: true +# +# +# + res2DX1: + histogramTitle: 'track dx/dz vs. residuals (x) - cluster size 1' + variable: 'pathX/pathZ:localPosX-rhLocalPosX' + baseCuts: 'tof<12.5&&hasRecHit>0&&clusterSize==1' + xTitle: '#Delta x [cm]' + xNbins: 300 + xMax: 0.075 + xMin: -0.075 + yNbins: 21 + yMin: -1.05 + yMax: 1.05 + logZ: true +# +# +# + res2DX2: + histogramTitle: 'track dz/dx vs. residuals (x) - cluster size 2' + variable: 'pathX/pathZ:localPosX-rhLocalPosX' + baseCuts: 'tof<12.5&&hasRecHit>0&&clusterSize==2' + xTitle: '#Delta x [cm]' + xNbins: 800 + xMax: 0.040 + xMin: -0.040 + yNbins: 21 + yMin: -1.05 + yMax: 1.05 + logZ: true +# +# +# + resY: + histogramTitle: 'pulls (y)' + variable: '(localPosY-rhLocalPosY)/rhLocalErrY' + baseCuts: 'tof<12.5&&hasRecHit>0' + xTitle: 'pull local y' + xNbins: 200 + xMax: 5 + xMin: -5 + logY: true +# +# +# diff --git a/DrawHits/drawHitsTmp.yaml b/DrawHits/drawHitsTmp.yaml index 67873c4..3915a9f 100644 --- a/DrawHits/drawHitsTmp.yaml +++ b/DrawHits/drawHitsTmp.yaml @@ -1,7 +1,7 @@ # # # -effAlphaLoose: +effAlphaTight: baseCuts: 'tof<12.5&&pabs>0.3' effCuts: 'hasRecHit>0&&abs(localPos.x()-rhLocalPos.x())<0.0075' histogramTitle: 'Efficiency vs. alpha(xz) tight' @@ -16,7 +16,7 @@ effAlphaLoose: # # # -effAlphaTight: +effAlphaLoose: baseCuts: 'tof<12.5&&pabs>0.3' effCuts: 'hasRecHit>0&&abs(localPos.x()-rhLocalPos.x())<0.1' histogramTitle: 'Efficiency vs. alpha(xz) loose' @@ -421,3 +421,244 @@ widthVsPath: yNbins: 25 yMin: 0.0 yMax: 25 + +res3DXVsDxDzModX: + histogramTitle: 'residuals 3D (x)' + baseCuts: 'tof<12.5&&hasRecHit>0' + xTitle: '#Delta x [#mum]' + xNbins: 100 + xMin: -500 + xMax: 500 + zTitle: 'local dx/dz' + zNbins: 21 + zMin: -1.05 + zMax: 1.05 + mType23: + variable: 'path.x()/path.z():floatMod(10000*localPos.x(),100):10000*(localPos.x()-rhLocalPos.x())' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 20 + yMin: 0 + yMax: 100 + mType24: + variable: 'path.x()/path.z():floatMod(10000*localPos.x(),100):10000*(localPos.x()-rhLocalPos.x())' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 20 + yMin: 0 + yMax: 100 + mType25: + variable: 'path.x()/path.z():floatMod(10000*localPos.x(),90):10000*(localPos.x()-rhLocalPos.x())' + yTitle: 'x modulo 90 #mum [#mum]' + yNbins: 18 + yMin: 0 + yMax: 90 + +res3DXVsDxDz: + histogramTitle: 'residuals 3D (x)' + baseCuts: 'tof<12.5&&hasRecHit>0' + xTitle: '#Delta x [#mum]' + xNbins: 100 + xMin: -200 + xMax: 200 + zTitle: 'local dx/dz' + zNbins: 21 + zMin: -1.05 + zMax: 1.05 + mType23: + variable: 'path.x()/path.z():floatMod(10000*localPos.x(),100):10000*(localPos.x()-rhLocalPos.x())' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType24: + variable: 'path.x()/path.z():floatMod(10000*localPos.x(),100):10000*(localPos.x()-rhLocalPos.x())' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType25: + variable: 'path.x()/path.z():floatMod(10000*localPos.x(),90):10000*(localPos.x()-rhLocalPos.x())' + yTitle: 'x modulo 90 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 90 + +res3DXVsDxDzW1: + histogramTitle: 'residuals 3D (x)' + baseCuts: 'tof<12.5&&hasRecHit>0&&clusterSize==1' + xTitle: '#Delta x [#mum]' + xNbins: 100 + xMin: -200 + xMax: 200 + zTitle: 'local dx/dz' + zNbins: 21 + zMin: -1.05 + zMax: 1.05 + mType23: + variable: 'path.x()/path.z():floatMod(10000*localPos.x(),100):10000*(localPos.x()-rhLocalPos.x())' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType24: + variable: 'path.x()/path.z():floatMod(10000*localPos.x(),100):10000*(localPos.x()-rhLocalPos.x())' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType25: + variable: 'path.x()/path.z():floatMod(10000*localPos.x(),90):10000*(localPos.x()-rhLocalPos.x())' + yTitle: 'x modulo 90 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 90 + +res3DXVsDxDzW2: + histogramTitle: 'residuals 3D (x)' + baseCuts: 'tof<12.5&&hasRecHit>0&&clusterSize==2' + xTitle: '#Delta x [#mum]' + xNbins: 100 + xMin: -200 + xMax: 200 + zTitle: 'local dx/dz' + zNbins: 21 + zMin: -1.05 + zMax: 1.05 + mType23: + variable: 'path.x()/path.z():floatMod(10000*localPos.x(),100):10000*(localPos.x()-rhLocalPos.x())' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType24: + variable: 'path.x()/path.z():floatMod(10000*localPos.x(),100):10000*(localPos.x()-rhLocalPos.x())' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType25: + variable: 'path.x()/path.z():floatMod(10000*localPos.x(),90):10000*(localPos.x()-rhLocalPos.x())' + yTitle: 'x modulo 90 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 90 + +pull3DXVsDxDz: + histogramTitle: 'pulls 3D (x)' + baseCuts: 'tof<12.5&&hasRecHit>0' + xTitle: '#Delta x [#mum]' + xNbins: 100 + xMin: -10 + xMax: 10 + zTitle: 'local dx/dz' + zNbins: 21 + zMin: -1.05 + zMax: 1.05 + mType23: + variable: 'path.x()/path.z():floatMod(10000*localPos.x(),100):(localPos.x()-rhLocalPos.x())/rhLocalErr.x()' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType24: + variable: 'path.x()/path.z():floatMod(10000*localPos.x(),100):(localPos.x()-rhLocalPos.x())/rhLocalErr.x()' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType25: + variable: 'path.x()/path.z():floatMod(10000*localPos.x(),90):(localPos.x()-rhLocalPos.x())/rhLocalErr.x()' + yTitle: 'x modulo 90 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 90 + +pull3DXVsDxDzW1: + histogramTitle: 'pulls 3D (x)' + baseCuts: 'tof<12.5&&hasRecHit>0&&clusterSize==1' + xTitle: '#Delta x [#mum]' + xNbins: 100 + xMin: -10 + xMax: 10 + zTitle: 'local dx/dz' + zNbins: 21 + zMin: -1.05 + zMax: 1.05 + mType23: + variable: 'path.x()/path.z():floatMod(10000*localPos.x(),100):(localPos.x()-rhLocalPos.x())/rhLocalErr.x()' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType24: + variable: 'path.x()/path.z():floatMod(10000*localPos.x(),100):(localPos.x()-rhLocalPos.x())/rhLocalErr.x()' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType25: + variable: 'path.x()/path.z():floatMod(10000*localPos.x(),90):(localPos.x()-rhLocalPos.x())/rhLocalErr.x()' + yTitle: 'x modulo 90 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 90 + +pull3DXVsDxDzW2: + histogramTitle: 'pulls 3D (x)' + baseCuts: 'tof<12.5&&hasRecHit>0&&clusterSize==2' + xTitle: '#Delta x [#mum]' + xNbins: 100 + xMin: -10 + xMax: 10 + zTitle: 'local dx/dz' + zNbins: 21 + zMin: -1.05 + zMax: 1.05 + mType23: + variable: 'path.x()/path.z():floatMod(10000*localPos.x(),100):(localPos.x()-rhLocalPos.x())/rhLocalErr.x()' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType24: + variable: 'path.x()/path.z():floatMod(10000*localPos.x(),100):(localPos.x()-rhLocalPos.x())/rhLocalErr.x()' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType25: + variable: 'path.x()/path.z():floatMod(10000*localPos.x(),90):(localPos.x()-rhLocalPos.x())/rhLocalErr.x()' + yTitle: 'x modulo 90 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 90 + +effDxDz: + variable: 'path.x()/path.z()' + histogramTitle: 'efficiency (dx/dz)' + baseCuts: 'tof<12.5' + effCuts: 'hasRecHit>0' + xTitle: 'local dx/dz' + xNbins: 21 + xMin: -1.05 + xMax: 1.05 + +effDxDzW1: + variable: 'path.x()/path.z()' + histogramTitle: 'efficiency, 1-strip clusters (dx/dz)' + baseCuts: 'tof<12.5' + effCuts: 'hasRecHit>0&&clusterSize==1' + xTitle: 'local dx/dz' + xNbins: 21 + xMin: -1.05 + xMax: 1.05 + +effDxDzW2: + variable: 'path.x()/path.z()' + histogramTitle: 'efficiency, 2-strip clusters (dx/dz)' + baseCuts: 'tof<12.5' + effCuts: 'hasRecHit>0&&clusterSize==2' + xTitle: 'local dx/dz' + xNbins: 21 + xMin: -1.05 + xMax: 1.05 + diff --git a/DrawHits/drawHitsTmpRDF.yaml b/DrawHits/drawHitsTmpRDF.yaml new file mode 100644 index 0000000..22f9909 --- /dev/null +++ b/DrawHits/drawHitsTmpRDF.yaml @@ -0,0 +1,955 @@ +# +# sensor dimensions: +# +# 2S: 10 x 10 cm2; 2 x 1016 strips at 90 μm x 5 cm +# PS: 5 x 10 cm2; 2 x 960 strips at 100 μm x 2.5 cm +# PSp: 5 x 10 cm2; 32 x 960 pixel at 100 μm x 1.5 mm + +# +# definitions of new variables +# +variables: + pathX: 'path.fCoordinates.fX' + pathY: 'path.fCoordinates.fY' + pathZ: 'path.fCoordinates.fZ' + dxdz: 'pathX/pathZ' + alpha: 'atan2(abs(pathX),abs(pathZ))/3.1415*180' + localPosX: 'localPos.fCoordinates.fX' + localPosY: 'localPos.fCoordinates.fY' + rhLocalPosX: 'rhLocalPos.fCoordinates.fX' + rhLocalPosY: 'rhLocalPos.fCoordinates.fY' + rhLocalErrX: 'rhLocalErr.fCoordinates.fX' + rhLocalErrY: 'rhLocalErr.fCoordinates.fY' + localPosXMod90: 'floatMod(10000*localPos.fCoordinates.fX,90)' + localPosXMod100: 'floatMod(10000*localPos.fCoordinates.fX,100)' + localPosYMod500: 'floatMod(10*localPosY,50)/10.' + localPosYMod250: 'floatMod(10*localPosY,25)/10.' + localPosYMod15: 'floatMod(100*localPosY,15)/10' + globRho2: 'globalDir.fCoordinates.fX*globalDir.fCoordinates.fX+globalDir.fCoordinates.fY*globalDir.fCoordinates.fY' + globR2: 'globRho2+globalDir.fCoordinates.fZ*globalDir.fCoordinates.fZ' + globPt: 'pabs*sqrt(sqrt(globRho2/globR2))' + detNormalT: 'detNormal.fCoordinates.fX*detNormal.fCoordinates.fX+detNormal.fCoordinates.fY*detNormal.fCoordinates.fY' +# +# definitions of canvases and histograms +# +histograms: +# +# +# + effAlphaTight: + baseCuts: 'tof<12.5&&pabs>0.3' + effCuts: 'hasRecHit>0&&abs(localPosX-rhLocalPosX)<0.0075' + histogramTitle: 'Efficiency vs. alpha(xz) tight' + variable: 'alpha' + xMax: 100.0 + xMin: 0 + xNbins: 100 + xTitle: 'alpha(xz) [deg]' + yMax: 1.05 + yMin: 0.5 + yTitle: 'efficiency' +# +# +# + effAlphaLoose: + baseCuts: 'tof<12.5&&pabs>0.3' + effCuts: 'hasRecHit>0&&abs(localPosX-rhLocalPosX)<0.1' + histogramTitle: 'Efficiency vs. alpha(xz) loose' + variable: 'alpha' + xMax: 100.0 + xMin: 0 + xNbins: 100 + xTitle: 'alpha(xz) [deg]' + yMax: 1.05 + yMin: 0.5 + yTitle: 'efficiency' +# +# +# + effVsModX: + baseCuts: 'tof<12.5&&pabs>0.3' + effCuts: 'hasRecHit>0&&abs(localPosX-rhLocalPosX)<0.0075' + yMin: 0 + yMax: 1.05 + mType23: + histogramTitle: 'Efficiency vs. x % 100 #mu m' + variable: 'localPosXMod100' + xTitle: 'x modulo 100 #mu m [#mu m]' + xNbins: 220 + xMax: 105 + xMin: -5 + mType24: + histogramTitle: 'Efficiency vs. x % 100 #mu m' + variable: 'localPosXMod100' + xTitle: 'x modulo 100 #mu m [#mu m]' + xNbins: 220 + xMax: 105 + xMin: -5 + mType25: + histogramTitle: 'Efficiency vs. x % 90 #mu m' + variable: 'localPosXMod90' + xTitle: 'x modulo 90 #mu m [#mu m]' + xNbins: 200 + xMax: 95 + xMin: -5 + + effVsModX1: + baseCuts: 'tof<12.5&&pabs>0.3&&globPt>0.5&&hasRecHit>0' + effCuts: 'clusterSize==1' + yMin: 0 + yMax: 1.05 + mType23: + histogramTitle: 'Efficiency vs. x % 100 #mu m' + variable: 'localPosXMod100' + xTitle: 'x modulo 100 #mu m [#mu m]' + xNbins: 220 + xMax: 105 + xMin: -5 + mType24: + histogramTitle: 'Efficiency vs. x % 100 #mu m' + variable: 'localPosXMod100' + xTitle: 'x modulo 100 #mu m [#mu m]' + xNbins: 220 + xMax: 105 + xMin: -5 + mType25: + histogramTitle: 'Efficiency vs. x % 90 #mu m' + variable: 'localPosXMod90' + xTitle: 'x modulo 90 #mu m [#mu m]' + xNbins: 200 + xMax: 95 + xMin: -5 + + effVsModX2: + baseCuts: 'tof<12.5&&pabs>0.3&&globPt>0.5&&hasRecHit>0' + effCuts: 'clusterSize==2' + yMin: 0 + yMax: 1.05 + mType23: + histogramTitle: 'Efficiency vs. x % 100 #mu m' + variable: 'localPosXMod100' + xTitle: 'x modulo 100 #mu m [#mu m]' + xNbins: 220 + xMax: 105 + xMin: -5 + mType24: + histogramTitle: 'Efficiency vs. x % 100 #mu m' + variable: 'localPosXMod100' + xTitle: 'x modulo 100 #mu m [#mu m]' + xNbins: 220 + xMax: 105 + xMin: -5 + mType25: + histogramTitle: 'Efficiency vs. x % 90 #mu m' + variable: 'localPosXMod90' + xTitle: 'x modulo 90 #mu m [#mu m]' + xNbins: 200 + xMax: 95 + xMin: -5 + + effVsModX3p: + baseCuts: 'tof<12.5&&pabs>0.3&&globPt>0.5&&hasRecHit>0' + effCuts: 'clusterSize>2' + yMin: 0 + yMax: 1.05 + mType23: + histogramTitle: 'Efficiency vs. x % 100 #mu m' + variable: 'localPosXMod100' + xTitle: 'x modulo 100 #mu m [#mu m]' + xNbins: 220 + xMax: 105 + xMin: -5 + mType24: + histogramTitle: 'Efficiency vs. x % 100 #mu m' + variable: 'localPosXMod100' + xTitle: 'x modulo 100 #mu m [#mu m]' + xNbins: 220 + xMax: 105 + xMin: -5 + mType25: + histogramTitle: 'Efficiency vs. x % 90 #mu m' + variable: 'localPosXMod90' + xTitle: 'x modulo 90 #mu m [#mu m]' + xNbins: 200 + xMax: 95 + xMin: -5 +# +# +# + effVsModY: + baseCuts: 'tof<12.5&&pabs>0.3' + effCuts: 'hasRecHit>0&&abs(localPosX-rhLocalPosX)<0.0075' + yMin: 0 + yMax: 1.05 + mType23: + histogramTitle: 'Efficiency vs. y % 1.5 mm' + variable: 'localPosYMod15' + xTitle: 'y modulo 1.5mm [mm]' + xNbins: 400 + xMax: 1.75 + xMin: -0.25 + mType24: + histogramTitle: 'Efficiency vs. y % 2.5 cm' + variable: 'localPosYMod250' + xTitle: 'y modulo 2.5cm [cm]' + xNbins: 600 + xMax: 2.75 + xMin: -0.25 + mType25: + histogramTitle: 'Efficiency vs. y % 5 cm' + variable: 'localPosYMod500' + xTitle: 'y modulo 5cm [cm]' + xNbins: 550 + xMax: 5.25 + xMin: -0.25 +# +# +# + effX: + baseCuts: 'tof<12.5&&pabs>0.3' + effCuts: 'hasRecHit>0&&abs(localPosX-rhLocalPosX)<0.0075' + histogramTitle: 'Efficiency 1D' + variable: 'localPosX' + xTitle: 'local x [cm]' + yMax: 1.05 + yMin: 0.8 + yTitle: 'efficiency' + mType23: + xMax: 5.5 + xMin: -5.5 + xNbins: 550 + mType24: + xMax: 5.0 + xMin: -5.0 + xNbins: 250 + mType25: + xMax: 5 + xMin: -5 + xNbins: 250 +# +# +# + effY: + variable: 'localPosY' + baseCuts: 'tof<12.5&&pabs>0.3' + effCuts: 'hasRecHit>0&&abs(localPosX-rhLocalPosX)<0.0075' + histogramTitle: 'Efficiency Y 1D' + xTitle: 'local y [cm]' + xNbins: 550 + xMax: 5.5 + xMin: -5.5 + yTitle: 'efficiency' + yMax: 1.05 + yMin: 0.8 +# +# +# + pathX: + histogramTitle: 'path (x)' + variable: '10000*abs(pathX)' + baseCuts: 'tof<12.5&&pabs>0.3' + xTitle: 'SimHit path (x) [#mum]' + xNbins: 200 + xMax: 1000 + xMin: 0 +# +# +# + pathY: + histogramTitle: 'path (y)' + variable: '10000*abs(pathY)' + baseCuts: 'tof<12.5&&pabs>0.3' + xTitle: 'SimHit path (y [#mum])' + xNbins: 200 + xMax: 2000 + xMin: 0 +# +# +# + pullX: + histogramTitle: 'pulls (x)' + variable: '(localPosX-rhLocalPosX)/rhLocalErrX' + baseCuts: 'tof<12.5&&hasRecHit>0' + xTitle: 'pull local x' + xNbins: 200 + xMax: 5 + xMin: -5 + logY: true +# +# +# + resX: + histogramTitle: 'residuals (x)' + variable: 'localPosX-rhLocalPosX' + baseCuts: 'tof<12.5&&hasRecHit>0' + xTitle: '#Delta x [cm]' + xNbins: 300 + xMax: 0.075 + xMin: -0.075 + logY: true +# +# +# + resX1: + histogramTitle: 'residuals (x) - cluster size 1' + variable: 'localPosX-rhLocalPosX' + baseCuts: 'tof<12.5&&hasRecHit>0&&clusterSize==1' + xTitle: '#Delta x [cm]' + xNbins: 300 + xMax: 0.075 + xMin: -0.075 + logY: true +# +# +# + resX2: + histogramTitle: 'residuals (x) - cluster size 2' + variable: 'localPosX-rhLocalPosX' + baseCuts: 'tof<12.5&&hasRecHit>0&&clusterSize==2' + xTitle: '#Delta x [cm]' + xNbins: 300 + xMax: 0.075 + xMin: -0.075 + logY: true +# +# +# + resX3P: + histogramTitle: 'residuals (x) - cluster size >2' + variable: 'localPosX-rhLocalPosX' + baseCuts: 'tof<12.5&&hasRecHit>0&&clusterSize>2' + xTitle: '#Delta x [cm]' + xNbins: 300 + xMax: 0.075 + xMin: -0.075 + logY: true +# +# +# + resY: + histogramTitle: 'pulls (y)' + variable: '(localPosY-rhLocalPosY)/rhLocalErrY' + baseCuts: 'tof<12.5&&hasRecHit>0' + xTitle: 'pull local y' + xNbins: 200 + xMax: 5 + xMin: -5 + logY: true +# +# +# + stdRes2D: + histogramTitle: 'my resolution' + variable: 'abs(pathX):localPosX-rhLocalPosX' + baseCuts: 'tof<12.5&&hasRecHit>0' + xTitle: 'residual x [cm]' + xNbins: 200 + xMin: -0.1 + xMax: 0.1 + yTitle: 'dx [cm]' + yNbins: 20 + yMin: 0.0 + yMax: 0.1 +# +# +# + widthVsMod2DX: + baseCuts: 'tof<12.5&&pabs>0.3' + yTitle: 'cluster width' + yNbins: 10 + yMin: 0 + yMax: 10 + mType23: + histogramTitle: 'width vs x modula 100 #mu m' + variable: 'clusterSize:localPosXMod100' + xTitle: 'x modulo 100 #mu m [#mu m]' + xNbins: 110 + xMin: -5 + xMax: 105 + mType24: + histogramTitle: 'width vs x modula 100 #mu m' + variable: 'clusterSize:localPosXMod100' + xTitle: 'x modulo 100 #mu m [#mu m]' + xNbins: 110 + xMin: -5 + xMax: 105 + mType25: + histogramTitle: 'width vs x modula 90 #mu m' + variable: 'clusterSize:localPosXMod90' + xTitle: 'x modulo 90 #mu m [#mu m]' + xNbins: 100 + xMin: -5 + xMax: 95 +# +# +# + widthVsModX: + baseCuts: 'tof<12.5&&pabs>0.3' + profile: true + yTitle: 'cluster width' + yMin: 0 + yMax: 25 + mType23: + histogramTitle: 'width vs x modula 100 #mu m' + variable: 'clusterSize:localPosXMod100' + xTitle: 'x modulo 100 #mu m [#mu m]' + xNbins: 110 + xMin: -5 + xMax: 105 + mType24: + histogramTitle: 'width vs x modula 100 #mu m' + variable: 'clusterSize:localPosXMod100' + xTitle: 'x modulo 100 #mu m [#mu m]' + xNbins: 110 + xMin: -5 + xMax: 105 + mType25: + histogramTitle: 'width vs x modula 90 #mu m' + variable: 'clusterSize:localPosXMod90' + xTitle: 'x modulo 90 #mu m [#mu m]' + xNbins: 100 + xMin: -5 + xMax: 95 +# +# +# + widthVsModY: + baseCuts: 'tof<12.5&&pabs>0.3' + yTitle: 'cluster width' + yMin: 0 + yMax: 25 + profile: true + mType23: + histogramTitle: 'Width vs. y % 1.5 mm' + variable: 'clusterSize:localPosYMod15' + xTitle: 'y modulo 1.5mm [mm]' + xNbins: 200 + xMin: -0.25 + xMax: 1.75 + mType24: + histogramTitle: 'Width vs. y % 2.5 cm' + variable: 'clusterSize:localPosYMod250' + xTitle: 'y modulo 2.5cm [cm]' + xNbins: 300 + xMin: -0.25 + xMax: 2.75 + mType25: + histogramTitle: 'Width vs. y % 5 cm' + variable: 'clusterSize:localPosYMod500' + xTitle: 'y modulo 5cm [cm]' + xNbins: 275 + xMin: -0.25 + xMax: 5.25 +# +# +# + widthVsPath: + histogramTitle: 'width vs path' + variable: 'clusterSize:abs(pathX)' + baseCuts: 'tof<12.5&&hasRecHit>0' + profile: true + xTitle: 'SimHit path (x)' + xNbins: 100 + xMin: 0 + xMax: 0.2 + yTitle: 'clustersize' + yNbins: 25 + yMin: 0.0 + yMax: 25 + + res3DXVsDxDzModX: + histogramTitle: 'residuals 3D (x)' + baseCuts: 'tof<12.5&&hasRecHit>0' + xTitle: '#Delta x [#mum]' + xNbins: 100 + xMin: -500 + xMax: 500 + zTitle: 'local dx/dz' + zNbins: 21 + zMin: -1.05 + zMax: 1.05 + mType23: + variable: 'dxdz:localPosXMod100:10000*(localPosX-rhLocalPosX)' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 20 + yMin: 0 + yMax: 100 + mType24: + variable: 'dxdz:localPosXMod100:10000*(localPosX-rhLocalPosX)' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 20 + yMin: 0 + yMax: 100 + mType25: + variable: 'dxdz:localPosXMod90:10000*(localPosX-rhLocalPosX)' + yTitle: 'x modulo 90 #mum [#mum]' + yNbins: 18 + yMin: 0 + yMax: 90 + + res3DXVsDxDz: + histogramTitle: 'residuals 3D (x)' + baseCuts: 'tof<12.5&&hasRecHit>0' + xTitle: '#Delta x [#mum]' + xNbins: 100 + xMin: -200 + xMax: 200 + zTitle: 'local dx/dz' + zNbins: 21 + zMin: -1.05 + zMax: 1.05 + mType23: + variable: 'dxdz:localPosXMod100:10000*(localPosX-rhLocalPosX)' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType24: + variable: 'dxdz:localPosXMod100:10000*(localPosX-rhLocalPosX)' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType25: + variable: 'dxdz:localPosXMod90:10000*(localPosX-rhLocalPosX)' + yTitle: 'x modulo 90 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 90 + + res3DXVsDxDzW1: + histogramTitle: 'residuals 3D (x)' + baseCuts: 'tof<12.5&&hasRecHit>0&&clusterSize==1' + xTitle: '#Delta x [#mum]' + xNbins: 100 + xMin: -200 + xMax: 200 + zTitle: 'local dx/dz' + zNbins: 21 + zMin: -1.05 + zMax: 1.05 + mType23: + variable: 'dxdz:localPosXMod100:10000*(localPosX-rhLocalPosX)' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType24: + variable: 'dxdz:localPosXMod100:10000*(localPosX-rhLocalPosX)' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType25: + variable: 'dxdz:localPosXMod90:10000*(localPosX-rhLocalPosX)' + yTitle: 'x modulo 90 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 90 + + res3DXVsDxDzW2: + histogramTitle: 'residuals 3D (x)' + baseCuts: 'tof<12.5&&hasRecHit>0&&clusterSize==2' + xTitle: '#Delta x [#mum]' + xNbins: 100 + xMin: -200 + xMax: 200 + zTitle: 'local dx/dz' + zNbins: 21 + zMin: -1.05 + zMax: 1.05 + mType23: + variable: 'dxdz:localPosXMod100:10000*(localPosX-rhLocalPosX)' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType24: + variable: 'dxdz:localPosXMod100:10000*(localPosX-rhLocalPosX)' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType25: + variable: 'dxdz:localPosXMod90:10000*(localPosX-rhLocalPosX)' + yTitle: 'x modulo 90 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 90 + + pull3DXVsDxDz: + histogramTitle: 'pulls 3D (x)' + baseCuts: 'tof<12.5&&hasRecHit>0' + xTitle: '#Delta x [#mum]' + xNbins: 100 + xMin: -10 + xMax: 10 + zTitle: 'local dx/dz' + zNbins: 21 + zMin: -1.05 + zMax: 1.05 + mType23: + variable: 'dxdz:localPosXMod100:(localPosX-rhLocalPosX)/rhLocalErrX' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType24: + variable: 'dxdz:localPosXMod100:(localPosX-rhLocalPosX)/rhLocalErrX' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType25: + variable: 'dxdz:localPosXMod90:(localPosX-rhLocalPosX)/rhLocalErrX' + yTitle: 'x modulo 90 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 90 + + pull3DXVsDxDzW1: + histogramTitle: 'pulls 3D (x)' + baseCuts: 'tof<12.5&&hasRecHit>0&&clusterSize==1' + xTitle: '#Delta x [#mum]' + xNbins: 100 + xMin: -10 + xMax: 10 + zTitle: 'local dx/dz' + zNbins: 21 + zMin: -1.05 + zMax: 1.05 + mType23: + variable: 'dxdz:localPosXMod100:(localPosX-rhLocalPosX)/rhLocalErrX' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType24: + variable: 'dxdz:localPosXMod100:(localPosX-rhLocalPosX)/rhLocalErrX' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType25: + variable: 'dxdz:localPosXMod90:(localPosX-rhLocalPosX)/rhLocalErrX' + yTitle: 'x modulo 90 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 90 + + pull3DXVsDxDzW2: + histogramTitle: 'pulls 3D (x)' + baseCuts: 'tof<12.5&&hasRecHit>0&&clusterSize==2' + xTitle: '#Delta x [#mum]' + xNbins: 100 + xMin: -10 + xMax: 10 + zTitle: 'local dx/dz' + zNbins: 21 + zMin: -1.05 + zMax: 1.05 + mType23: + variable: 'dxdz:localPosXMod100:(localPosX-rhLocalPosX)/rhLocalErrX' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType24: + variable: 'dxdz:localPosXMod100:(localPosX-rhLocalPosX)/rhLocalErrX' + yTitle: 'x modulo 100 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 100 + mType25: + variable: 'dxdz:localPosXMod90:(localPosX-rhLocalPosX)/rhLocalErrX' + yTitle: 'x modulo 90 #mum [#mum]' + yNbins: 1 + yMin: 0 + yMax: 90 + + effDxDz: + variable: 'dxdz' + histogramTitle: 'efficiency (dx/dz)' + baseCuts: 'tof<12.5' + effCuts: 'hasRecHit>0' + xTitle: 'local dx/dz' + xNbins: 21 + xMin: -1.05 + xMax: 1.05 + + effDxDzW1: + variable: 'dxdz' + histogramTitle: 'efficiency, 1-strip clusters (dx/dz)' + baseCuts: 'tof<12.5' + effCuts: 'hasRecHit>0&&clusterSize==1' + xTitle: 'local dx/dz' + xNbins: 21 + xMin: -1.05 + xMax: 1.05 + + effDxDzW2: + variable: 'dxdz' + histogramTitle: 'efficiency, 2-strip clusters (dx/dz)' + baseCuts: 'tof<12.5' + effCuts: 'hasRecHit>0&&clusterSize==2' + xTitle: 'local dx/dz' + xNbins: 21 + xMin: -1.05 + xMax: 1.05 + + nrhVsModX: + baseCuts: 'tof<12.5&&pabs>0.3&&globPt>0.5' + yTitle: 'nr of rechits' + yNbins: 10 + yMin: 0 + yMax: 10 + mType23: + histogramTitle: 'nr of rechits vs. x % 100 #mu m' + variable: 'rhNMatched:localPosXMod100' + xTitle: 'x modulo 100 #mu m [#mu m]' + xNbins: 220 + xMax: 105 + xMin: -5 + mType24: + histogramTitle: 'nr of rechits vs. x % 100 #mu m' + variable: 'rhNMatched:localPosXMod100' + xTitle: 'x modulo 100 #mu m [#mu m]' + xNbins: 220 + xMax: 105 + xMin: -5 + mType25: + histogramTitle: 'nr of rechits vs. x % 90 #mu m' + variable: 'rhNMatched:localPosXMod90' + xTitle: 'x modulo 90 #mu m [#mu m]' + xNbins: 200 + xMax: 95 + xMin: -5 + + nrhEffVsModX1: + baseCuts: 'tof<12.5&&pabs>0.3&&globPt>0.5' + effCuts: 'rhNMatched==1' + yTitle: 'Efficieny (nrh==1)' + yMin: 0 + yMax: 1.05 + mType23: + histogramTitle: 'Efficiency (nrh=1) vs. x % 100 #mu m' + variable: 'localPosXMod100' + xTitle: 'x modulo 100 #mu m [#mu m]' + xNbins: 220 + xMax: 105 + xMin: -5 + mType24: + histogramTitle: 'Efficiency (nrh=1) vs. x % 100 #mu m' + variable: 'localPosXMod100' + xTitle: 'x modulo 100 #mu m [#mu m]' + xNbins: 220 + xMax: 105 + xMin: -5 + mType25: + histogramTitle: 'Efficiency (nrh=1) vs. x % 90 #mu m' + variable: 'localPosXMod90' + xTitle: 'x modulo 90 #mu m [#mu m]' + xNbins: 200 + xMax: 95 + xMin: -5 + + nrhEffVsModX2: + baseCuts: 'tof<12.5&&pabs>0.3&&globPt>0.5' + effCuts: 'rhNMatched==2' + yTitle: 'Efficieny (nrh==2)' + yMin: 0 + yMax: 1.05 + mType23: + histogramTitle: 'Efficiency (nrh=2) vs. x % 100 #mu m' + variable: 'localPosXMod100' + xTitle: 'x modulo 100 #mu m [#mu m]' + xNbins: 220 + xMax: 105 + xMin: -5 + mType24: + histogramTitle: 'Efficiency (nrh=2) vs. x % 100 #mu m' + variable: 'localPosXMod100' + xTitle: 'x modulo 100 #mu m [#mu m]' + xNbins: 220 + xMax: 105 + xMin: -5 + mType25: + histogramTitle: 'Efficiency (nrh=2) vs. x % 90 #mu m' + variable: 'localPosXMod90' + xTitle: 'x modulo 90 #mu m [#mu m]' + xNbins: 200 + xMax: 95 + xMin: -5 + + nrhEffVsModX3p: + baseCuts: 'tof<12.5&&pabs>0.3&&globPt>0.5' + effCuts: 'rhNMatched>2' + yTitle: 'Efficieny (nrh>2)' + yMin: 0 + yMax: 1.05 + mType23: + histogramTitle: 'Efficiency (nrh>2) vs. x % 100 #mu m' + variable: 'localPosXMod100' + xTitle: 'x modulo 100 #mu m [#mu m]' + xNbins: 220 + xMax: 105 + xMin: -5 + mType24: + histogramTitle: 'Efficiency (nrh>2) vs. x % 100 #mu m' + variable: 'localPosXMod100' + xTitle: 'x modulo 100 #mu m [#mu m]' + xNbins: 220 + xMax: 105 + xMin: -5 + mType25: + histogramTitle: 'Efficiency (nrh>2) vs. x % 90 #mu m' + variable: 'localPosXMod90' + xTitle: 'x modulo 90 #mu m [#mu m]' + xNbins: 200 + xMax: 95 + xMin: -5 + + nrhVsModY: + baseCuts: 'tof<12.5&&pabs>0.3&&globPt>0.5' + yTitle: 'nr of rechits' + yNbins: 10 + yMin: 0 + yMax: 10 + variable: 'rhNMatched:localPosY' + histogramTitle: 'nr of rechits vs. local y' + mType23: + histogramTitle: 'nr of rechits vs. y % 1.5 mm' + variable: 'rhNMatched:localPosYMod15' + xTitle: 'y modulo 1.5mm [mm]' + xNbins: 400 + xMax: 1.75 + xMin: -0.25 + mType24: + histogramTitle: 'nr of rechits vs. y % 2.5 cm' + variable: 'rhNMatched:localPosYMod250' + xTitle: 'y modulo 2.5cm [cm]' + xNbins: 600 + xMax: 2.75 + xMin: -0.25 + mType25: + histogramTitle: 'nr of rechits vs. y % 5 cm' + variable: 'rhNMatched:localPosYMod500' + xTitle: 'y modulo 5cm [cm]' + xNbins: 550 + xMax: 5.25 + xMin: -0.25 + + nrhEffVsModY1: + baseCuts: 'tof<12.5&&pabs>0.3&&globPt>0.5' + effCuts: 'rhNMatched==1' + yTitle: 'Efficieny (nrh==1)' + yMin: 0 + yMax: 1.05 + mType23: + histogramTitle: 'Efficiency (nrh=1) vs. y % 1.5 mm' + variable: 'localPosYMod15' + xTitle: 'y modulo 1.5mm [mm]' + xNbins: 400 + xMax: 1.75 + xMin: -0.25 + mType24: + histogramTitle: 'Efficiency (nrh=1) vs. y % 2.5 cm' + variable: 'localPosYMod250' + xTitle: 'y modulo 2.5cm [cm]' + xNbins: 600 + xMax: 2.75 + xMin: -0.25 + mType25: + histogramTitle: 'Efficiency (nrh=1) vs. y % 5 cm' + variable: 'localPosYMod500' + xTitle: 'y modulo 5cm [cm]' + xNbins: 550 + xMax: 5.25 + xMin: -0.25 + + nrhEffVsModY2: + baseCuts: 'tof<12.5&&pabs>0.3&&globPt>0.5' + effCuts: 'rhNMatched==2' + yTitle: 'Efficiency (nrh==2)' + yMin: 0 + yMax: 1.05 + mType23: + histogramTitle: 'Efficiency (nrh=2) vs. y % 1.5 mm' + variable: 'localPosYMod15' + xTitle: 'y modulo 1.5mm [mm]' + xNbins: 400 + xMax: 1.75 + xMin: -0.25 + mType24: + histogramTitle: 'Efficiency (nrh=2) vs. y % 2.5 cm' + variable: 'localPosYMod250' + xTitle: 'y modulo 2.5cm [cm]' + xNbins: 600 + xMax: 2.75 + xMin: -0.25 + mType25: + histogramTitle: 'Efficiency (nrh=2) vs. y % 5 cm' + variable: 'localPosYMod500' + xTitle: 'y modulo 5cm [cm]' + xNbins: 550 + xMax: 5.25 + xMin: -0.25 + + nrhEffVsModY3p: + baseCuts: 'tof<12.5&&pabs>0.3&&globPt>0.5' + effCuts: 'rhNMatched>2' + yTitle: 'Efficiency (nrh>2)' + yMin: 0 + yMax: 1.05 + mType23: + histogramTitle: 'Efficiency (nrh>2) vs. y % 1.5 mm' + variable: 'localPosYMod15' + xTitle: 'y modulo 1.5mm [mm]' + xNbins: 400 + xMax: 1.75 + xMin: -0.25 + mType24: + histogramTitle: 'Efficiency (nrh>2) vs. y % 2.5 cm' + variable: 'localPosYMod250' + xTitle: 'y modulo 2.5cm [cm]' + xNbins: 600 + xMax: 2.75 + xMin: -0.25 + mType25: + histogramTitle: 'Efficiency (nrh>2) vs. y % 5 cm' + variable: 'localPosYMod500' + xTitle: 'y modulo 5cm [cm]' + xNbins: 550 + xMax: 5.25 + xMin: -0.25 + + nrRecHitsAll: + baseCuts: 'tof<12.' + histogramTitle: 'Number of RecHits (all)' + variable: 'rhNMatched' + xTitle: 'number of matched RecHits' + xNbins: 100 + xMax: 99.5 + xMin: -0.5 + logY: true + + nrRecHitsPabs: + baseCuts: 'tof<12.&&pabs>3.' + histogramTitle: 'Number of RecHits (all)' + variable: 'rhNMatched' + xTitle: 'number of matched RecHits' + xNbins: 100 + xMax: 99.5 + xMin: -0.5 + logY: true + + nrRecHitsPt: + baseCuts: 'tof<12.&&pabs>0.3&&globPt>0.5' + histogramTitle: 'Number of RecHits (all)' + variable: 'rhNMatched' + xTitle: 'number of matched RecHits' + xNbins: 100 + xMax: 99.5 + xMin: -0.5 + logY: true + diff --git a/DrawHits/fitRes.py b/DrawHits/fitRes.py new file mode 100644 index 0000000..3f43967 --- /dev/null +++ b/DrawHits/fitRes.py @@ -0,0 +1,779 @@ +import sys,os +from math import sqrt,log +import ROOT +import argparse + +class FitHistogram: + + def interpolate(y,y1,y2,x1=0.,x2=1.): + ''' Linear interpolation between two points (x1,y1) and (x2,y2) + to yield x corresponding to the target value y. + ''' + return x1 + (y-y1)/(y2-y1)*(x2-x1) + + def __init__(self,histogram): + self.hist = histogram + self.graph_ = None + self.cumulativeGraph_ = None + self.cumNormGraph_ = None + self.cumNormTGraph_ = None + self.cumNormSpline_ = None # spline corresponding to normalized cumulative graph + + def graph(self): + ''' Histogram converted to list of (x,y) coordinates with x = bin center and y = bin contents. + Ignores under- / overflow. + ''' + if self.graph_==None: + self.graph_ = [ ] + + for i in range(self.hist.GetNbinsX()): + xbin = self.hist.GetBinLowEdge(i+1) + self.hist.GetBinWidth(i+1)/2. + self.graph_.append( (xbin,self.hist.GetBinContent(i+1)) ) + + return self.graph_ + + def cumulativeGraph(self): + ''' Histogram converted to list of (x,y) coordinates with x = bin center and y = cumulated bin contents. + Ignores under- / overflow. + ''' + if self.cumulativeGraph_==None: + self.cumulativeGraph_ = [ ] + sum = 0. + for i in range(self.hist.GetNbinsX()): + xbin = self.hist.GetBinLowEdge(i+1) + self.hist.GetBinWidth(i+1)/2. + sum += self.hist.GetBinContent(i+1) + self.cumulativeGraph_.append( (xbin,sum) ) + + return self.cumulativeGraph_ + + def cumulativeNormGraph(self): + ''' Histogram converted to list of (x,y) coordinates with x = bin center and y = cumulated bin contents. + Normalized to total histogram contents (including under- / overflow). + ''' + if self.cumNormGraph_==None: + cg = self.cumulativeGraph() + sum = self.hist.GetSumOfWeights() + self.cumNormGraph_ = [ ( x,y/sum ) for x,y in cg ] + + return self.cumNormGraph_ + + def cumulativeNormSpline(self): + ''' Create (TGraph and) TSpline3 from cumulativeNormGraph + ''' + if self.cumNormTGraph_==None: + self.cumNormTGraph_ = ROOT.TGraph() + for x,y in self.cumulativeNormGraph(): + self.cumNormTGraph_.SetPoint(self.cumNormTGraph_.GetN(),x,y) + + if self.cumNormSpline_==None: + self.cumNormSpline_ = ROOT.TSpline3(self.hist.GetName()+"-spline",self.cumNormTGraph_) + + return self.cumNormSpline_ + + + def intersects(self,value,cumulative=False,norm=False,direction=0): + ''' Calculate x-coordinates for intersection(s) of a graph defined by a list of (x,y) points sorted in x + with y==value. Uses linear interpolation. + Arguments: + value ....... target value + cumulative .. if true, use cumulative graph + norm ........ if true, use normalized cumulative graph + direction ... three possible values: 0 = any intersection, +1/-1 = consider only segments + with positive / negative slope. + ''' + # + # Get graph and do basic check + # + graph = None + if cumulative: + graph = self.cumulativeNormGraph() if norm else self.cumulativeGraph() + else: + assert norm==False + graph = self.graph() + + result = [ ] + if len(graph)<2: + return result + # + # loop over adjacent pairs of points + # + x1 = None + y1 = None + for x2,y2 in graph: + # + # start checking at 2nd point + # + if x1!=None: + assert x2>x1 + # + # value in interval? + # + if value>=min(y1,y2) and value<=max(y1,y2): + # check if dy is positive or negative, and compare with required sign of direction + if direction==0 or direction*(y2-y1)>0: + result.append(FitHistogram.interpolate(value,y2,y1,x2,x1)) + # + # move to next point + # + x1 = x2 + y1 = y2 + + return result + + def fwhm(self): + ''' Return lowest / highest x corresponding to ymax/2, and ymax/2 + ''' + # + # target y value (1/2 maximum) + # + y = self.hist.GetBinContent(self.hist.GetMaximumBin())/2. + # + # use first intersection with upward slope + # + xups = self.intersects(y,cumulative=False,direction=1) + print("xups",xups) + xlow = xups[0] if xups else None + # + # use last intersection with downward slope + # + xdowns = self.intersects(y,cumulative=False,direction=-1) + print("xdowns",xdowns) + xhigh = xdowns[-1] if xdowns else None + #dx = self.hist.GetBinWidth(1) + ### + ## preset result (low / high x values) + ## + #xlow = None + #xhigh = None + ## + ## loop over pairs of adjacent histogram bins + ## + #x1 = self.hist.GetBinLowEdge(1) + #y1 = self.hist.GetBinContent(1) + #for i in range(2,self.hist.GetNbinsX()): + # # check if target value between contents of neighbouring bins + # x2 = self.hist.GetBinLowEdge(i) + # y2 = self.hist.GetBinContent(i) + # first = None + # if y>=y1 and yy2: + # first = False + # + # if first!=None: + # # calculate interpolated x value + # x = FitHistogram.interpolate(y,y1,y2,x1,x2) + # # store lowest and highest x corresponding to y + # if first and xlow==None: + # xlow = x + # elif not first: + # xhigh = x + # # move to next bin + # x1 = x2 + # y1 = y2 + + # + # require result ( assumes that first and last bins are < ymax/2 ) and + # correct for bin width / 2 ( assumes that bin value corresponds to center of bin ) + assert xlow!=None and xhigh!=None + return (xlow,xhigh,y) + #return (xlow+dx,xhigh+dx,y) + + def quantile(self,prob): + ''' Return x-value to quantile q. + ''' + result = self.intersects(prob,cumulative=True,norm=True,direction=1) + #print(prob,result) + assert len(result)<2 + + return result[0] if result else None + + def findRootSpline(self,value,eps=0.001): + ''' Find position where spline derived from cumulative NormGraph= value. Assumes that the spline is \ + monotonously increasing. Tolerance is eps*value or cGraph[-1][1]=ymax: + print("findRootSpline: last value <= first value") + return None + # + # Check if inside values spanned by spline + # + if valueymax: + print("findRootSpline: required value outside range") + return None + # + # tolerance for distance between points + # + dxmax = eps*(ymax-ymin) + # + # loop with cutoff in case of non-convergence + # + found = False + for i in range(1000): + if (xh-xl)prob1 + result = ( None, None ) + xmin = self.hist.GetXaxis().GetXmin() + if prob1>0.: + xmin = self.findRootSpline(prob1) + if xmin==None: + return result + xmax = self.hist.GetXaxis().GetXmax() + if prob2<1.: + xmax = self.findRootSpline(prob2) + if xmin==None: + return result + + fitFunc = ROOT.TF1("mygaus","gaus(0)",xmin,xmax) + fitFunc.SetParameter(0,self.hist.GetMaximum()) + fitFunc.SetParameter(1,0.) + fitFunc.SetParLimits(1,-10.,10.) + fitFunc.SetParameter(2,self.hist.GetRMS()) + + fitPtr = self.hist.Fit(fitFunc,"S0") + #fitPtr = self.hist.Fit("gaus","S0","",xmin,xmax) + if ( not fitPtr.IsValid() ) or fitPtr.IsEmpty(): + return result + + #return fitPtr,self.hist.GetFunction("gaus") + return fitPtr,fitFunc + + +class FitCanvas: + + def __init__(self,canvas,mtype): + self.canvas = canvas + self.mtype = mtype + self.pad = None + self.fhist = None + + padName = canvas.GetName()+str(mtype-22) + for o in canvas.GetListOfPrimitives(): + if o.InheritsFrom(ROOT.TPad.Class()) and o.GetName()==padName: + assert self.pad==None + self.pad = o + assert self.pad!=None + + histNameBase = "h"+canvas.GetName()[1:]+str(mtype) + self.fhist = None + for o in self.pad.GetListOfPrimitives(): + if o.InheritsFrom(ROOT.TH1.Class()) and o.GetName().startswith(histNameBase): + assert self.fhist==None + self.fhist = FitHistogram(o) + assert self.fhist!=None + + def histogram(self): + #print("fhist",self.fhist,self.fhist.hist) + return self.fhist.hist + + def fwhm(self): + return self.fhist.fwhm() + +def saveCanvas(canvas,inputName,outputDir,formats): + canvasName = canvas.GetName() + outputBaseName = os.path.splitext(os.path.basename(inputName))[0] + outputName = os.path.join(outputDir,outputBaseName) + if canvasName.split("_")[0]==outputBaseName.split("_")[0]: + canvasName = "_".join(canvasName.split("_")[1:]) + elif canvasName.startswith("c"): + canvasName = canvasName[1].lower() + canvasName[2:] + outputName += "_" + canvasName + for fmt in formats: + print("Saving canvas as",outputName+"."+fmt) + canvas.SaveAs(outputName+"."+fmt) + +def setHistogramMinMax(histo,limits,margins=0.05): + ''' Set minimum / maximum for histogram considering values within a range. + Arguments: + histo ...... histogram (TH1 or TH2) + limits ..... lower/upper limit of the range defining values to be considered + margins .... add this fraction of the min-max range below and above + Return value: + min / max value found within range + ''' + # + # min / max value can't be outside the range + # + vmin = limits[1] + vmax = limits[0] + # + # loop over all bins and updated vmin/vmax + # + for ibx in range(histo.GetNbinsX()): + for iby in range(histo.GetNbinsY()): + c = histo.GetBinContent(ibx+1,iby+1) + if c>=limits[0] and c<=limits[1]: + vmin = min(vmin,c) + vmax = max(vmax,c) + # + # apply min / max and return result + # + if (vmax-vmin)/(limits[1]-limits[0])<0.0001: + vmin,vmax = limits + histo.SetMinimum(vmin-margins*(vmax-vmin)) + histo.SetMaximum(vmax+margins*(vmax-vmin)) + + return (vmin,vmax) + +def defineMeanWidthHistograms(h,nby,ymin,ymax,nbz,zmin,zmax,isigma=None,sigma=None): + if isigma!=None: + y1 = ROOT.TMath.Freq(-sigma) + y2 = ROOT.TMath.Freq(sigma) + hcName = h.GetName()+"Qc"+str(isigma) + hcTitle = h.GetTitle()+" (median)" + hwName = h.GetName()+"Qhw"+str(isigma) + hwTitle = h.GetTitle()+" (#sigma from quantiles)".format(y1,y2) + hwAxis1Title = "median [#mum]" + hwAxis2Title = "#sigma from {:4.1%}/{:4.1%} quantiles [#mum]".format(y1,y2) + else: + y1 = None + y2 = None + hcName = h.GetName()+"Qc fit" + hcTitle = h.GetTitle()+" (fitted mean)" + hwName = h.GetName()+"Qhw fit" + hwTitle = h.GetTitle()+" (fitted sigma)" + hwAxis1Title = "mean [#mum]" + hwAxis2Title = "#sigma from fit [#mum]".format(y1,y2) + if nby==1: + hcentre = ROOT.TH1F(hcName,hcTitle,nbz,zmin,zmax) + hcentre.GetXaxis().SetTitle(h.GetZaxis().GetTitle()) + hcentre.GetYaxis().SetTitle(hwAxis1Title) + hhalfwidth = ROOT.TH1F(hwName,hwTitle,nbz,zmin,zmax) + hhalfwidth.GetXaxis().SetTitle(h.GetZaxis().GetTitle()) + hhalfwidth.GetYaxis().SetTitle(hwAxis2Title) + elif nbz==1: + hcentre = ROOT.TH1F(hcName,hcTitle,nby,ymin,ymax) + hcentre.GetXaxis().SetTitle(h.GetYaxis().GetTitle()) + hcentre.GetYaxis().SetTitle(hwAxis1Title) + hhalfwidth = ROOT.TH2F(hwName,hwTitle,nby,ymin,ymax) + hhalfwidth.GetXaxis().SetTitle(h.GetYaxis().GetTitle()) + hhalfwidth.GetYaxis().SetTitle(hwAxis2Title) + else: + hcentre = ROOT.TH2F(hcName,hctitle,nby,ymin,ymax,nbz,zmin,zmax) + hcentre.GetXaxis().SetTitle(h.GetYaxis().GetTitle()) + hcentre.GetYaxis().SetTitle(h.GetZaxis().GetTitle()) + hcentre.GetZaxis().SetTitle(hwAxis1Title) + hhalfwidth = ROOT.TH2F(hwName,hwTitle,nby,ymin,ymax,nbz,zmin,zmax) + hhalfwidth.GetXaxis().SetTitle(h.GetYaxis().GetTitle()) + hhalfwidth.GetYaxis().SetTitle(h.GetZaxis().GetTitle()) + hhalfwidth.GetZaxis().SetTitle(hwAxis2Title) + hhalfwidth.SetMinimum(0.) + + return hcentre,hhalfwidth + + +parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter) +parser.add_argument('--moduleType', '-m', help='module type', type=int, choices=[23,24,25], default=23) +parser.add_argument('--fwhm', help='determining FWHM', action='store_true', default=False) +parser.add_argument('--slices', '-s', help='fit in slices in y', action='store_true', default=False) +parser.add_argument('--quantile', '-q', help='calculate quantiles corresponding to +- x sigma (can be repeated)', \ + action='append', type=float, default=[]) +parser.add_argument('--dbgQuantiles', help='modX,dxdz pairs defining individual bins for the quantile determination', \ + action='append', type=str, default=[ ]) +parser.add_argument('--dbgAllQuantiles', help='show residuals distributions for all bins', \ + action='store_true', default=False) +parser.add_argument('--output', '-o', help='output directory', type=str, default=None) +parser.add_argument('--formats', help='comma-separated list of extensions for output files', \ + type=str, default="pdf,png") +parser.add_argument('file', help='input file', type=str, nargs=1, default=None) +args = parser.parse_args() + +tf = ROOT.TFile(args.file[0]) +tf.ls() +mainCnv = None +for k in tf.GetListOfKeys(): + o = k.ReadObj() + if o.IsA()==ROOT.TCanvas.Class(): + assert mainCnv==None + mainCnv = k.ReadObj() +assert mainCnv!=None + +dbgQuantPoints = [ ] +for qp in args.dbgQuantiles: + dbgQuantPoints.append( tuple([float(x.strip()) for x in qp.strip().split(",")]) ) + +if args.output!=None: + assert os.path.isdir(args.output) + outputFormats = [ x.strip() for x in args.formats.strip().split(",") ] + +canvases = [ ] +slices = [ ] + +mainCnv.Draw() +fitCanvas = FitCanvas(mainCnv,args.moduleType) +fwhmArrow = ROOT.TArrow() +fwhmArrow.SetLineColor(2) +fwhmArrow.SetLineWidth(2) +quantLines = [ ] +for i in range(len(args.quantile)): + quantLine = ROOT.TLine() + quantLine.SetLineColor(4) + #quantLine.SetLineWidth(2) + quantLine.SetLineStyle(i+1) + quantLines.append(quantLine) +#quantLine.SetLineWidth(2) +if fitCanvas.histogram().GetDimension()==1 and ( not args.slices ): + assert fitCanvas.histogram().GetDimension()==1 + #print(fitCanvas.fwhm()) + x1,x2,y = fitCanvas.fwhm() + print(" = {:6.1f}um, dx = {:6.1f}um, sig = {:6.1f}um ( interval {:6.4f} - {:6.4f}cm )".format( \ + 10000*(x1+x2)/2.,10000*(x2-x1),10000*(x2-x1)/2/sqrt(2*log(2)),x1,x2)) + + fitCanvas.pad.cd() + fwhmArrow.DrawArrow(x1,y,x2,y,0.005,"<>") + + quantiles = [ ] + qlmax = slices[-1].hist.GetMaximum()/10. + for isigma,sigma in enumerate(args.quantile): + qs = [ ] + for sgn in [-1,0,1]: + p = ROOT.TMath.Freq(sgn*sigma) + #qs.append(slices[-1].quantile(p)) + qs.append(slices[-1].findRootSpline(p)) + #print(sigma,sgn,p,qs[-1]) + quantiles.append(qs) + if not ( None in qs ): + quantLines[isigma].DrawLine(qs[0],0.,qs[0],qlmax) + quantLines[0].DrawLine(qs[1],0.,qs[1],qlmax) + quantLines[isigma].DrawLine(qs[2],0.,qs[2],qlmax) + + fitCanvas.pad.Update() + +elif fitCanvas.histogram().GetDimension()==2 and args.slices: + hist2D = fitCanvas.histogram() + assert hist2D.GetDimension()==2 + yaxis = hist2D.GetYaxis() + nby = hist2D.GetNbinsY() + ncol = int(sqrt(nby)) + while ncol*ncol") + + if hProj.GetMaximum()>100: + quantiles = [ ] + qlmax = slices[-1].hist.GetMaximum()/10. + for sigma in args.quantile: + qs = [ ] + for sgn in [-1,0,1]: + p = ROOT.TMath.Freq(sgn*sigma) + #qs.append(slices[-1].quantile(p)) + qs.append(slices[-1].findRootSpline(p)) + #print(sigma,sgn,p,qs[-1]) + quantiles.append(qs) + if not ( None in qs ): + quantLine.DrawLine(qs[0],0.,qs[0],qlmax) + quantLine.DrawLine(qs[1],0.,qs[1],qlmax) + quantLine.DrawLine(qs[2],0.,qs[2],qlmax) + + ROOT.gPad.Update() + +elif fitCanvas.histogram().GetDimension()==3: + # + # calculate summary of quantiles in x + # + h = fitCanvas.histogram() + nby = h.GetNbinsY() + ymin = h.GetYaxis().GetXmin() + ymax = h.GetYaxis().GetXmax() + nbz = h.GetNbinsZ() + zmin = h.GetZaxis().GetXmin() + zmax = h.GetZaxis().GetXmax() + # + # define summary histograms + # + summaries = [ ] + for isigma,sigma in enumerate(args.quantile): +#!# y1 = ROOT.TMath.Freq(-sigma) +#!# y2 = ROOT.TMath.Freq(sigma) +#!# if nby==1: +#!# hcentre = ROOT.TH1F(h.GetName()+"Qc"+str(isigma), \ +#!# h.GetTitle()+" (median)", \ +#!# nbz,zmin,zmax) +#!# hcentre.GetXaxis().SetTitle(h.GetZaxis().GetTitle()) +#!# hcentre.GetYaxis().SetTitle("median [#mum]") +#!# hhalfwidth = ROOT.TH1F(h.GetName()+"Qhw"+str(isigma), \ +#!# h.GetTitle()+" (#sigma from quantiles)".format(y1,y2), \ +#!# nbz,zmin,zmax) +#!# hhalfwidth.GetXaxis().SetTitle(h.GetZaxis().GetTitle()) +#!# hhalfwidth.GetYaxis().SetTitle("#sigma from {:4.1%}/{:4.1%} quantiles [#mum]".format(y1,y2)) +#!# elif nbz==1: +#!# hcentre = ROOT.TH1F(h.GetName()+"Qc"+str(isigma), \ +#!# h.GetTitle()+" (median)", \ +#!# nby,ymin,ymax) +#!# hcentre.GetXaxis().SetTitle(h.GetYaxis().GetTitle()) +#!# hcentre.GetYaxis().SetTitle("median [#mum]") +#!# hhalfwidth = ROOT.TH2F(h.GetName()+"Qhw"+str(isigma), \ +#!# h.GetTitle()+" (#sigma from quantiles)".format(y1,y2), \ +#!# nby,ymin,ymax) +#!# hhalfwidth.GetXaxis().SetTitle(h.GetYaxis().GetTitle()) +#!# hhalfwidth.GetYaxis().SetTitle("#sigma from {:4.1%}/{:4.1%} quantiles [#mum]".format(y1,y2)) +#!# else: +#!# hcentre = ROOT.TH2F(h.GetName()+"Qc"+str(isigma), \ +#!# h.GetTitle()+" (median)", \ +#!# nby,ymin,ymax,nbz,zmin,zmax) +#!# hcentre.GetXaxis().SetTitle(h.GetYaxis().GetTitle()) +#!# hcentre.GetYaxis().SetTitle(h.GetZaxis().GetTitle()) +#!# hcentre.GetZaxis().SetTitle("median [#mum]") + #!# hhalfwidth = ROOT.TH2F(h.GetName()+"Qhw"+str(isigma), \ +#!# h.GetTitle()+" (#sigma from quantiles)".format(y1,y2), \ +#!# nby,ymin,ymax,nbz,zmin,zmax) +#!# hhalfwidth.GetXaxis().SetTitle(h.GetYaxis().GetTitle()) +#!# hhalfwidth.GetYaxis().SetTitle(h.GetZaxis().GetTitle()) +#!# hhalfwidth.GetZaxis().SetTitle("#sigma from {:4.1%}/{:4.1%} quantiles [#mum]".format(y1,y2)) + #!# hhalfwidth.SetMinimum(0.) + hcentre,hhalfwidth = defineMeanWidthHistograms(h,nby,ymin,ymax,nbz,zmin,zmax,isigma,sigma) + summaries.append((hcentre,hhalfwidth)) + hcentre,hhalfwidth = defineMeanWidthHistograms(h,nby,ymin,ymax,nbz,zmin,zmax,None,None) + summaries.append((hcentre,hhalfwidth)) + + dbgBins = set() + if args.dbgAllQuantiles: + # get histograms for all bins + for iy in range(nby): + for iz in range(nbz): + dbgBins.add((iy+1,iz+1)) + else: + for y,z in dbgQuantPoints: + dbgBins.add( (h.GetYaxis().FindBin(y), h.GetZaxis().FindBin(z)) ) + dbgBins = sorted(dbgBins) + + allDbgObjects = { x:[ ] for x in range(len(args.quantile)) } + allFitObjects = [ ] + nDbg = len(dbgBins) + nrow = ncol = int(sqrt(nDbg)) + while nrow*ncol100: + # ensure sufficient statistics + p1 = ROOT.TMath.Freq(-sigma) + if htmp.GetSumOfWeights()>3/p1: + #q1 = fhtmp.quantile(ROOT.TMath.Freq(-sigma)) + #q2 = fhtmp.quantile(ROOT.TMath.Freq(sigma)) + qm1 = fhtmp.findRootSpline(ROOT.TMath.Freq(-sigma)) + q0 = fhtmp.findRootSpline(ROOT.TMath.Freq(0)) + qp1 = fhtmp.findRootSpline(ROOT.TMath.Freq(sigma)) + else: + qm1 = None + q0 = None + qp1 = None + quantiles.append((qm1,q0,qp1)) + hsum = summaries[isigma][0] + if nby==1: + ibin = hsum.GetBin(iz+1) + elif nbz==1: + ibin = hsum.GetBin(iy+1) + else: + ibin = hsum.GetBin(iy+1,iz+1) + if q0!=None: + hsum.SetBinContent(ibin,q0) + else: + print("No median for bins",iy+1,iz+1,"of",htmp.GetName()) + hsum.SetBinContent(ibin,-999999) + if htmp.GetName()=="hPull3DXVsDxDzW223_1_px" and (iy+1)==1 and (iz+1)==3: + print("...",isigma,sigma,qm1,q0,qp1,htmp.GetSumOfWeights()) + hsum = summaries[isigma][1] + if qm1!=None and qp1!=None: + hsum.SetBinContent(ibin,(qp1-qm1)/2./sigma) + + #allDbgObjects[isigma] = [ ] + if (iy+1,iz+1) in dbgBins: + y1 = htmp.GetYaxis().GetBinLowEdge(iy+1) + y2 = y1 + htmp.GetYaxis().GetBinWidth(iy+1) + z1 = htmp.GetZaxis().GetBinLowEdge(iz+1) + z2 = z1 + htmp.GetZaxis().GetBinWidth(iz+1) + + hResName = "hResDbg_mT"+str(args.moduleType)+"_y{:02d}_z{:02d}".format(iy+1,iz+1) + hResTitle = "residual mType "+str(args.moduleType) + if nby>1: + hResTitle += " {:6.3f}1: + hResTitle += " {:6.3f}100: + fitPtr,fitFunc = fhtmp.fitGaus(ROOT.TMath.Freq(-2.),ROOT.TMath.Freq(2.)) + print(fitPtr,fitFunc) + if fitPtr!=None: + fitFunc2 = fitFunc.Clone() + fitFunc2.SetParameter(0,fitFunc.GetParameter(0)/hResDbgMax) + fitFunc2.Draw("same") + allFitObjects.append((fitPtr,fitFunc2)) + dbgObjects[0].Update() + hsum = summaries[-1][0] + if nby==1: + ibin = hsum.GetBin(iz+1) + elif nbz==1: + ibin = hsum.GetBin(iy+1) + else: + ibin = hsum.GetBin(iy+1,iz+1) + hsum.SetBinContent(ibin,fitPtr.Parameter(1)) + hsum.SetBinError(ibin,fitPtr.Error(1)) + hsum = summaries[-1][1] + if nby==1: + ibin = hsum.GetBin(iz+1) + elif nbz==1: + ibin = hsum.GetBin(iy+1) + else: + ibin = hsum.GetBin(iy+1,iz+1) + hsum.SetBinContent(ibin,fitPtr.Parameter(2)) + hsum.SetBinError(ibin,fitPtr.Error(2)) + else: + print("No median for bins",iy+1,iz+1,"of",htmp.GetName()) + hsum.SetBinContent(ibin,-999999) + + + allDbgObjects['canvas'].Update() + + if args.output: + saveCanvas(allDbgObjects['canvas'],args.file[0],args.output,outputFormats) + savePad.cd() + + nsigma = len(args.quantile) + # create from histogram name (remove version number last "_", split of last two digits (mod. type) + cName = "c"+"_".join(h.GetName()[1:].split("_")[:-1]) + cName = cName[:-2] + "_mT" + cName[-2:] + "_MeanWidth" + cName = "cResMeanWidth_mT"+str(args.moduleType) + + c = ROOT.TCanvas(cName,h.GetTitle()+" (mean/width summary)",500*nsigma,1000) + c.Divide(nsigma+1,2) + for i in range(nsigma+1): + hopt = "HIST" if i #include +#include +using RVecF = ROOT::VecOps::RVec; float floatMod(float a, float b) { // @@ -10,3 +12,15 @@ float floatMod(float a, float b) { if ( result < 0. ) result += b; return result; } + + +RVecF floatMod(const RVecF& a, float b) { + // + // apply floatMod to RVec of floats + // + RVecF result = RVecF(a.size()); + for ( unsigned int i=0; i0.3 + alpha_loose: tof<12.5 && pabs>0.3 & hasRecHit>0 & abs(localPosX-rhLocalPosX)<0.1 + alpha_tight: tof<12.5 && pabs>0.3 & hasRecHit>0 & abs(localPosX-rhLocalPosX)<0.0075 + variables: + - alpha + 2d_plots: + effVsModXM23: + selections: + modXM23_base: tof<12.5 && pabs>0.3 && moduleType==23 + modXM23_tight: tof<12.5 && pabs>0.3 && hasRecHit>0 && abs(localPosX-rhLocalPosX)<0.0075 && moduleType==23 + variables: + - localPosXMod100 + 2d_plots: + effVsModXM24: + selections: + modXM24_base: tof<12.5 && pabs>0.3 && moduleType==24 + modXM24_tight: tof<12.5 && pabs>0.3 && hasRecHit>0 && abs(localPosX-rhLocalPosX)<0.0075 && moduleType==24 + variables: + - localPosXMod100 + 2d_plots: + effVsModXM25: + selections: + modXM25_base: tof<12.5 && pabs>0.3 && moduleType==25 + modXM25_tight: tof<12.5 && pabs>0.3 && hasRecHit>0 && abs(localPosX-rhLocalPosX)<0.0075 && moduleType==25 + variables: + - localPosXMod90 + 2d_plots: + effVsModX1M23: + selections: + modX1M23_base: tof<12.5 && pabs>0.3 && pabs*globPt>0.5 && hasRecHit>0 + modX1M23_1: clusterSize==1 + variables: + - localPosXMod100 + 2d_plots: + effVsModX2M23: + selections: + modX2M23_base: tof<12.5 && pabs>0.3 && pabs*globPt>0.5 && hasRecHit>0 + modX2M23_2: clusterSize==2 + variables: + - localPosXMod100 + 2d_plots: + effVsModX3pM23: + selections: + modX3pM23_base: tof<12.5 && pabs>0.3 && pabs*globPt>0.5 && hasRecHit>0 + modX3pM23_3p: clusterSize>2 + variables: + - localPosXMod100 + 2d_plots: + +# effVsModXM24: +# selections: +# modXM24_base: tof<12.5 && pabs>0.3 +# modXM24_tight: tof<12.5 && pabs>0.3 & hasRecHit>0 & abs(localPosX-rhLocalPosX)<0.0075 & moduleType=24 +# variables: +# - localPosXMod100 +# 2d_plots: +# effVsModXM25: +# selections: +# modXM25_base: tof<12.5 && pabs>0.3 +# modXM25_tight: tof<12.5 && pabs>0.3 & hasRecHit>0 & abs(localPosX-rhLocalPosX)<0.0075 & moduleType=25 +# variables: +# - localPosXMod90 +# 2d_plots: + effRecHits: + selections: + dxdz_all: tof<12.5 + dxdz_recHit: tof<12.5 & hasRecHit>0 + variables: + - dxdz + 2d_plots: + + +### PLOT SETTING +# This section defines the name, titles, binning of histograms +# usage: variable_name: [name,title,number of bins,range mininum,range maximum] +plot_setting: + dxdz: [dxdz,'dx/dz (all):local dx/dz:entries / bin',21,-1.05,1.05] + alpha: [alpha,'Efficiency vs. alpha(xz): alpha(xz) [deg]:entries / bin',100,0.,100.] + localPosXMod100: [localPosXMod100,'Efficiency vs. x % 100 #mu m: x modulo 100 #mu m [#mu m]:entries/bin',220,-5,105] + localPosXMod90: [localPosXMod90,'Efficiency vs. x % 90 #mu m: x modulo 90 #mu m [#mu m]:entries/bin',220,-5,95] diff --git a/DrawHits/showHitsConfiguration.py b/DrawHits/showHitsConfiguration.py new file mode 100644 index 0000000..925a7c7 --- /dev/null +++ b/DrawHits/showHitsConfiguration.py @@ -0,0 +1,50 @@ +# +# Show definitions for histograms from a configuration file +# +import sys +from drawHitsConfiguration import * +import argparse + +def printDict(d,indent=0): + for k,v in d.items(): + if type(v)==dict: + print(indent*" "+k+":") + printDict(v,indent+2) + else: + print(indent*" "+k+":",v) + +if __name__=="__main__": + parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter) + parser.add_argument('--isEff', help='Select efficiency plots', action='store_true', default=False) + parser.add_argument('--isNoEff', help='Veto efficiency plots', action='store_true', default=False) + parser.add_argument('--isProfile', help='Select profile histograms', action='store_true', default=False) + parser.add_argument('--isNoProfile', help='Veto profile histograms', action='store_true', default=False) + parser.add_argument('--is1D', help='Select 1D histograms', action='store_true', default=False) + parser.add_argument('--is2D', help='Select 2D histograms', action='store_true', default=False) + parser.add_argument('--is3D', help='Select 3D histograms', action='store_true', default=False) + parser.add_argument('--verbose', '-v', help='Show details for selected histograms', \ + action='store_true', default=False) + parser.add_argument('file', help='yaml defining efficiency histograms', type=str, nargs=1, default=None) + args = parser.parse_args() + + allVDefs,allHDefs = loadConfiguration(args.file[0],['*']) + print() + + for hDef in allHDefs.allDefinitions.values(): + skip = True + if args.is1D and hDef('yNbins')==None: + skip = False + if args.is2D and hDef('yNbins')!=None and hDef('zNbins')==None: + skip = False + if args.is3D and hDef('zNbins')!=None: + skip = False + if skip: + continue + if ( args.isEff and hDef('effCuts')==None ) or ( args.isNoEff and hDef('effCuts')!=None ) : + continue + if ( args.isProfile and hDef('profile')==None ) or ( args.isNoProfile and hDef('profile')!=None ) : + continue + print(hDef.name) + if args.verbose: + printDict(hDef.parameters,2) + diff --git a/HitAnalyzer/interface/RecHitInfo.h b/HitAnalyzer/interface/RecHitInfo.h index 1c1e561..8c95146 100644 --- a/HitAnalyzer/interface/RecHitInfo.h +++ b/HitAnalyzer/interface/RecHitInfo.h @@ -59,6 +59,7 @@ class RecHitInfo { std::vector moduleType; std::vector detNormal; std::vector trackId; + std::vector nSimTracks; }; RecHitData recHitData; diff --git a/HitAnalyzer/src/RecHitInfo.cc b/HitAnalyzer/src/RecHitInfo.cc index d386dfb..9816d69 100644 --- a/HitAnalyzer/src/RecHitInfo.cc +++ b/HitAnalyzer/src/RecHitInfo.cc @@ -41,7 +41,8 @@ void RecHitInfo::setBranches(TTree& tree) { tree.Branch("layer", &recHitData.layer); tree.Branch("moduleType", &recHitData.moduleType); tree.Branch("detNormal", &recHitData.detNormal); - tree.Branch("trackId", &recHitData.trackId); + tree.Branch("trackId", &recHitData.trackId); + tree.Branch("nSimTracks", &recHitData.nSimTracks); }; std::vector RecHitInfo::getSimTrackId( @@ -124,7 +125,7 @@ void RecHitInfo::fillRecHitInfo(const Phase2TrackerRecHit1D& recHit, unsigned in const Phase2TrackerCluster1D* clustIt = &*recHit.cluster(); // Get all the simTracks that form the cluster - std::vector clusterSimTrackIds; + std::set clusterSimTrackIds; for (unsigned int i(0); i < clustIt->size(); ++i) { unsigned int channel(Phase2TrackerDigi::pixelToChannel(clustIt->firstRow() + i, clustIt->column())); std::vector simTrackIds_unselected(getSimTrackId(*pixelSimLinks, detId, channel)); @@ -134,18 +135,9 @@ void RecHitInfo::fillRecHitInfo(const Phase2TrackerRecHit1D& recHit, unsigned in if (istfind != simTracks.end()) simTrackIds.push_back(istId); } - for (unsigned int i = 0; i < simTrackIds.size(); ++i) { - bool add = true; - for (unsigned int j = 0; j < clusterSimTrackIds.size(); ++j) { - // only save simtrackids that are not present yet - if (simTrackIds.at(i) == clusterSimTrackIds.at(j)) - add = false; - } - if (add) - clusterSimTrackIds.push_back(simTrackIds.at(i)); - } + clusterSimTrackIds.insert(simTrackIds.begin(),simTrackIds.end()); } - + // find the closest simhit // this is needed because otherwise you get cases with simhits and clusters being swapped // when there are more than 1 cluster with common simtrackids @@ -158,9 +150,9 @@ void RecHitInfo::fillRecHitInfo(const Phase2TrackerRecHit1D& recHit, unsigned in ++simhitIt) { // check SimHit detId is the same with the RecHit if (rawid == simhitIt->detUnitId()) { - auto it = std::lower_bound(clusterSimTrackIds.begin(), clusterSimTrackIds.end(), simhitIt->trackId()); - // check SimHit track id is included in the cluster - if (it != clusterSimTrackIds.end() && *it == simhitIt->trackId()) { + // auto it = std::lower_bound(clusterSimTrackIds.begin(), clusterSimTrackIds.end(), simhitIt->trackId()); + // // check SimHit track id is included in the cluster + if ( clusterSimTrackIds.find(simhitIt->trackId()) != clusterSimTrackIds.end() ) { trackSimHits.push_back(&*simhitIt); if (!simhit || fabs(simhitIt->localPosition().x() - localPosClu.x()) < minx) { minx = fabs(simhitIt->localPosition().x() - localPosClu.x()); @@ -255,13 +247,14 @@ void RecHitInfo::fillRecHitInfo(const Phase2TrackerRecHit1D& recHit, unsigned in recHitData.Hit_cluster_closestSimHit_local_z.push_back(simhit->localPosition().z()); fillSimHitInfo(*simhit); } + recHitData.nSimTracks.push_back(clusterSimTrackIds.size()); recHitData.Hit_det_rawid.push_back(rawid); recHitData.Hit_cluster_firstStrip.push_back(clustIt->firstStrip()); recHitData.Hit_cluster_firstRow.push_back(clustIt->firstRow()); recHitData.Hit_cluster_column.push_back(clustIt->column()); recHitData.Hit_cluster_edge.push_back(clustIt->edge()); recHitData.Hit_cluster_threshold.push_back(clustIt->threshold()); - + }; void RecHitInfo::clear() { @@ -301,4 +294,5 @@ void RecHitInfo::clear() { recHitData.moduleType.clear(); recHitData.detNormal.clear(); recHitData.trackId.clear(); + recHitData.nSimTracks.clear(); }; diff --git a/HitAnalyzer/src/SimHitInfo.cc b/HitAnalyzer/src/SimHitInfo.cc index 99dfa18..6a24cbd 100644 --- a/HitAnalyzer/src/SimHitInfo.cc +++ b/HitAnalyzer/src/SimHitInfo.cc @@ -148,16 +148,8 @@ SimHitInfo::matchRecHitOnDet(const PSimHit* simHit, const DetId& detId, // find SimTracks contributing to the channel unsigned int channel(Phase2TrackerDigi::pixelToChannel(cluster.firstRow() + i, cluster.column())); std::vector simTrackIds = getSimTrackId(detId,channel); - // std::cout << "SimHitInfo " << simTrackIds.size() << " simTracks" << std::endl; - for ( std::vector::const_iterator ist=simTrackIds.begin(); ist!=simTrackIds.end(); ++ist ) { - // std::cout << "SimHitInfo track id " << *ist << " (simhit track id " << (*simHit).trackId() << " ) " - // << std::endl; - // compare to track id of the SimHit - if ( (*ist)==(*simHit).trackId() ) { - matched = true; - break; - } - } + // // std::cout << "SimHitInfo " << simTrackIds.size() << " simTracks" << std::endl; + matched = std::find(simTrackIds.begin(),simTrackIds.end(),(*simHit).trackId())!=simTrackIds.end(); // if match was found: consider RecHit // std::cout << "SimHitInfo matched : " << matched << std::endl; if ( matched ) break; @@ -258,6 +250,24 @@ void SimHitInfo::fillSimHitInfo(const PSimHit& simHit) { // find all rechits with contributions from the SimTrack related to simHit // RecHitDistancePairs matchedRecHits = matchRecHitOnDet(&simHit,detId,ivrh->second); + // //// + // if ( matchedRecHits.size()>50 ) { + // std::cout << "SimHit with " << matchedRecHits.size() << " matched hits" << std::endl; + // std::cout << " " << "at " << &simHit << " detid " << simHit.detUnitId() << " loc pos " + // << simHit.localPosition().x() << " / " << simHit.localPosition().y() << std::endl; + // std::cout << " " << "entry " << simHit.entryPoint().x() << " / " << simHit.entryPoint().y() + // << " ; exit " << simHit.exitPoint().x() << " / " << simHit.exitPoint().y() << std::endl; + // std::cout << " " << "trackId " << simHit.trackId() << " ; pType " << simHit.particleType() + // << " ; pabs " << simHit.pabs() << std::endl; + // for ( SimHitInfo::RecHitDistancePairs::const_iterator irrr=matchedRecHits.begin(); + // irrr!=matchedRecHits.end(); ++irrr ) { + // std::cout << " RecHit on detid " << irrr->first->geographicalId() << " loc pos " + // << irrr->first->localPosition().x() << " / " << irrr->first->localPosition().y() << std::endl; + // const Phase2TrackerCluster1D& cluster = *(irrr->first->cluster()); + // std::cout << " first row / column " << cluster.firstRow() << " / " << cluster.column() << std::endl; + // } + // } + // //// // store best match (if any) nMatched = matchedRecHits.size(); if ( nMatched>0 ) rechit = matchedRecHits[0].first; diff --git a/HitAnalyzer/test/floatMod.C b/HitAnalyzer/test/floatMod.C deleted file mode 100644 index 1950020..0000000 --- a/HitAnalyzer/test/floatMod.C +++ /dev/null @@ -1,12 +0,0 @@ -#include -#include - -float floatMod(float a, float b) { - // - // change fmod behaviour for negative numbers - // (fmod will return -fmod(|a|,b) for a<0) - // - float result = std::fmod(a,b); - if ( result < 0. ) result += b; - return result; -} diff --git a/Plotter/autoplotter.py b/Plotter/autoplotter.py index 31d33af..d4cb4c9 100644 --- a/Plotter/autoplotter.py +++ b/Plotter/autoplotter.py @@ -2,7 +2,7 @@ import math import uuid import shutil -import Phase2Tracking.Plotter.plotter as p +import Phase2Tracking.Plotter.python.plotter as p import argparse @@ -45,6 +45,8 @@ def getFileList(inputdir): help='The number of jobs to submit for each sample. Could be a list of a single value.') parser.add_argument('--nfiles', type=int, default=-1, help='The number of files per job to submit for each sample. Could be a list of a single value.') +parser.add_argument('--loadMacro', type=str, nargs='*', default=[ ], \ + help='File name for C++ macro(s) to be loaded (can be repeated)') args = parser.parse_args() @@ -137,6 +139,6 @@ def getFileList(inputdir): else: if not os.path.exists(args.filelist): raise Exception("file list not given for plotting") - plotter = p.Plotter(name=args.name,treeName=str(args.treeName),outputDir=args.output,input_filelist=args.filelist,config=args.config,isData=args.data,postfix=args.postfix) + plotter = p.Plotter(name=args.name,treeName=str(args.treeName),outputDir=args.output,input_filelist=args.filelist,config=args.config,isData=args.data,postfix=args.postfix,macros=args.loadMacro) plotter.makeHistFiles() diff --git a/Plotter/filelist.txt b/Plotter/filelist.txt index 6d254cb..d487306 100644 --- a/Plotter/filelist.txt +++ b/Plotter/filelist.txt @@ -1 +1 @@ -root://eos.grid.vbc.ac.at//eos/vbc/experiments/cms/store/user/lian/RelValSingleMuPt10/SingleMuPt10_noPU_2_lian/250916_143902/0000/rechits_tree_1.root \ No newline at end of file +../HitAnalyzer/test/rechits_tree_tt_3.root diff --git a/Plotter/python/plotter.py b/Plotter/python/plotter.py index c299674..8e66213 100644 --- a/Plotter/python/plotter.py +++ b/Plotter/python/plotter.py @@ -11,7 +11,8 @@ ROOT.gStyle.SetOptStat(0) class Plotter: - def __init__(self,name,treeName,outputDir="./",input_filelist=None,config="",isData=False,postfix=""): + def __init__(self,name,treeName,outputDir="./",input_filelist=None,config="",isData=False,postfix="", \ + macros=[]): self.name = name self.treeName = treeName self.outputDir = outputDir @@ -21,6 +22,8 @@ def __init__(self,name,treeName,outputDir="./",input_filelist=None,config="",isD with open(config, "r") as f_cfg: cfg = yaml.load(f_cfg, Loader=yaml.FullLoader) self.cfg = cfg + for m in macros: + ROOT.gROOT.ProcessLine(".L "+m+"+") def getFileList(self): self.filelist = [] @@ -35,11 +38,13 @@ def getFileList(self): print("No files provided as input!") def AddVars(self,d): + print(type(self.cfg)) vars_to_define = ['new_variables'] for v in vars_to_define: if (not v in self.cfg) or (self.cfg[v] is None): continue if self.cfg[v] is not None: + print(v,type(self.cfg[v])) for newvar in self.cfg[v]: if isinstance(self.cfg[v][newvar],list): formatstr = [self.cfg[self.cfg[v][newvar][i]] for i in range(1,len(self.cfg[v][newvar]))] @@ -47,6 +52,7 @@ def AddVars(self,d): elif isinstance(self.cfg[v][newvar],str): var_define = self.cfg[v][newvar] d = d.Define(newvar,var_define) + print("*** Define AddVars",newvar,var_define) return d def AddVarsWithSelection(self,d): @@ -59,8 +65,10 @@ def AddVarsWithSelection(self,d): for v in variables: if selections[sel]: d = d.Define(v+sel,"{0}[{1}]".format(v,selections[sel])) + print("*** Define AddVarsWithSelection1",sel,v+sel,"{0}[{1}]".format(v,selections[sel])) else: d = d.Define(v+sel,"{0}".format(v)) + print("*** Define AddVarsWithSelection",sel,v+sel,"{0}".format(v)) if ('nm1' in self.cfg['objects'][obj]) and (self.cfg['objects'][obj]['nm1']): nm1s = self.cfg['objects'][obj]['nm1'] cutstr_objsel = "" @@ -75,16 +83,19 @@ def AddVarsWithSelection(self,d): cutstr = "&&".join(cutstrs) cutstr = cutstr_objsel + "({})".format(cutstr) d = d.Define(nm1s[i][0]+sel+'_nm1',"{0}[{1}]".format(nm1s[i][0],cutstr)) + print("*** Define AddVarsWithSelection2",nm1s[i][0]+sel+'_nm1',"{0}[{1}]".format(nm1s[i][0],cutstr)) #print("define {}: {}".format(nm1s[i][0]+sel+'_nm1',"{0}[{1}]".format(nm1s[i][0],cutstr))) return d def FilterEvents(self,d): d_filter = d.Filter(self.presel) + print("*** Filter",self.presel) return d_filter def AddWeights(self,d,weight): d = d.Define("evt_weight","{0}".format(weight)) + print("*** Define AddWeights","evt_weight","{0}".format(weight)) return d def getRDF(self): @@ -99,6 +110,7 @@ def getRDF(self): d = self.AddVarsWithSelection(d) if self.cfg['presel'] is not None: d = d.Filter(self.cfg['presel']) + print("*** Filter getRDF",self.cfg['presel']) xsec_weights = 1 d = self.AddWeights(d,xsec_weights) return d,xsec_weights @@ -111,14 +123,18 @@ def getplots(self,d,weight,plots_1d,plots_2d,plots_nm1,varlabel): plots_2d = [] if plots_nm1 is None: plots_nm1 = [] - +#!# + self.isData = True +#!# for plt in plots_1d: if not plt in self.cfg['plot_setting']: print("{} not registered in plot setting!".format(plt)) if self.isData: h = d.Histo1D(tuple(self.cfg['plot_setting'][plt]),plt+varlabel) + print("*** Histo1D getplots",plt,self.isData,tuple(self.cfg['plot_setting'][plt]),plt+varlabel) else: h = d.Histo1D(tuple(self.cfg['plot_setting'][plt]),plt+varlabel,weight) + print("*** Histo1D getplots",plt,self.isData,tuple(self.cfg['plot_setting'][plt]),plt+varlabel,weight) hs.append(h) for plt in plots_nm1: @@ -179,6 +195,7 @@ def makeHistFiles(self): d_sr = d if self.cfg['regions'][sr] is not None: d_sr = d_sr.Filter(self.cfg['regions'][sr]) + print("*** Filter makeHistFiles1",sr,self.cfg['regions'][sr]) newd_evt = fout.mkdir("{}_evt".format(sr)) hs = self.getplots(d_sr,weight="evt_weight",plots_1d=self.cfg['event_variables'],plots_2d=self.cfg['event_2d_plots'],plots_nm1=self.cfg.get('event_nm1'),varlabel="")