Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-11-28 23:11:14

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