python/PlottingTools/plotResidualsPerLayer.py

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