From a61a72d7856829c89aac3d6c11f9e9acc1770ee1 Mon Sep 17 00:00:00 2001 From: Wolfgang Adam Date: Tue, 25 Nov 2025 14:46:44 +0100 Subject: [PATCH 01/16] remove files that are now in DrawHits --- HitAnalyzer/test/floatMod.C | 12 -- HitAnalyzer/test/histogramDefinition.py | 206 ------------------------ 2 files changed, 218 deletions(-) delete mode 100644 HitAnalyzer/test/floatMod.C delete mode 100644 HitAnalyzer/test/histogramDefinition.py 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/HitAnalyzer/test/histogramDefinition.py b/HitAnalyzer/test/histogramDefinition.py deleted file mode 100644 index 8bfc321..0000000 --- a/HitAnalyzer/test/histogramDefinition.py +++ /dev/null @@ -1,206 +0,0 @@ -# -# Definitions for histograms from a configuration file -# -from fnmatch import fnmatch - -class HistogramDefinition: - ''' A single histogram definition. - ''' - reqGenFields = [ 'canvasName', 'histogramName', 'histogramTitle', 'variable', 'baseCuts' ] - 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 - - def __init__(self,name,inputDict): - ''' Define histogram and drawing parameters from dictionary. - Arguments: - name ....... name of the histogram definition - inputDict .. dictionary with values - ''' - # - assert name.isalnum() - self.name = name - # - # read parameters - # - # default is None for all parameters - # - self.parameters = { x:None for x in HistogramDefinition.allFields } - # - # loop over main dictionary - # - for k,v in inputDict.items(): - # - # skip standard python entries - # - if k.startswith('__'): - continue - # - # standard entry in main dictionary - store value - # - if k in HistogramDefinition.allFields: - # - # general variable - # - self.parameters[k] = v - # - # nested dictionary with parameters for a specific module type - # - elif k.startswith("mType"): - # - # mType-specific histogram parameters - # - assert type(v)==dict and ( not k in self.parameters ) and len(k)>5 and k[5:].isdigit() - self.parameters[k] = { x:None for x in HistogramDefinition.allHistFields } - for kh,vh in v.items(): - if kh in HistogramDefinition.allHistFields: - self.parameters[k][kh] = vh - else: - print("Warning: key",kh, \ - "is not a standard histogram field name - ignoring the entry in", \ - self.name) - # - # unknown key - # - else: - print("Warning: key",k,"is not a standard field name - ignoring the entry in",self.name) - # - # set some parameters from other inputs if unspecified - # - if self.parameters['canvasName']==None: - self.parameters['canvasName'] = "c" + self.name[0].upper() + self.name[1:] - if self.parameters['histogramName']==None: - self.parameters['histogramName'] = "h" + self.name[0].upper() + self.name[1:] - # - # make sure all required fields are present - # - for f in HistogramDefinition.requiredFields: - assert ( f in self.parameters ) and self.parameters[f]!=None - - #def __getitem__(self,field): - # if field in self.parameters: - # return self.parameters[field] - # return None - - def getParameter(self,name,mType=None): - ''' Retrieve a single parameter (optionally a specific one for a given module type) - ''' - result = None - # - # give priority to parameter specific to a module type - # - mTName = "mType"+str(mType) if mType!=None else None - if ( mTName in self.parameters ) and ( name in self.parameters[mTName] ): - result = self.parameters[mTName][name] - if result!=None: - return self.parameters[mTName][name] - # - # not found: use general parameter - # - if name in self.parameters: - return self.parameters[name] - return None - - def vetoMType(self,mType): - ''' Check for an mType entry with display = False - ''' - # - # try to get 'display' parameter - # - name = 'display' - mTName = "mType"+str(mType) if mType!=None else None - if ( mTName in self.parameters ) and ( name in self.parameters[mTName] ): - if self.parameters[mTName][name]!=None: - return not self.parameters[mTName][name] - return False - -class HistogramDefinitions: - - def __init__(self): - self.allDefinitions = { } - self.allHistoNames = set() - self.allCanvases = set() - self.byCanvas = { } - - def add(self,hdef): - assert not hdef.name in self.allDefinitions - - hName = hdef.getParameter('histogramName') - assert not hName in self.allHistoNames - self.allHistoNames.add(hName) - - cName = hdef.getParameter('canvasName') - assert not cName in self.allCanvases - self.allDefinitions[hdef.name] = hdef - self.allCanvases.add(cName) - - if not cName in self.byCanvas: - self.byCanvas[cName] = { } - self.byCanvas[cName][hdef.name] = hdef - - def canvasNames(self): - return self.allCanvases - - def __getitem__(self,name): - if name in self.allDefinitions: - return self.allDefinitions[name] - return None - -def loadHistogramDefinitions(configName,selectedNames=[],vetoedNames=[]): - ''' Load histogram definitions from a configuration file. The configuration file - 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) - ''' - # load histogram definitions - # - result = HistogramDefinitions() - if configName==None: - return result - - 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() - assert type(hDict)==dict - # - # check if in list of histograms to be displayed - # - selFlg = False - for p in selectedNames: - if fnmatch(n,p): - selFlg = True - break - if not selFlg: - continue - # - # check if in list of histograms to be vetoed - # - selFlg = True - for p in vetoedNames: - if fnmatch(n,p): - selFlg = False - break - if not selFlg: - continue - # - # add histogram - # - hDef = HistogramDefinition(n,hDict) - result.add(hDef) - print("Added",hDef.getParameter('canvasName')) - - return result - From 56e97dd779f96d0e53f06707bf36af49fe81656e Mon Sep 17 00:00:00 2001 From: Wolfgang Adam Date: Mon, 1 Dec 2025 18:25:57 +0100 Subject: [PATCH 02/16] logZ option --- DrawHits/drawHits.py | 10 +++++++--- DrawHits/histogramDefinition.py | 2 +- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/DrawHits/drawHits.py b/DrawHits/drawHits.py index 3748b9f..e7ce1a8 100755 --- a/DrawHits/drawHits.py +++ b/DrawHits/drawHits.py @@ -316,7 +316,7 @@ def fillHistoByDef(tree,hDef,extraCuts): 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 @@ -405,6 +405,8 @@ def drawHistoByDef(histos,hDef,logY=False,same=False): histos[mType][0].Draw("ZCOL") if logY or hDef.getParameter('logY',mType): ROOT.gPad.SetLogy(1) + if not is1D and ( logZ or hDef.getParameter('logZ',mType) ): + ROOT.gPad.SetLogz(1) ROOT.gPad.Update() #cnv.cd() @@ -484,7 +486,8 @@ 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) @@ -591,7 +594,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'] diff --git a/DrawHits/histogramDefinition.py b/DrawHits/histogramDefinition.py index 580f3f2..0685f28 100644 --- a/DrawHits/histogramDefinition.py +++ b/DrawHits/histogramDefinition.py @@ -13,7 +13,7 @@ class HistogramDefinition: reqHistFields = [ ] requiredFields = list(set(reqGenFields + reqHistFields)) # optional fields in the general section - optGenFields = [ 'variable', 'baseCuts', 'effCuts', 'logY' ] + optGenFields = [ 'variable', 'baseCuts', 'effCuts', 'logY', 'logZ' ] # optional fields in the histogram section optHistFields = [ 'variable', 'histogramName', 'histogramTitle', \ 'xNbins', 'xMin', 'xMax', 'xTitle', 'yTitle', 'yNbins', 'yMin', 'yMax', \ From 7f0de31edaf074c93fcbd57c3427a400f25f1827 Mon Sep 17 00:00:00 2001 From: Wolfgang Adam Date: Tue, 2 Dec 2025 12:06:25 +0100 Subject: [PATCH 03/16] utility for histograms --- DrawHits/fitRes.py | 178 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 178 insertions(+) create mode 100644 DrawHits/fitRes.py diff --git a/DrawHits/fitRes.py b/DrawHits/fitRes.py new file mode 100644 index 0000000..1669f5e --- /dev/null +++ b/DrawHits/fitRes.py @@ -0,0 +1,178 @@ +import sys +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.intGraph = None + + def fwhm(self): + ''' Return lowest / highest x corresponding to ymax/2, and ymax/2 + ''' + # + # target y value (1/2 maximum) and bin width (assume constant bin width) + # + y = self.hist.GetBinContent(self.hist.GetMaximumBin())/2. + 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+dx,xhigh+dx,y) + + def quantile(self,q): + +class FitCanvas: + + def __init__(self,canvas,mtype): + self.canvas = canvas + self.mtype = mtype + self.pad = None + self.fhist = None + + padName = mainCnv.GetName()+str(mtype-22) + for o in mainCnv.GetListOfPrimitives(): + if o.InheritsFrom(ROOT.TPad.Class()) and o.GetName()==padName: + assert self.pad==None + self.pad = o + assert self.pad!=None + + histNameBase = "h"+mainCnv.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() + +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('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 + +canvases = [ ] +slices = [ ] + +mainCnv.Draw() +fitCanvas = FitCanvas(mainCnv,args.moduleType) +if 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 = ROOT.TArrow() + fwhmArrow.SetLineColor(2) + fwhmArrow.SetLineWidth(2) + fwhmArrow.DrawArrow(x1,y,x2,y,0.005,"<>") + fitCanvas.pad.Update() + +else: + hist2D = fitCanvas.histogram() + assert hist2D.GetDimension()==2 + yaxis = hist2D.GetYaxis() + nby = hist2D.GetNbinsY() + ncol = int(sqrt(nby)) + while ncol*ncol") + ROOT.gPad.Update() + #cnv.Update() + From 2b8128b41999e29259eb9086c3cbe46d3402a1d3 Mon Sep 17 00:00:00 2001 From: Wolfgang Adam Date: Tue, 2 Dec 2025 15:32:46 +0100 Subject: [PATCH 04/16] some more options for fits --- DrawHits/drawHits.py | 35 +++- DrawHits/drawHitsTmp.yaml | 30 ++++ DrawHits/fitRes.py | 281 +++++++++++++++++++++++++++----- DrawHits/histogramDefinition.py | 2 +- 4 files changed, 301 insertions(+), 47 deletions(-) diff --git a/DrawHits/drawHits.py b/DrawHits/drawHits.py index e7ce1a8..80820c2 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) @@ -290,7 +292,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,6 +312,20 @@ 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() @@ -344,6 +360,8 @@ def drawHistoByDef(histos,hDef,logY=False,logZ=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 +380,7 @@ def drawHistoByDef(histos,hDef,logY=False,logZ=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 @@ -387,7 +406,7 @@ def drawHistoByDef(histos,hDef,logY=False,logZ=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,9 +422,19 @@ def drawHistoByDef(histos,hDef,logY=False,logZ=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 not is1D and ( logZ or hDef.getParameter('logZ',mType) ): + if is2D and ( logZ or hDef.getParameter('logZ',mType) ): ROOT.gPad.SetLogz(1) ROOT.gPad.Update() diff --git a/DrawHits/drawHitsTmp.yaml b/DrawHits/drawHitsTmp.yaml index 67873c4..11e05f9 100644 --- a/DrawHits/drawHitsTmp.yaml +++ b/DrawHits/drawHitsTmp.yaml @@ -421,3 +421,33 @@ 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 diff --git a/DrawHits/fitRes.py b/DrawHits/fitRes.py index 1669f5e..554869e 100644 --- a/DrawHits/fitRes.py +++ b/DrawHits/fitRes.py @@ -13,56 +13,165 @@ def interpolate(y,y1,y2,x1=0.,x2=1.): def __init__(self,histogram): self.hist = histogram - self.graph = None - self.intGraph = None + self.graph_ = None + self.cumulativeGraph_ = None + self.cumNormGraph_ = None + + 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 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) and bin width (assume constant bin width) + # target y value (1/2 maximum) # y = self.hist.GetBinContent(self.hist.GetMaximumBin())/2. - dx = self.hist.GetBinWidth(1) # - # preset result (low / high x values) + # use first intersection with upward slope # - xlow = None - xhigh = None + xups = self.intersects(y,cumulative=False,direction=1) + print("xups",xups) + xlow = xups[0] if xups else None # - # loop over pairs of adjacent histogram bins + # use last intersection with downward slope # - 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 + 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+dx,xhigh+dx,y) + 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 - def quantile(self,q): + return result[0] if result else None class FitCanvas: @@ -88,7 +197,7 @@ def __init__(self,canvas,mtype): assert self.fhist!=None def histogram(self): - print("fhist",self.fhist,self.fhist.hist) + #print("fhist",self.fhist,self.fhist.hist) return self.fhist.hist def fwhm(self): @@ -98,6 +207,8 @@ def fwhm(self): 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('file', help='input file', type=str, nargs=1, default=None) args = parser.parse_args() @@ -116,7 +227,13 @@ def fwhm(self): mainCnv.Draw() fitCanvas = FitCanvas(mainCnv,args.moduleType) -if not args.slices: +fwhmArrow = ROOT.TArrow() +fwhmArrow.SetLineColor(2) +fwhmArrow.SetLineWidth(2) +quantLine = ROOT.TLine() +quantLine.SetLineColor(4) +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() @@ -124,13 +241,24 @@ def fwhm(self): 10000*(x1+x2)/2.,10000*(x2-x1),10000*(x2-x1)/2/sqrt(2*log(2)),x1,x2)) fitCanvas.pad.cd() - fwhmArrow = ROOT.TArrow() - fwhmArrow.SetLineColor(2) - fwhmArrow.SetLineWidth(2) fwhmArrow.DrawArrow(x1,y,x2,y,0.005,"<>") + + quantiles = [ ] + qlmax = slices[-1].hist.GetMaximum()/10. + for sigma in args.quantile: + qs = [ ] + for sgn in [-1,1]: + p = ROOT.TMath.Freq(sgn*sigma) + qs.append(slices[-1].quantile(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) + fitCanvas.pad.Update() -else: +elif fitCanvas.histogram().GetDimension()==2 and args.slices: hist2D = fitCanvas.histogram() assert hist2D.GetDimension()==2 yaxis = hist2D.GetYaxis() @@ -142,9 +270,6 @@ def fwhm(self): #cnv = ROOT.TCanvas("cslice","cslice",1200,1200) #cnv.Divide(ncol,ncol) #slices = [ ] - fwhmArrow = ROOT.TArrow() - fwhmArrow.SetLineColor(2) - fwhmArrow.SetLineWidth(2) hNameBase = hist2D.GetName() hTitleBase = hist2D.GetTitle() for iby in range(nby): @@ -173,6 +298,76 @@ def fwhm(self): print(fmt.format(yaxis.GetBinLowEdge(iby+1),yaxis.GetBinUpEdge(iby+1), \ 10000*(x1+x2)/2.,10000*(x2-x1),10000*(x2-x1)/2/sqrt(2*log(2)),x1,x2)) fwhmArrow.DrawArrow(x1,y,x2,y,0.005,"<>") + + quantiles = [ ] + qlmax = slices[-1].hist.GetMaximum()/10. + for sigma in args.quantile: + qs = [ ] + for sgn in [-1,1]: + p = ROOT.TMath.Freq(sgn*sigma) + qs.append(slices[-1].quantile(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) + + 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() + print(nby,ymin,ymax,nbz,zmin,zmax) + # + # define summary histograms + # + summaries = [ ] + for isigma,sigma in enumerate(args.quantile): + hcentre = ROOT.TH2F(h.GetName()+"Qc"+str(isigma), \ + h.GetTitle()+" (center {:3.1f}#sigma interval)".format(sigma), \ + nby,ymin,ymax,nbz,zmin,zmax) + hcentre.GetXaxis().SetTitle(h.GetYaxis().GetTitle()) + hcentre.GetYaxis().SetTitle(h.GetZaxis().GetTitle()) + hhalfwidth = ROOT.TH2F(h.GetName()+"Qhw"+str(isigma), \ + h.GetTitle()+" (half-width {:3.1f}#sigma interval)".format(sigma), \ + nby,ymin,ymax,nbz,zmin,zmax) + hhalfwidth.GetXaxis().SetTitle(h.GetYaxis().GetTitle()) + hhalfwidth.GetYaxis().SetTitle(h.GetZaxis().GetTitle()) + summaries.append((hcentre,hhalfwidth)) + + for iy in range(nby): + for iz in range(nbz): + htmp = h.ProjectionX(iymin=iy+1,iymax=iy+1,izmin=iz+1,izmax=iz+1) + if htmp.GetSumOfWeights()<100: + continue + fhtmp = FitHistogram(htmp) + for isigma,sigma in enumerate(args.quantile): + q1 = fhtmp.quantile(ROOT.TMath.Freq(-sigma)) + q2 = fhtmp.quantile(ROOT.TMath.Freq(sigma)) + if q1!=None and q2!=None: + hsum = summaries[isigma][0] + ibin = hsum.GetBin(iy+1,iz+1) + hsum.SetBinContent(ibin,(q1+q2)/2.) + hsum = summaries[isigma][1] + hsum.SetBinContent(ibin,(q2-q1)/2.) + + nsigma = len(args.quantile) + c = ROOT.TCanvas("cSum","cSum",500*nsigma,1000) + c.Divide(nsigma,2) + for i in range(nsigma): + c.cd(i+1) + summaries[i][0].Draw("ZCOL") + ROOT.gPad.Update() + c.cd(nsigma+i+1) + summaries[i][1].Draw("ZCOL") ROOT.gPad.Update() - #cnv.Update() + c.Update() diff --git a/DrawHits/histogramDefinition.py b/DrawHits/histogramDefinition.py index 0685f28..78904ac 100644 --- a/DrawHits/histogramDefinition.py +++ b/DrawHits/histogramDefinition.py @@ -17,7 +17,7 @@ class HistogramDefinition: # optional fields in the histogram section optHistFields = [ 'variable', 'histogramName', 'histogramTitle', \ 'xNbins', 'xMin', 'xMax', 'xTitle', 'yTitle', 'yNbins', 'yMin', 'yMax', \ - 'zMin', 'zMax', 'display', 'profile' ] + 'zNbins', 'zMin', 'zMax', 'display', 'profile' ] optionalFields = list(set(optGenFields + optHistFields)) allFields = list(set(requiredFields + optionalFields)) allHistFields = list(set(reqHistFields + optHistFields)) From 567be6225a7d9224ab68cded26e80645fcf485ad Mon Sep 17 00:00:00 2001 From: Wolfgang Adam Date: Tue, 2 Dec 2025 15:55:37 +0100 Subject: [PATCH 05/16] cosmetics for res mean/width in 2d --- DrawHits/drawHits.py | 18 ++++++++++++++++++ DrawHits/histogramDefinition.py | 3 ++- 2 files changed, 20 insertions(+), 1 deletion(-) diff --git a/DrawHits/drawHits.py b/DrawHits/drawHits.py index 80820c2..7b04884 100755 --- a/DrawHits/drawHits.py +++ b/DrawHits/drawHits.py @@ -520,6 +520,8 @@ def addHistogram(varString,cuts,effCuts=None,name='userHist'): 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 = [ ] @@ -530,6 +532,20 @@ 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? +# +match args.zone: + case "barrel": + zoneCuts = "detNormal.Rho()>0.99" + case "endcap": + zoneCuts = "detNormal.Rho()<0.01" + case "tilted": + zoneCuts = "detNormal.Rho()>0.05&&detNormal.Rho()<0.095" + case _: + zoneCuts = "" +args.cuts = cutString(args.cuts,zoneCuts) + # # load histogram definitions # @@ -632,6 +648,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/DrawHits/histogramDefinition.py b/DrawHits/histogramDefinition.py index 78904ac..2d585c0 100644 --- a/DrawHits/histogramDefinition.py +++ b/DrawHits/histogramDefinition.py @@ -16,7 +16,8 @@ class HistogramDefinition: optGenFields = [ 'variable', 'baseCuts', 'effCuts', 'logY', 'logZ' ] # optional fields in the histogram section optHistFields = [ 'variable', 'histogramName', 'histogramTitle', \ - 'xNbins', 'xMin', 'xMax', 'xTitle', 'yTitle', 'yNbins', 'yMin', 'yMax', \ + 'xNbins', 'xMin', 'xMax', 'xTitle', 'yTitle', 'zTitle', \ + 'yNbins', 'yMin', 'yMax', \ 'zNbins', 'zMin', 'zMax', 'display', 'profile' ] optionalFields = list(set(optGenFields + optHistFields)) allFields = list(set(requiredFields + optionalFields)) From 14b0a3444c89e295934f1254ca16328f800f5d00 Mon Sep 17 00:00:00 2001 From: Wolfgang Adam Date: Tue, 2 Dec 2025 16:19:31 +0100 Subject: [PATCH 06/16] change match to classical if elif statement to stay compatible with older py verions --- DrawHits/drawHits.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/DrawHits/drawHits.py b/DrawHits/drawHits.py index 7b04884..a64113a 100755 --- a/DrawHits/drawHits.py +++ b/DrawHits/drawHits.py @@ -535,15 +535,14 @@ def addHistogram(varString,cuts,effCuts=None,name='userHist'): # # add cut for zone definition? # -match args.zone: - case "barrel": - zoneCuts = "detNormal.Rho()>0.99" - case "endcap": - zoneCuts = "detNormal.Rho()<0.01" - case "tilted": - zoneCuts = "detNormal.Rho()>0.05&&detNormal.Rho()<0.095" - case _: - zoneCuts = "" +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) # From 5fc9c8828f454cb1776eaea63c242877d7c46269 Mon Sep 17 00:00:00 2001 From: Wolfgang Adam Date: Wed, 3 Dec 2025 12:11:49 +0100 Subject: [PATCH 07/16] cosmetics --- DrawHits/drawHitsTmp.yaml | 30 +++++++++++ DrawHits/fitRes.py | 107 +++++++++++++++++++++++++++++++++----- 2 files changed, 123 insertions(+), 14 deletions(-) diff --git a/DrawHits/drawHitsTmp.yaml b/DrawHits/drawHitsTmp.yaml index 11e05f9..77125ef 100644 --- a/DrawHits/drawHitsTmp.yaml +++ b/DrawHits/drawHitsTmp.yaml @@ -451,3 +451,33 @@ res3DXVsDxDzModX: 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 diff --git a/DrawHits/fitRes.py b/DrawHits/fitRes.py index 554869e..6d3af34 100644 --- a/DrawHits/fitRes.py +++ b/DrawHits/fitRes.py @@ -209,6 +209,8 @@ def fwhm(self): 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('file', help='input file', type=str, nargs=1, default=None) args = parser.parse_args() @@ -222,6 +224,10 @@ def fwhm(self): mainCnv = k.ReadObj() assert mainCnv!=None +dbgQuantPoints = [ ] +for qp in args.dbgQuantiles: + dbgQuantPoints.append( tuple([float(x.strip()) for x in qp.strip().split(",")]) ) + canvases = [ ] slices = [ ] @@ -331,18 +337,50 @@ def fwhm(self): # summaries = [ ] for isigma,sigma in enumerate(args.quantile): - hcentre = ROOT.TH2F(h.GetName()+"Qc"+str(isigma), \ - h.GetTitle()+" (center {:3.1f}#sigma interval)".format(sigma), \ - nby,ymin,ymax,nbz,zmin,zmax) - hcentre.GetXaxis().SetTitle(h.GetYaxis().GetTitle()) - hcentre.GetYaxis().SetTitle(h.GetZaxis().GetTitle()) - hhalfwidth = ROOT.TH2F(h.GetName()+"Qhw"+str(isigma), \ - h.GetTitle()+" (half-width {:3.1f}#sigma interval)".format(sigma), \ - nby,ymin,ymax,nbz,zmin,zmax) - hhalfwidth.GetXaxis().SetTitle(h.GetYaxis().GetTitle()) - hhalfwidth.GetYaxis().SetTitle(h.GetZaxis().GetTitle()) + y1 = ROOT.TMath.Freq(-sigma) + y2 = ROOT.TMath.Freq(sigma) + if nby==1: + hcentre = ROOT.TH1F(h.GetName()+"Qc"+str(isigma), \ + h.GetTitle()+" (mean from quantiles)".format(y1,y2), \ + nbz,zmin,zmax) + hcentre.GetXaxis().SetTitle(h.GetZaxis().GetTitle()) + hcentre.GetYaxis().SetTitle("mean of {:4.1%}/{:4.1%} quantiles".format(y1,y2)) + 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".format(y1,y2)) + elif nbz==1: + hcentre = ROOT.TH1F(h.GetName()+"Qc"+str(isigma), \ + h.GetTitle()+" (mean from quantiles)".format(y1,y2), \ + nby,ymin,ymax) + hcentre.GetXaxis().SetTitle(h.GetYaxis().GetTitle()) + hcentre.GetYaxis().SetTitle("mean of {:4.1%}/{:4.1%} quantiles".format(y1,y2)) + 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".format(y1,y2)) + else: + hcentre = ROOT.TH2F(h.GetName()+"Qc"+str(isigma), \ + h.GetTitle()+" (mean from quantiles)".format(y1,y2), \ + nby,ymin,ymax,nbz,zmin,zmax) + hcentre.GetXaxis().SetTitle(h.GetYaxis().GetTitle()) + hcentre.GetYaxis().SetTitle(h.GetZaxis().GetTitle()) + hcentre.GetZaxis().SetTitle("mean of {:4.1%}/{:4.1%} quantiles".format(y1,y2)) + 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".format(y1,y2)) summaries.append((hcentre,hhalfwidth)) + dbgBins = [ ] + for y,z in dbgQuantPoints: + dbgBins.append( (h.GetYaxis().FindBin(y), h.GetZaxis().FindBin(z)) ) + + allDbgObjects = { x:[ ] for x in range(len(args.quantile)) } for iy in range(nby): for iz in range(nbz): htmp = h.ProjectionX(iymin=iy+1,iymax=iy+1,izmin=iz+1,izmax=iz+1) @@ -354,20 +392,61 @@ def fwhm(self): q2 = fhtmp.quantile(ROOT.TMath.Freq(sigma)) if q1!=None and q2!=None: hsum = summaries[isigma][0] - ibin = hsum.GetBin(iy+1,iz+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,(q1+q2)/2.) hsum = summaries[isigma][1] - hsum.SetBinContent(ibin,(q2-q1)/2.) + hsum.SetBinContent(ibin,(q2-q1)/2./sigma) + + #allDbgObjects[isigma] = [ ] + if (iy+1,iz+1) in dbgBins: + dbgObjects = [ ] + hResDbg =htmp.Clone("hResDbg_"+str(isigma)+"_"+str(len(allDbgObjects[isigma]))) + hResDbg.SetTitle("hResDbg "+str(isigma)+" y"+str(iy+1)+" z"+str(iz+1)) + hResDbg.GetXaxis().SetTitle(h.GetXaxis().GetTitle()) + savePad = ROOT.gPad + dbgObjects.append(ROOT.TCanvas("c"+hResDbg.GetName()[1:],"c"+hResDbg.GetName()[1:],500,500)) + dbgObjects[0].SetGridx(1) + dbgObjects[0].SetGridy(1) + dbgObjects.append(hResDbg) + hResDbg.Scale(1./hResDbg.GetMaximum()) + hResDbg.Draw("hist") + if q1!=None and q2!=None: + quantLine.DrawLine(q1,0.,q1,hResDbg.GetMaximum()/1.) + quantLine.DrawLine(q2,0.,q2,hResDbg.GetMaximum()/1.) + g = ROOT.TGraph() + g.SetLineColor(ROOT.kMagenta) + g.SetLineStyle(2) + #g.SetMarkerStyle(20) + for x,y in fhtmp.cumulativeNormGraph(): + g.SetPoint(g.GetN(),x,y) + g.Draw("L") + dbgObjects.append(g) + + dbgObjects[0].Update() + allDbgObjects[isigma].append(dbgObjects) + savePad.cd() nsigma = len(args.quantile) c = ROOT.TCanvas("cSum","cSum",500*nsigma,1000) c.Divide(nsigma,2) for i in range(nsigma): c.cd(i+1) - summaries[i][0].Draw("ZCOL") + if summaries[i][0].GetDimension()==2: + ROOT.gPad.SetRightMargin(0.15) + summaries[i][0].Draw("ZCOL") + else: + summaries[i][0].Draw("HIST") ROOT.gPad.Update() c.cd(nsigma+i+1) - summaries[i][1].Draw("ZCOL") + if summaries[i][1].GetDimension()==2: + summaries[i][1].Draw("ZCOL") + else: + summaries[i][1].Draw("HIST") ROOT.gPad.Update() c.Update() From 2fbcf3e1b42c8764c8cf8906ebad2bb1e886da07 Mon Sep 17 00:00:00 2001 From: Wolfgang Adam Date: Wed, 3 Dec 2025 15:35:54 +0100 Subject: [PATCH 08/16] add histograms residuals (vs. pos on strip) vs. angle / cluster size --- DrawHits/drawHitsTmp.yaml | 61 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) diff --git a/DrawHits/drawHitsTmp.yaml b/DrawHits/drawHitsTmp.yaml index 77125ef..5bd60c5 100644 --- a/DrawHits/drawHitsTmp.yaml +++ b/DrawHits/drawHitsTmp.yaml @@ -481,3 +481,64 @@ res3DXVsDxDz: 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 + From 331c60aff5834f78bf580a663bad1fae7293e3cb Mon Sep 17 00:00:00 2001 From: Wolfgang Adam Date: Wed, 3 Dec 2025 15:36:32 +0100 Subject: [PATCH 09/16] define quantiles from cumulative histo of residuals via spline --- DrawHits/fitRes.py | 83 ++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 77 insertions(+), 6 deletions(-) diff --git a/DrawHits/fitRes.py b/DrawHits/fitRes.py index 6d3af34..0333051 100644 --- a/DrawHits/fitRes.py +++ b/DrawHits/fitRes.py @@ -16,6 +16,8 @@ def __init__(self,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. @@ -54,7 +56,20 @@ def cumulativeNormGraph(self): 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 @@ -172,7 +187,58 @@ def quantile(self,prob): 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*=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) Date: Thu, 4 Dec 2025 12:49:18 +0100 Subject: [PATCH 10/16] mainly cosmetics --- DrawHits/fitRes.py | 278 ++++++++++++++++++++++++++++++++------------- 1 file changed, 198 insertions(+), 80 deletions(-) diff --git a/DrawHits/fitRes.py b/DrawHits/fitRes.py index 0333051..2835cdc 100644 --- a/DrawHits/fitRes.py +++ b/DrawHits/fitRes.py @@ -1,4 +1,4 @@ -import sys +import sys,os from math import sqrt,log import ROOT import argparse @@ -247,14 +247,14 @@ def __init__(self,canvas,mtype): self.pad = None self.fhist = None - padName = mainCnv.GetName()+str(mtype-22) - for o in mainCnv.GetListOfPrimitives(): + 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"+mainCnv.GetName()[1:]+str(mtype) + 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): @@ -268,7 +268,53 @@ def histogram(self): 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) + 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) @@ -277,6 +323,11 @@ def fwhm(self): 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() @@ -294,6 +345,10 @@ def fwhm(self): 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 = [ ] @@ -319,7 +374,7 @@ def fwhm(self): qlmax = slices[-1].hist.GetMaximum()/10. for sigma in args.quantile: qs = [ ] - for sgn in [-1,1]: + for sgn in [-1,0,1]: p = ROOT.TMath.Freq(sgn*sigma) #qs.append(slices[-1].quantile(p)) qs.append(slices[-1].findRootSpline(p)) @@ -328,6 +383,7 @@ def fwhm(self): 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) fitCanvas.pad.Update() @@ -372,19 +428,21 @@ def fwhm(self): 10000*(x1+x2)/2.,10000*(x2-x1),10000*(x2-x1)/2/sqrt(2*log(2)),x1,x2)) fwhmArrow.DrawArrow(x1,y,x2,y,0.005,"<>") - quantiles = [ ] - qlmax = slices[-1].hist.GetMaximum()/10. - for sigma in args.quantile: - qs = [ ] - for sgn in [-1,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) + 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() @@ -399,7 +457,6 @@ def fwhm(self): nbz = h.GetNbinsZ() zmin = h.GetZaxis().GetXmin() zmax = h.GetZaxis().GetXmax() - print(nby,ymin,ymax,nbz,zmin,zmax) # # define summary histograms # @@ -409,108 +466,167 @@ def fwhm(self): y2 = ROOT.TMath.Freq(sigma) if nby==1: hcentre = ROOT.TH1F(h.GetName()+"Qc"+str(isigma), \ - h.GetTitle()+" (mean from quantiles)".format(y1,y2), \ + h.GetTitle()+" (median)", \ nbz,zmin,zmax) hcentre.GetXaxis().SetTitle(h.GetZaxis().GetTitle()) - hcentre.GetYaxis().SetTitle("mean of {:4.1%}/{:4.1%} quantiles".format(y1,y2)) + 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".format(y1,y2)) + 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()+" (mean from quantiles)".format(y1,y2), \ + h.GetTitle()+" (median)", \ nby,ymin,ymax) hcentre.GetXaxis().SetTitle(h.GetYaxis().GetTitle()) - hcentre.GetYaxis().SetTitle("mean of {:4.1%}/{:4.1%} quantiles".format(y1,y2)) + 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".format(y1,y2)) + 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()+" (mean from quantiles)".format(y1,y2), \ + h.GetTitle()+" (median)", \ nby,ymin,ymax,nbz,zmin,zmax) hcentre.GetXaxis().SetTitle(h.GetYaxis().GetTitle()) hcentre.GetYaxis().SetTitle(h.GetZaxis().GetTitle()) - hcentre.GetZaxis().SetTitle("mean of {:4.1%}/{:4.1%} quantiles".format(y1,y2)) + 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".format(y1,y2)) + hhalfwidth.GetZaxis().SetTitle("#sigma from {:4.1%}/{:4.1%} quantiles [#mum]".format(y1,y2)) + hhalfwidth.SetMinimum(0.) summaries.append((hcentre,hhalfwidth)) - dbgBins = [ ] - for y,z in dbgQuantPoints: - dbgBins.append( (h.GetYaxis().FindBin(y), h.GetZaxis().FindBin(z)) ) + 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)) } + nDbg = len(dbgBins) + nrow = ncol = int(sqrt(nDbg)) + while nrow*ncol100: + #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 + 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) + 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} Date: Thu, 4 Dec 2025 14:40:35 +0100 Subject: [PATCH 11/16] style --- DrawHits/fitRes.py | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/DrawHits/fitRes.py b/DrawHits/fitRes.py index 2835cdc..6132fcb 100644 --- a/DrawHits/fitRes.py +++ b/DrawHits/fitRes.py @@ -357,9 +357,14 @@ def setHistogramMinMax(histo,limits,margins=0.05): fwhmArrow = ROOT.TArrow() fwhmArrow.SetLineColor(2) fwhmArrow.SetLineWidth(2) -quantLine = ROOT.TLine() -quantLine.SetLineColor(4) -quantLine.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()) @@ -372,7 +377,7 @@ def setHistogramMinMax(histo,limits,margins=0.05): quantiles = [ ] qlmax = slices[-1].hist.GetMaximum()/10. - for sigma in args.quantile: + for isigma,sigma in enumerate(args.quantile): qs = [ ] for sgn in [-1,0,1]: p = ROOT.TMath.Freq(sgn*sigma) @@ -381,9 +386,9 @@ def setHistogramMinMax(histo,limits,margins=0.05): #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) + 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() @@ -575,6 +580,8 @@ def setHistogramMinMax(histo,limits,margins=0.05): hResDbg.SetTitle(hResTitle) hResDbg.GetXaxis().SetTitle(h.GetXaxis().GetTitle()) hResDbg.GetYaxis().SetTitle("fraction per bin / cumulatitive fraction") + hResDbg.SetFillStyle(1001) + hResDbg.SetFillColor(ROOT.kYellow) #dbgObjects.append(ROOT.TCanvas("c"+hResDbg.GetName()[1:],"c"+hResDbg.GetName()[1:],500,500)) iDbgPad += 1 dbgObjects.append(allDbgObjects['canvas'].cd(iDbgPad)) @@ -588,11 +595,11 @@ def setHistogramMinMax(histo,limits,margins=0.05): q0 = fhtmp.findRootSpline(ROOT.TMath.Freq(0.)) qp1 = fhtmp.findRootSpline(ROOT.TMath.Freq(sigma)) if qm1!=None: - quantLine.DrawLine(qm1,0.,qm1,hResDbg.GetMaximum()/1.) + quantLines[isigma].DrawLine(qm1,0.,qm1,hResDbg.GetMaximum()/1.) if q0!=None: - quantLine.DrawLine(q0,0.,q0,hResDbg.GetMaximum()/1.) + quantLines[0].DrawLine(q0,0.,q0,hResDbg.GetMaximum()/1.) if qp1!=None: - quantLine.DrawLine(qp1,0.,qp1,hResDbg.GetMaximum()/1.) + quantLines[isigma].DrawLine(qp1,0.,qp1,hResDbg.GetMaximum()/1.) g = ROOT.TGraph() g.SetLineColor(ROOT.kMagenta) g.SetLineStyle(2) From 2e2288ea98f41e21a1d9114c7da52b39c7f64674 Mon Sep 17 00:00:00 2001 From: Wolfgang Adam Date: Fri, 5 Dec 2025 13:48:27 +0100 Subject: [PATCH 12/16] tuning of quantile / width calculation --- DrawHits/drawHits.py | 3 + DrawHits/drawHitsTmp.yaml | 120 ++++++++++++++++++++++++++++++++ DrawHits/fitRes.py | 24 +++++-- DrawHits/histogramDefinition.py | 2 +- 4 files changed, 143 insertions(+), 6 deletions(-) diff --git a/DrawHits/drawHits.py b/DrawHits/drawHits.py index a64113a..bac27e7 100755 --- a/DrawHits/drawHits.py +++ b/DrawHits/drawHits.py @@ -398,6 +398,9 @@ def drawHistoByDef(histos,hDef,logY=False,logZ=False,same=False): 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) diff --git a/DrawHits/drawHitsTmp.yaml b/DrawHits/drawHitsTmp.yaml index 5bd60c5..2ba0371 100644 --- a/DrawHits/drawHitsTmp.yaml +++ b/DrawHits/drawHitsTmp.yaml @@ -542,3 +542,123 @@ res3DXVsDxDzW2: 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==2' + 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/fitRes.py b/DrawHits/fitRes.py index 6132fcb..cf205d2 100644 --- a/DrawHits/fitRes.py +++ b/DrawHits/fitRes.py @@ -197,6 +197,12 @@ def findRootSpline(self,value,eps=0.001): # spline = self.cumulativeNormSpline() # + # No result if first / last point is above / below required value + # + cGraph = self.cumulativeGraph() + if len(cGraph)==0 or cGraph[0][1]>value or cGraph[-1][1]100: + #if htmp.GetSumOfWeights()>100: + # 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)) @@ -543,6 +553,7 @@ def setHistogramMinMax(histo,limits,margins=0.05): qm1 = None q0 = None qp1 = None + quantiles.append((qm1,q0,qp1)) hsum = summaries[isigma][0] if nby==1: ibin = hsum.GetBin(iz+1) @@ -555,6 +566,8 @@ def setHistogramMinMax(histo,limits,margins=0.05): 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) @@ -581,7 +594,7 @@ def setHistogramMinMax(histo,limits,margins=0.05): hResDbg.GetXaxis().SetTitle(h.GetXaxis().GetTitle()) hResDbg.GetYaxis().SetTitle("fraction per bin / cumulatitive fraction") hResDbg.SetFillStyle(1001) - hResDbg.SetFillColor(ROOT.kYellow) + hResDbg.SetFillColor(ROOT.kOrange) #dbgObjects.append(ROOT.TCanvas("c"+hResDbg.GetName()[1:],"c"+hResDbg.GetName()[1:],500,500)) iDbgPad += 1 dbgObjects.append(allDbgObjects['canvas'].cd(iDbgPad)) @@ -591,9 +604,10 @@ def setHistogramMinMax(histo,limits,margins=0.05): hResDbg.Scale(1./hResDbg.GetMaximum()) hResDbg.Draw("hist") for isigma,sigma in enumerate(args.quantile): - qm1 = fhtmp.findRootSpline(ROOT.TMath.Freq(-sigma)) - q0 = fhtmp.findRootSpline(ROOT.TMath.Freq(0.)) - qp1 = fhtmp.findRootSpline(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)) + qm1,q0,qp1 = quantiles[isigma] if qm1!=None: quantLines[isigma].DrawLine(qm1,0.,qm1,hResDbg.GetMaximum()/1.) if q0!=None: diff --git a/DrawHits/histogramDefinition.py b/DrawHits/histogramDefinition.py index 2d585c0..0e8c451 100644 --- a/DrawHits/histogramDefinition.py +++ b/DrawHits/histogramDefinition.py @@ -15,7 +15,7 @@ class HistogramDefinition: # optional fields in the general section optGenFields = [ 'variable', 'baseCuts', 'effCuts', 'logY', 'logZ' ] # optional fields in the histogram section - optHistFields = [ 'variable', 'histogramName', 'histogramTitle', \ + optHistFields = [ 'variable', 'histogramName', 'histogramTitle', 'fit', \ 'xNbins', 'xMin', 'xMax', 'xTitle', 'yTitle', 'zTitle', \ 'yNbins', 'yMin', 'yMax', \ 'zNbins', 'zMin', 'zMax', 'display', 'profile' ] From 446a242522ed009e23928c600827fb057705c403 Mon Sep 17 00:00:00 2001 From: Wolfgang Adam Date: Fri, 5 Dec 2025 16:51:50 +0100 Subject: [PATCH 13/16] add efficiency as f(dx/dz) --- DrawHits/drawHits.py | 3 ++- DrawHits/drawHitsTmp.yaml | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/DrawHits/drawHits.py b/DrawHits/drawHits.py index bac27e7..a140072 100755 --- a/DrawHits/drawHits.py +++ b/DrawHits/drawHits.py @@ -280,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] @@ -391,7 +392,7 @@ def drawHistoByDef(histos,hDef,logY=False,logZ=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) diff --git a/DrawHits/drawHitsTmp.yaml b/DrawHits/drawHitsTmp.yaml index 2ba0371..7957666 100644 --- a/DrawHits/drawHitsTmp.yaml +++ b/DrawHits/drawHitsTmp.yaml @@ -646,7 +646,7 @@ effDxDzW1: variable: 'path.x()/path.z()' histogramTitle: 'efficiency, 1-strip clusters (dx/dz)' baseCuts: 'tof<12.5' - effCuts: 'hasRecHit>0&&clusterSize==2' + effCuts: 'hasRecHit>0&&clusterSize==1' xTitle: 'local dx/dz' xNbins: 21 xMin: -1.05 From e6572f96a7a57c44a036837263268db39e9dee9a Mon Sep 17 00:00:00 2001 From: Wolfgang Adam Date: Thu, 11 Dec 2025 17:11:13 +0100 Subject: [PATCH 14/16] adding result of fit --- DrawHits/fitRes.py | 168 +++++++++++++++++++++++++++++++++------------ 1 file changed, 126 insertions(+), 42 deletions(-) diff --git a/DrawHits/fitRes.py b/DrawHits/fitRes.py index cf205d2..5bdbdd6 100644 --- a/DrawHits/fitRes.py +++ b/DrawHits/fitRes.py @@ -244,7 +244,27 @@ def findRootSpline(self,value,eps=0.001): return (xl+xh)/2. + def fitGaus(self,prob1=0.,prob2=1.): + assert prob2>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 + + fitPtr = self.hist.Fit("gaus","S0","",xmin,xmax) + if ( not fitPtr.IsValid() ) or fitPtr.IsEmpty(): + return result + return fitPtr,self.hist.GetFunction("gaus") + + class FitCanvas: def __init__(self,canvas,mtype): @@ -321,6 +341,53 @@ def setHistogramMinMax(histo,limits,margins=0.05): 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) @@ -473,45 +540,48 @@ def setHistogramMinMax(histo,limits,margins=0.05): # 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.) +#!# 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: @@ -525,6 +595,7 @@ def setHistogramMinMax(histo,limits,margins=0.05): dbgBins = sorted(dbgBins) allDbgObjects = { x:[ ] for x in range(len(args.quantile)) } + allFitObjects = [ ] nDbg = len(dbgBins) nrow = ncol = int(sqrt(nDbg)) while nrow*ncol100: + 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() allDbgObjects['canvas'].Update() if args.output: @@ -639,8 +722,8 @@ def setHistogramMinMax(histo,limits,margins=0.05): cName = "cResMeanWidth_mT"+str(args.moduleType) c = ROOT.TCanvas(cName,h.GetTitle()+" (mean/width summary)",500*nsigma,1000) - c.Divide(nsigma,2) - for i in range(nsigma): + c.Divide(nsigma+1,2) + for i in range(nsigma+1): c.cd(i+1) if summaries[i][0].GetDimension()==2: ROOT.gPad.SetRightMargin(0.15) @@ -650,12 +733,13 @@ def setHistogramMinMax(histo,limits,margins=0.05): setHistogramMinMax(summaries[i][0],(-250,250),margins=0.05) summaries[i][0].Draw("HIST") ROOT.gPad.Update() - c.cd(nsigma+i+1) + c.cd(nsigma+1+i+1) if summaries[i][1].GetDimension()==2: summaries[i][1].Draw("ZCOL") else: summaries[i][1].Draw("HIST") ROOT.gPad.Update() + c.Update() if args.output!=None: saveCanvas(c,args.file[0],args.output,outputFormats) From 53950fae086ec048ae3c14d8efb2bd36d08b8b1d Mon Sep 17 00:00:00 2001 From: Wolfgang Adam Date: Thu, 11 Dec 2025 17:17:15 +0100 Subject: [PATCH 15/16] adding result of fit --- DrawHits/fitRes.py | 28 ++++++++++++++++++++++++++-- 1 file changed, 26 insertions(+), 2 deletions(-) diff --git a/DrawHits/fitRes.py b/DrawHits/fitRes.py index 5bdbdd6..9578c3a 100644 --- a/DrawHits/fitRes.py +++ b/DrawHits/fitRes.py @@ -709,6 +709,29 @@ def defineMeanWidthHistograms(h,nby,ymin,ymax,nbz,zmin,zmax,isigma=None,sigma=No 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: @@ -724,6 +747,7 @@ def defineMeanWidthHistograms(h,nby,ymin,ymax,nbz,zmin,zmax,isigma=None,sigma=No 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 Date: Fri, 9 Jan 2026 12:59:48 +0100 Subject: [PATCH 16/16] update to fitRes --- DrawHits/fitRes.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/DrawHits/fitRes.py b/DrawHits/fitRes.py index 9578c3a..3f43967 100644 --- a/DrawHits/fitRes.py +++ b/DrawHits/fitRes.py @@ -258,11 +258,19 @@ def fitGaus(self,prob1=0.,prob2=1.): if xmin==None: return result - fitPtr = self.hist.Fit("gaus","S0","",xmin,xmax) + 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,self.hist.GetFunction("gaus") + return fitPtr,fitFunc class FitCanvas: @@ -700,6 +708,7 @@ def defineMeanWidthHistograms(h,nby,ymin,ymax,nbz,zmin,zmax,isigma=None,sigma=No allDbgObjects[isigma].append(dbgObjects) fitPtr = None fitFunc = None + print(htmp.GetName(),"Sum of weights",htmp.GetSumOfWeights()) if htmp.GetSumOfWeights()>100: fitPtr,fitFunc = fhtmp.fitGaus(ROOT.TMath.Freq(-2.),ROOT.TMath.Freq(2.)) print(fitPtr,fitFunc)