python/PlottingTools/plotResidualsPerLayer.py

0001 import ROOT
0002 from .fitResidual import fitResidual
0003 from .drawHistoAllChambers import drawHisto
0004
0005 layerCorrectionFactors = {'SL1':(1.17,1.16,1.15,1.14),
0006                           'SL2':(1.83,1.20,1.20,1.83),
0007                           'SL3':(1.14,1.15,1.16,1.17)}
0008
0009 def plotResLayer(fileName,sl,layer,
0010                           dir='DQMData/Run 1/DT/Run summary/DTCalibValidation',
0011                           option="HISTOPE1",draw=True):
0012
0013     mean_ymin = -0.02
0014     mean_ymax =  0.02
0015     sig_ymin = 0.
0016     sig_ymax = 0.1
0017
0018     slType = sl
0019     slStr = "SL%d" % slType
0020     layerType = layer
0021     layerStr = "Layer%d" % layerType
0022     verbose = False
0023     nSigmas = 2
0024
0025     ROOT.TH1.AddDirectory(False)
0026
0027     file = ROOT.TFile(fileName,'read')
0028
0029     wheels = (-2,-1,0,1,2)
0030     stations = (1,2,3,4)
0031
0032     # (Wh-2 MB1 Sec1 ... Wh-2 MB1 Sec12 ... Wh-1 MB1 Sec1 ... Wh-1 MB1 Sec12 ...)
0033     # (Wh-2 MB2 Sec1 ... Wh-2 MB2 Sec12 ... Wh-1 MB2 Sec1 ... Wh-1 MB1 Sec12 ...) ...
0034     nBins = 250
0035     if slType == 2: nBins = 180
0036     histoMean = ROOT.TH1F("h_ResMeanAll_%s_%s" % (slStr,layerStr),"Mean of residuals",nBins,0,nBins)
0037     histoSigma = ROOT.TH1F("h_ResSigmaAll_%s_%s" % (slStr,layerStr),"Sigma of residuals",nBins,0,nBins)
0038     for st in stations:
0039         nSectors = 12
0040         if st == 4: nSectors = 14
0041         if st == 4 and slType == 2: continue
0042         if verbose: print("Station",st)
0043         for wh in wheels:
0044             if verbose: print("Wheel",wh)
0045             for sec in range(1,nSectors+1):
0046                 if verbose: print("Sector",sec)
0047                 # Get histogram
0048                 histoName = "%s/Wheel%d/Station%d/Sector%d/%s/hResDist_STEP3_W%d_St%d_Sec%d_%s_%s" % (dir,wh,st,sec,slStr,wh,st,sec,slStr,layerStr)
0049                 print("Accessing",histoName)
0050                 histo = file.Get(histoName)
0051                 (histo,fitFunc) = fitResidual(histo,nSigmas,verbose)
0052                 fitMean = fitFunc.GetParameter(1)
0053                 fitMeanErr = fitFunc.GetParError(1)
0054                 fitSigma = fitFunc.GetParameter(2)
0055                 fitSigmaErr = fitFunc.GetParError(2)
0056
0057                 layerIdx = (layer - 1)
0058                 corrFactor = layerCorrectionFactors[slStr][layerIdx]
0059
0060                 binHistoNew = (st - 1)*60 + (wh + 2)*nSectors + sec
0061                 if verbose: print("Bin in summary histo",binHistoNew)
0062                 histoMean.SetBinContent(binHistoNew,fitMean)
0063                 histoMean.SetBinError(binHistoNew,fitMeanErr)
0064                 histoSigma.SetBinContent(binHistoNew,fitSigma*corrFactor)
0065                 histoSigma.SetBinError(binHistoNew,fitSigmaErr*corrFactor)
0066
0067                 if sec == 1:
0068                     label = "Wheel %d" % wh
0069                     if wh == -2: label += " MB%d" % st
0070                     histoMean.GetXaxis().SetBinLabel(binHistoNew,label)
0071                     histoSigma.GetXaxis().SetBinLabel(binHistoNew,label)
0072
0073     objectsMean = drawHisto(histoMean,title="Mean of residuals (cm)",
0074                                       ymin=mean_ymin,ymax=mean_ymax,option=option,draw=draw)
0075     objectsSigma = drawHisto(histoSigma,title="Resolution (cm)",
0076                                         ymin=sig_ymin,ymax=sig_ymax,option=option,draw=draw)
0077
0078     return (objectsMean,objectsSigma)
0079
0080 def plot(fileName,sl,
0081                   dir='DQMData/Run 1/DT/Run summary/DTCalibValidation',type='mean',option='HISTOPE1'):
0082     colors = (2,4,12,44,55,38,27,46)
0083     markers = (24,25,26,27,28,30,32,5)
0084     labels=['Layer 1','Layer 2','Layer 3','Layer 4']
0085     idx_type = None
0086     if type == 'mean': idx_type = 0
0087     elif type == 'sigma': idx_type = 1
0088     else: raise RuntimeError("Wrong option: %s" % type)
0089
0090     idx = 0
0091     canvas = None
0092     objects = None
0093     histos = []
0094     for layer in range(1,5):
0095         draw = False
0096         if not idx: draw = True
0097
0098         objs = plotResLayer(fileName,sl,layer,dir,option,draw)
0099         histos.append(objs[idx_type][1])
0100         histos[-1].SetName( "%s_%d" % (histos[-1].GetName(),idx) )
0101         if not idx:
0102             canvas = objs[idx_type][0]
0103             objects = objs[idx_type][2]
0104
0105         canvas.cd()
0106         if idx:
0107             histos[-1].SetLineColor(colors[ (idx - 1) % len(colors) ])
0108             histos[-1].SetMarkerColor(colors[ (idx - 1) % len(colors) ])
0109             histos[-1].SetMarkerStyle(markers[ (idx - 1) % len(markers) ])
0110
0111             histos[-1].Draw(option + "SAME")
0112
0113         idx += 1
0114
0115     legend = ROOT.TLegend(0.4,0.7,0.95,0.8)
0116     for idx in range( len(histos) ):
0117         histo = histos[idx]
0118         label = histo.GetName()
0119         if len(labels): label = labels[idx]
0120         legend.AddEntry(histo,label,"LP")
0121
0122         idx += 1
0123
0124     canvas.cd()
0125     legend.SetFillColor( canvas.GetFillColor() )
0126     legend.Draw("SAME")
0127     objects.append(legend)
0128
0129     # Compute averages
0130     # (Wh-2 MB1 Sec1 ... Wh-2 MB1 Sec12 ... Wh-1 MB1 Sec1 ... Wh-1 MB1 Sec12 ...)
0131     # (Wh-2 MB2 Sec1 ... Wh-2 MB2 Sec12 ... Wh-1 MB2 Sec1 ... Wh-1 MB1 Sec12 ...) ...
0132     import math
0133     wheels = (-2,-1,0,1,2)
0134     stations = (1,2,3,4)
0135     slType = sl
0136     slStr = "SL%d" % slType
0137
0138     nBinsAve = len(stations)*len(wheels)
0139     histoAverage = ROOT.TH1F("h_AverageAll_" + slStr,"",nBinsAve,0,nBinsAve)
0140     averages = {}
0141     averagesErr = {}
0142     averagesSumw = {}
0143     print("Averages:")
0144     for st in stations:
0145         nSectors = 12
0146         if st == 4: nSectors = 14
0147         if st == 4 and slType == 2: continue
0148         for wh in wheels:
0149             binHistoAve = (st - 1)*5 + (wh + 2) + 1
0150             label = "Wheel %d" % wh
0151             if wh == -2: label += " MB%d" % st
0152             histoAverage.GetXaxis().SetBinLabel(binHistoAve,label)
0153
0154             averages[(st,wh)] = 0.
0155             averagesSumw[(st,wh)] = 0.
0156             for sec in range(1,nSectors+1):
0157                 binHisto = (st - 1)*60 + (wh + 2)*nSectors + sec
0158                 for idx in range( len(histos) ):
0159                     histo = histos[idx]
0160                     value = histo.GetBinContent( binHisto )
0161                     error = histo.GetBinError( binHisto )
0162                     averages[(st,wh)]     += value/( error*error )
0163                     averagesSumw[(st,wh)] += 1./( error*error )
0164             # Average per (st,wh)
0165             averages[(st,wh)] = averages[(st,wh)]/averagesSumw[(st,wh)]
0166             averagesErr[(st,wh)] = math.sqrt( 1./averagesSumw[(st,wh)] )
0167             histoAverage.SetBinContent(binHistoAve,averages[(st,wh)])
0168             histoAverage.SetBinError(binHistoAve,averagesErr[(st,wh)])
0169             print("Station %d, Wheel %d: %.4f +/- %.6f" % (st,wh,averages[(st,wh)],averagesErr[(st,wh)]))
0170
0171     canvasAverage = ROOT.TCanvas("c_" + histoAverage.GetName())
0172     canvasAverage.SetGridx()
0173     canvasAverage.SetGridy()
0174     canvasAverage.SetFillColor( 0 )
0175     canvasAverage.cd()
0176     mean_ymin = -0.02
0177     mean_ymax =  0.02
0178     sig_ymin = 0.
0179     sig_ymax = 0.1
0180     if type == 'mean':
0181         histoAverage.GetYaxis().SetTitle("Mean of residuals (cm)")
0182         histoAverage.GetYaxis().SetRangeUser(mean_ymin,mean_ymax)
0183     elif type == 'sigma':
0184         histoAverage.GetYaxis().SetTitle("Resolution (cm)")
0185         histoAverage.GetYaxis().SetRangeUser(sig_ymin,sig_ymax)
0186
0187     histoAverage.SetStats(0)
0188     histoAverage.SetLineWidth(2)
0189     histoAverage.SetMarkerStyle( 27 )
0190     histoAverage.SetMarkerSize( 1.5 )
0191     histoAverage.LabelsOption("d","X")
0192     histoAverage.Draw("E2")
0193
0194     return ( (canvas,canvasAverage),(histos,histoAverage),objects )
0195
0196 def plotMean(fileName,sl,dir='DQMData/Run 1/DT/Run summary/DTCalibValidation',option='HISTOPE1'):
0197     type = 'mean'
0198     objs = plot(fileName,sl,dir,type,option)
0199     return objs
0200
0201 def plotSigma(fileName,sl,dir='DQMData/Run 1/DT/Run summary/DTCalibValidation',option='HISTOPE1'):
0202     type = 'sigma'
0203     objs = plot(fileName,sl,dir,type,option)
0204     return objs
0205
0206 def plotSigmaAll(fileName,dir='DQMData/Run 1/DT/Run summary/DTCalibValidation',option='HISTOPE1',outputFileName=''):
0207     colors = (2,4,12,44,55,38,27,46)
0208     markers = (24,25,26,27,28,30,32,5)
0209
0210     slList = (1,2,3)
0211     labels = ('R-#phi SL1','R-z SL2','R-#phi SL3')
0212     canvas = None
0213     objects = None
0214     histos = []
0215     idx = 0
0216     for sl in slList:
0217         draw = False
0218         if not idx: draw = True
0219
0220         objs = plotSigma(fileName,sl,dir,option)
0221         histos.append(objs[1][1])
0222         histos[-1].SetName( "%s_%d" % (histos[-1].GetName(),idx) )
0223         if not idx:
0224             canvas = objs[0][1]
0225             #objects = objs[2][1]
0226
0227         canvas.cd()
0228         if idx:
0229             histos[-1].SetLineColor(colors[ (idx - 1) % len(colors) ])
0230             histos[-1].SetMarkerColor(colors[ (idx - 1) % len(colors) ])
0231             histos[-1].SetMarkerStyle(markers[ (idx - 1) % len(markers) ])
0232
0233             histos[-1].Draw(option + "SAME")
0234
0235         idx += 1
0236
0237     legend = ROOT.TLegend(0.4,0.7,0.95,0.8)
0238     for idx in range( len(histos) ):
0239         histo = histos[idx]
0240         label = histo.GetName()
0241         if len(labels): label = labels[idx]
0242         legend.AddEntry(histo,label,"LP")
0243
0244         idx += 1
0245
0246     canvas.cd()
0247     legend.SetFillColor( canvas.GetFillColor() )
0248     legend.Draw("SAME")
0249     if not objects: objects = [legend]
0250     else:           objects.append(legend)
0251
0252     if outputFileName:
0253         outputFile = ROOT.TFile(outputFileName,'recreate')
0254         outputFile.cd()
0255         for histo in histos: histo.Write()
0256         outputFile.Close()
0257         return 0
0258     else:
0259         return (canvas,histos,objects)
0260
0261 #def plotDataVsMCFromFile(fileNameData,fileNameMC,labels=[]):
0262 def plotFromFile(fileNames,labels=[]):
0263
0264     AddDirectoryStatus_ = ROOT.TH1.AddDirectoryStatus()
0265     ROOT.TH1.AddDirectory(False)
0266
0267     #fileData = ROOT.TFile(fileNameData,'read')
0268     #fileMC = ROOT.TFile(fileNameMC,'read')
0269     rootFiles = []
0270     for file in fileNames: rootFiles.append( ROOT.TFile(file,'read') )
0271
0272     variables = ['h_AverageAll_SL1_0',
0273                  'h_AverageAll_SL2_1',
0274                  'h_AverageAll_SL3_2']
0275
0276     colors = (1,2,4,12,44,55,38,27,46)
0277     markers = (20,24,25,26,27,28,30,32,5)
0278     objects = None
0279     canvases = []
0280     legends = []
0281     histos = []
0282     idx_var = 0
0283     for var in variables:
0284         print("Accessing",var)
0285         #histoData = fileData.Get(var)
0286         #histoData.SetName(histoData.GetName() + "_data")
0287         #histoMC = fileMC.Get(var)
0288         #histoMC.SetName(histoMC.GetName() + "_mc")
0289         #histoData.SetLineColor(1)
0290         #histoData.SetMarkerStyle(20)
0291         #histoData.SetMarkerSize(1.4)
0292         #histoData.SetMarkerColor(1)
0293
0294         #histoMC.SetLineColor(2)
0295         #histoMC.SetMarkerStyle(24)
0296         #histoMC.SetMarkerSize(1.4)
0297         #histoMC.SetMarkerColor(2)
0298
0299         histos_tmp = []
0300         idx = 0
0301         for file in rootFiles:
0302             histos_tmp.append( file.Get(var) )
0303             histos_tmp[-1].SetName( "%s_%d" % (histos_tmp[-1].GetName(),idx) )
0304             print("Created",histos_tmp[-1].GetName())
0305             histos_tmp[-1].SetLineColor(colors[ idx % len(colors) ])
0306             histos_tmp[-1].SetMarkerColor(colors[ idx % len(colors) ])
0307             histos_tmp[-1].SetMarkerStyle(markers[ idx % len(markers) ])
0308             histos_tmp[-1].SetMarkerSize(1.4)
0309             idx += 1
0310         histos.append( histos_tmp )
0311
0312         canvases.append( ROOT.TCanvas("c_" + var,var) )
0313         canvases[-1].SetGridx()
0314         canvases[-1].SetGridy()
0315         canvases[-1].SetFillColor(0)
0316         canvases[-1].cd()
0317         #histoData.Draw()
0318         #histoMC.Draw("SAME")
0319         #histos.append( (histoData,histoMC) )
0320         histos[-1][0].Draw()
0321         for histo in histos[-1][1:]: histo.Draw("SAME")
0322
0323         if len(labels):
0324             #labelData = labels[0]
0325             #labelMC = labels[1]
0326             legends.append( ROOT.TLegend(0.4,0.7,0.95,0.8) )
0327             idx = 0
0328             for histo in histos[-1]:
0329                 legends[-1].AddEntry(histo,labels[idx],"LP")
0330                 idx += 1
0331
0332             legends[-1].SetFillColor( canvases[-1].GetFillColor() )
0333             legends[-1].Draw("SAME")
0334
0335         idx_var += 1
0336
0337     if not objects: objects = [legends]
0338     else:           objects.append(legends)
0339
0340     ROOT.TH1.AddDirectory(AddDirectoryStatus_)
0341
0342     return (canvases,histos,objects)