File indexing completed on 2024-11-28 23:11:15
0001
0002
0003
0004
0005 import sys
0006 import numpy as np
0007 import pandas as pd
0008 import matplotlib.pyplot as plt
0009
0010 from array import array
0011 oldargv = sys.argv[:]
0012 sys.argv = [ '-b-' ]
0013 from ROOT import TCanvas, TLegend, TPaveText, THStack, TFile, TLatex, TPaveLabel
0014 from ROOT import TProfile, TProfile2D, TH1D, TH2D, TH2F, TPaletteAxis, TH1, TH1F, TColor, TExec, TLine
0015 from ROOT import kBlack, kWhite, kOrange, kAzure, kBlue
0016 from ROOT import gROOT, gStyle
0017 gROOT.SetBatch(True)
0018 sys.argv = oldargv
0019
0020 from Validation.Geometry.plot_utils import setTDRStyle, Plot_params, drawEtaValues
0021
0022 from Validation.Geometry.plot_hgcal_utils import plots, COMPOUNDS, DETECTORS, sDETS, hist_label_to_num, TwikiPrintout, acustompalette, dEdx, MatXo, drawHalfEtaValues
0023 from collections import namedtuple, OrderedDict
0024 import sys, os
0025 import argparse
0026
0027 def paramsGood_(detector, plot):
0028 """Check the validity of the arguments.
0029
0030 Common function to check the validity of the parameters passed
0031 in. It returns a tuple composed by a bool and a string. The
0032 bool indicates if all checks are ok, the string the name of the
0033 appropriate ROOT file to open (empty string in case the any
0034 check failed)
0035
0036 """
0037
0038 if plot not in plots.keys():
0039 print("Error, unknown plot %s" % plot)
0040 return (False, '')
0041
0042 if detector not in DETECTORS and detector not in COMPOUNDS.keys():
0043 print('Error, unknown detector: %s' % detector)
0044 return (False, '')
0045
0046 theDetectorFilename = ''
0047 if detector in DETECTORS:
0048 theDetectorFilename = 'matbdg_%s.root' % detector
0049 else:
0050 theDetectorFilename = 'matbdg_%s.root' % COMPOUNDS[detector][0]
0051
0052 if not checkFile_(theDetectorFilename):
0053 print("Error, missing file %s" % theDetectorFilename)
0054 raise RuntimeError
0055 return (True, theDetectorFilename)
0056
0057 def checkFile_(filename):
0058 return os.path.exists(filename)
0059
0060 def setColorIfExists_(histos, h, color):
0061 if h in histos.keys():
0062 histos[h].SetFillColor(color)
0063
0064 def assignOrAddIfExists_(h, p):
0065 """Assign the projection of p to h.
0066
0067 Function to assign the projection of p to h, in the case in
0068 which h is None, otherwise add the projection to the already
0069 valid h object
0070
0071 """
0072
0073 if not h:
0074 h = p.ProjectionX()
0075 else:
0076 h.Add(p.ProjectionX("B_%s" % h.GetName()), +1.000)
0077 return h
0078
0079 def createPlots_(plot, compounddetectorname):
0080 """Cumulative material budget from simulation.
0081
0082 Internal function that will produce a cumulative profile of the
0083 material budget inferred from the simulation starting from the
0084 single detectors that compose the tracker. It will iterate over
0085 all existing detectors contained in the DETECTORS
0086 dictionary. The function will automatically skip non-existent
0087 detectors.
0088
0089 """
0090
0091 theDirname = "Figures"
0092
0093 hist_X0_detectors = OrderedDict()
0094 if plot not in plots.keys():
0095 print("Error: chosen plot name not known %s" % plot)
0096 return
0097
0098
0099
0100 subDetectorFiles = []
0101
0102 hist_X0_elements = OrderedDict()
0103 prof_X0_elements = OrderedDict()
0104 for subDetector,color in DETECTORS.items():
0105 subDetectorFilename = "matbdg_%s.root" % subDetector
0106 if not checkFile_(subDetectorFilename):
0107 print("Error opening file: %s" % subDetectorFilename)
0108 continue
0109
0110 subDetectorFiles.append(TFile(subDetectorFilename))
0111 subDetectorFile = subDetectorFiles[-1]
0112 print ("Opening file: %s" % subDetectorFilename)
0113 prof_X0_XXX = subDetectorFile.Get("%d" % plots[plot].plotNumber)
0114
0115 hist_X0_detectors[subDetector] = prof_X0_XXX.ProjectionX()
0116
0117
0118 for label, [num, color, leg] in hist_label_to_num.items():
0119 prof_X0_elements[label] = subDetectorFile.Get("%d" % (num + plots[plot].plotNumber))
0120 hist_X0_elements[label] = assignOrAddIfExists_(hist_X0_elements.setdefault(label, None),
0121 prof_X0_elements[label])
0122
0123 cumulative_matbdg = TH1D("CumulativeSimulMatBdg",
0124 "CumulativeSimulMatBdg",
0125 hist_X0_detectors["BeamPipe"].GetNbinsX(),
0126 hist_X0_detectors["BeamPipe"].GetXaxis().GetXmin(),
0127 hist_X0_detectors["BeamPipe"].GetXaxis().GetXmax())
0128 cumulative_matbdg.SetDirectory(0)
0129
0130
0131 for det, color in DETECTORS.items():
0132 setColorIfExists_(hist_X0_detectors, det, color)
0133
0134 for label, [num, color, leg] in hist_label_to_num.items():
0135 hist_X0_elements[label].SetFillColor(color)
0136
0137
0138
0139 stackTitle_SubDetectors = "Material Budget;%s;%s" % (
0140 plots[plot].abscissa,plots[plot].ordinate)
0141 stack_X0_SubDetectors = THStack("stack_X0",stackTitle_SubDetectors)
0142 for det, histo in hist_X0_detectors.items():
0143 stack_X0_SubDetectors.Add(histo)
0144 cumulative_matbdg.Add(histo, 1)
0145
0146
0147 can_SubDetectors = TCanvas("can_SubDetectors","can_SubDetectors",800,800)
0148
0149 can_SubDetectors.SetFillColor(kWhite)
0150
0151
0152 stack_X0_SubDetectors.SetMinimum(plots[plot].ymin)
0153 stack_X0_SubDetectors.SetMaximum(plots[plot].ymax)
0154 stack_X0_SubDetectors.Draw("HIST")
0155
0156
0157
0158
0159 theLegend_SubDetectors = TLegend(0.130,0.7,0.93,0.90)
0160 theLegend_SubDetectors.SetNColumns(2)
0161 theLegend_SubDetectors.SetFillColor(0)
0162 theLegend_SubDetectors.SetFillStyle(0)
0163 theLegend_SubDetectors.SetBorderSize(0)
0164
0165 for det, histo in hist_X0_detectors.items():
0166 theLegend_SubDetectors.AddEntry(histo, det, "f")
0167
0168 theLegend_SubDetectors.Draw()
0169
0170
0171 text_SubDetectors = TPaveText(0.130,0.627,0.352,0.687,"NDC")
0172 text_SubDetectors.SetFillColor(0)
0173 text_SubDetectors.SetBorderSize(0)
0174 text_SubDetectors.AddText("CMS Simulation")
0175 text_SubDetectors.SetTextAlign(11)
0176 text_SubDetectors.Draw()
0177
0178
0179 can_SubDetectors.Update()
0180 if not checkFile_(theDirname):
0181 os.mkdir(theDirname)
0182 can_SubDetectors.SaveAs("%s/MaterialBdg_%s_%s.pdf" % (theDirname, compounddetectorname,plot))
0183 can_SubDetectors.SaveAs("%s/MaterialBdg_%s_%s.png" % (theDirname, compounddetectorname,plot))
0184 can_SubDetectors.SaveAs("%s/MaterialBdg_%s_%s.root" % (theDirname, compounddetectorname,plot))
0185
0186 if plot == "x_vs_eta" or plot == "l_vs_eta":
0187 canname = "MBCan_1D_%s_%s_total" % (compounddetectorname, plot)
0188 can2 = TCanvas(canname, canname, 800, 800)
0189 can2.Range(0,0,25,25)
0190 can2.SetFillColor(kWhite)
0191 gStyle.SetOptStat(0)
0192 gStyle.SetOptTitle(0);
0193
0194 stack_X0_SubDetectors.GetStack().Last().SetMarkerStyle(34)
0195 stack_X0_SubDetectors.GetStack().Last().GetXaxis().SetRangeUser( 1.0, 3.5)
0196 stack_X0_SubDetectors.GetStack().Last().Draw();
0197 stack_X0_SubDetectors.GetYaxis().SetTitleOffset(1.15);
0198 can2.Update()
0199 can2.Modified()
0200 can2.SaveAs("%s/%s_%s_total_Zplus.pdf" % (theDirname, compounddetectorname, plot))
0201 can2.SaveAs("%s/%s_%s_total_Zplus.png" % (theDirname, compounddetectorname, plot))
0202 stack_X0_SubDetectors.GetStack().Last().GetXaxis().SetRangeUser( -3.5, -1.0)
0203 stack_X0_SubDetectors.GetStack().Last().Draw();
0204 stack_X0_SubDetectors.GetYaxis().SetTitleOffset(1.15);
0205 can2.Update()
0206 can2.Modified()
0207 can2.SaveAs("%s/%s_%s_total_Zminus.pdf" % (theDirname, compounddetectorname, plot))
0208 can2.SaveAs("%s/%s_%s_total_Zminus.png" % (theDirname, compounddetectorname, plot))
0209
0210
0211 etavalues = []
0212 matbudginX0 = []
0213 matbudginIntLen = []
0214 for binx in range(0,stack_X0_SubDetectors.GetStack().Last().GetXaxis().GetNbins()):
0215 bincontent = stack_X0_SubDetectors.GetStack().Last().GetBinContent(binx)
0216 if bincontent == 0: continue
0217 etavalues.append( stack_X0_SubDetectors.GetStack().Last().GetBinCenter(binx) )
0218 if plot == "x_vs_eta":
0219 matbudginX0.append( bincontent )
0220 d1 = {'Eta': etavalues, 'MatBudInX0': matbudginX0}
0221 df1 = pd.DataFrame(data=d1).round(2)
0222 df1.to_csv(r'/afs/cern.ch/work/a/apsallid/CMS/PFCalStudies/CMS-HGCAL/matbudV10fromVertexToBackofHGCal/CMSSW_11_0_X_2019-06-04-2300/src/Validation/Geometry/test/EtavsMatBudinXo.txt',sep=' ', index=False, header=False)
0223
0224 if plot == "l_vs_eta":
0225 matbudginIntLen.append( bincontent )
0226 d2 = {'Eta': etavalues, 'MatBudInIntLen': matbudginIntLen}
0227 df2 = pd.DataFrame(data=d2).round(2)
0228 df2.to_csv(r'/afs/cern.ch/work/a/apsallid/CMS/PFCalStudies/CMS-HGCAL/matbudV10fromVertexToBackofHGCal/CMSSW_11_0_X_2019-06-04-2300/src/Validation/Geometry/test/EtavsMatBudInIntLen.txt',sep=' ', index=False, header=False)
0229
0230
0231 return cumulative_matbdg
0232
0233 def createPlots2D_(plot, compounddetectorname):
0234 """2D material budget map to know exactly what we are adding.
0235 """
0236
0237
0238 theDirname = "Figures"
0239
0240 hist_X0_detectors = OrderedDict()
0241 if plot not in plots.keys():
0242 print("Error: chosen plot name not known %s" % plot)
0243 return
0244
0245
0246
0247 subDetectorFiles = []
0248
0249 hist_X0_elements = OrderedDict()
0250 prof_X0_elements = OrderedDict()
0251
0252 for subDetector,color in DETECTORS.items():
0253 subDetectorFilename = "matbdg_%s.root" % subDetector
0254 if not checkFile_(subDetectorFilename):
0255 print("Error opening file: %s" % subDetectorFilename)
0256 continue
0257
0258 subDetectorFiles.append(TFile(subDetectorFilename))
0259 subDetectorFile = subDetectorFiles[-1]
0260 print ("Opening file: %s" % subDetectorFilename)
0261 prof_X0_XXX = subDetectorFile.Get("%d" % plots[plot].plotNumber)
0262
0263
0264 hist_X0_detectors[subDetector] = prof_X0_XXX.ProjectionXY("_pxy","B")
0265 print(subDetector)
0266
0267
0268
0269
0270 minX = 1.03*hist_X0_detectors["BeamPipe"].GetXaxis().GetXmin()
0271 maxX = 1.03*hist_X0_detectors["BeamPipe"].GetXaxis().GetXmax()
0272 minY = 1.03*hist_X0_detectors["BeamPipe"].GetYaxis().GetXmin()
0273 maxY = 1.03*hist_X0_detectors["BeamPipe"].GetYaxis().GetXmax()
0274
0275 frame = TH2F("frame", "", 10, minX, maxX, 10, minY, maxY);
0276 frame.SetMinimum(0.1)
0277 frame.SetMaximum(10.)
0278 frame.GetXaxis().SetTickLength(frame.GetXaxis().GetTickLength()*0.50)
0279 frame.GetYaxis().SetTickLength(frame.GetXaxis().GetTickLength()/4.)
0280
0281 hist2d_X0_total = hist_X0_detectors["BeamPipe"]
0282
0283
0284 hist2dTitle = ('%s %s;%s;%s;%s' % (plots[plot].quotaName,
0285 "All detectors",
0286 plots[plot].abscissa,
0287 plots[plot].ordinate,
0288 plots[plot].quotaName))
0289
0290 hist2d_X0_total.SetTitle(hist2dTitle)
0291 frame.SetTitle(hist2dTitle)
0292 frame.SetTitleOffset(0.5,"Y")
0293
0294
0295
0296 if plots[plot].histoMin != -1.:
0297 hist2d_X0_total.SetMinimum(plots[plot].histoMin)
0298 if plots[plot].histoMax != -1.:
0299 hist2d_X0_total.SetMaximum(plots[plot].histoMax)
0300
0301
0302
0303 can_SubDetectors = TCanvas("can_SubDetectors","can_SubDetectors",2480+248, 580+58+58)
0304 can_SubDetectors.SetTopMargin(0.1)
0305 can_SubDetectors.SetBottomMargin(0.1)
0306 can_SubDetectors.SetLeftMargin(0.04)
0307 can_SubDetectors.SetRightMargin(0.06)
0308 can_SubDetectors.SetFillColor(kWhite)
0309 gStyle.SetOptStat(0)
0310 gStyle.SetTitleFillColor(0)
0311 gStyle.SetTitleBorderSize(0)
0312 gStyle.SetOptTitle(0)
0313
0314 hist2d_X0_total.GetYaxis().SetTickLength(hist2d_X0_total.GetXaxis().GetTickLength()/4.)
0315 hist2d_X0_total.GetYaxis().SetTickLength(hist2d_X0_total.GetXaxis().GetTickLength()/4.)
0316 hist2d_X0_total.SetTitleOffset(0.5,"Y")
0317 hist2d_X0_total.GetYaxis().SetTitleOffset(0.50);
0318
0319
0320
0321
0322
0323
0324 for det, color in DETECTORS.items():
0325 hist_X0_detectors[det].SetMarkerColor(color)
0326 hist_X0_detectors[det].SetFillColor(color)
0327
0328 for det, histo in hist_X0_detectors.items():
0329 print(det)
0330 histo.Draw("same")
0331
0332
0333 theLegend_SubDetectors = TLegend(0.100,0.7,0.90,0.90)
0334 theLegend_SubDetectors.SetNColumns(3)
0335 theLegend_SubDetectors.SetFillColor(0)
0336 theLegend_SubDetectors.SetFillStyle(0)
0337 theLegend_SubDetectors.SetBorderSize(0)
0338
0339 for det, histo in hist_X0_detectors.items():
0340 theLegend_SubDetectors.AddEntry(histo, det, "f")
0341
0342
0343 theLegend_SubDetectors.Draw()
0344
0345
0346
0347 text_SubDetectors = TPaveText(0.100,0.627,0.322,0.687,"NDC")
0348 text_SubDetectors.SetFillColor(0)
0349 text_SubDetectors.SetBorderSize(0)
0350 text_SubDetectors.AddText("CMS Simulation")
0351 text_SubDetectors.SetTextAlign(11)
0352 text_SubDetectors.Draw()
0353
0354
0355
0356 keep_alive = []
0357 if plots[plot].iDrawEta:
0358 keep_alive.extend(drawEtaValues())
0359
0360
0361 can_SubDetectors.Update()
0362 if not checkFile_(theDirname):
0363 os.mkdir(theDirname)
0364 can_SubDetectors.SaveAs("%s/MaterialBdg_%s_%s.png" % (theDirname, compounddetectorname, plot))
0365
0366
0367
0368
0369 hist2d_X0_total.GetXaxis().SetRangeUser( 0., 7000.)
0370
0371 keep_alive = []
0372 keep_alive.extend(drawHalfEtaValues())
0373
0374 can_SubDetectors.Update()
0375 can_SubDetectors.Modified()
0376 can_SubDetectors.SaveAs("%s/MaterialBdg_%s_%s_Zpluszoom.png" % (theDirname, compounddetectorname, plot))
0377
0378 def createPlotsReco_(reco_file, label, debug=False):
0379 """Cumulative material budget from reconstruction.
0380
0381
0382 Internal function that will produce a cumulative profile of the
0383 material budget in the reconstruction starting from the single
0384 detectors that compose the tracker. It will iterate over all
0385 existing detectors contained in the sDETS dictionary. The
0386 function will automatically stop everytime it encounters a
0387 non-existent detector, until no more detectors are left to
0388 try. For this reason the keys in the sDETS dictionary can be as
0389 inclusive as possible.
0390
0391 """
0392
0393 cumulative_matbdg = None
0394 sPREF = ["Original_RadLen_vs_Eta_", "RadLen_vs_Eta_"]
0395
0396 c = TCanvas("c", "c", 1024, 1024);
0397 diffs = []
0398 if not checkFile_(reco_file):
0399 print("Error: missing file %s" % reco_file)
0400 raise RuntimeError
0401 file = TFile(reco_file)
0402 prefix = "/DQMData/Run 1/Tracking/Run summary/RecoMaterial/"
0403 for s in sPREF:
0404 hs = THStack("hs","");
0405 histos = []
0406 for det, color in sDETS.items():
0407 layer_number = 0
0408 while True:
0409 layer_number += 1
0410 name = "%s%s%s%d" % (prefix, s, det, layer_number)
0411 prof = file.Get(name)
0412
0413
0414
0415
0416 if not prof:
0417 if debug:
0418 print("Missing profile %s" % name)
0419 break
0420 else:
0421 histos.append(prof.ProjectionX("_px", "hist"));
0422 diffs.append(histos[-1]);
0423 histos[-1].SetFillColor(color + layer_number);
0424 histos[-1].SetLineColor(color + layer_number + 1);
0425
0426 name = "CumulativeRecoMatBdg_%s" % s
0427 if s == "RadLen_vs_Eta_":
0428 cumulative_matbdg = TH1D(name, name,
0429 histos[0].GetNbinsX(),
0430 histos[0].GetXaxis().GetXmin(),
0431 histos[0].GetXaxis().GetXmax())
0432 cumulative_matbdg.SetDirectory(0)
0433 for h in histos:
0434 hs.Add(h)
0435 if cumulative_matbdg:
0436 cumulative_matbdg.Add(h, 1.)
0437 hs.Draw()
0438 hs.GetYaxis().SetTitle("RadLen")
0439 c.Update()
0440 c.Modified()
0441 c.SaveAs("%sstacked_%s.png" % (s, label))
0442 hs = THStack("diff","")
0443 for d in range(0,len(diffs)/2):
0444 diffs[d+len(diffs)/2].Add(diffs[d], -1.)
0445 hs.Add(diffs[d+len(diffs)/2]);
0446 hs.Draw()
0447 hs.GetYaxis().SetTitle("RadLen")
0448 c.Update()
0449 c.Modified()
0450 c.SaveAs("RadLen_difference_%s.png" % label)
0451 return cumulative_matbdg
0452
0453 def materialBudget_Simul_vs_Reco(reco_file, label, debug=False):
0454 """Plot reco vs simulation material budget.
0455
0456 Function are produces a direct comparison of the material
0457 budget as extracted from the reconstruction geometry and
0458 inferred from the simulation one.
0459
0460 """
0461
0462 setTDRStyle()
0463
0464
0465 cumulative_matbdg_sim = createPlots_("x_vs_eta")
0466 cumulative_matbdg_rec = createPlotsReco_(reco_file, label, debug=False)
0467
0468 cc = TCanvas("cc", "cc", 1024, 1024)
0469 cumulative_matbdg_sim.SetMinimum(0.)
0470 cumulative_matbdg_sim.SetMaximum(3.5)
0471 cumulative_matbdg_sim.GetXaxis().SetRangeUser(-3.0, 3.0)
0472 cumulative_matbdg_sim.SetLineColor(kOrange)
0473 cumulative_matbdg_rec.SetMinimum(0.)
0474 cumulative_matbdg_rec.SetMaximum(3.)
0475 cumulative_matbdg_rec.SetLineColor(kAzure+1)
0476 l = TLegend(0.18, 0.8, 0.95, 0.92)
0477 l.AddEntry(cumulative_matbdg_sim, "Sim Material", "f")
0478 l.AddEntry(cumulative_matbdg_rec, "Reco Material", "f")
0479 cumulative_matbdg_sim.Draw("HIST")
0480 cumulative_matbdg_rec.Draw("HIST SAME")
0481 l.Draw()
0482 filename = "MaterialBdg_Reco_vs_Simul_%s.png" % label
0483 cc.SaveAs(filename)
0484
0485 def createCompoundPlots(detector, plot):
0486 """Produce the requested plot for the specified detector.
0487
0488 Function that will plot the requested @plot for the specified
0489 @detector. The specified detector could either be a real
0490 detector or a compound one. The list of available plots are the
0491 keys of plots dictionary (imported from plot_utils.
0492
0493 """
0494
0495 theDirname = 'Images'
0496 if not checkFile_(theDirname):
0497 os.mkdir(theDirname)
0498
0499 goodToGo, theDetectorFilename = paramsGood_(detector, plot)
0500 if not goodToGo:
0501 return
0502
0503 theDetectorFile = TFile(theDetectorFilename)
0504
0505
0506
0507 prof_X0_elements = OrderedDict()
0508 hist_X0_elements = OrderedDict()
0509 for label, [num, color, leg] in hist_label_to_num.items():
0510
0511 prof_X0_elements[label] = theDetectorFile.Get("%d" % (num + plots[plot].plotNumber))
0512 hist_X0_elements[label] = prof_X0_elements[label].ProjectionX()
0513 hist_X0_elements[label].SetFillColor(color)
0514 hist_X0_elements[label].SetLineColor(kBlack)
0515
0516 files = []
0517 if detector in COMPOUNDS.keys():
0518 for subDetector in COMPOUNDS[detector][1:]:
0519 subDetectorFilename = "matbdg_%s.root" % subDetector
0520
0521
0522 if not checkFile_(subDetectorFilename):
0523 continue
0524
0525 subDetectorFile = TFile(subDetectorFilename)
0526 files.append(subDetectorFile)
0527 print("*** Open file... %s" % subDetectorFilename)
0528
0529
0530 for label, [num, color, leg] in hist_label_to_num.items():
0531 prof_X0_elements[label] = subDetectorFile.Get("%d" % (num + plots[plot].plotNumber))
0532 hist_X0_elements[label].Add(prof_X0_elements[label].ProjectionX("B_%s" % prof_X0_elements[label].GetName())
0533 , +1.000)
0534
0535
0536 stackTitle = "Material Budget %s;%s;%s" % (detector,
0537 plots[plot].abscissa,
0538 plots[plot].ordinate)
0539 stack_X0 = THStack("stack_X0", stackTitle);
0540 for label, [num, color, leg] in hist_label_to_num.items():
0541 stack_X0.Add(hist_X0_elements[label])
0542
0543
0544 canname = "MBCan_1D_%s_%s" % (detector, plot)
0545 can = TCanvas(canname, canname, 800, 800)
0546 can.Range(0,0,25,25)
0547 can.SetFillColor(kWhite)
0548 gStyle.SetOptStat(0)
0549 gStyle.SetOptTitle(1);
0550
0551
0552 stack_X0.Draw("HIST");
0553 stack_X0.GetYaxis().SetTitleOffset(1.15);
0554
0555
0556 theLegend = TLegend(0.40, 0.65, 0.60, 0.89)
0557 if plot == "x_vs_phi" or plot == "l_vs_phi": theLegend = TLegend(0.65, 0.30, 0.89, 0.70)
0558 if plot == "x_vs_R" or plot == "l_vs_R": theLegend = TLegend(0.75, 0.60, 0.95, 0.90)
0559
0560 for label, [num, color, leg] in hist_label_to_num.items():
0561 theLegend.AddEntry(hist_X0_elements[label], leg, "f")
0562 theLegend.Draw();
0563
0564
0565 can.Update();
0566 can.SaveAs( "%s/%s_%s.pdf" % (theDirname, detector, plot))
0567 can.SaveAs( "%s/%s_%s.png" % (theDirname, detector, plot))
0568
0569
0570
0571 if plot == "x_vs_eta" or plot == "l_vs_eta":
0572 canname = "MBCan_1D_%s_%s_total" % (detector, plot)
0573 can2 = TCanvas(canname, canname, 800, 800)
0574 can2.Range(0,0,25,25)
0575 can2.SetFillColor(kWhite)
0576 gStyle.SetOptStat(0)
0577 gStyle.SetOptTitle(0);
0578
0579 stack_X0.GetStack().Last().GetXaxis().SetRangeUser( 0., 3.)
0580 stack_X0.GetStack().Last().Draw();
0581 stack_X0.GetYaxis().SetTitleOffset(1.15);
0582 can2.Update()
0583 can2.Modified()
0584 can2.SaveAs("%s/%s_%s_total_Zplus.pdf" % (theDirname, detector, plot))
0585 can2.SaveAs("%s/%s_%s_total_Zplus.png" % (theDirname, detector, plot))
0586 stack_X0.GetStack().Last().GetXaxis().SetRangeUser( -3., 0.)
0587 stack_X0.GetStack().Last().Draw();
0588 stack_X0.GetYaxis().SetTitleOffset(1.15);
0589 can2.Update()
0590 can2.Modified()
0591 can2.SaveAs("%s/%s_%s_total_Zminus.pdf" % (theDirname, detector, plot))
0592 can2.SaveAs("%s/%s_%s_total_Zminus.png" % (theDirname, detector, plot))
0593
0594
0595 def create2DPlots(detector, plot, plotnum, plotmat, dosingledetector = True):
0596 """Produce the requested plot for the specified detector.
0597
0598 Function that will plot the requested 2D-@plot for the
0599 specified @detector. The specified detector could either be a
0600 real detector or a compound one. The list of available plots
0601 are the keys of plots dictionary imported from plot_utils.
0602 """
0603
0604
0605
0606 plotmat = plotmat.replace(" ", "_")
0607
0608
0609 if plotmat != "":
0610 theDirname = ('Images/%s' % plotmat).replace(" ", "")
0611 else:
0612 theDirname = 'Images'
0613
0614 if not checkFile_(theDirname):
0615 os.mkdir(theDirname)
0616 if not os.path.isdir(('Images/%s/ZPlusZoom' % plotmat).replace(" ", "")):
0617 os.mkdir( ('Images/%s/ZPlusZoom' % plotmat).replace(" ", "") )
0618 if not os.path.isdir(('Images/%s/ZMinusZoom' % plotmat).replace(" ", "")):
0619 os.mkdir( ('Images/%s/ZMinusZoom' % plotmat).replace(" ", "") )
0620
0621 goodToGo, theDetectorFilename = paramsGood_(detector, plot)
0622 if not goodToGo:
0623 return
0624
0625 theDetectorFile = TFile(theDetectorFilename)
0626
0627 prof2d_X0_det_total = TProfile2D()
0628 prof2d_X0_det_total.Reset()
0629
0630
0631
0632 prof2d_X0_det_total = theDetectorFile.Get('%s' % plotnum)
0633 print("==================================================================")
0634 print(plotnum)
0635
0636
0637 prof2d_X0_det_total.__class__ = TProfile2D
0638 hist_X0_total = prof2d_X0_det_total.ProjectionXY()
0639
0640
0641 files = []
0642 if detector in COMPOUNDS.keys() and not dosingledetector:
0643
0644
0645
0646
0647
0648 for subDetector in COMPOUNDS[detector]:
0649
0650 subDetectorFilename = "matbdg_%s.root" % subDetector
0651
0652
0653 if not checkFile_(subDetectorFilename):
0654 print("Error, missing file %s" % subDetectorFilename)
0655 continue
0656
0657 subDetectorFile = TFile(subDetectorFilename)
0658 files.append(subDetectorFile)
0659 print("*** Open file... %s" % subDetectorFilename)
0660
0661
0662 prof2d_X0_det_total = subDetectorFile.Get('%s' % plots[plot].plotNumber)
0663 prof2d_X0_det_total.__class__ = TProfile2D
0664
0665
0666 hist_X0_total.Add(prof2d_X0_det_total.ProjectionXY("B_%s" % prof2d_X0_det_total.GetName()), +1.000 )
0667
0668
0669
0670 gStyle.SetStripDecimals(False)
0671
0672
0673
0674 minX = 1.03*prof2d_X0_det_total.GetXaxis().GetXmin()
0675 maxX = 1.03*prof2d_X0_det_total.GetXaxis().GetXmax()
0676 minY = 1.03*prof2d_X0_det_total.GetYaxis().GetXmin()
0677 maxY = 1.03*prof2d_X0_det_total.GetYaxis().GetXmax()
0678
0679 frame = TH2F("frame", "", 10, minX, maxX, 10, minY, maxY);
0680 frame.SetMinimum(0.1)
0681 frame.SetMaximum(10.)
0682 frame.GetXaxis().SetTickLength(frame.GetXaxis().GetTickLength()*0.50)
0683 frame.GetYaxis().SetTickLength(frame.GetXaxis().GetTickLength()/4.)
0684
0685
0686 if plots[plot].iRebin:
0687 prof2d_X0_det_total.Rebin2D()
0688
0689
0690 hist2dTitle = ('%s %s;%s;%s;%s' % (plots[plot].quotaName,
0691 detector,
0692 plots[plot].abscissa,
0693 plots[plot].ordinate,
0694 plots[plot].quotaName))
0695
0696 if dosingledetector:
0697 hist2d_X0_total = prof2d_X0_det_total
0698 else :
0699 hist2d_X0_total = hist_X0_total
0700 hist2d_X0_total.SetTitle(hist2dTitle)
0701 frame.SetTitle(hist2dTitle)
0702 frame.SetTitleOffset(0.5,"Y")
0703
0704
0705
0706 if plots[plot].histoMin != -1.:
0707 hist2d_X0_total.SetMinimum(plots[plot].histoMin)
0708 if plots[plot].histoMax != -1.:
0709 hist2d_X0_total.SetMaximum(plots[plot].histoMax)
0710
0711
0712 can2name = "MBCan_2D_%s_%s_%s" % (detector, plot, plotmat)
0713 can2 = TCanvas(can2name, can2name, 2480+248, 580+58+58)
0714 can2.SetTopMargin(0.1)
0715 can2.SetBottomMargin(0.1)
0716 can2.SetLeftMargin(0.04)
0717 can2.SetRightMargin(0.06)
0718 can2.SetFillColor(kWhite)
0719 gStyle.SetOptStat(0)
0720 gStyle.SetTitleFillColor(0)
0721 gStyle.SetTitleBorderSize(0)
0722
0723
0724
0725
0726 acustompalette()
0727 ex1 = TExec("ex1","acustompalette();");
0728 ex1.Draw();
0729
0730
0731
0732
0733
0734
0735 can2.SetLogz(plots[plot].zLog)
0736
0737
0738
0739
0740 hist2d_X0_total.Draw("COLZ")
0741
0742
0743 can2.Update()
0744
0745
0746 palette = hist2d_X0_total.GetListOfFunctions().FindObject("palette")
0747 if palette:
0748 palette.__class__ = TPaletteAxis
0749 palette.SetX1NDC(0.945)
0750 palette.SetX2NDC(0.96)
0751 palette.SetY1NDC(0.1)
0752 palette.SetY2NDC(0.9)
0753 palette.GetAxis().SetTickSize(.01)
0754 palette.GetAxis().SetTitle("")
0755 if plots[plot].zLog:
0756 palette.GetAxis().SetLabelOffset(-0.01)
0757 paletteTitle = TLatex(1.12*maxX, maxY, plots[plot].quotaName)
0758 paletteTitle.SetTextAngle(90.)
0759 paletteTitle.SetTextSize(0.05)
0760 paletteTitle.SetTextAlign(31)
0761 paletteTitle.Draw()
0762 hist2d_X0_total.GetYaxis().SetTickLength(hist2d_X0_total.GetXaxis().GetTickLength()/4.)
0763 hist2d_X0_total.GetYaxis().SetTickLength(hist2d_X0_total.GetXaxis().GetTickLength()/4.)
0764 hist2d_X0_total.SetTitleOffset(0.5,"Y")
0765 hist2d_X0_total.GetYaxis().SetTitleOffset(0.45);
0766
0767
0768
0769
0770
0771 keep_alive = []
0772 if plots[plot].iDrawEta:
0773 keep_alive.extend(drawEtaValues())
0774
0775 can2.Modified()
0776 hist2d_X0_total.SetContour(255)
0777
0778
0779 can2.Update()
0780 can2.Modified()
0781
0782 can2.SaveAs( "%s/%s_%s%s.pdf" % (theDirname, detector, plot, plotmat))
0783 can2.SaveAs( "%s/%s_%s%s.png" % (theDirname, detector, plot, plotmat))
0784
0785
0786
0787 if plot == "x_vs_z_vs_Rsum" or plot == "l_vs_z_vs_Rsum" or plot == "x_vs_z_vs_Rsumcos" or plot == "l_vs_z_vs_Rsumcos" or plot == "x_vs_z_vs_Rloc" or plot == "l_vs_z_vs_Rloc" or plot == "x_vs_z_vs_Rloccos" or plot == "l_vs_z_vs_Rloccos":
0788
0789
0790 if dosingledetector: hist2d_X0_total.GetXaxis().SetRangeUser( 3100., 5400.)
0791 else : hist2d_X0_total.GetXaxis().SetRangeUser( 0., 7000.)
0792
0793 keep_alive = []
0794
0795 can2.Update()
0796 can2.Modified()
0797 can2.SaveAs( "%s/%s/%s_%s%s_ZplusZoom.pdf" % (theDirname, "ZPlusZoom", detector, plot, plotmat))
0798 can2.SaveAs( "%s/%s/%s_%s%s_ZplusZoom.png" % (theDirname, "ZPlusZoom", detector, plot, plotmat))
0799
0800
0801 if dosingledetector: hist2d_X0_total.GetXaxis().SetRangeUser( -5400., -3100.)
0802 else : hist2d_X0_total.GetXaxis().SetRangeUser( 0., -7000.)
0803
0804 keep_alive = []
0805
0806 can2.Update()
0807 can2.Modified()
0808 can2.SaveAs( "%s/%s/%s_%s%s_ZminusZoom.pdf" % (theDirname, "ZMinusZoom", detector, plot, plotmat))
0809 can2.SaveAs( "%s/%s/%s_%s%s_ZminusZoom.png" % (theDirname, "ZMinusZoom", detector, plot, plotmat))
0810
0811
0812 gStyle.SetStripDecimals(True)
0813
0814 def createRatioPlots(detector, plot):
0815 """Create ratio plots.
0816
0817 Function that will make the ratio between the radiation length
0818 and interaction length, for the specified detector. The
0819 specified detector could either be a real detector or a
0820 compound one.
0821
0822 """
0823
0824 goodToGo, theDetectorFilename = paramsGood_(detector, plot)
0825 if not goodToGo:
0826 return
0827
0828 theDirname = 'Images'
0829 if not os.path.exists(theDirname):
0830 os.mkdir(theDirname)
0831
0832 theDetectorFile = TFile(theDetectorFilename)
0833
0834 prof_x0_det_total = theDetectorFile.Get('%d' % plots[plot].plotNumber)
0835 prof_l0_det_total = theDetectorFile.Get('%d' % (10000+plots[plot].plotNumber))
0836
0837
0838 hist_x0_total = prof_x0_det_total.ProjectionX()
0839 hist_l0_total = prof_l0_det_total.ProjectionX()
0840
0841 if detector in COMPOUNDS.keys():
0842 for subDetector in COMPOUNDS[detector][1:]:
0843
0844
0845 subDetectorFilename = "matbdg_%s.root" % subDetector
0846
0847
0848 if not checkFile_(subDetectorFilename):
0849 print("Error, missing file %s" % subDetectorFilename)
0850 continue
0851
0852 subDetectorFile = TFile(subDetectorFilename)
0853
0854
0855 prof_x0_det_total = subDetectorFile.Get('%d' % plots[plot].plotNumber)
0856 prof_l0_det_total = subDetectorFile.Get('%d' % (10000+plots[plot].plotNumber))
0857
0858 hist_x0_total.Add(prof_x0_det_total.ProjectionX("B_%s" % prof_x0_det_total.GetName()), +1.000 )
0859 hist_l0_total.Add(prof_l0_det_total.ProjectionX("B_%s" % prof_l0_det_total.GetName()), +1.000 )
0860
0861 hist_x0_over_l0_total = hist_x0_total
0862 hist_x0_over_l0_total.Divide(hist_l0_total)
0863 histTitle = "Material Budget %s;%s;%s" % (detector,
0864 plots[plot].abscissa,
0865 plots[plot].ordinate)
0866 hist_x0_over_l0_total.SetTitle(histTitle)
0867
0868 hist_x0_over_l0_total.SetMarkerStyle(1)
0869 hist_x0_over_l0_total.SetMarkerSize(3)
0870 hist_x0_over_l0_total.SetMarkerColor(kBlue)
0871
0872
0873 canRname = "MBRatio_%s_%s" % (detector, plot)
0874 canR = TCanvas(canRname,canRname,800,800)
0875 canR.Range(0,0,25,25)
0876 canR.SetFillColor(kWhite)
0877 gStyle.SetOptStat(0)
0878
0879
0880 hist_x0_over_l0_total.Draw("E1")
0881
0882
0883 canR.Update()
0884 canR.SaveAs("%s/%s_%s.pdf" % (theDirname, detector, plot))
0885 canR.SaveAs("%s/%s_%s.png" % (theDirname, detector, plot))
0886
0887
0888 def GetSiliconZValuesFromXML():
0889 """Someone can use the xml files of the geometry. Just do
0890 cd $CMSSW_BASE/src/Geometry/HGCalCommonData/test
0891 cmsRun testHGCalParameters_cfg.py
0892 The scintillator position in the HEB is very close to the silicon but they have a
0893 1.95 mm difference. So, we take the silicon position.
0894
0895 This returns a np.array of the z position of the silicon values in mm.
0896 """
0897
0898
0899
0900
0901 layersmaxZ = ( 319.05, 319.815, 320.725, 322.255, 323.165, 324.695, 325.605, 327.135, 328.045, 329.575, 330.485, 332.015, 332.925, 334.455, 335.365, 336.895, 337.805, 339.335, 340.245, 341.775, 342.685, 344.215, 345.125, 346.655, 347.565, 349.095, 350.005, 351.535, 352.445, 357.715, 362.615, 367.515, 372.415, 377.315, 382.215, 387.115, 392.015, 397.105, 402.195, 407.285, 412.375, 420.765, 429.155, 437.545, 445.935, 454.325, 462.715, 471.105, 479.495, 487.885, 496.275, 504.665, 513.055)
0902
0903 layersmaxZ = np.asarray(layersmaxZ)
0904 layersmaxZ = 10. * layersmaxZ
0905
0906 print("Total number of layers from XML " , layersmaxZ.size - 1)
0907
0908 return layersmaxZ
0909
0910
0911 if __name__ == '__main__':
0912 parser = argparse.ArgumentParser(description='Generic Material Plotter',
0913 formatter_class=argparse.ArgumentDefaultsHelpFormatter)
0914 parser.add_argument('-r', '--reco',
0915 help='Input Reconstruction Material file, DQM format')
0916 parser.add_argument('-l', '--label',
0917 help='Label to use in naming the plots')
0918 parser.add_argument('-c', '--compare',
0919 help='Compare simulation and reco materials',
0920 action='store_true',
0921 default=False)
0922 parser.add_argument('-m', '--multi',
0923 help='Combining multiple detectors',
0924 action='store_true',
0925 default=False)
0926 parser.add_argument('-s', '--single',
0927 help='Material budget for single detector from simulation',
0928 action='store_true',
0929 default=False)
0930 parser.add_argument('-d', '--detector',
0931 help='Detector for which you want to compute the material budget',
0932 type=str,)
0933 args = parser.parse_args()
0934
0935 if args.compare and args.single:
0936 print("Error, too many actions required")
0937 raise RuntimeError
0938
0939 if args.compare:
0940 if args.reco is None:
0941 print("Error, missing inpur reco file")
0942 raise RuntimeError
0943 if args.label is None:
0944 print("Error, missing label")
0945 raise RuntimeError
0946 materialBudget_Simul_vs_Reco(args.reco, args.label, debug=False)
0947
0948 if args.multi:
0949
0950 cumulative_matbdg_sim = createPlots_("x_vs_eta",args.detector)
0951 cumulative_matbdg_sim = createPlots_("l_vs_eta",args.detector)
0952 cumulative_matbdg_sim = createPlots2D_("x_vs_z_vs_Rsum",args.detector)
0953 cumulative_matbdg_sim = createPlots2D_("l_vs_z_vs_Rsum",args.detector)
0954
0955 required_2DplotsForalldetectors = ["x_vs_z_vs_Rsum", "l_vs_z_vs_Rsum"]
0956
0957 for p in required_2DplotsForalldetectors:
0958
0959 create2DPlots(args.detector, p, plots[p].plotNumber, "", not args.multi)
0960
0961
0962
0963
0964
0965 if args.single:
0966 if args.detector is None:
0967 print("Error, missing detector")
0968 raise RuntimeError
0969 required_2Dplots = ["x_vs_eta_vs_phi", "l_vs_eta_vs_phi", "x_vs_z_vs_Rsum", "l_vs_z_vs_Rsum", "x_vs_z_vs_Rsumcos", "l_vs_z_vs_Rsumcos", "x_vs_z_vs_Rloc", "l_vs_z_vs_Rloc", "x_vs_z_vs_Rloccos", "l_vs_z_vs_Rloccos"]
0970
0971
0972
0973
0974 required_plots = ["x_vs_eta", "x_vs_phi", "x_vs_R", "l_vs_eta", "l_vs_phi", "l_vs_R"]
0975
0976 required_ratio_plots = ["x_over_l_vs_eta", "x_over_l_vs_phi"]
0977
0978
0979
0980
0981 print("---+++ Results: Plots for individual material in all HGCal")
0982 for label, [num, color, leg] in hist_label_to_num.items():
0983 for p in ["x_vs_z_vs_Rsum", "l_vs_z_vs_Rsum", "x_vs_z_vs_Rsumcos", "l_vs_z_vs_Rsumcos", "x_vs_z_vs_Rloc", "l_vs_z_vs_Rloc"]:
0984 TwikiPrintout(p, leg, "all")
0985
0986 print("---+++ Results: Plots for individual material in Z+ Endcap")
0987 for label, [num, color, leg] in hist_label_to_num.items():
0988 for p in ["x_vs_z_vs_Rsum", "l_vs_z_vs_Rsum", "x_vs_z_vs_Rsumcos", "l_vs_z_vs_Rsumcos", "x_vs_z_vs_Rloc", "l_vs_z_vs_Rloc"]:
0989 TwikiPrintout(p, leg, "zplus")
0990
0991 print("---+++ Results: Plots for individual material in Z- Endcap")
0992 for label, [num, color, leg] in hist_label_to_num.items():
0993 for p in ["x_vs_z_vs_Rsum", "l_vs_z_vs_Rsum", "x_vs_z_vs_Rsumcos", "l_vs_z_vs_Rsumcos", "x_vs_z_vs_Rloc", "l_vs_z_vs_Rloc"]:
0994 TwikiPrintout(p, leg, "zminus")
0995
0996
0997
0998 print(GetSiliconZValuesFromXML())
0999
1000 for p in required_2Dplots:
1001
1002 create2DPlots(args.detector, p, plots[p].plotNumber, "")
1003
1004 for label, [num, color, leg] in hist_label_to_num.items():
1005
1006 create2DPlots(args.detector, p, num + plots[p].plotNumber, leg)
1007
1008
1009 for p in required_plots:
1010 createCompoundPlots(args.detector, p)
1011 for p in required_ratio_plots:
1012 createRatioPlots(args.detector, p)