Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-11-23 04:29:23

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