File indexing completed on 2024-11-28 23:11:14
0001
0002
0003
0004
0005 import sys
0006 oldargv = sys.argv[:]
0007 sys.argv = [ '-b-' ]
0008 from ROOT import TCanvas, TPad, TGaxis, TLegend, TPaveText, THStack, TFile, TLatex
0009 from ROOT import TProfile, TProfile2D, TH1D, TH2F, TPaletteAxis, TStyle, TColor
0010 from ROOT import kBlack, kWhite, kOrange, kAzure, kBlue, kRed, kGreen
0011 from ROOT import kGreyScale, kTemperatureMap
0012 from ROOT import kTRUE, kFALSE
0013 from ROOT import gROOT, gStyle, gPad
0014 gROOT.SetBatch(True)
0015 sys.argv = oldargv
0016
0017 from Validation.Geometry.plot_utils import setTDRStyle, Plot_params, plots, COMPOUNDS, DETECTORS, sDETS, hist_label_to_num, drawEtaValues
0018 from collections import namedtuple, OrderedDict
0019 import sys, os, copy
0020 import argparse
0021
0022 def paramsGood_(detector, plot, geometryOld = '', geometryNew = ''):
0023 """Check the validity of the arguments.
0024
0025 Common function to check the validity of the parameters passed
0026 in. It returns a tuple composed by a bool and a string. The
0027 bool indicates if all checks are ok, the string the appropriate
0028 ROOT filename to open (empty string in case any check failed)
0029 If geometry comparison is being made, a list of strings is
0030 returned instead.
0031
0032 """
0033
0034 theFiles = []
0035
0036 if plot not in plots.keys():
0037 print("Error, unknown plot %s" % plot)
0038 return (False, '')
0039
0040 if detector not in DETECTORS and detector not in COMPOUNDS.keys():
0041 print('Error, unknown detector: %s' % detector)
0042 return (False, '')
0043
0044 if detector not in DETECTORS:
0045 detector = COMPOUNDS[detector][0]
0046
0047 if geometryNew:
0048 oldgeoFilename = 'matbdg_%s_%s.root' % (detector,geometryOld)
0049 theFiles.append(oldgeoFilename)
0050 newgeoFilename = 'matbdg_%s_%s.root' % (detector,geometryNew)
0051 theFiles.append(newgeoFilename)
0052 else:
0053 theFiles.append('matbdg_%s_%s.root' % (detector,geometryOld))
0054
0055 for thisFile in theFiles:
0056 if not checkFile_(thisFile):
0057 print("Error, missing file %s" % thisFile)
0058 raise RuntimeError
0059
0060 if len(theFiles) > 1:
0061 return (True, theFiles)
0062 else:
0063 return (True, theFiles[0])
0064
0065 def checkFile_(filename):
0066 return os.path.exists(filename)
0067
0068 def setColorIfExists_(histos, h, color):
0069 if h in histos.keys():
0070 histos[h].SetFillColor(color)
0071
0072 def assignOrAddIfExists_(h1, h2):
0073 """Assign the projection of h2 to h1.
0074
0075 Function to assign the h2 to h1 in the case in
0076 which h1 is None, otherwise add h2 to the already
0077 valid h1 object
0078
0079 """
0080
0081 if not h1:
0082 h1 = h2
0083 else:
0084 h1.Add(h2, +1.000)
0085 return h1
0086
0087 def get1DHisto_(detector,plotNumber,geometry):
0088 """
0089 This function opens the appropiate ROOT file,
0090 extracts the TProfile and turns it into a Histogram,
0091 if it is a compound detector, this function
0092 takes care of the subdetectors' addition unless the
0093 detector's ROOT file is present in which case no addition
0094 is performed and the detector ROOT file is used.
0095 """
0096 histo = None
0097 rootFile = TFile()
0098
0099 detectorFilename = 'matbdg_%s_%s.root'%(detector,geometry)
0100 if detector not in COMPOUNDS.keys() or checkFile_(detectorFilename):
0101 if not checkFile_(detectorFilename):
0102 print('Warning: %s not found' % detectorFilename)
0103 return 0
0104 print('Reading from: %s File' % detectorFilename)
0105 rootFile = TFile.Open(detectorFilename,'READ')
0106 prof = rootFile.Get("%d" % plotNumber)
0107 if not prof: return 0
0108
0109 prof.SetName('%u_%s_%s' %(plotNumber,detector,geometry))
0110 histo = prof.ProjectionX()
0111 else:
0112 theFiles = []
0113 histos = OrderedDict()
0114 for subDetector in COMPOUNDS[detector]:
0115 subDetectorFilename = 'matbdg_%s_%s.root' % (subDetector,geometry)
0116 if not checkFile_(subDetectorFilename):
0117 print('Warning: %s not found'%subDetectorFilename)
0118 continue
0119 print('Reading from: %s File' % subDetectorFilename)
0120 subDetectorFile = TFile.Open(subDetectorFilename,'READ')
0121 theFiles.append(subDetectorFile)
0122 prof = subDetectorFile.Get('%d'%(plotNumber))
0123 if not prof: return 0
0124 prof.__class__ = TProfile
0125 histo = assignOrAddIfExists_(histo,prof.ProjectionX())
0126
0127 return copy.deepcopy(histo)
0128
0129 def get2DHisto_(detector,plotNumber,geometry):
0130 """
0131 This function opens the appropiate ROOT file,
0132 extracts the TProfile2D and turns it into a Histogram,
0133 if it is a compound detector, this function
0134 takes care of the subdetectors' addition.
0135
0136 Note that it takes plotNumber as opposed to plot
0137 """
0138 histo = None
0139 rootFile = TFile()
0140
0141 detectorFilename = 'matbdg_%s_%s.root'%(detector,geometry)
0142 if detector not in COMPOUNDS.keys() or checkFile_(detectorFilename):
0143 if not checkFile_(detectorFilename):
0144 print('Warning: %s not found' % detectorFilename)
0145 return 0
0146 rootFile = TFile.Open(detectorFilename,'READ')
0147 prof = rootFile.Get("%d" % plotNumber)
0148 if not prof: return 0
0149
0150 prof.SetName('%u_%s_%s' %(plotNumber,detector,geometry))
0151 prof.__class__ = TProfile2D
0152 histo = prof.ProjectionXY()
0153 else:
0154 histos = OrderedDict()
0155 theFiles = []
0156 for subDetector in COMPOUNDS[detector]:
0157 subDetectorFilename = 'matbdg_%s_%s.root' % (subDetector,geometry)
0158 if not checkFile_(subDetectorFilename):
0159 print('Warning: %s not found'%subDetectorFilename)
0160 continue
0161 subDetectorFile = TFile.Open(subDetectorFilename,'READ')
0162 theFiles.append(subDetectorFile)
0163 print('*** Open file... %s' % subDetectorFilename)
0164 prof = subDetectorFile.Get('%d'%plotNumber)
0165 if not prof: return 0
0166 prof.__class__ = TProfile2D
0167 if not histo:
0168 histo = prof.ProjectionXY('B_%s' % prof.GetName())
0169 else:
0170 histo.Add(prof.ProjectionXY('B_%s' % prof.GetName()))
0171
0172 return copy.deepcopy(histo)
0173
0174 def createCompoundPlotsGeometryComparison(detector, plot, geometryOld,
0175 geometryNew):
0176
0177 setTDRStyle()
0178
0179 goodToGo, theFiles = paramsGood_(detector,plot,
0180 geometryOld,geometryNew)
0181
0182 if not goodToGo:
0183 return
0184
0185 oldHistos = OrderedDict()
0186 newHistos = OrderedDict()
0187 ratioHistos = OrderedDict()
0188 diffHistos = OrderedDict()
0189
0190 def setUpCanvas(canvas):
0191
0192 gStyle.SetOptStat(False)
0193
0194 mainPadTop = [
0195 TPad("mainPadTop"+str(i)+'_'+canvas.GetName(),
0196 "mainPad"+str(i),
0197 i*0.25, 0.60, (i+1)*0.25, 1.0)
0198 for i in range(4)
0199 ]
0200
0201 subPadTop = [
0202 TPad("subPadTop"+str(i)+'_'+canvas.GetName(),
0203 "subPad"+str(i),
0204 i*0.25, 0.50, (i+1)*0.25, 0.6)
0205 for i in range(4)
0206 ]
0207
0208 mainPadBottom = [
0209 TPad("mainPadBottom"+str(i)+'_'+canvas.GetName(),
0210 "subPad"+str(i),
0211 i*0.25, 0.10, (i+1)*0.25, 0.5)
0212 for i in range(4)
0213 ]
0214
0215 subPadBottom = [
0216 TPad("subPadBottom"+str(i)+'_'+canvas.GetName(),
0217 "subPad"+str(i),
0218 i*0.25, 0.00, (i+1)*0.25, 0.1)
0219 for i in range(4)
0220 ]
0221
0222 mainPad = mainPadTop + mainPadBottom
0223 subPad = subPadTop + subPadBottom
0224
0225 leftMargin = 0.12
0226 rightMargin = 0.12
0227 topMargin = 0.12
0228 bottomMargin = 0.3
0229 for i in range(8):
0230 mainPad[i].SetLeftMargin(leftMargin)
0231 mainPad[i].SetRightMargin(rightMargin)
0232 mainPad[i].SetTopMargin(topMargin)
0233 mainPad[i].SetBottomMargin(1e-3)
0234 mainPad[i].Draw()
0235 subPad[i].SetLeftMargin(leftMargin)
0236 subPad[i].SetRightMargin(rightMargin)
0237 subPad[i].SetTopMargin(1e-3)
0238 subPad[i].SetBottomMargin(bottomMargin)
0239 subPad[i].Draw()
0240
0241 return mainPad, subPad
0242
0243 canComparison = TCanvas("canComparison","canComparison",2400,1200)
0244 mainPad, subPad = setUpCanvas(canComparison)
0245
0246
0247 def setStyleHistoSubPad(histo):
0248 histo.SetTitle('')
0249 histo.SetMarkerColor(kBlack)
0250 histo.SetMarkerStyle(20)
0251 histo.SetMarkerSize(.5)
0252 histo.SetLineWidth(1)
0253
0254 histo.GetYaxis().SetTitleSize(14)
0255 histo.GetYaxis().SetTitleFont(43)
0256 histo.GetYaxis().SetLabelSize(0.17)
0257 histo.GetYaxis().SetTitleOffset(5.0)
0258 histo.GetYaxis().SetNdivisions(6,3,0)
0259
0260 histo.GetXaxis().SetTitleSize(25)
0261 histo.GetXaxis().SetTitleFont(43)
0262 histo.GetXaxis().SetTitleOffset(6.0)
0263 histo.GetXaxis().SetLabelSize(0.17)
0264
0265 return histo
0266
0267 def makeRatio(histoX,histoY):
0268
0269 histoXOverY = copy.deepcopy(histoX)
0270 histoXOverY.Divide(histoY)
0271 histoXOverY.GetYaxis().SetTitle('#frac{%s}{%s}' % (geometryNew,geometryOld))
0272
0273 return histoXOverY
0274
0275 def makeDiff(histoNew,histoOld):
0276
0277 diff = copy.deepcopy(histoNew)
0278 diff.Add(histoOld,-1.0)
0279 diff.GetYaxis().SetTitle(geometryNew
0280 + " - "
0281 + geometryOld)
0282 diff.GetYaxis().SetNdivisions(6,3,0)
0283
0284 diff.GetXaxis().SetTitleSize(25)
0285 diff.GetXaxis().SetTitleFont(43)
0286 diff.GetXaxis().SetTitleOffset(3.5)
0287 diff.GetXaxis().SetLabelSize(0.17)
0288
0289 return diff
0290
0291
0292
0293
0294 def setUpTitle(detector,label,plot):
0295 title = 'Material Budget %s [%s];%s;%s' % (detector,label,
0296 plots[plot].abscissa,
0297 plots[plot].ordinate)
0298 return title
0299
0300 def setUpLegend(gOld,gNew,label):
0301 legend = TLegend(0.4,0.7,0.7,0.85)
0302 legend.AddEntry(gOld,"%s %s [%s]"%(detector,geometryOld,label),"F")
0303 legend.AddEntry(gNew,"%s %s [%s]"%(detector,geometryNew,label),"P")
0304 legend.SetTextFont(42)
0305 legend.SetTextSize(0.03)
0306 return legend
0307
0308 def setRanges(h):
0309 legendSpace = 1. + 0.3
0310 minX = h.GetXaxis().GetXmin()
0311 maxX = h.GetXaxis().GetXmax()
0312 minY = h.GetYaxis().GetXmin()
0313 maxY = h.GetBinContent(h.GetMaximumBin()) * legendSpace
0314 h.GetYaxis().SetRangeUser(minY, maxY)
0315 h.GetXaxis().SetRangeUser(minX, maxX)
0316
0317
0318
0319
0320 counter = 0
0321 legends = OrderedDict()
0322 for label, [num, color, leg] in hist_label_to_num.items():
0323
0324 mainPad[counter].cd()
0325 oldHistos[label] = get1DHisto_(detector,
0326 num+plots[plot].plotNumber
0327 ,geometryOld)
0328 oldHistos[label].SetTitle(setUpTitle(detector,leg,plot))
0329 oldHistos[label].SetFillColor(color)
0330 oldHistos[label].SetLineColor(kBlack)
0331 oldHistos[label].SetLineWidth(1)
0332 setRanges(oldHistos[label])
0333 oldHistos[label].Draw("HIST")
0334
0335 newHistos[label] = get1DHisto_(detector,
0336 num+plots[plot].plotNumber
0337 ,geometryNew)
0338 newHistos[label].SetMarkerSize(.5)
0339 newHistos[label].SetMarkerStyle(20)
0340 newHistos[label].Draw('SAME P')
0341
0342 legends[label]= setUpLegend(oldHistos[label],newHistos[label],
0343 leg);
0344 legends[label].Draw()
0345
0346
0347 subPad[counter].cd()
0348 ratioHistos[label] = makeRatio( newHistos[label],oldHistos[label] )
0349 ratioHistos[label] = setStyleHistoSubPad(ratioHistos[label])
0350 ratioHistos[label].Draw("HIST P")
0351
0352 counter += 1
0353
0354 theDirname = "Images"
0355
0356 if not checkFile_(theDirname):
0357 os.mkdir(theDirname)
0358
0359 canComparison.SaveAs( "%s/%s_ComparisonRatio_%s_%s_vs_%s.png"
0360 % (theDirname,detector,plot,geometryOld,geometryNew) )
0361
0362
0363
0364 canDiff = TCanvas("canDiff","canDiff",2400,1200)
0365
0366 mainPadDiff, subPadDiff = setUpCanvas(canDiff)
0367
0368 counter = 0
0369 for label, [num, color, leg] in hist_label_to_num.items():
0370 mainPadDiff[counter].cd()
0371 oldHistos[label].SetTitle(setUpTitle(detector,leg,plot))
0372 oldHistos[label].Draw("HIST")
0373 newHistos[label].Draw('SAME P')
0374
0375 legends[label].Draw()
0376
0377
0378 subPadDiff[counter].cd()
0379 diffHistos[label] = makeDiff( newHistos[label],oldHistos[label] )
0380 diffHistos[label] = setStyleHistoSubPad(diffHistos[label])
0381 diffHistos[label].SetTitle('')
0382 diffHistos[label].SetFillColor(color+1)
0383 diffHistos[label].Draw("HIST")
0384 counter +=1
0385
0386 canDiff.SaveAs( "%s/%s_ComparisonDifference_%s_%s_vs_%s.png"
0387 % (theDirname,detector,plot,geometryOld,geometryNew) )
0388
0389
0390 def setUpPalette(histo2D, plot) :
0391
0392
0393
0394 minX = 1.03*histo2D.GetXaxis().GetXmin();
0395 maxX = 1.03*histo2D.GetXaxis().GetXmax();
0396 minY = 1.03*histo2D.GetYaxis().GetXmin();
0397 maxY = 1.03*histo2D.GetYaxis().GetXmax();
0398
0399 palette = histo2D.GetListOfFunctions().FindObject("palette")
0400 if palette:
0401 palette.__class__ = TPaletteAxis
0402 palette.SetX1NDC(0.945)
0403 palette.SetY1NDC(gPad.GetBottomMargin())
0404 palette.SetX2NDC(0.96)
0405 palette.SetY2NDC(1-gPad.GetTopMargin())
0406 palette.GetAxis().SetTickSize(.01)
0407 palette.GetAxis().SetTitle("")
0408 if plots[plot].zLog:
0409 palette.GetAxis().SetLabelOffset(-0.01)
0410 if histo2D.GetMaximum()/histo2D.GetMinimum() < 1e3 :
0411 palette.GetAxis().SetMoreLogLabels(True)
0412 palette.GetAxis().SetNoExponent(True)
0413
0414 paletteTitle = TLatex(1.12*maxX, maxY, plots[plot].quotaName)
0415 paletteTitle.SetTextAngle(90.)
0416 paletteTitle.SetTextSize(0.05)
0417 paletteTitle.SetTextAlign(31)
0418 paletteTitle.Draw()
0419
0420 histo2D.GetXaxis().SetTickLength(histo2D.GetXaxis().GetTickLength()/4.)
0421 histo2D.GetYaxis().SetTickLength(histo2D.GetYaxis().GetTickLength()/4.)
0422 histo2D.SetTitleOffset(0.5,'Y')
0423 histo2D.GetXaxis().SetNoExponent(True)
0424 histo2D.GetYaxis().SetNoExponent(True)
0425
0426 def create2DPlotsGeometryComparison(detector, plot,
0427 geometryOld, geometryNew):
0428
0429 setTDRStyle()
0430
0431 print('Extracting plot: %s.'%(plot))
0432 goodToGo, theFiles = paramsGood_(detector,plot,
0433 geometryOld,geometryNew)
0434
0435 if not goodToGo:
0436 return
0437
0438 gStyle.SetOptStat(False)
0439
0440 old2DHisto = get2DHisto_(detector,plots[plot].plotNumber,geometryOld)
0441 new2DHisto = get2DHisto_(detector,plots[plot].plotNumber,geometryNew)
0442
0443 if plots[plot].iRebin:
0444 old2DHisto.Rebin2D()
0445 new2DHisto.Rebin2D()
0446
0447 def setRanges(h):
0448 h.GetXaxis().SetRangeUser(plots[plot].xmin, plots[plot].xmax)
0449 h.GetYaxis().SetRangeUser(plots[plot].ymin, plots[plot].ymax)
0450 if plots[plot].histoMin != -1.:
0451 h.SetMinimum(plots[plot].histoMin)
0452 if plots[plot].histoMax != -1.:
0453 h.SetMaximum(plots[plot].histoMax)
0454
0455 ratio2DHisto = copy.deepcopy(new2DHisto)
0456 ratio2DHisto.Divide(old2DHisto)
0457
0458
0459
0460 ratio2DHisto.SetMinimum(0.2)
0461 ratio2DHisto.SetMaximum(1.8)
0462 setRanges(ratio2DHisto)
0463
0464 diff2DHisto = copy.deepcopy(new2DHisto)
0465 diff2DHisto.Add(old2DHisto,-1.0)
0466 setRanges(diff2DHisto)
0467
0468
0469 def setPadStyle():
0470 gPad.SetLeftMargin(0.05)
0471 gPad.SetRightMargin(0.08)
0472 gPad.SetTopMargin(0.10)
0473 gPad.SetBottomMargin(0.10)
0474 gPad.SetLogz(plots[plot].zLog)
0475 gPad.SetFillColor(kWhite)
0476 gPad.SetBorderMode(0)
0477
0478 can = TCanvas('can','can',
0479 2724,1336)
0480 can.Divide(1,2)
0481 can.cd(1)
0482 setPadStyle()
0483 gPad.SetLogz(plots[plot].zLog)
0484
0485 gStyle.SetOptStat(0)
0486 gStyle.SetFillColor(kWhite)
0487 gStyle.SetPalette(kTemperatureMap)
0488
0489 ratio2DHisto.SetTitle("%s, Ratio: %s/%s;%s;%s"
0490 %(plots[plot].quotaName,
0491 geometryOld, geometryNew,
0492 plots[plot].abscissa,
0493 plots[plot].ordinate))
0494 ratio2DHisto.Draw('COLZ')
0495
0496 can.Update()
0497
0498 setUpPalette(ratio2DHisto,plot)
0499
0500 etasTop = []
0501 if plots[plot].iDrawEta:
0502 etasTop.extend(drawEtaValues())
0503
0504 can.cd(2)
0505
0506 diff2DHisto.SetTitle('%s, Difference: %s - %s %s;%s;%s'
0507 %(plots[plot].quotaName,geometryNew,geometryOld,detector,
0508 plots[plot].abscissa,plots[plot].ordinate))
0509 setPadStyle()
0510 diff2DHisto.Draw("COLZ")
0511 can.Update()
0512 setUpPalette(diff2DHisto,plot)
0513
0514 etasBottom = []
0515 if plots[plot].iDrawEta:
0516 etasBottom.extend(drawEtaValues())
0517
0518 can.Modified()
0519
0520 theDirname = "Images"
0521
0522 if not checkFile_(theDirname):
0523 os.mkdir(theDirname)
0524
0525 can.SaveAs( "%s/%s_Comparison_%s_%s_vs_%s.png"
0526 % (theDirname,detector,plot,geometryOld,geometryNew) )
0527 gStyle.SetStripDecimals(True)
0528
0529 def createPlots_(plot, geometry):
0530 """Cumulative material budget from simulation.
0531
0532 Internal function that will produce a cumulative profile of the
0533 material budget inferred from the simulation starting from the
0534 single detectors that compose the tracker. It will iterate over
0535 all existing detectors contained in the DETECTORS
0536 dictionary. The function will automatically skip non-existent
0537 detectors.
0538
0539 """
0540
0541 IBs = ["InnerServices", "Phase2PixelBarrel", "TIB", "TIDF", "TIDB"]
0542 theDirname = "Figures"
0543
0544 if plot not in plots.keys():
0545 print("Error: chosen plot name not known %s" % plot)
0546 return
0547
0548 hist_X0_detectors = OrderedDict()
0549 hist_X0_IB = None
0550 hist_X0_elements = OrderedDict()
0551
0552 for subDetector,color in DETECTORS.items():
0553 h = get1DHisto_(subDetector,plots[plot].plotNumber,geometry)
0554 if not h:
0555 print('Warning: Skipping %s'%subDetector)
0556 continue
0557 hist_X0_detectors[subDetector] = h
0558
0559
0560
0561 if subDetector in IBs:
0562 hist_X0_IB = assignOrAddIfExists_(
0563 hist_X0_IB,
0564 hist_X0_detectors[subDetector]
0565 )
0566
0567
0568 for label, [num, color, leg] in hist_label_to_num.items():
0569 if label == 'SUM': continue
0570 hist_label = get1DHisto_(subDetector, num + plots[plot].plotNumber, geometry)
0571 hist_X0_elements[label] = assignOrAddIfExists_(
0572 hist_X0_elements.setdefault(label,None),
0573 hist_label,
0574 )
0575 hist_X0_elements[label].SetFillColor(color)
0576
0577
0578 cumulative_matbdg = TH1D("CumulativeSimulMatBdg",
0579 "CumulativeSimulMatBdg",
0580 hist_X0_IB.GetNbinsX(),
0581 hist_X0_IB.GetXaxis().GetXmin(),
0582 hist_X0_IB.GetXaxis().GetXmax())
0583 cumulative_matbdg.SetDirectory(0)
0584
0585
0586 for det, color in DETECTORS.items():
0587 setColorIfExists_(hist_X0_detectors, det, color)
0588
0589
0590
0591 stackTitle_SubDetectors = "Tracker Material Budget;%s;%s" % (
0592 plots[plot].abscissa,plots[plot].ordinate)
0593 stack_X0_SubDetectors = THStack("stack_X0",stackTitle_SubDetectors)
0594 for det, histo in hist_X0_detectors.items():
0595 stack_X0_SubDetectors.Add(histo)
0596 cumulative_matbdg.Add(histo, 1)
0597
0598
0599 can_SubDetectors = TCanvas("can_SubDetectors","can_SubDetectors",800,800)
0600 can_SubDetectors.Range(0,0,25,25)
0601 can_SubDetectors.SetFillColor(kWhite)
0602
0603
0604 stack_X0_SubDetectors.SetMinimum(plots[plot].ymin)
0605 stack_X0_SubDetectors.SetMaximum(plots[plot].ymax)
0606 stack_X0_SubDetectors.Draw("HIST")
0607 stack_X0_SubDetectors.GetXaxis().SetLimits(plots[plot].xmin, plots[plot].xmax)
0608
0609
0610
0611 theLegend_SubDetectors = TLegend(0.180,0.8,0.98,0.92)
0612 theLegend_SubDetectors.SetNColumns(3)
0613 theLegend_SubDetectors.SetFillColor(0)
0614 theLegend_SubDetectors.SetFillStyle(0)
0615 theLegend_SubDetectors.SetBorderSize(0)
0616
0617 for det, histo in hist_X0_detectors.items():
0618 theLegend_SubDetectors.AddEntry(histo, det, "f")
0619
0620 theLegend_SubDetectors.Draw()
0621
0622
0623 text_SubDetectors = TPaveText(0.180,0.727,0.402,0.787,"NDC")
0624 text_SubDetectors.SetFillColor(0)
0625 text_SubDetectors.SetBorderSize(0)
0626 text_SubDetectors.AddText("CMS Simulation")
0627 text_SubDetectors.SetTextAlign(11)
0628 text_SubDetectors.Draw()
0629
0630
0631 can_SubDetectors.Update()
0632 if not checkFile_(theDirname):
0633 os.mkdir(theDirname)
0634 can_SubDetectors.SaveAs("%s/Tracker_SubDetectors_%s.pdf" % (theDirname, plot))
0635 can_SubDetectors.SaveAs("%s/Tracker_SubDetectors_%s.root" % (theDirname, plot))
0636
0637
0638
0639
0640 stackTitle_Materials = "Tracker Material Budget;%s;%s" % (plots[plot].abscissa,
0641 plots[plot].ordinate)
0642 stack_X0_Materials = THStack("stack_X0",stackTitle_Materials)
0643 stack_X0_Materials.Add(hist_X0_detectors["BeamPipe"])
0644 for label, [num, color, leg] in hist_label_to_num.items():
0645 if label == 'SUM':
0646 continue
0647 stack_X0_Materials.Add(hist_X0_elements[label])
0648
0649
0650 can_Materials = TCanvas("can_Materials","can_Materials",800,800)
0651 can_Materials.Range(0,0,25,25)
0652 can_Materials.SetFillColor(kWhite)
0653
0654
0655 stack_X0_Materials.SetMinimum(plots[plot].ymin)
0656 stack_X0_Materials.SetMaximum(plots[plot].ymax)
0657 stack_X0_Materials.Draw("HIST")
0658 stack_X0_Materials.GetXaxis().SetLimits(plots[plot].xmin, plots[plot].xmax)
0659
0660
0661 theLegend_Materials = TLegend(0.180,0.8,0.95,0.92)
0662 theLegend_Materials.SetNColumns(3)
0663 theLegend_Materials.SetFillColor(0)
0664 theLegend_Materials.SetBorderSize(0)
0665
0666 theLegend_Materials.AddEntry(hist_X0_detectors["BeamPipe"], "Beam Pipe", "f")
0667 for label, [num, color, leg] in hist_label_to_num.items():
0668 if label == 'SUM':
0669 continue
0670 theLegend_Materials.AddEntry(hist_X0_elements[label], leg, "f")
0671 theLegend_Materials.Draw()
0672
0673
0674 text_Materials = TPaveText(0.180,0.727,0.402,0.787,"NDC")
0675 text_Materials.SetFillColor(0)
0676 text_Materials.SetBorderSize(0)
0677 text_Materials.AddText("CMS Simulation")
0678 text_Materials.SetTextAlign(11)
0679 text_Materials.Draw()
0680
0681
0682 can_Materials.Update()
0683 can_Materials.SaveAs("%s/Tracker_Materials_%s.pdf" % (theDirname, plot))
0684 can_Materials.SaveAs("%s/Tracker_Materials_%s.root" % (theDirname, plot))
0685
0686 return cumulative_matbdg
0687
0688 def createPlotsReco_(reco_file, label, debug=False):
0689 """Cumulative material budget from reconstruction.
0690
0691
0692 Internal function that will produce a cumulative profile of the
0693 material budget in the reconstruction starting from the single
0694 detectors that compose the tracker. It will iterate over all
0695 existing detectors contained in the sDETS dictionary. The
0696 function will automatically stop everytime it encounters a
0697 non-existent detector, until no more detectors are left to
0698 try. For this reason the keys in the sDETS dictionary can be as
0699 inclusive as possible.
0700
0701 """
0702
0703 cumulative_matbdg = None
0704 sPREF = ["Original_RadLen_vs_Eta_", "RadLen_vs_Eta_"]
0705
0706 c = TCanvas("c", "c", 1024, 1024);
0707 diffs = []
0708 if not checkFile_(reco_file):
0709 print("Error: missing file %s" % reco_file)
0710 raise RuntimeError
0711 file = TFile(reco_file)
0712 prefix = "/DQMData/Run 1/Tracking/Run summary/RecoMaterial/"
0713 for s in sPREF:
0714 hs = THStack("hs","");
0715 histos = []
0716 for det, color in sDETS.items():
0717 layer_number = 0
0718 while True:
0719 layer_number += 1
0720 name = "%s%s%s%d" % (prefix, s, det, layer_number)
0721 prof = file.Get(name)
0722
0723
0724
0725
0726 if not prof:
0727 if debug:
0728 print("Missing profile %s" % name)
0729 break
0730 else:
0731 histos.append(prof.ProjectionX("_px", "hist"));
0732 diffs.append(histos[-1]);
0733 histos[-1].SetFillColor(color + layer_number);
0734 histos[-1].SetLineColor(color + layer_number + 1);
0735
0736 name = "CumulativeRecoMatBdg_%s" % s
0737 if s == "RadLen_vs_Eta_":
0738 cumulative_matbdg = TH1D(name, name,
0739 histos[0].GetNbinsX(),
0740 histos[0].GetXaxis().GetXmin(),
0741 histos[0].GetXaxis().GetXmax())
0742 cumulative_matbdg.SetDirectory(0)
0743 for h in histos:
0744 hs.Add(h)
0745 if cumulative_matbdg:
0746 cumulative_matbdg.Add(h, 1.)
0747 hs.Draw()
0748 hs.GetYaxis().SetTitle("RadLen")
0749 c.Update()
0750 c.Modified()
0751 c.SaveAs("%sstacked_%s.png" % (s, label))
0752 hs = THStack("diff","")
0753 for d in range(0,len(diffs)/2):
0754 diffs[d+len(diffs)/2].Add(diffs[d], -1.)
0755 hs.Add(diffs[d+len(diffs)/2]);
0756 hs.Draw()
0757 hs.GetYaxis().SetTitle("RadLen")
0758 c.Update()
0759 c.Modified()
0760 c.SaveAs("RadLen_difference_%s.png" % label)
0761 return cumulative_matbdg
0762
0763 def materialBudget_Simul_vs_Reco(reco_file, label, geometry, debug=False):
0764 """Plot reco vs simulation material budget.
0765
0766 Function are produces a direct comparison of the material
0767 budget as extracted from the reconstruction geometry and
0768 inferred from the simulation one.
0769
0770 """
0771
0772 setTDRStyle()
0773
0774
0775 cumulative_matbdg_sim = createPlots_("x_vs_eta", geometry)
0776 cumulative_matbdg_rec = createPlotsReco_(reco_file, label, debug=False)
0777
0778 cc = TCanvas("cc", "cc", 1024, 1024)
0779 cumulative_matbdg_sim.SetMinimum(0.)
0780 cumulative_matbdg_sim.SetMaximum(3.5)
0781 cumulative_matbdg_sim.GetXaxis().SetRangeUser(-3.0, 3.0)
0782 cumulative_matbdg_sim.SetLineColor(kOrange)
0783 cumulative_matbdg_rec.SetMinimum(0.)
0784 cumulative_matbdg_rec.SetMaximum(3.)
0785 cumulative_matbdg_rec.SetLineColor(kAzure+1)
0786 l = TLegend(0.18, 0.8, 0.95, 0.92)
0787 l.AddEntry(cumulative_matbdg_sim, "Sim Material", "f")
0788 l.AddEntry(cumulative_matbdg_rec, "Reco Material", "f")
0789 cumulative_matbdg_sim.Draw("HIST")
0790 cumulative_matbdg_rec.Draw("HIST SAME")
0791 l.Draw()
0792 filename = "MaterialBdg_Reco_vs_Simul_%s.png" % label
0793 cc.SaveAs(filename)
0794
0795 def createCompoundPlots(detector, plot, geometry):
0796 """Produce the requested plot for the specified detector.
0797
0798 Function that will plot the requested @plot for the specified
0799 @detector. The specified detector could either be a real
0800 detector or a compound one. The list of available plots are the
0801 keys of plots dictionary (imported from plot_utils.
0802
0803 """
0804 setTDRStyle()
0805
0806 theDirname = 'Images'
0807 if not checkFile_(theDirname):
0808 os.mkdir(theDirname)
0809
0810 goodToGo, theDetectorFilename = paramsGood_(detector, plot, geometry)
0811 if not goodToGo:
0812 return
0813
0814 hist_X0_elements = OrderedDict()
0815
0816
0817 stackTitle = "%s;%s;%s" % (detector,
0818 plots[plot].abscissa,
0819 plots[plot].ordinate)
0820 stack_X0 = THStack("stack_X0", stackTitle);
0821 theLegend = TLegend(0.50, 0.70, 0.70, 0.90);
0822
0823 def setRanges(h):
0824 legendSpace = 1. + 0.3
0825 minY = h.GetYaxis().GetXmin()
0826 maxY = h.GetBinContent(h.GetMaximumBin()) * legendSpace
0827 h.GetYaxis().SetRangeUser(minY, maxY)
0828
0829 for label, [num, color, leg] in hist_label_to_num.items():
0830
0831 if label == 'SUM':
0832 continue
0833 hist_X0_elements[label] = get1DHisto_(detector,
0834 num + plots[plot].plotNumber,
0835 geometry)
0836 hist_X0_elements[label].SetFillColor(color)
0837 hist_X0_elements[label].SetLineColor(kBlack)
0838 stack_X0.Add(hist_X0_elements[label])
0839 if hist_X0_elements[label].Integral() > 0.: theLegend.AddEntry(hist_X0_elements[label], leg, "f")
0840
0841
0842 canname = "MBCan_1D_%s_%s" % (detector, plot)
0843 can = TCanvas(canname, canname, 800, 800)
0844 can.Range(0,0,25,25)
0845 gStyle.SetOptTitle(0)
0846
0847
0848 setRanges(stack_X0.GetStack().Last())
0849 stack_X0.Draw("HIST");
0850 stack_X0.GetXaxis().SetLabelSize(0.035)
0851 stack_X0.GetYaxis().SetLabelSize(0.035)
0852 theLegend.Draw();
0853
0854 cmsMark = TLatex()
0855 cmsMark.SetNDC();
0856 cmsMark.SetTextAngle(0);
0857 cmsMark.SetTextColor(kBlack);
0858 cmsMark.SetTextFont(61)
0859 cmsMark.SetTextSize(5e-2)
0860 cmsMark.SetTextAlign(11)
0861 cmsMark.DrawLatex(0.16,0.86,"CMS")
0862
0863 simuMark = TLatex()
0864 simuMark.SetNDC();
0865 simuMark.SetTextAngle(0);
0866 simuMark.SetTextColor(kBlack);
0867 simuMark.SetTextSize(3e-2)
0868 simuMark.SetTextAlign(11)
0869 simuMark.DrawLatex(0.16,0.82,"#font[52]{Simulation Internal}")
0870
0871
0872 can.Update();
0873 can.SaveAs( "%s/%s_%s_%s.pdf"
0874 % (theDirname, detector, plot, geometry))
0875 can.SaveAs( "%s/%s_%s_%s.png"
0876 % (theDirname, detector, plot, geometry))
0877
0878
0879 def create2DPlots(detector, plot, geometry):
0880 """Produce the requested plot for the specified detector.
0881
0882 Function that will plot the requested 2D-@plot for the
0883 specified @detector. The specified detector could either be a
0884 real detector or a compound one. The list of available plots
0885 are the keys of plots dictionary (imported from plot_utils).
0886
0887 """
0888
0889 theDirname = 'Images'
0890 if not checkFile_(theDirname):
0891 os.mkdir(theDirname)
0892
0893 goodToGo, theDetectorFilename = paramsGood_(detector, plot, geometry)
0894 if not goodToGo:
0895 return
0896
0897 theDetectorFile = TFile(theDetectorFilename)
0898
0899 hist_X0_total = get2DHisto_(detector,plots[plot].plotNumber,geometry)
0900
0901
0902 gStyle.SetStripDecimals(False)
0903
0904
0905 if plots[plot].iRebin:
0906 hist_X0_total.Rebin2D()
0907
0908
0909 hist2dTitle = ('%s %s;%s;%s;%s' % (plots[plot].quotaName,
0910 detector,
0911 plots[plot].abscissa,
0912 plots[plot].ordinate,
0913 plots[plot].quotaName))
0914
0915 hist_X0_total.SetTitle(hist2dTitle)
0916 hist_X0_total.SetTitleOffset(0.5,"Y")
0917
0918 if plots[plot].histoMin != -1.:
0919 hist_X0_total.SetMinimum(plots[plot].histoMin)
0920 if plots[plot].histoMax != -1.:
0921 hist_X0_total.SetMaximum(plots[plot].histoMax)
0922
0923 can2name = "MBCan_2D_%s_%s" % (detector, plot)
0924 can2 = TCanvas(can2name, can2name, 2480+248, 580+58+58)
0925 can2.SetTopMargin(0.1)
0926 can2.SetBottomMargin(0.1)
0927 can2.SetLeftMargin(0.04)
0928 can2.SetRightMargin(0.06)
0929 can2.SetFillColor(kWhite)
0930 gStyle.SetOptStat(0)
0931 gStyle.SetTitleFillColor(0)
0932 gStyle.SetTitleBorderSize(0)
0933
0934
0935 gStyle.SetPalette(kGreyScale)
0936
0937
0938 can2.SetLogz(plots[plot].zLog)
0939
0940
0941 hist_X0_total.Draw("COLZ")
0942
0943
0944 can2.Update()
0945
0946
0947 setUpPalette(hist_X0_total,plot)
0948
0949
0950 keep_alive = []
0951 if plots[plot].iDrawEta:
0952 keep_alive.extend(drawEtaValues())
0953
0954 can2.Modified()
0955 hist_X0_total.SetContour(255)
0956
0957
0958 can2.Update()
0959 can2.Modified()
0960
0961 can2.SaveAs( "%s/%s_%s_%s_bw.pdf"
0962 % (theDirname, detector, plot, geometry))
0963 can2.SaveAs( "%s/%s_%s_%s_bw.png"
0964 % (theDirname, detector, plot, geometry))
0965 gStyle.SetStripDecimals(True)
0966
0967 def createRatioPlots(detector, plot, geometry):
0968 """Create ratio plots.
0969
0970 Function that will make the ratio between the radiation length
0971 and interaction length, for the specified detector. The
0972 specified detector could either be a real detector or a
0973 compound one.
0974
0975 """
0976
0977 goodToGo, theDetectorFilename = paramsGood_(detector, plot, geometry)
0978 if not goodToGo:
0979 return
0980
0981 theDirname = 'Images'
0982 if not os.path.exists(theDirname):
0983 os.mkdir(theDirname)
0984
0985 theDetectorFile = TFile(theDetectorFilename)
0986
0987 prof_x0_det_total = theDetectorFile.Get('%d' % plots[plot].plotNumber)
0988 prof_l0_det_total = theDetectorFile.Get('%d' % (1000+plots[plot].plotNumber))
0989
0990
0991 hist_x0_total = get1DHisto_(detector,plots[plot].plotNumber,geometry)
0992 hist_l0_total = get1DHisto_(detector,1000+plots[plot].plotNumber,geometry)
0993
0994 hist_x0_over_l0_total = hist_x0_total
0995 hist_x0_over_l0_total.Divide(hist_l0_total)
0996
0997 histTitle = "Material Budget %s;%s;%s" % (detector,
0998 plots[plot].abscissa,
0999 plots[plot].ordinate)
1000 hist_x0_over_l0_total.SetTitle(histTitle)
1001
1002 hist_x0_over_l0_total.SetMarkerStyle(1)
1003 hist_x0_over_l0_total.SetMarkerSize(3)
1004 hist_x0_over_l0_total.SetMarkerColor(kBlue)
1005
1006
1007 canRname = "MBRatio_%s_%s" % (detector, plot)
1008 canR = TCanvas(canRname,canRname,800,800)
1009 canR.Range(0,0,25,25)
1010 canR.SetFillColor(kWhite)
1011 gStyle.SetOptStat(0)
1012
1013
1014 hist_x0_over_l0_total.Draw("E1")
1015
1016
1017 canR.Update()
1018 canR.SaveAs("%s/%s_%s_%s.pdf"
1019 % (theDirname, detector, plot, geometry))
1020 canR.SaveAs("%s/%s_%s_%s.png"
1021 % (theDirname, detector, plot, geometry))
1022
1023 if __name__ == '__main__':
1024 parser = argparse.ArgumentParser(description='Generic Material Plotter',
1025 formatter_class=argparse.ArgumentDefaultsHelpFormatter)
1026 parser.add_argument('-r', '--reco',
1027 help='Input Reconstruction Material file, DQM format')
1028 parser.add_argument('-l', '--label',
1029 help='Label to use in naming the plots')
1030 parser.add_argument('-c', '--compare',
1031 help='Compare simulation and reco materials',
1032 action='store_true',
1033 default=False)
1034 parser.add_argument('-sw','--subdetector-wise',
1035 help="Subdetector-wise categorization. Individual ROOT " \
1036 "files for each subdetector are required.",
1037 action="store_true",
1038 default=False)
1039 parser.add_argument('-s', '--single',
1040 help='Material budget for single detector from simulation',
1041 action='store_true',
1042 default=False)
1043 parser.add_argument('-d', '--detector',
1044 help='Detector for which you want to compute the material budget',
1045 type=str,)
1046 parser.add_argument('-g', '--geometry',
1047 help='Geometry, used to determine filenames',
1048 type=str)
1049 parser.add_argument('-gc', '--geometry-comparison',
1050 help='Compare the material budget for two different geometries'
1051 +'-g should be specied',
1052 type=str)
1053 args = parser.parse_args()
1054
1055
1056 if args.geometry is None:
1057 print("Error, missing geometry. -g is Required")
1058 raise RuntimeError
1059
1060 if args.geometry_comparison and args.geometry is None:
1061 print("Error, geometry comparison requires two geometries.")
1062 print("use -gc option")
1063 raise RuntimeError
1064
1065 if args.geometry_comparison and args.geometry:
1066
1067
1068
1069
1070 required_plots = ["x_vs_eta","x_vs_phi","x_vs_R",
1071 "l_vs_eta","l_vs_phi","l_vs_R"]
1072 required_2Dplots = ["x_vs_eta_vs_phi",
1073 "l_vs_eta_vs_phi",
1074 "x_vs_z_vs_R",
1075 "l_vs_z_vs_R_geocomp",
1076 "x_vs_z_vs_Rsum",
1077 "l_vs_z_vs_Rsum"]
1078
1079 for p in required_plots:
1080 createCompoundPlotsGeometryComparison(args.detector, p, args.geometry,
1081 args.geometry_comparison)
1082 for p in required_2Dplots:
1083 create2DPlotsGeometryComparison(args.detector, p, args.geometry,
1084 args.geometry_comparison)
1085
1086 if args.compare and args.single:
1087 print("Error, too many actions required")
1088 raise RuntimeError
1089
1090 if args.compare:
1091 if args.reco is None:
1092 print("Error, missing input reco file")
1093 raise RuntimeError
1094 if args.label is None:
1095 print("Error, missing label")
1096 raise RuntimeError
1097 materialBudget_Simul_vs_Reco(args.reco, args.label, args.geometry, debug=False)
1098
1099 if args.single and not args.geometry_comparison:
1100 if args.detector is None:
1101 print("Error, missing detector")
1102 raise RuntimeError
1103 required_2Dplots = ["x_vs_eta_vs_phi", "l_vs_eta_vs_phi", "x_vs_z_vs_R",
1104 "l_vs_z_vs_R", "x_vs_z_vs_Rsum", "l_vs_z_vs_Rsum"]
1105 required_plots = ["x_vs_eta", "x_vs_phi", "l_vs_eta", "l_vs_phi"]
1106
1107 required_ratio_plots = ["x_over_l_vs_eta", "x_over_l_vs_phi"]
1108
1109 for p in required_2Dplots:
1110 create2DPlots(args.detector, p, args.geometry)
1111 for p in required_plots:
1112 createCompoundPlots(args.detector, p, args.geometry)
1113 for p in required_ratio_plots:
1114 createRatioPlots(args.detector, p, args.geometry)
1115
1116 if args.subdetector_wise and args.geometry:
1117 required_plots = ["x_vs_eta","l_vs_eta"]
1118
1119 for p in required_plots:
1120 createPlots_(p, args.geometry)
1121