CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
plotResidualsPerLayer.py
Go to the documentation of this file.
1 from __future__ import print_function
2 from __future__ import absolute_import
3 import ROOT
4 from .fitResidual import fitResidual
5 from .drawHistoAllChambers import drawHisto
6 
7 layerCorrectionFactors = {'SL1':(1.17,1.16,1.15,1.14),
8  'SL2':(1.83,1.20,1.20,1.83),
9  'SL3':(1.14,1.15,1.16,1.17)}
10 
11 def plotResLayer(fileName,sl,layer,
12  dir='DQMData/Run 1/DT/Run summary/DTCalibValidation',
13  option="HISTOPE1",draw=True):
14 
15  mean_ymin = -0.02
16  mean_ymax = 0.02
17  sig_ymin = 0.
18  sig_ymax = 0.1
19 
20  slType = sl
21  slStr = "SL%d" % slType
22  layerType = layer
23  layerStr = "Layer%d" % layerType
24  verbose = False
25  nSigmas = 2
26 
27  ROOT.TH1.AddDirectory(False)
28 
29  file = ROOT.TFile(fileName,'read')
30 
31  wheels = (-2,-1,0,1,2)
32  stations = (1,2,3,4)
33 
34  # (Wh-2 MB1 Sec1 ... Wh-2 MB1 Sec12 ... Wh-1 MB1 Sec1 ... Wh-1 MB1 Sec12 ...)
35  # (Wh-2 MB2 Sec1 ... Wh-2 MB2 Sec12 ... Wh-1 MB2 Sec1 ... Wh-1 MB1 Sec12 ...) ...
36  nBins = 250
37  if slType == 2: nBins = 180
38  histoMean = ROOT.TH1F("h_ResMeanAll_%s_%s" % (slStr,layerStr),"Mean of residuals",nBins,0,nBins)
39  histoSigma = ROOT.TH1F("h_ResSigmaAll_%s_%s" % (slStr,layerStr),"Sigma of residuals",nBins,0,nBins)
40  for st in stations:
41  nSectors = 12
42  if st == 4: nSectors = 14
43  if st == 4 and slType == 2: continue
44  if verbose: print("Station",st)
45  for wh in wheels:
46  if verbose: print("Wheel",wh)
47  for sec in range(1,nSectors+1):
48  if verbose: print("Sector",sec)
49  # Get histogram
50  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)
51  print("Accessing",histoName)
52  histo = file.Get(histoName)
53  (histo,fitFunc) = fitResidual(histo,nSigmas,verbose)
54  fitMean = fitFunc.GetParameter(1)
55  fitMeanErr = fitFunc.GetParError(1)
56  fitSigma = fitFunc.GetParameter(2)
57  fitSigmaErr = fitFunc.GetParError(2)
58 
59  layerIdx = (layer - 1)
60  corrFactor = layerCorrectionFactors[slStr][layerIdx]
61 
62  binHistoNew = (st - 1)*60 + (wh + 2)*nSectors + sec
63  if verbose: print("Bin in summary histo",binHistoNew)
64  histoMean.SetBinContent(binHistoNew,fitMean)
65  histoMean.SetBinError(binHistoNew,fitMeanErr)
66  histoSigma.SetBinContent(binHistoNew,fitSigma*corrFactor)
67  histoSigma.SetBinError(binHistoNew,fitSigmaErr*corrFactor)
68 
69  if sec == 1:
70  label = "Wheel %d" % wh
71  if wh == -2: label += " MB%d" % st
72  histoMean.GetXaxis().SetBinLabel(binHistoNew,label)
73  histoSigma.GetXaxis().SetBinLabel(binHistoNew,label)
74 
75  objectsMean = drawHisto(histoMean,title="Mean of residuals (cm)",
76  ymin=mean_ymin,ymax=mean_ymax,option=option,draw=draw)
77  objectsSigma = drawHisto(histoSigma,title="Resolution (cm)",
78  ymin=sig_ymin,ymax=sig_ymax,option=option,draw=draw)
79 
80  return (objectsMean,objectsSigma)
81 
82 def plot(fileName,sl,
83  dir='DQMData/Run 1/DT/Run summary/DTCalibValidation',type='mean',option='HISTOPE1'):
84  colors = (2,4,12,44,55,38,27,46)
85  markers = (24,25,26,27,28,30,32,5)
86  labels=['Layer 1','Layer 2','Layer 3','Layer 4']
87  idx_type = None
88  if type == 'mean': idx_type = 0
89  elif type == 'sigma': idx_type = 1
90  else: raise RuntimeError("Wrong option: %s" % type)
91 
92  idx = 0
93  canvas = None
94  objects = None
95  histos = []
96  for layer in range(1,5):
97  draw = False
98  if not idx: draw = True
99 
100  objs = plotResLayer(fileName,sl,layer,dir,option,draw)
101  histos.append(objs[idx_type][1])
102  histos[-1].SetName( "%s_%d" % (histos[-1].GetName(),idx) )
103  if not idx:
104  canvas = objs[idx_type][0]
105  objects = objs[idx_type][2]
106 
107  canvas.cd()
108  if idx:
109  histos[-1].SetLineColor(colors[ (idx - 1) % len(colors) ])
110  histos[-1].SetMarkerColor(colors[ (idx - 1) % len(colors) ])
111  histos[-1].SetMarkerStyle(markers[ (idx - 1) % len(markers) ])
112 
113  histos[-1].Draw(option + "SAME")
114 
115  idx += 1
116 
117  legend = ROOT.TLegend(0.4,0.7,0.95,0.8)
118  for idx in range( len(histos) ):
119  histo = histos[idx]
120  label = histo.GetName()
121  if len(labels): label = labels[idx]
122  legend.AddEntry(histo,label,"LP")
123 
124  idx += 1
125 
126  canvas.cd()
127  legend.SetFillColor( canvas.GetFillColor() )
128  legend.Draw("SAME")
129  objects.append(legend)
130 
131  # Compute averages
132  # (Wh-2 MB1 Sec1 ... Wh-2 MB1 Sec12 ... Wh-1 MB1 Sec1 ... Wh-1 MB1 Sec12 ...)
133  # (Wh-2 MB2 Sec1 ... Wh-2 MB2 Sec12 ... Wh-1 MB2 Sec1 ... Wh-1 MB1 Sec12 ...) ...
134  import math
135  wheels = (-2,-1,0,1,2)
136  stations = (1,2,3,4)
137  slType = sl
138  slStr = "SL%d" % slType
139 
140  nBinsAve = len(stations)*len(wheels)
141  histoAverage = ROOT.TH1F("h_AverageAll_" + slStr,"",nBinsAve,0,nBinsAve)
142  averages = {}
143  averagesErr = {}
144  averagesSumw = {}
145  print("Averages:")
146  for st in stations:
147  nSectors = 12
148  if st == 4: nSectors = 14
149  if st == 4 and slType == 2: continue
150  for wh in wheels:
151  binHistoAve = (st - 1)*5 + (wh + 2) + 1
152  label = "Wheel %d" % wh
153  if wh == -2: label += " MB%d" % st
154  histoAverage.GetXaxis().SetBinLabel(binHistoAve,label)
155 
156  averages[(st,wh)] = 0.
157  averagesSumw[(st,wh)] = 0.
158  for sec in range(1,nSectors+1):
159  binHisto = (st - 1)*60 + (wh + 2)*nSectors + sec
160  for idx in range( len(histos) ):
161  histo = histos[idx]
162  value = histo.GetBinContent( binHisto )
163  error = histo.GetBinError( binHisto )
164  averages[(st,wh)] += value/( error*error )
165  averagesSumw[(st,wh)] += 1./( error*error )
166  # Average per (st,wh)
167  averages[(st,wh)] = averages[(st,wh)]/averagesSumw[(st,wh)]
168  averagesErr[(st,wh)] = math.sqrt( 1./averagesSumw[(st,wh)] )
169  histoAverage.SetBinContent(binHistoAve,averages[(st,wh)])
170  histoAverage.SetBinError(binHistoAve,averagesErr[(st,wh)])
171  print("Station %d, Wheel %d: %.4f +/- %.6f" % (st,wh,averages[(st,wh)],averagesErr[(st,wh)]))
172 
173  canvasAverage = ROOT.TCanvas("c_" + histoAverage.GetName())
174  canvasAverage.SetGridx()
175  canvasAverage.SetGridy()
176  canvasAverage.SetFillColor( 0 )
177  canvasAverage.cd()
178  mean_ymin = -0.02
179  mean_ymax = 0.02
180  sig_ymin = 0.
181  sig_ymax = 0.1
182  if type == 'mean':
183  histoAverage.GetYaxis().SetTitle("Mean of residuals (cm)")
184  histoAverage.GetYaxis().SetRangeUser(mean_ymin,mean_ymax)
185  elif type == 'sigma':
186  histoAverage.GetYaxis().SetTitle("Resolution (cm)")
187  histoAverage.GetYaxis().SetRangeUser(sig_ymin,sig_ymax)
188 
189  histoAverage.SetStats(0)
190  histoAverage.SetLineWidth(2)
191  histoAverage.SetMarkerStyle( 27 )
192  histoAverage.SetMarkerSize( 1.5 )
193  histoAverage.LabelsOption("d","X")
194  histoAverage.Draw("E2")
195 
196  return ( (canvas,canvasAverage),(histos,histoAverage),objects )
197 
198 def plotMean(fileName,sl,dir='DQMData/Run 1/DT/Run summary/DTCalibValidation',option='HISTOPE1'):
199  type = 'mean'
200  objs = plot(fileName,sl,dir,type,option)
201  return objs
202 
203 def plotSigma(fileName,sl,dir='DQMData/Run 1/DT/Run summary/DTCalibValidation',option='HISTOPE1'):
204  type = 'sigma'
205  objs = plot(fileName,sl,dir,type,option)
206  return objs
207 
208 def plotSigmaAll(fileName,dir='DQMData/Run 1/DT/Run summary/DTCalibValidation',option='HISTOPE1',outputFileName=''):
209  colors = (2,4,12,44,55,38,27,46)
210  markers = (24,25,26,27,28,30,32,5)
211 
212  slList = (1,2,3)
213  labels = ('R-#phi SL1','R-z SL2','R-#phi SL3')
214  canvas = None
215  objects = None
216  histos = []
217  idx = 0
218  for sl in slList:
219  draw = False
220  if not idx: draw = True
221 
222  objs = plotSigma(fileName,sl,dir,option)
223  histos.append(objs[1][1])
224  histos[-1].SetName( "%s_%d" % (histos[-1].GetName(),idx) )
225  if not idx:
226  canvas = objs[0][1]
227  #objects = objs[2][1]
228 
229  canvas.cd()
230  if idx:
231  histos[-1].SetLineColor(colors[ (idx - 1) % len(colors) ])
232  histos[-1].SetMarkerColor(colors[ (idx - 1) % len(colors) ])
233  histos[-1].SetMarkerStyle(markers[ (idx - 1) % len(markers) ])
234 
235  histos[-1].Draw(option + "SAME")
236 
237  idx += 1
238 
239  legend = ROOT.TLegend(0.4,0.7,0.95,0.8)
240  for idx in range( len(histos) ):
241  histo = histos[idx]
242  label = histo.GetName()
243  if len(labels): label = labels[idx]
244  legend.AddEntry(histo,label,"LP")
245 
246  idx += 1
247 
248  canvas.cd()
249  legend.SetFillColor( canvas.GetFillColor() )
250  legend.Draw("SAME")
251  if not objects: objects = [legend]
252  else: objects.append(legend)
253 
254  if outputFileName:
255  outputFile = ROOT.TFile(outputFileName,'recreate')
256  outputFile.cd()
257  for histo in histos: histo.Write()
258  outputFile.Close()
259  return 0
260  else:
261  return (canvas,histos,objects)
262 
263 #def plotDataVsMCFromFile(fileNameData,fileNameMC,labels=[]):
264 def plotFromFile(fileNames,labels=[]):
265 
266  AddDirectoryStatus_ = ROOT.TH1.AddDirectoryStatus()
267  ROOT.TH1.AddDirectory(False)
268 
269  #fileData = ROOT.TFile(fileNameData,'read')
270  #fileMC = ROOT.TFile(fileNameMC,'read')
271  rootFiles = []
272  for file in fileNames: rootFiles.append( ROOT.TFile(file,'read') )
273 
274  variables = ['h_AverageAll_SL1_0',
275  'h_AverageAll_SL2_1',
276  'h_AverageAll_SL3_2']
277 
278  colors = (1,2,4,12,44,55,38,27,46)
279  markers = (20,24,25,26,27,28,30,32,5)
280  objects = None
281  canvases = []
282  legends = []
283  histos = []
284  idx_var = 0
285  for var in variables:
286  print("Accessing",var)
287  #histoData = fileData.Get(var)
288  #histoData.SetName(histoData.GetName() + "_data")
289  #histoMC = fileMC.Get(var)
290  #histoMC.SetName(histoMC.GetName() + "_mc")
291  #histoData.SetLineColor(1)
292  #histoData.SetMarkerStyle(20)
293  #histoData.SetMarkerSize(1.4)
294  #histoData.SetMarkerColor(1)
295 
296  #histoMC.SetLineColor(2)
297  #histoMC.SetMarkerStyle(24)
298  #histoMC.SetMarkerSize(1.4)
299  #histoMC.SetMarkerColor(2)
300 
301  histos_tmp = []
302  idx = 0
303  for file in rootFiles:
304  histos_tmp.append( file.Get(var) )
305  histos_tmp[-1].SetName( "%s_%d" % (histos_tmp[-1].GetName(),idx) )
306  print("Created",histos_tmp[-1].GetName())
307  histos_tmp[-1].SetLineColor(colors[ idx % len(colors) ])
308  histos_tmp[-1].SetMarkerColor(colors[ idx % len(colors) ])
309  histos_tmp[-1].SetMarkerStyle(markers[ idx % len(markers) ])
310  histos_tmp[-1].SetMarkerSize(1.4)
311  idx += 1
312  histos.append( histos_tmp )
313 
314  canvases.append( ROOT.TCanvas("c_" + var,var) )
315  canvases[-1].SetGridx()
316  canvases[-1].SetGridy()
317  canvases[-1].SetFillColor(0)
318  canvases[-1].cd()
319  #histoData.Draw()
320  #histoMC.Draw("SAME")
321  #histos.append( (histoData,histoMC) )
322  histos[-1][0].Draw()
323  for histo in histos[-1][1:]: histo.Draw("SAME")
324 
325  if len(labels):
326  #labelData = labels[0]
327  #labelMC = labels[1]
328  legends.append( ROOT.TLegend(0.4,0.7,0.95,0.8) )
329  idx = 0
330  for histo in histos[-1]:
331  legends[-1].AddEntry(histo,labels[idx],"LP")
332  idx += 1
333 
334  legends[-1].SetFillColor( canvases[-1].GetFillColor() )
335  legends[-1].Draw("SAME")
336 
337  idx_var += 1
338 
339  if not objects: objects = [legends]
340  else: objects.append(legends)
341 
342  ROOT.TH1.AddDirectory(AddDirectoryStatus_)
343 
344  return (canvases,histos,objects)
sistrip::SpyUtilities::range
const uint16_t range(const Frame &aFrame)
Definition: SiStripSpyUtilities.cc:76
print
void print(TMatrixD &m, const char *label=nullptr, bool mathematicaFormat=false)
Definition: Utilities.cc:47
Generated for CMSSW Reference Manual by doxygen 1.8.5