Geometry/test/MaterialBudgetHGCal.py

0001 #! /usr/bin/env python3
0002
0003 # Pure trick to start ROOT in batch mode, pass this only option to it
0004 # and the rest of the command line options to this code.
0005 from __future__ import print_function
0006 import sys
0007 import numpy as np
0008 import pandas as pd
0009 import matplotlib.pyplot as plt
0010
0011 from array import array
0012 oldargv = sys.argv[:]
0013 sys.argv = [ '-b-' ]
0014 from ROOT import TCanvas, TLegend, TPaveText, THStack, TFile, TLatex, TPaveLabel
0015 from ROOT import TProfile, TProfile2D, TH1D, TH2D, TH2F, TPaletteAxis, TH1, TH1F, TColor, TExec, TLine
0016 from ROOT import kBlack, kWhite, kOrange, kAzure, kBlue
0017 from ROOT import gROOT, gStyle
0018 gROOT.SetBatch(True)
0019 sys.argv = oldargv
0020
0021 from Validation.Geometry.plot_utils import setTDRStyle, Plot_params, drawEtaValues
0022
0023 from Validation.Geometry.plot_hgcal_utils import plots, COMPOUNDS, DETECTORS, sDETS, hist_label_to_num, TwikiPrintout, acustompalette, dEdx, MatXo, drawHalfEtaValues
0024 from collections import namedtuple, OrderedDict
0025 import sys, os
0026 import argparse
0027
0028 def paramsGood_(detector, plot):
0029     """Check the validity of the arguments.
0030
0031        Common function to check the validity of the parameters passed
0032        in. It returns a tuple composed by a bool and a string. The
0033        bool indicates if all checks are ok, the string the name of the
0034        appropriate ROOT file to open (empty string in case the any
0035        check failed)
0036
0037     """
0038
0039     if plot not in plots.keys():
0040         print("Error, unknown plot %s" % plot)
0041         return (False, '')
0042
0043     if detector not in DETECTORS and detector not in COMPOUNDS.keys():
0044         print('Error, unknown detector: %s' % detector)
0045         return (False, '')
0046
0047     theDetectorFilename = ''
0048     if detector in DETECTORS:
0049         theDetectorFilename = 'matbdg_%s.root' % detector
0050     else:
0051         theDetectorFilename = 'matbdg_%s.root' % COMPOUNDS[detector][0]
0052
0053     if not checkFile_(theDetectorFilename):
0054         print("Error, missing file %s" % theDetectorFilename)
0055         raise RuntimeError
0056     return (True, theDetectorFilename)
0057
0058 def checkFile_(filename):
0059     return os.path.exists(filename)
0060
0061 def setColorIfExists_(histos, h, color):
0062     if h in histos.keys():
0063         histos[h].SetFillColor(color)
0064
0065 def assignOrAddIfExists_(h, p):
0066     """Assign the projection of p to h.
0067
0068        Function to assign the projection of p to h, in the case in
0069        which h is None, otherwise add the projection to the already
0070        valid h object
0071
0072     """
0073
0074     if not h:
0075         h = p.ProjectionX()
0076     else:
0077         h.Add(p.ProjectionX("B_%s" % h.GetName()), +1.000)
0078     return h
0079
0080 def createPlots_(plot, compounddetectorname):
0081     """Cumulative material budget from simulation.
0082
0083        Internal function that will produce a cumulative profile of the
0084        material budget inferred from the simulation starting from the
0085        single detectors that compose the tracker. It will iterate over
0086        all existing detectors contained in the DETECTORS
0087        dictionary. The function will automatically skip non-existent
0088        detectors.
0089
0090     """
0091
0092     theDirname = "Figures"
0093
0094     hist_X0_detectors = OrderedDict()
0095     if plot not in plots.keys():
0096         print("Error: chosen plot name not known %s" % plot)
0097         return
0098
0099     # We need to keep the file content alive for the lifetime of the
0100     # full function....
0101     subDetectorFiles = []
0102
0103     hist_X0_elements = OrderedDict()
0104     prof_X0_elements = OrderedDict()
0105     for subDetector,color in DETECTORS.items():
0106         subDetectorFilename = "matbdg_%s.root" % subDetector
0107         if not checkFile_(subDetectorFilename):
0108             print("Error opening file: %s" % subDetectorFilename)
0109             continue
0110
0111         subDetectorFiles.append(TFile(subDetectorFilename))
0112         subDetectorFile = subDetectorFiles[-1]
0113         print ("Opening file: %s" % subDetectorFilename)
0114         prof_X0_XXX = subDetectorFile.Get("%d" % plots[plot].plotNumber)
0115
0116         hist_X0_detectors[subDetector] = prof_X0_XXX.ProjectionX()
0117
0118         # category profiles
0119         for label, [num, color, leg] in hist_label_to_num.items():
0120             prof_X0_elements[label] = subDetectorFile.Get("%d" % (num + plots[plot].plotNumber))
0121             hist_X0_elements[label] = assignOrAddIfExists_(hist_X0_elements.setdefault(label, None),
0122                                                           prof_X0_elements[label])
0123
0124     cumulative_matbdg = TH1D("CumulativeSimulMatBdg",
0125                              "CumulativeSimulMatBdg",
0126                              hist_X0_detectors["BeamPipe"].GetNbinsX(),
0127                              hist_X0_detectors["BeamPipe"].GetXaxis().GetXmin(),
0128                              hist_X0_detectors["BeamPipe"].GetXaxis().GetXmax())
0129     cumulative_matbdg.SetDirectory(0)
0130
0131     # colors
0132     for det, color in DETECTORS.items():
0133         setColorIfExists_(hist_X0_detectors, det, color)
0134
0135     for label, [num, color, leg] in hist_label_to_num.items():
0136         hist_X0_elements[label].SetFillColor(color)
0137
0138     # First Plot: BeamPipe + Tracker + ECAL + HCal + HGCal + MB + MGNT
0139     # stack
0140     stackTitle_SubDetectors = "Material Budget;%s;%s" % (
0141         plots[plot].abscissa,plots[plot].ordinate)
0142     stack_X0_SubDetectors = THStack("stack_X0",stackTitle_SubDetectors)
0143     for det, histo in hist_X0_detectors.items():
0144         stack_X0_SubDetectors.Add(histo)
0145         cumulative_matbdg.Add(histo, 1)
0146
0147     # canvas
0148     can_SubDetectors = TCanvas("can_SubDetectors","can_SubDetectors",800,800)
0149     #can_SubDetectors.Range(0,0,25,25)
0150     can_SubDetectors.SetFillColor(kWhite)
0151
0152     # Draw
0153     stack_X0_SubDetectors.SetMinimum(plots[plot].ymin)
0154     stack_X0_SubDetectors.SetMaximum(plots[plot].ymax)
0155     stack_X0_SubDetectors.Draw("HIST")
0156     #stack_X0_SubDetectors.GetXaxis().SetLimits(plots[plot].xmin, plots[plot].xmax)
0157
0158
0159     # Legenda
0160     theLegend_SubDetectors = TLegend(0.130,0.7,0.93,0.90) #(0.180,0.8,0.98,0.90)
0161     theLegend_SubDetectors.SetNColumns(2)
0162     theLegend_SubDetectors.SetFillColor(0)
0163     theLegend_SubDetectors.SetFillStyle(0)
0164     theLegend_SubDetectors.SetBorderSize(0)
0165
0166     for det, histo in hist_X0_detectors.items():
0167         theLegend_SubDetectors.AddEntry(histo, det,  "f")
0168
0169     theLegend_SubDetectors.Draw()
0170
0171     # text
0172     text_SubDetectors = TPaveText(0.130,0.627,0.352,0.687,"NDC") #(0.180,0.727,0.402,0.787,"NDC")
0173     text_SubDetectors.SetFillColor(0)
0174     text_SubDetectors.SetBorderSize(0)
0175     text_SubDetectors.AddText("CMS Simulation")
0176     text_SubDetectors.SetTextAlign(11)
0177     text_SubDetectors.Draw()
0178
0179     # Store
0180     can_SubDetectors.Update()
0181     if not checkFile_(theDirname):
0182         os.mkdir(theDirname)
0183     can_SubDetectors.SaveAs("%s/MaterialBdg_%s_%s.pdf" % (theDirname, compounddetectorname,plot))
0184     can_SubDetectors.SaveAs("%s/MaterialBdg_%s_%s.png" % (theDirname, compounddetectorname,plot))
0185     can_SubDetectors.SaveAs("%s/MaterialBdg_%s_%s.root" % (theDirname, compounddetectorname,plot))
0186
0187     if plot == "x_vs_eta" or plot == "l_vs_eta":
0188         canname = "MBCan_1D_%s_%s_total"  % (compounddetectorname, plot)
0189         can2 = TCanvas(canname, canname, 800, 800)
0190         can2.Range(0,0,25,25)
0191         can2.SetFillColor(kWhite)
0192         gStyle.SetOptStat(0)
0193         gStyle.SetOptTitle(0);
0194         #title = TPaveLabel(.11,.95,.35,.99,"Total accumulated material budget","brndc")
0195         stack_X0_SubDetectors.GetStack().Last().SetMarkerStyle(34)
0196         stack_X0_SubDetectors.GetStack().Last().GetXaxis().SetRangeUser( 1.0, 3.5)
0197         stack_X0_SubDetectors.GetStack().Last().Draw();
0198         stack_X0_SubDetectors.GetYaxis().SetTitleOffset(1.15);
0199         can2.Update()
0200         can2.Modified()
0201         can2.SaveAs("%s/%s_%s_total_Zplus.pdf" % (theDirname, compounddetectorname, plot))
0202         can2.SaveAs("%s/%s_%s_total_Zplus.png" % (theDirname, compounddetectorname, plot))
0203         stack_X0_SubDetectors.GetStack().Last().GetXaxis().SetRangeUser( -3.5, -1.0)
0204         stack_X0_SubDetectors.GetStack().Last().Draw();
0205         stack_X0_SubDetectors.GetYaxis().SetTitleOffset(1.15);
0206         can2.Update()
0207         can2.Modified()
0208         can2.SaveAs("%s/%s_%s_total_Zminus.pdf" % (theDirname, compounddetectorname, plot))
0209         can2.SaveAs("%s/%s_%s_total_Zminus.png" % (theDirname, compounddetectorname, plot))
0210
0211         #Also print them to give them exact numbers
0212         etavalues = []
0213         matbudginX0 = []
0214         matbudginIntLen = []
0215         for binx in range(0,stack_X0_SubDetectors.GetStack().Last().GetXaxis().GetNbins()):
0216             bincontent = stack_X0_SubDetectors.GetStack().Last().GetBinContent(binx)
0217             if bincontent == 0: continue
0218             etavalues.append( stack_X0_SubDetectors.GetStack().Last().GetBinCenter(binx)  )
0219             if plot == "x_vs_eta":
0220                 matbudginX0.append( bincontent  )
0221                 d1 = {'Eta': etavalues, 'MatBudInX0': matbudginX0}
0222                 df1 = pd.DataFrame(data=d1).round(2)
0223                 df1.to_csv(r'/afs/cern.ch/work/a/apsallid/CMS/PFCalStudies/CMS-HGCAL/matbudV10fromVertexToBackofHGCal/CMSSW_11_0_X_2019-06-04-2300/src/Validation/Geometry/test/EtavsMatBudinXo.txt',sep=' ', index=False, header=False)
0224                 #print df1
0225             if plot == "l_vs_eta":
0226                 matbudginIntLen.append( bincontent )
0227                 d2 = {'Eta': etavalues, 'MatBudInIntLen': matbudginIntLen}
0228                 df2 = pd.DataFrame(data=d2).round(2)
0229                 df2.to_csv(r'/afs/cern.ch/work/a/apsallid/CMS/PFCalStudies/CMS-HGCAL/matbudV10fromVertexToBackofHGCal/CMSSW_11_0_X_2019-06-04-2300/src/Validation/Geometry/test/EtavsMatBudInIntLen.txt',sep=' ', index=False, header=False)
0230                 #print df2
0231
0232     return cumulative_matbdg
0233
0234 def createPlots2D_(plot, compounddetectorname):
0235     """2D material budget map to know exactly what we are adding.
0236     """
0237
0238     #IBs = ["InnerServices", "Phase2PixelBarrel", "TIB", "TIDF", "TIDB"]
0239     theDirname = "Figures"
0240
0241     hist_X0_detectors = OrderedDict()
0242     if plot not in plots.keys():
0243         print("Error: chosen plot name not known %s" % plot)
0244         return
0245
0246     # We need to keep the file content alive for the lifetime of the
0247     # full function....
0248     subDetectorFiles = []
0249
0250     hist_X0_elements = OrderedDict()
0251     prof_X0_elements = OrderedDict()
0252
0253     for subDetector,color in DETECTORS.items():
0254         subDetectorFilename = "matbdg_%s.root" % subDetector
0255         if not checkFile_(subDetectorFilename):
0256             print("Error opening file: %s" % subDetectorFilename)
0257             continue
0258
0259         subDetectorFiles.append(TFile(subDetectorFilename))
0260         subDetectorFile = subDetectorFiles[-1]
0261         print ("Opening file: %s" % subDetectorFilename)
0262         prof_X0_XXX = subDetectorFile.Get("%d" % plots[plot].plotNumber)
0263
0264         #hist_X0_detectors[subDetector] = prof_X0_XXX
0265         hist_X0_detectors[subDetector] = prof_X0_XXX.ProjectionXY("_pxy","B")
0266         print(subDetector)
0267
0268     # First Plot: BeamPipe + Tracker + ECAL + HCal + HGCal + MB + MGNT
0269
0270     # Create "null" histo
0271     minX = 1.03*hist_X0_detectors["BeamPipe"].GetXaxis().GetXmin()
0272     maxX = 1.03*hist_X0_detectors["BeamPipe"].GetXaxis().GetXmax()
0273     minY = 1.03*hist_X0_detectors["BeamPipe"].GetYaxis().GetXmin()
0274     maxY = 1.03*hist_X0_detectors["BeamPipe"].GetYaxis().GetXmax()
0275
0276     frame = TH2F("frame", "", 10, minX, maxX, 10, minY, maxY);
0277     frame.SetMinimum(0.1)
0278     frame.SetMaximum(10.)
0279     frame.GetXaxis().SetTickLength(frame.GetXaxis().GetTickLength()*0.50)
0280     frame.GetYaxis().SetTickLength(frame.GetXaxis().GetTickLength()/4.)
0281
0282     hist2d_X0_total = hist_X0_detectors["BeamPipe"]
0283
0284     # stack
0285     hist2dTitle = ('%s %s;%s;%s;%s' % (plots[plot].quotaName,
0286                                        "All detectors",
0287                                        plots[plot].abscissa,
0288                                        plots[plot].ordinate,
0289                                        plots[plot].quotaName))
0290
0291     hist2d_X0_total.SetTitle(hist2dTitle)
0292     frame.SetTitle(hist2dTitle)
0293     frame.SetTitleOffset(0.5,"Y")
0294
0295     #If here you put different histomin,histomaxin plot_utils you won't see anything
0296     #for the material plots.
0297     if plots[plot].histoMin != -1.:
0298         hist2d_X0_total.SetMinimum(plots[plot].histoMin)
0299     if plots[plot].histoMax != -1.:
0300         hist2d_X0_total.SetMaximum(plots[plot].histoMax)
0301
0302     #
0303     # canvas
0304     can_SubDetectors = TCanvas("can_SubDetectors","can_SubDetectors",2480+248, 580+58+58)
0305     can_SubDetectors.SetTopMargin(0.1)
0306     can_SubDetectors.SetBottomMargin(0.1)
0307     can_SubDetectors.SetLeftMargin(0.04)
0308     can_SubDetectors.SetRightMargin(0.06)
0309     can_SubDetectors.SetFillColor(kWhite)
0310     gStyle.SetOptStat(0)
0311     gStyle.SetTitleFillColor(0)
0312     gStyle.SetTitleBorderSize(0)
0313     gStyle.SetOptTitle(0)
0314
0315     hist2d_X0_total.GetYaxis().SetTickLength(hist2d_X0_total.GetXaxis().GetTickLength()/4.)
0316     hist2d_X0_total.GetYaxis().SetTickLength(hist2d_X0_total.GetXaxis().GetTickLength()/4.)
0317     hist2d_X0_total.SetTitleOffset(0.5,"Y")
0318     hist2d_X0_total.GetYaxis().SetTitleOffset(0.50);
0319     #hist2d_X0_total.GetXaxis().SetTitleOffset(1.15);
0320     #hist2d_X0_total.GetXaxis().SetNoExponent(True)
0321     #hist2d_X0_total.GetYaxis().SetNoExponent(True)
0322
0323
0324     # colors
0325     for det, color in DETECTORS.items():
0326         hist_X0_detectors[det].SetMarkerColor(color)
0327         hist_X0_detectors[det].SetFillColor(color)
0328
0329     for det, histo in hist_X0_detectors.items():
0330         print(det)
0331         histo.Draw("same")
0332
0333     # Legenda
0334     theLegend_SubDetectors = TLegend(0.100,0.7,0.90,0.90)#(0.180,0.8,0.98,0.90)
0335     theLegend_SubDetectors.SetNColumns(3)
0336     theLegend_SubDetectors.SetFillColor(0)
0337     theLegend_SubDetectors.SetFillStyle(0)
0338     theLegend_SubDetectors.SetBorderSize(0)
0339
0340     for det, histo in hist_X0_detectors.items():
0341         theLegend_SubDetectors.AddEntry(histo, det,  "f")
0342     #theLegend_SubDetectors.AddEntry(hgbound1, "HGCal Eta Boundaries [1.3, 3.0]",  "l")
0343
0344     theLegend_SubDetectors.Draw()
0345
0346
0347     # text
0348     text_SubDetectors = TPaveText(0.100,0.627,0.322,0.687,"NDC")#(0.180,0.727,0.402,0.787,"NDC")
0349     text_SubDetectors.SetFillColor(0)
0350     text_SubDetectors.SetBorderSize(0)
0351     text_SubDetectors.AddText("CMS Simulation")
0352     text_SubDetectors.SetTextAlign(11)
0353     text_SubDetectors.Draw()
0354
0355
0356     #Add eta labels
0357     keep_alive = []
0358     if plots[plot].iDrawEta:
0359         keep_alive.extend(drawEtaValues())
0360
0361     # Store
0362     can_SubDetectors.Update()
0363     if not checkFile_(theDirname):
0364         os.mkdir(theDirname)
0365     can_SubDetectors.SaveAs("%s/MaterialBdg_%s_%s.png" % (theDirname, compounddetectorname, plot))
0366     #It seems that it is too heavy to create .pdf and .root
0367     #can_SubDetectors.SaveAs("%s/MaterialBdg_FromVertexToEndofHGCal_%s.pdf" % (theDirname, plot))
0368     #can_SubDetectors.SaveAs("%s/MaterialBdg_FromVertexToEndofHGCal_%s.root" % (theDirname, plot))
0369
0370     hist2d_X0_total.GetXaxis().SetRangeUser( 0., 7000.)
0371     #Draw eta values in the zoom case
0372     keep_alive = []
0373     keep_alive.extend(drawHalfEtaValues())
0374     #hist2d_X0_total.Draw("COLZ")
0375     can_SubDetectors.Update()
0376     can_SubDetectors.Modified()
0377     can_SubDetectors.SaveAs("%s/MaterialBdg_%s_%s_Zpluszoom.png" % (theDirname, compounddetectorname, plot))
0378
0379 def createPlotsReco_(reco_file, label, debug=False):
0380     """Cumulative material budget from reconstruction.
0381
0382
0383        Internal function that will produce a cumulative profile of the
0384        material budget in the reconstruction starting from the single
0385        detectors that compose the tracker. It will iterate over all
0386        existing detectors contained in the sDETS dictionary. The
0387        function will automatically stop everytime it encounters a
0388        non-existent detector, until no more detectors are left to
0389        try. For this reason the keys in the sDETS dictionary can be as
0390        inclusive as possible.
0391
0392     """
0393
0394     cumulative_matbdg = None
0395     sPREF = ["Original_RadLen_vs_Eta_", "RadLen_vs_Eta_"]
0396
0397     c = TCanvas("c", "c", 1024, 1024);
0398     diffs = []
0399     if not checkFile_(reco_file):
0400         print("Error: missing file %s" % reco_file)
0401         raise RuntimeError
0402     file = TFile(reco_file)
0403     prefix = "/DQMData/Run 1/Tracking/Run summary/RecoMaterial/"
0404     for s in sPREF:
0405         hs = THStack("hs","");
0406         histos = []
0407         for det, color in sDETS.items():
0408             layer_number = 0
0409             while True:
0410                 layer_number += 1
0411                 name = "%s%s%s%d" % (prefix, s, det, layer_number)
0412                 prof = file.Get(name)
0413                 # If we miss an object, since we are incrementally
0414                 # searching for consecutive layers, we may safely
0415                 # assume that there are no additional layers and skip
0416                 # to the next detector.
0417                 if not prof:
0418                     if debug:
0419                         print("Missing profile %s" % name)
0420                     break
0421                 else:
0422                     histos.append(prof.ProjectionX("_px", "hist"));
0423                     diffs.append(histos[-1]);
0424                     histos[-1].SetFillColor(color + layer_number);
0425                     histos[-1].SetLineColor(color + layer_number + 1);
0426
0427         name = "CumulativeRecoMatBdg_%s" % s
0428         if s == "RadLen_vs_Eta_":
0429             cumulative_matbdg = TH1D(name, name,
0430                                      histos[0].GetNbinsX(),
0431                                      histos[0].GetXaxis().GetXmin(),
0432                                      histos[0].GetXaxis().GetXmax())
0433             cumulative_matbdg.SetDirectory(0)
0434         for h in histos:
0435             hs.Add(h)
0436             if cumulative_matbdg:
0437                 cumulative_matbdg.Add(h, 1.)
0438         hs.Draw()
0439         hs.GetYaxis().SetTitle("RadLen")
0440         c.Update()
0441         c.Modified()
0442         c.SaveAs("%sstacked_%s.png" % (s, label))
0443     hs = THStack("diff","")
0444     for d in range(0,len(diffs)/2):
0445         diffs[d+len(diffs)/2].Add(diffs[d], -1.)
0446         hs.Add(diffs[d+len(diffs)/2]);
0447     hs.Draw()
0448     hs.GetYaxis().SetTitle("RadLen")
0449     c.Update()
0450     c.Modified()
0451     c.SaveAs("RadLen_difference_%s.png" % label)
0452     return cumulative_matbdg
0453
0454 def materialBudget_Simul_vs_Reco(reco_file, label, debug=False):
0455     """Plot reco vs simulation material budget.
0456
0457        Function are produces a direct comparison of the material
0458        budget as extracted from the reconstruction geometry and
0459        inferred from the simulation one.
0460
0461     """
0462
0463     setTDRStyle()
0464
0465     # plots
0466     cumulative_matbdg_sim = createPlots_("x_vs_eta")
0467     cumulative_matbdg_rec = createPlotsReco_(reco_file, label, debug=False)
0468
0469     cc = TCanvas("cc", "cc", 1024, 1024)
0470     cumulative_matbdg_sim.SetMinimum(0.)
0471     cumulative_matbdg_sim.SetMaximum(3.5)
0472     cumulative_matbdg_sim.GetXaxis().SetRangeUser(-3.0, 3.0)
0473     cumulative_matbdg_sim.SetLineColor(kOrange)
0474     cumulative_matbdg_rec.SetMinimum(0.)
0475     cumulative_matbdg_rec.SetMaximum(3.)
0476     cumulative_matbdg_rec.SetLineColor(kAzure+1)
0477     l = TLegend(0.18, 0.8, 0.95, 0.92)
0478     l.AddEntry(cumulative_matbdg_sim, "Sim Material", "f")
0479     l.AddEntry(cumulative_matbdg_rec, "Reco Material", "f")
0480     cumulative_matbdg_sim.Draw("HIST")
0481     cumulative_matbdg_rec.Draw("HIST SAME")
0482     l.Draw()
0483     filename = "MaterialBdg_Reco_vs_Simul_%s.png" % label
0484     cc.SaveAs(filename)
0485
0486 def createCompoundPlots(detector, plot):
0487     """Produce the requested plot for the specified detector.
0488
0489        Function that will plot the requested @plot for the specified
0490        @detector. The specified detector could either be a real
0491        detector or a compound one. The list of available plots are the
0492        keys of plots dictionary (imported from plot_utils.
0493
0494     """
0495
0496     theDirname = 'Images'
0497     if not checkFile_(theDirname):
0498         os.mkdir(theDirname)
0499
0500     goodToGo, theDetectorFilename = paramsGood_(detector, plot)
0501     if not goodToGo:
0502         return
0503
0504     theDetectorFile = TFile(theDetectorFilename)
0505     #
0506
0507     # get TProfiles
0508     prof_X0_elements = OrderedDict()
0509     hist_X0_elements = OrderedDict()
0510     for label, [num, color, leg] in hist_label_to_num.items():
0511         #print label, num, color, leg
0512         prof_X0_elements[label] = theDetectorFile.Get("%d" % (num + plots[plot].plotNumber))
0513         hist_X0_elements[label] = prof_X0_elements[label].ProjectionX()
0514         hist_X0_elements[label].SetFillColor(color)
0515         hist_X0_elements[label].SetLineColor(kBlack)
0516
0517     files = []
0518     if detector in COMPOUNDS.keys():
0519         for subDetector in COMPOUNDS[detector][1:]:
0520             subDetectorFilename = "matbdg_%s.root" % subDetector
0521
0522             # open file
0523             if not checkFile_(subDetectorFilename):
0524                 continue
0525
0526             subDetectorFile = TFile(subDetectorFilename)
0527             files.append(subDetectorFile)
0528             print("*** Open file... %s" %  subDetectorFilename)
0529
0530             # subdetector profiles
0531             for label, [num, color, leg] in hist_label_to_num.items():
0532                 prof_X0_elements[label] = subDetectorFile.Get("%d" % (num + plots[plot].plotNumber))
0533                 hist_X0_elements[label].Add(prof_X0_elements[label].ProjectionX("B_%s" % prof_X0_elements[label].GetName())
0534                                             , +1.000)
0535
0536     # stack
0537     stackTitle = "Material Budget %s;%s;%s" % (detector,
0538                                                plots[plot].abscissa,
0539                                                plots[plot].ordinate)
0540     stack_X0 = THStack("stack_X0", stackTitle);
0541     for label, [num, color, leg] in hist_label_to_num.items():
0542         stack_X0.Add(hist_X0_elements[label])
0543
0544     # canvas
0545     canname = "MBCan_1D_%s_%s"  % (detector, plot)
0546     can = TCanvas(canname, canname, 800, 800)
0547     can.Range(0,0,25,25)
0548     can.SetFillColor(kWhite)
0549     gStyle.SetOptStat(0)
0550     gStyle.SetOptTitle(1);
0551
0552     # Draw
0553     stack_X0.Draw("HIST");
0554     stack_X0.GetYaxis().SetTitleOffset(1.15);
0555
0556     # Legenda
0557     theLegend = TLegend(0.40, 0.65, 0.60, 0.89)
0558     if plot == "x_vs_phi" or plot == "l_vs_phi": theLegend = TLegend(0.65, 0.30, 0.89, 0.70)
0559     if plot == "x_vs_R" or plot == "l_vs_R": theLegend = TLegend(0.75, 0.60, 0.95, 0.90)
0560
0561     for label, [num, color, leg] in hist_label_to_num.items():
0562         theLegend.AddEntry(hist_X0_elements[label], leg, "f")
0563     theLegend.Draw();
0564
0565     # Store
0566     can.Update();
0567     can.SaveAs( "%s/%s_%s.pdf" % (theDirname, detector, plot))
0568     can.SaveAs( "%s/%s_%s.png" % (theDirname, detector, plot))
0569
0570     #Let's also save the total accumulated budget vs eta since muon id relies
0571     #on adequate calorimeter thickness
0572     if plot == "x_vs_eta" or plot == "l_vs_eta":
0573         canname = "MBCan_1D_%s_%s_total"  % (detector, plot)
0574         can2 = TCanvas(canname, canname, 800, 800)
0575         can2.Range(0,0,25,25)
0576         can2.SetFillColor(kWhite)
0577         gStyle.SetOptStat(0)
0578         gStyle.SetOptTitle(0);
0579         #title = TPaveLabel(.11,.95,.35,.99,"Total accumulated material budget","brndc")
0580         stack_X0.GetStack().Last().GetXaxis().SetRangeUser( 0., 3.)
0581         stack_X0.GetStack().Last().Draw();
0582         stack_X0.GetYaxis().SetTitleOffset(1.15);
0583         can2.Update()
0584         can2.Modified()
0585         can2.SaveAs("%s/%s_%s_total_Zplus.pdf" % (theDirname, detector, plot))
0586         can2.SaveAs("%s/%s_%s_total_Zplus.png" % (theDirname, detector, plot))
0587         stack_X0.GetStack().Last().GetXaxis().SetRangeUser( -3., 0.)
0588         stack_X0.GetStack().Last().Draw();
0589         stack_X0.GetYaxis().SetTitleOffset(1.15);
0590         can2.Update()
0591         can2.Modified()
0592         can2.SaveAs("%s/%s_%s_total_Zminus.pdf" % (theDirname, detector, plot))
0593         can2.SaveAs("%s/%s_%s_total_Zminus.png" % (theDirname, detector, plot))
0594
0595
0596 def create2DPlots(detector, plot, plotnum, plotmat, dosingledetector = True):
0597     """Produce the requested plot for the specified detector.
0598
0599        Function that will plot the requested 2D-@plot for the
0600        specified @detector. The specified detector could either be a
0601        real detector or a compound one. The list of available plots
0602        are the keys of plots dictionary imported from plot_utils.
0603     """
0604
0605     #gStyle.Reset()
0606     #Better to use an underscore.
0607     plotmat = plotmat.replace(" ", "_")
0608
0609
0610     if plotmat != "":
0611         theDirname = ('Images/%s' % plotmat).replace(" ", "")
0612     else:
0613         theDirname = 'Images'
0614
0615     if not checkFile_(theDirname):
0616         os.mkdir(theDirname)
0617     if not os.path.isdir(('Images/%s/ZPlusZoom' % plotmat).replace(" ", "")):
0618         os.mkdir( ('Images/%s/ZPlusZoom' % plotmat).replace(" ", "") )
0619     if not os.path.isdir(('Images/%s/ZMinusZoom' % plotmat).replace(" ", "")):
0620         os.mkdir( ('Images/%s/ZMinusZoom' % plotmat).replace(" ", "") )
0621
0622     goodToGo, theDetectorFilename = paramsGood_(detector, plot)
0623     if not goodToGo:
0624         return
0625
0626     theDetectorFile = TFile(theDetectorFilename)
0627
0628     prof2d_X0_det_total = TProfile2D()
0629     prof2d_X0_det_total.Reset()
0630
0631     # get TProfiles
0632     #prof2d_X0_det_total = theDetectorFile.Get('%s' % plots[plot].plotNumber)
0633     prof2d_X0_det_total = theDetectorFile.Get('%s' % plotnum)
0634     print("==================================================================")
0635     print(plotnum)
0636
0637     # histos
0638     prof2d_X0_det_total.__class__ = TProfile2D
0639     hist_X0_total = prof2d_X0_det_total.ProjectionXY()
0640
0641     # keep files live forever
0642     files = []
0643     if detector in COMPOUNDS.keys() and not dosingledetector:
0644         #When the loop was:
0645         #for subDetector in COMPOUNDS[detector][1:]:
0646         #and the detector was single it never went in the loop and read the single file
0647         #from above. I alter this to COMPOUNDS[detector] to do the multi material budget plot.
0648         #This won't effect the single detector due to the alter in the if above
0649         for subDetector in COMPOUNDS[detector]:
0650             # filenames of single components
0651             subDetectorFilename = "matbdg_%s.root" % subDetector
0652
0653             # open file
0654             if not checkFile_(subDetectorFilename):
0655                 print("Error, missing file %s" % subDetectorFilename)
0656                 continue
0657
0658             subDetectorFile = TFile(subDetectorFilename)
0659             files.append(subDetectorFile)
0660             print("*** Open file... %s" %  subDetectorFilename)
0661
0662             # subdetector profiles
0663             prof2d_X0_det_total = subDetectorFile.Get('%s' % plots[plot].plotNumber)
0664             prof2d_X0_det_total.__class__ = TProfile2D
0665
0666             # add to summary histogram
0667             hist_X0_total.Add(prof2d_X0_det_total.ProjectionXY("B_%s" % prof2d_X0_det_total.GetName()), +1.000 )
0668
0669     # # properties
0670     #gStyle.SetPalette(1)
0671     gStyle.SetStripDecimals(False)
0672     # #
0673
0674     # Create "null" histo
0675     minX = 1.03*prof2d_X0_det_total.GetXaxis().GetXmin()
0676     maxX = 1.03*prof2d_X0_det_total.GetXaxis().GetXmax()
0677     minY = 1.03*prof2d_X0_det_total.GetYaxis().GetXmin()
0678     maxY = 1.03*prof2d_X0_det_total.GetYaxis().GetXmax()
0679
0680     frame = TH2F("frame", "", 10, minX, maxX, 10, minY, maxY);
0681     frame.SetMinimum(0.1)
0682     frame.SetMaximum(10.)
0683     frame.GetXaxis().SetTickLength(frame.GetXaxis().GetTickLength()*0.50)
0684     frame.GetYaxis().SetTickLength(frame.GetXaxis().GetTickLength()/4.)
0685
0686     # Ratio
0687     if plots[plot].iRebin:
0688         prof2d_X0_det_total.Rebin2D()
0689
0690     # stack
0691     hist2dTitle = ('%s %s;%s;%s;%s' % (plots[plot].quotaName,
0692                                        detector,
0693                                        plots[plot].abscissa,
0694                                        plots[plot].ordinate,
0695                                        plots[plot].quotaName))
0696
0697     if dosingledetector:
0698         hist2d_X0_total = prof2d_X0_det_total
0699     else :
0700         hist2d_X0_total = hist_X0_total
0701     hist2d_X0_total.SetTitle(hist2dTitle)
0702     frame.SetTitle(hist2dTitle)
0703     frame.SetTitleOffset(0.5,"Y")
0704
0705     #If here you put different histomin,histomaxin plot_utils you won't see anything
0706     #for the material plots.
0707     if plots[plot].histoMin != -1.:
0708         hist2d_X0_total.SetMinimum(plots[plot].histoMin)
0709     if plots[plot].histoMax != -1.:
0710         hist2d_X0_total.SetMaximum(plots[plot].histoMax)
0711
0712     #
0713     can2name = "MBCan_2D_%s_%s_%s" % (detector, plot, plotmat)
0714     can2 = TCanvas(can2name, can2name, 2480+248, 580+58+58)
0715     can2.SetTopMargin(0.1)
0716     can2.SetBottomMargin(0.1)
0717     can2.SetLeftMargin(0.04)
0718     can2.SetRightMargin(0.06)
0719     can2.SetFillColor(kWhite)
0720     gStyle.SetOptStat(0)
0721     gStyle.SetTitleFillColor(0)
0722     gStyle.SetTitleBorderSize(0)
0723
0724     #hist2d_X0_total.SetMaximum(hist2d_X0_total.GetMaximum())
0725     # Color palette
0726     # gStyle.SetPalette()#1
0727     acustompalette()
0728     ex1 = TExec("ex1","acustompalette();");
0729     ex1.Draw();
0730
0731     #for i in range(100): MyPaletteArray.append(i+1)
0732
0733     #gStyle.SetPalette(first_color_number);
0734
0735     # Log?
0736     can2.SetLogz(plots[plot].zLog)
0737
0738     # Draw in colors
0739     #frame.Draw()
0740     #hist2d_X0_total.Draw("COLZsame") #Dummy draw to create the palette object
0741     hist2d_X0_total.Draw("COLZ") #Dummy draw to create the palette object
0742
0743     # Store
0744     can2.Update()
0745
0746     #Aesthetic
0747     palette = hist2d_X0_total.GetListOfFunctions().FindObject("palette")
0748     if palette:
0749         palette.__class__ = TPaletteAxis
0750         palette.SetX1NDC(0.945)
0751         palette.SetX2NDC(0.96)
0752         palette.SetY1NDC(0.1)
0753         palette.SetY2NDC(0.9)
0754         palette.GetAxis().SetTickSize(.01)
0755         palette.GetAxis().SetTitle("")
0756         if plots[plot].zLog:
0757             palette.GetAxis().SetLabelOffset(-0.01)
0758     paletteTitle = TLatex(1.12*maxX, maxY, plots[plot].quotaName)
0759     paletteTitle.SetTextAngle(90.)
0760     paletteTitle.SetTextSize(0.05)
0761     paletteTitle.SetTextAlign(31)
0762     paletteTitle.Draw()
0763     hist2d_X0_total.GetYaxis().SetTickLength(hist2d_X0_total.GetXaxis().GetTickLength()/4.)
0764     hist2d_X0_total.GetYaxis().SetTickLength(hist2d_X0_total.GetXaxis().GetTickLength()/4.)
0765     hist2d_X0_total.SetTitleOffset(0.5,"Y")
0766     hist2d_X0_total.GetYaxis().SetTitleOffset(0.45);
0767     #hist2d_X0_total.GetXaxis().SetTitleOffset(1.15);
0768     #hist2d_X0_total.GetXaxis().SetNoExponent(True)
0769     #hist2d_X0_total.GetYaxis().SetNoExponent(True)
0770
0771     #Add eta labels
0772     keep_alive = []
0773     if plots[plot].iDrawEta:
0774         keep_alive.extend(drawEtaValues())
0775
0776     can2.Modified()
0777     hist2d_X0_total.SetContour(255)
0778
0779     # Store
0780     can2.Update()
0781     can2.Modified()
0782
0783     can2.SaveAs( "%s/%s_%s%s.pdf" % (theDirname, detector, plot, plotmat))
0784     can2.SaveAs( "%s/%s_%s%s.png" % (theDirname, detector, plot, plotmat))
0785     #can2.SaveAs( "%s/%s_%s%s.root" % (theDirname, detector, plot, plotmat))
0786
0787     #Zoom in a little bit
0788     if plot == "x_vs_z_vs_Rsum" or plot == "l_vs_z_vs_Rsum" or plot == "x_vs_z_vs_Rsumcos" or plot == "l_vs_z_vs_Rsumcos" or plot == "x_vs_z_vs_Rloc" or plot == "l_vs_z_vs_Rloc" or  plot == "x_vs_z_vs_Rloccos" or plot == "l_vs_z_vs_Rloccos":
0789         #Z+
0790         #hist2d_X0_total.GetXaxis().SetLimits( 3100., 5200.)
0791         if dosingledetector: hist2d_X0_total.GetXaxis().SetRangeUser( 3100., 5400.)
0792         else : hist2d_X0_total.GetXaxis().SetRangeUser( 0., 7000.)
0793         #Do not draw eta values in the zoom case
0794         keep_alive = []
0795         #hist2d_X0_total.Draw("COLZ")
0796         can2.Update()
0797         can2.Modified()
0798         can2.SaveAs( "%s/%s/%s_%s%s_ZplusZoom.pdf" % (theDirname, "ZPlusZoom", detector, plot, plotmat))
0799         can2.SaveAs( "%s/%s/%s_%s%s_ZplusZoom.png" % (theDirname, "ZPlusZoom", detector, plot, plotmat))
0800         #Z-
0801         #hist2d_X0_total.GetXaxis().SetLimits( 3100., 5200.)
0802         if dosingledetector: hist2d_X0_total.GetXaxis().SetRangeUser( -5400., -3100.)
0803         else : hist2d_X0_total.GetXaxis().SetRangeUser( 0., -7000.)
0804         #Do not draw eta values in the zoom case
0805         keep_alive = []
0806         #hist2d_X0_total.Draw("COLZ")
0807         can2.Update()
0808         can2.Modified()
0809         can2.SaveAs( "%s/%s/%s_%s%s_ZminusZoom.pdf" % (theDirname, "ZMinusZoom", detector, plot, plotmat))
0810         can2.SaveAs( "%s/%s/%s_%s%s_ZminusZoom.png" % (theDirname, "ZMinusZoom", detector, plot, plotmat))
0811
0812
0813     gStyle.SetStripDecimals(True)
0814
0815 def createRatioPlots(detector, plot):
0816     """Create ratio plots.
0817
0818        Function that will make the ratio between the radiation length
0819        and interaction length, for the specified detector. The
0820        specified detector could either be a real detector or a
0821        compound one.
0822
0823     """
0824
0825     goodToGo, theDetectorFilename = paramsGood_(detector, plot)
0826     if not goodToGo:
0827         return
0828
0829     theDirname = 'Images'
0830     if not os.path.exists(theDirname):
0831         os.mkdir(theDirname)
0832
0833     theDetectorFile = TFile(theDetectorFilename)
0834     # get TProfiles
0835     prof_x0_det_total = theDetectorFile.Get('%d' % plots[plot].plotNumber)
0836     prof_l0_det_total = theDetectorFile.Get('%d' % (10000+plots[plot].plotNumber))
0837
0838     # histos
0839     hist_x0_total = prof_x0_det_total.ProjectionX()
0840     hist_l0_total = prof_l0_det_total.ProjectionX()
0841
0842     if detector in COMPOUNDS.keys():
0843         for subDetector in COMPOUNDS[detector][1:]:
0844
0845             # file name
0846             subDetectorFilename = "matbdg_%s.root" % subDetector
0847
0848             # open file
0849             if not checkFile_(subDetectorFilename):
0850                 print("Error, missing file %s" % subDetectorFilename)
0851                 continue
0852
0853             subDetectorFile = TFile(subDetectorFilename)
0854
0855             # subdetector profiles
0856             prof_x0_det_total = subDetectorFile.Get('%d' % plots[plot].plotNumber)
0857             prof_l0_det_total = subDetectorFile.Get('%d' % (10000+plots[plot].plotNumber))
0858             # add to summary histogram
0859             hist_x0_total.Add(prof_x0_det_total.ProjectionX("B_%s" % prof_x0_det_total.GetName()), +1.000 )
0860             hist_l0_total.Add(prof_l0_det_total.ProjectionX("B_%s" % prof_l0_det_total.GetName()), +1.000 )
0861     #
0862     hist_x0_over_l0_total = hist_x0_total
0863     hist_x0_over_l0_total.Divide(hist_l0_total)
0864     histTitle = "Material Budget %s;%s;%s" % (detector,
0865                                               plots[plot].abscissa,
0866                                               plots[plot].ordinate)
0867     hist_x0_over_l0_total.SetTitle(histTitle)
0868     # properties
0869     hist_x0_over_l0_total.SetMarkerStyle(1)
0870     hist_x0_over_l0_total.SetMarkerSize(3)
0871     hist_x0_over_l0_total.SetMarkerColor(kBlue)
0872
0873     # canvas
0874     canRname = "MBRatio_%s_%s" % (detector, plot)
0875     canR = TCanvas(canRname,canRname,800,800)
0876     canR.Range(0,0,25,25)
0877     canR.SetFillColor(kWhite)
0878     gStyle.SetOptStat(0)
0879
0880     # Draw
0881     hist_x0_over_l0_total.Draw("E1")
0882
0883     # Store
0884     canR.Update()
0885     canR.SaveAs("%s/%s_%s.pdf" % (theDirname, detector, plot))
0886     canR.SaveAs("%s/%s_%s.png" % (theDirname, detector, plot))
0887
0888
0889 def GetSiliconZValuesFromXML():
0890     """Someone can use the xml files of the geometry. Just do
0891     cd $CMSSW_BASE/src/Geometry/HGCalCommonData/test
0892     cmsRun testHGCalParameters_cfg.py
0893     The scintillator position in the HEB is very close to the silicon but they have a
0894     1.95 mm difference. So, we take the silicon position.
0895
0896     This returns a np.array of the z position of the silicon values in mm.
0897     """
0898
0899     #First as they come out of the testHGCalParameters_cfg.py as cm
0900     #plus either the thermal screen front (297.0 cm)
0901     #or the EE front (319.05 cm)
0902     layersmaxZ = ( 319.05, 319.815, 320.725, 322.255, 323.165, 324.695, 325.605, 327.135, 328.045, 329.575, 330.485, 332.015, 332.925, 334.455, 335.365, 336.895, 337.805, 339.335, 340.245, 341.775, 342.685, 344.215, 345.125, 346.655, 347.565, 349.095, 350.005, 351.535, 352.445, 357.715, 362.615, 367.515, 372.415, 377.315, 382.215, 387.115, 392.015, 397.105, 402.195, 407.285, 412.375, 420.765, 429.155, 437.545, 445.935, 454.325, 462.715, 471.105, 479.495, 487.885, 496.275, 504.665, 513.055)
0903
0904     layersmaxZ = np.asarray(layersmaxZ)
0905     layersmaxZ = 10. * layersmaxZ # in mm
0906
0907     print("Total number of layers from XML " , layersmaxZ.size - 1) # Minus 1 for the EE front input by hand.
0908
0909     return layersmaxZ
0910
0911
0912 if __name__ == '__main__':
0913     parser = argparse.ArgumentParser(description='Generic Material Plotter',
0914                                      formatter_class=argparse.ArgumentDefaultsHelpFormatter)
0915     parser.add_argument('-r', '--reco',
0916                         help='Input Reconstruction Material file, DQM format')
0917     parser.add_argument('-l', '--label',
0918                         help='Label to use in naming the plots')
0919     parser.add_argument('-c', '--compare',
0920                         help='Compare simulation and reco materials',
0921                         action='store_true',
0922                         default=False)
0923     parser.add_argument('-m', '--multi',
0924                         help='Combining multiple detectors',
0925                         action='store_true',
0926                         default=False)
0927     parser.add_argument('-s', '--single',
0928                         help='Material budget for single detector from simulation',
0929                         action='store_true',
0930                         default=False)
0931     parser.add_argument('-d', '--detector',
0932                         help='Detector for which you want to compute the material budget',
0933                         type=str,)
0934     args = parser.parse_args()
0935
0936     if args.compare and args.single:
0937         print("Error, too many actions required")
0938         raise RuntimeError
0939
0940     if args.compare:
0941         if args.reco is None:
0942             print("Error, missing inpur reco file")
0943             raise RuntimeError
0944         if args.label is None:
0945             print("Error, missing label")
0946             raise RuntimeError
0947         materialBudget_Simul_vs_Reco(args.reco, args.label, debug=False)
0948
0949     if args.multi:
0950         #Material Budget plot from Vertex to end of HGCal
0951         cumulative_matbdg_sim = createPlots_("x_vs_eta",args.detector)
0952         cumulative_matbdg_sim = createPlots_("l_vs_eta",args.detector)
0953         cumulative_matbdg_sim = createPlots2D_("x_vs_z_vs_Rsum",args.detector)
0954         cumulative_matbdg_sim = createPlots2D_("l_vs_z_vs_Rsum",args.detector)
0955
0956         required_2DplotsForalldetectors = ["x_vs_z_vs_Rsum", "l_vs_z_vs_Rsum"]
0957
0958         for p in required_2DplotsForalldetectors:
0959             #First the total
0960             create2DPlots(args.detector, p, plots[p].plotNumber, "", not args.multi)
0961             #Then, the rest
0962             #for label, [num, color, leg] in hist_label_to_num.iteritems():
0963                 #print label, num, color, leg
0964                 #create2DPlots(args.detector, p, num + plots[p].plotNumber, leg)
0965
0966     if args.single:
0967         if args.detector is None:
0968             print("Error, missing detector")
0969             raise RuntimeError
0970         required_2Dplots = ["x_vs_eta_vs_phi", "l_vs_eta_vs_phi", "x_vs_z_vs_Rsum", "l_vs_z_vs_Rsum", "x_vs_z_vs_Rsumcos", "l_vs_z_vs_Rsumcos", "x_vs_z_vs_Rloc", "l_vs_z_vs_Rloc", "x_vs_z_vs_Rloccos", "l_vs_z_vs_Rloccos"]
0971         # The plots below require z bin 1 cm to have more statistics for eta, phi
0972         # while the plots above are with 200000 events for z bin 1 mm.
0973         #required_2Dplots = ["x_vs_eta_vs_phi", "l_vs_eta_vs_phi", "x_vs_z_vs_Rsum", "l_vs_z_vs_Rsum", "x_vs_z_vs_Rsumcos"]
0974
0975         required_plots = ["x_vs_eta", "x_vs_phi", "x_vs_R", "l_vs_eta", "l_vs_phi", "l_vs_R"]
0976
0977         required_ratio_plots = ["x_over_l_vs_eta", "x_over_l_vs_phi"]
0978
0979         #Function to help filling the twiki with all these plots
0980         #First I loop through labels to put the hide button in twiki
0981         #All HGCal
0982         print("---+++ Results: Plots for individual material in all HGCal")
0983         for label, [num, color, leg] in hist_label_to_num.items():
0984             for p in ["x_vs_z_vs_Rsum", "l_vs_z_vs_Rsum", "x_vs_z_vs_Rsumcos", "l_vs_z_vs_Rsumcos", "x_vs_z_vs_Rloc", "l_vs_z_vs_Rloc"]:
0985                 TwikiPrintout(p, leg, "all")
0986         #Z+
0987         print("---+++ Results: Plots for individual material in Z+ Endcap")
0988         for label, [num, color, leg] in hist_label_to_num.items():
0989             for p in ["x_vs_z_vs_Rsum", "l_vs_z_vs_Rsum", "x_vs_z_vs_Rsumcos", "l_vs_z_vs_Rsumcos", "x_vs_z_vs_Rloc", "l_vs_z_vs_Rloc"]:
0990                 TwikiPrintout(p, leg, "zplus")
0991         #Z-
0992         print("---+++ Results: Plots for individual material in Z- Endcap")
0993         for label, [num, color, leg] in hist_label_to_num.items():
0994             for p in ["x_vs_z_vs_Rsum", "l_vs_z_vs_Rsum", "x_vs_z_vs_Rsumcos", "l_vs_z_vs_Rsumcos", "x_vs_z_vs_Rloc", "l_vs_z_vs_Rloc"]:
0995                 TwikiPrintout(p, leg, "zminus")
0996
0997         #Below is the silicon position from the xml geometry file
0998         #Should we put them on top of plots like the eta values?
0999         print(GetSiliconZValuesFromXML())
1000
1001         for p in required_2Dplots:
1002             #First the total
1003             create2DPlots(args.detector, p, plots[p].plotNumber, "")
1004             #Then, the rest
1005             for label, [num, color, leg] in hist_label_to_num.items():
1006                 #print label, num, color, leg
1007                 create2DPlots(args.detector, p, num + plots[p].plotNumber, leg)
1008
1009
1010         for p in required_plots:
1011             createCompoundPlots(args.detector, p)
1012         for p in required_ratio_plots:
1013             createRatioPlots(args.detector, p)