Geometry/test/MaterialBudget.py

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