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 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         # Prevent memory leaking by specifing a unique name
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         # Prevent memory leaking by specifing a unique name
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) # Circles
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         # return stylized ratio histoX/histoY
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         # Return stylized histoNew - histoOld
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     # Plotting the different categories
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") #(F)illed Box
0303         legend.AddEntry(gNew,"%s %s [%s]"%(detector,geometryNew,label),"P") #(P)olymarker
0304         legend.SetTextFont(42)
0305         legend.SetTextSize(0.03)
0306         return legend
0307
0308     def setRanges(h):
0309         legendSpace = 1. + 0.3 # 30%
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     ########### Ratio ###########
0319
0320     counter = 0
0321     legends = OrderedDict() #KeepAlive
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         # Ratio
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     ######## Difference ########
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         # Difference
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     # Configure Palette for 2D Histos
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     # Ratio and Difference have the same call
0458     # But different 'Palette' range so we are
0459     # setting the range only for the Ratio
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         # Merge together the "inner barrel detectors".
0561         if subDetector in IBs:
0562             hist_X0_IB = assignOrAddIfExists_(
0563                 hist_X0_IB,
0564                 hist_X0_detectors[subDetector]
0565                 )
0566
0567         # category profiles
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     # colors
0586     for det, color in DETECTORS.items():
0587         setColorIfExists_(hist_X0_detectors, det,  color)
0588
0589     # First Plot: BeamPipe + Pixel + TIB/TID + TOB + TEC + Outside
0590     # stack
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     # canvas
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     # Draw
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     # Legenda
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     # text
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     # Store
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     # Second Plot: BeamPipe + SEN + ELE + CAB + COL + SUP + OTH/AIR +
0639     # Outside stack
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     # canvas
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     # Draw
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     # Legenda
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     # text
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     # Store
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                 # If we miss an object, since we are incrementally
0723                 # searching for consecutive layers, we may safely
0724                 # assume that there are no additional layers and skip
0725                 # to the next detector.
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     # plots
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     # stack
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 # 30%
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         # We don't want the sum to be added as part of the stack
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     # canvas
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     # Draw
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     # Store
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     # # properties
0902     gStyle.SetStripDecimals(False)
0903
0904     # Ratio
0905     if plots[plot].iRebin:
0906         hist_X0_total.Rebin2D()
0907
0908     # stack
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     # Color palette
0935     gStyle.SetPalette(kGreyScale)
0936
0937     # Log?
0938     can2.SetLogz(plots[plot].zLog)
0939
0940     # Draw in colors
0941     hist_X0_total.Draw("COLZ")
0942
0943     # Store
0944     can2.Update()
0945
0946     #Aesthetic
0947     setUpPalette(hist_X0_total,plot)
0948
0949     #Add eta labels
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     # Store
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     # get TProfiles
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     # histos
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     # properties
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     # canvas
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     # Draw
1014     hist_x0_over_l0_total.Draw("E1")
1015
1016     # Store
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         # For the definition of the properties of these graphs
1068         # check plot_utils.py
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