Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:27:23 

0001 ///////////////////////////////////////////////////////////////////////////////
0002 //
0003 // etaPhiPlot(fileName, plot, ifirst, ilast, drawLeg, ifEta, maxEta, tag, debug)
0004 //      Make the plots of integrated interaction/radiation/step lengths as a
0005 //      function of eta or phi
0006 // fileName (TString)  Name of the input ROOT file ("matbdg_HCAL.root")
0007 // plot     (TString)  Type of plot: IntLen/RadLen/StepLen ("IntLen")
0008 // ifirst   (int)      First depth to be displayed (0)
0009 // ilast    (int)      Last depth to be displayed (21)
0010 // drawLeg  (int)      Flag to show the legend or not (1)
0011 // ifEta    (bool)     Draw as a function of eta or phi (true)
0012 // maxEta   (double)   Maximum value of x-axis: if -1 use default (-1)
0013 // tag      (string)   Tag to be added to the name of the canvas ("Run")
0014 // debug    (bool)     print or not the debug information (false)
0015 //
0016 // etaPhiPlotComp(fileName1, fileName2, plot, lay, ifEta, maxEta, tag1,
0017 //                tag2, debug)
0018 //      Make superimposed plors of integrated interaction/radiation/step
0019 //      lengths from 2 files produced through 2 sources of geometry as a
0020 //      function of eta or phi
0021 //
0022 // etaPhiPlotHO(fileName, plot, drawLeg, ifEta, maxEta, tag)
0023 //      Make the same plots as *etaPhiPlot* for HO with same parameter meanings
0024 //
0025 // etaPhiPlotEC(fileName, plot, drawLeg, ifEta, maxEta, tag)
0026 //      Make the same plots as *etaPhiPlot* till front of HCAL with same
0027 //      parameter meanings
0028 //
0029 // etaPhiPlotHC(fileName, plot, drawLeg, ifEta, maxEta, tag)
0030 //      Make the same plots as *etaPhiPlot* for HC with same parameter meanings
0031 //
0032 // etaPhi2DPlot(fileName, plot, ifirst, ilast, drawLeg, tag)
0033 //      Make the 2-D plots as a function of eta and phi with same parameter
0034 //      meanings as those of *etaPhiPlot*
0035 //
0036 // etaPhi2DPlot(nslice, kslice, fileName, plot, ifirst, ilast, drawLeg, tag)
0037 //      Plot depth slices from kslice to nslice with other parameters having
0038 //      the same meanings as those of *etaPhiPlot*
0039 //
0040 // etaSlicePlot(fileName, plot, ifirst, ilast, ietaRange, drawLeg, tag, debug);
0041 //      Plot phi distributions of integrated interaction/radiation/step
0042 //      length for a given ietaRange (0 -> 28; 1 --> 29; 2 --> 9)
0043 //
0044 // printTable(fileName, outputFileName, inputFileName)
0045 //      Print a table of integrated lengths (or difference wrt an earlier one
0046 //      if inputFileName is other than None) for each depth in outputFileName
0047 //      with information from fileName
0048 //
0049 // plotDiff(fileName, plot, tag)
0050 //      Plots number of interaction/radiation/step lengths as a function of
0051 //      depth with parameters having the same meaning as in *etaPhiPlot*
0052 //
0053 ///////////////////////////////////////////////////////////////////////////////
0054 
0055 // include files
0056 #include <iostream>
0057 #include <iomanip>
0058 #include <fstream>
0059 #include <cmath>
0060 #include <vector>
0061 
0062 #include "TCanvas.h"
0063 #include "TDirectory.h"
0064 #include "TFile.h"
0065 #include "TGraph.h"
0066 #include "TLegend.h"
0067 #include "TMultiGraph.h"
0068 #include "TPaveText.h"
0069 #include "TProfile.h"
0070 #include "TProfile2D.h"
0071 #include "TStyle.h"
0072 
0073 const int nlaymax = 25;
0074 const int nbinmax = 41;
0075 double mean[nlaymax][nbinmax], diff[nlaymax][nbinmax];
0076 
0077 double towLow[41] = {0.000, 0.087, 0.174, 0.261, 0.348, 0.435, 0.522, 0.609, 0.696, 0.783, 0.870, 0.957, 1.044, 1.131,
0078                      1.218, 1.305, 1.392, 1.479, 1.566, 1.653, 1.740, 1.830, 1.930, 2.043, 2.172, 2.322, 2.500, 2.650,
0079                      2.853, 2.964, 3.139, 3.314, 3.489, 3.664, 3.839, 4.013, 4.191, 4.363, 4.538, 4.716, 4.889};
0080 double towHigh[41] = {0.087, 0.174, 0.261, 0.348, 0.435, 0.522, 0.609, 0.696, 0.783, 0.870, 0.957, 1.044, 1.131, 1.218,
0081                       1.305, 1.392, 1.479, 1.566, 1.653, 1.740, 1.830, 1.930, 2.043, 2.172, 2.322, 2.500, 2.650, 3.000,
0082                       2.964, 3.139, 3.314, 3.489, 3.664, 3.839, 4.013, 4.191, 4.363, 4.538, 4.716, 4.889, 5.191};
0083 
0084 int colorLayer[25] = {2, 7, 9, 30, 34, 38, 14, 40, 41, 42, 45, 46, 8, 49, 37, 28, 4, 1, 48, 50, 3, 6, 5, 156, 159};
0085 
0086 void etaPhiPlot(TString fileName = "matbdg_HCAL.root",
0087                 TString plot = "IntLen",
0088                 int ifirst = 0,
0089                 int ilast = 21,
0090                 int drawLeg = 1,
0091                 bool ifEta = true,
0092                 double maxEta = -1,
0093                 std::string tag = "Run",
0094                 bool debug = false);
0095 void etaPhiPlotComp(TString fileName1 = "matbdgHCAL_run3.root",
0096                     TString fileName2 = "matbdgHCAL_run3_dd4hep.root",
0097                     std::string plot = "IntLen",
0098                     int lay = 2,
0099                     bool ifEta = true,
0100                     double maxEta = -1,
0101                     std::string tag1 = "DDD",
0102                     std::string tag2 = "DD4hep",
0103                     bool debug = false);
0104 void etaPhiPlotHO(TString fileName = "matbdg_HCAL.root",
0105                   TString plot = "IntLen",
0106                   int drawLeg = 1,
0107                   bool ifEta = true,
0108                   double maxEta = -1,
0109                   std::string tag = "Run");
0110 void etaPhiPlotEC(TString fileName = "matbdg_HCAL.root",
0111                   TString plot = "IntLen",
0112                   int drawLeg = 1,
0113                   bool ifEta = true,
0114                   double maxEta = -1,
0115                   std::string tag = "Run");
0116 void etaPhiPlotHC(TString fileName = "matbdg_HCAL.root",
0117                   TString plot = "IntLen",
0118                   int drawLeg = 1,
0119                   bool ifEta = true,
0120                   double maxEta = -1,
0121                   std::string tag = "Run");
0122 void etaPhi2DPlot(TString fileName = "matbdg_HCAL.root",
0123                   TString plot = "IntLen",
0124                   int ifirst = 0,
0125                   int ilast = 19,
0126                   int drawLeg = 1,
0127                   std::string tag = "Run");
0128 void etaSlicePlot(TString fileName = "matbdg_HCAL.root",
0129                   TString plot = "IntLen",
0130                   int ifirst = 0,
0131                   int ilast = 19,
0132                   int ietaRange = 0,
0133                   int drawLeg = 1,
0134                   std::string tag = "Run",
0135                   bool debug = false);
0136 void etaPhi2DPlot(int nslice,
0137                   int kslice,
0138                   TString fileName = "matbdg_HCAL.root",
0139                   TString plot = "IntLen",
0140                   int ifirst = 0,
0141                   int ilast = 21,
0142                   int drawLeg = 1,
0143                   std::string tag = "Run");
0144 void printTable(TString fileName = "matbdg_HCAL.root",
0145                 TString outputFileName = "hcal.txt",
0146                 TString inputFileName = "None");
0147 void plotDiff(TString fileName = "matbdg_HCAL.root", TString plot = "IntLen", std::string tag = "Run");
0148 void getDiff(TString fileName = "matbdg_HCAL.root", TString plot = "IntLen");
0149 void plotHE(int flag = 0, int logy = 0, std::string tag = "Run", int save = 0);
0150 void etaPhiCastorPlot(TString fileName = "matbdg_Castor.root",
0151                       TString plot = "IntLen",
0152                       TString type = "All",
0153                       bool etaPlus = true,
0154                       int drawLeg = 1,
0155                       bool ifEta = true,
0156                       std::string tag = "Run",
0157                       bool debug = false);
0158 void efficiencyPlot(TString fileName = "matbdg_HCAL.root",
0159                     TString type = "All",
0160                     bool ifEtaPhi = true,
0161                     double maxEta = -1,
0162                     std::string tag = "Run",
0163                     bool debug = false);
0164 void etaPhiFwdPlot(TString fileName = "matbdg_Fwd.root",
0165                    TString plot = "IntLen",
0166                    int first = 0,
0167                    int last = 9,
0168                    int drawLeg = 1,
0169                    std::string tag = "Run",
0170                    bool debug = false);
0171 void setStyle();
0172 
0173 void standardPlot(TString fileName = "matbdgHCAL_run3.root",
0174                   std::string tag = "Run3",
0175                   int maxEta = 4.8,
0176                   TString outputFileName = "hcal.txt") {
0177   etaPhiPlot(fileName, "IntLen", 0, 21, 1, true, maxEta, tag);
0178   etaPhiPlot(fileName, "IntLen", 0, 20, 1, false, -1.0, tag);
0179   etaPhiPlot(fileName, "RadLen", 0, 21, 1, true, maxEta, tag);
0180   etaPhiPlot(fileName, "RadLen", 0, 20, 1, false, -1.0, tag);
0181   etaPhiPlot(fileName, "StepLen", 0, 21, 1, true, -1.0, tag);
0182   etaPhiPlot(fileName, "StepLen", 0, 20, 1, false, -1.0, tag);
0183   etaPhiPlotHO(fileName, "IntLen", 1, true, 2.5, tag);
0184   etaPhiPlotEC(fileName, "IntLen", 1, true, 2.5, tag);
0185   plotDiff(fileName, "IntLen", tag);
0186   plotDiff(fileName, "RadLen", tag);
0187   printTable(fileName, outputFileName);
0188   etaPhi2DPlot(fileName, "IntLen", 0, 21, 1, tag);
0189   etaPhi2DPlot(fileName, "RadLen", 0, 21, 1, tag);
0190 }
0191 
0192 void etaPhiPlot(TString fileName,
0193                 TString plot,
0194                 int ifirst,
0195                 int ilast,
0196                 int drawLeg,
0197                 bool ifEta,
0198                 double maxEta,
0199                 std::string tag,
0200                 bool debug) {
0201   TFile *hcalFile = new TFile(fileName);
0202   hcalFile->cd("g4SimHits");
0203   setStyle();
0204 
0205   TString xtit = TString("#eta");
0206   TString ztit = TString("eta");
0207   TString ytit = "none";
0208   int ymin = 0, ymax = 20, istart = 200;
0209   double xh = 0.90;
0210   if (plot.CompareTo("RadLen") == 0) {
0211     ytit = TString("HCal Material Budget (X_{0})");
0212     ymin = 0;
0213     ymax = 200;
0214     istart = 100;
0215   } else if (plot.CompareTo("StepLen") == 0) {
0216     ytit = TString("HCal Material Budget (Step Length)");
0217     ymin = 0;
0218     ymax = 15000;
0219     istart = 300;
0220     xh = 0.70;
0221   } else {
0222     ytit = TString("HCal Material Budget (#lambda)");
0223     ymin = 0;
0224     ymax = 20;
0225     istart = 200;
0226   }
0227   if (!ifEta) {
0228     istart += 400;
0229     xtit = TString("#phi");
0230     ztit = TString("phi");
0231   }
0232 
0233   TLegend *leg = new TLegend(xh - 0.13, 0.60, xh, 0.90);
0234   leg->SetBorderSize(1);
0235   leg->SetFillColor(10);
0236   leg->SetMargin(0.25);
0237   leg->SetTextSize(0.018);
0238 
0239   int nplots = 0;
0240   TProfile *prof[nlaymax];
0241   for (int ii = ilast; ii >= ifirst; ii--) {
0242     char hname[10], title[50];
0243     sprintf(hname, "%i", istart + ii);
0244     gDirectory->GetObject(hname, prof[nplots]);
0245     prof[nplots]->GetXaxis()->SetTitle(xtit);
0246     prof[nplots]->GetYaxis()->SetTitle(ytit);
0247     prof[nplots]->GetYaxis()->SetRangeUser(ymin, ymax);
0248     prof[nplots]->SetLineColor(colorLayer[ii]);
0249     prof[nplots]->SetFillColor(colorLayer[ii]);
0250     if (ifEta && maxEta > 0)
0251       prof[nplots]->GetXaxis()->SetRangeUser(-maxEta, maxEta);
0252     if (xh < 0.8)
0253       prof[nplots]->GetYaxis()->SetTitleOffset(1.7);
0254     else
0255       prof[nplots]->GetYaxis()->SetTitleOffset(1.0);
0256     int lay = ii - 1;
0257     if (lay > 0 && lay < 20) {
0258       sprintf(title, "Layer %d", lay);
0259     } else if (lay == 0) {
0260       sprintf(title, "After Crystal");
0261     } else if (lay >= 20) {
0262       sprintf(title, "After HF");
0263     } else {
0264       sprintf(title, "Before Crystal");
0265     }
0266     leg->AddEntry(prof[nplots], title, "lf");
0267     nplots++;
0268     if (ii == ilast && debug) {
0269       int nbinX = prof[0]->GetNbinsX();
0270       double xmin = prof[0]->GetXaxis()->GetXmin();
0271       double xmax = prof[0]->GetXaxis()->GetXmax();
0272       double dx = (xmax - xmin) / nbinX;
0273       std::cout << "Hist " << ii;
0274       for (int ibx = 0; ibx < nbinX; ibx++) {
0275         double xx1 = xmin + ibx * dx;
0276         double cont = prof[0]->GetBinContent(ibx + 1);
0277         std::cout << " | " << ibx << "(" << xx1 << ":" << (xx1 + dx) << ") " << cont;
0278       }
0279       std::cout << "\n";
0280     }
0281   }
0282 
0283   TString cname = "c_" + plot + ztit + tag;
0284   TCanvas *cc1 = new TCanvas(cname, cname, 700, 600);
0285   if (xh < 0.8) {
0286     cc1->SetLeftMargin(0.15);
0287     cc1->SetRightMargin(0.05);
0288   }
0289 
0290   prof[0]->Draw("h");
0291   for (int i = 1; i < nplots; i++)
0292     prof[i]->Draw("h sames");
0293   if (drawLeg > 0)
0294     leg->Draw("sames");
0295 }
0296 
0297 void etaPhiPlotComp(TString fileName1,
0298                     TString fileName2,
0299                     std::string plot,
0300                     int lay,
0301                     bool ifEta,
0302                     double maxEta,
0303                     std::string tag1,
0304                     std::string tag2,
0305                     bool debug) {
0306   setStyle();
0307   TFile *file1 = new TFile(fileName1);
0308   TFile *file2 = new TFile(fileName2);
0309   if (lay < 0)
0310     lay = 0;
0311   else if (lay > 21)
0312     lay = 21;
0313   if ((file1 != nullptr) && (file2 != nullptr)) {
0314     TDirectory *dir1 = (TDirectory *)(file1->FindObjectAny("g4SimHits"));
0315     TDirectory *dir2 = (TDirectory *)(file2->FindObjectAny("g4SimHits"));
0316 
0317     TString xtit = TString("#eta");
0318     std::string ztit("eta");
0319     std::string ytit("none");
0320     int ymin = 0, ymax = 20, istart = 200;
0321     double xh = 0.90;
0322     std::string plot1("L");
0323     if (plot == "RadLen") {
0324       ytit = "Radiation Length";
0325       ymin = 0;
0326       ymax = 200;
0327       istart = 100;
0328       plot1 = "X";
0329     } else if (plot == "StepLen") {
0330       ytit = "Step Length";
0331       ymin = 0;
0332       ymax = 15000;
0333       istart = 300;
0334       xh = 0.70;
0335       plot1 = "S";
0336     } else {
0337       ytit = "Interaction Length";
0338       ymin = 0;
0339       ymax = 20;
0340       istart = 200;
0341     }
0342     if (!ifEta) {
0343       istart += 400;
0344       xtit = TString("#phi");
0345       ztit = "phi";
0346     }
0347     double fac[23] = {0.05, 0.1, 0.3, 0.3, 0.3, 0.3, 0.3, 0.4, 0.4, 0.4, 0.5, 0.5,
0348                       0.5,  0.5, 0.6, 0.6, 0.7, 0.8, 0.8, 0.9, 1.0, 1.0, 1.0};
0349     ymax *= fac[lay];
0350 
0351     TLegend *leg = new TLegend(xh - 0.13, 0.80, xh, 0.90);
0352     leg->SetFillColor(kWhite);
0353     leg->SetBorderSize(1);
0354     leg->SetMargin(0.25);
0355 
0356     TProfile *prof[2];
0357     int color[2] = {2, 4};
0358     char hname[10], title[50], cname[100];
0359     if (lay > 1 && lay < 21) {
0360       sprintf(title, "%s (Layer %d)", ytit.c_str(), lay - 1);
0361     } else if (lay == 1) {
0362       sprintf(title, "%s (After Crystal)", ytit.c_str());
0363     } else if (lay > 20) {
0364       sprintf(title, "%s (After HF)", ytit.c_str());
0365     } else {
0366       sprintf(title, "%s (Before Crystal)", ytit.c_str());
0367     }
0368     sprintf(hname, "%i", istart + lay);
0369     dir1->GetObject(hname, prof[0]);
0370     dir2->GetObject(hname, prof[1]);
0371     if ((prof[0] != nullptr) && (prof[1] != nullptr)) {
0372       for (int k = 0; k < 2; ++k) {
0373         prof[k]->GetXaxis()->SetTitle(xtit);
0374         prof[k]->GetYaxis()->SetTitle(title);
0375         prof[k]->GetYaxis()->SetRangeUser(ymin, ymax);
0376         prof[k]->SetLineColor(color[k]);
0377         prof[k]->SetLineWidth(2);
0378         if (ifEta && maxEta > 0)
0379           prof[k]->GetXaxis()->SetRangeUser(-maxEta, maxEta);
0380         if (xh < 0.8)
0381           prof[k]->GetYaxis()->SetTitleOffset(1.7);
0382         else
0383           prof[k]->GetYaxis()->SetTitleOffset(1.0);
0384         if (k == 0)
0385           leg->AddEntry(prof[k], tag1.c_str(), "lf");
0386         else
0387           leg->AddEntry(prof[k], tag2.c_str(), "lf");
0388       }
0389 
0390       sprintf(cname, "c_%sLay%d%s%s%s", plot1.c_str(), lay, ztit.c_str(), tag1.c_str(), tag2.c_str());
0391       TCanvas *cc1 = new TCanvas(cname, cname, 700, 600);
0392       if (xh < 0.8) {
0393         cc1->SetLeftMargin(0.15);
0394         cc1->SetRightMargin(0.05);
0395       }
0396 
0397       prof[0]->Draw("h");
0398       prof[1]->Draw("h sames");
0399       leg->Draw("sames");
0400 
0401       if (debug) {
0402         int nbinX = prof[0]->GetNbinsX();
0403         double xmin = prof[0]->GetXaxis()->GetXmin();
0404         double xmax = prof[0]->GetXaxis()->GetXmax();
0405         double dx = (xmax - xmin) / nbinX;
0406         for (int ibx = 0; ibx < nbinX; ibx++) {
0407           double xx1 = xmin + ibx * dx;
0408           if (xx1 >= -1.5 && xx1 < 1.5) {
0409             double cont1 = prof[0]->GetBinContent(ibx + 1);
0410             double cont2 = prof[1]->GetBinContent(ibx + 1);
0411             if ((cont1 + cont2) > 0.01) {
0412               double frac = 2.0 * (cont1 - cont2) / (cont1 + cont2);
0413               std::cout << ibx << "(" << xx1 << ":" << (xx1 + dx) << ") " << cont1 << ":" << cont2 << ":" << frac
0414                         << std::endl;
0415             }
0416           }
0417         }
0418       }
0419     }
0420   }
0421 }
0422 
0423 void etaPhiPlotHO(TString fileName, TString plot, int drawLeg, bool ifEta, double maxEta, std::string tag) {
0424   TFile *hcalFile = new TFile(fileName);
0425   hcalFile->cd("g4SimHits");
0426   setStyle();
0427 
0428   int ihid[3] = {2, 18, 20};
0429   TString xtit = TString("#eta");
0430   TString ztit = TString("eta");
0431   TString ytit = "none";
0432   int ymin = 0, ymax = 20, istart = 200;
0433   double xh = 0.90;
0434   if (plot.CompareTo("RadLen") == 0) {
0435     ytit = TString("Material Budget (X_{0})");
0436     ymin = 0;
0437     ymax = 200;
0438     istart = 100;
0439   } else if (plot.CompareTo("StepLen") == 0) {
0440     ytit = TString("Material Budget (Step Length)");
0441     ymin = 0;
0442     ymax = 15000;
0443     istart = 300;
0444     xh = 0.70;
0445   } else {
0446     ytit = TString("Material Budget (#lambda)");
0447     ymin = 0;
0448     ymax = 20;
0449     istart = 200;
0450   }
0451   if (!ifEta) {
0452     istart += 400;
0453     xtit = TString("#phi");
0454     ztit = TString("phi");
0455   }
0456 
0457   TLegend *leg = new TLegend(xh - 0.25, 0.60, xh, 0.90);
0458   leg->SetBorderSize(1);
0459   leg->SetFillColor(10);
0460   leg->SetMargin(0.45);
0461   leg->SetTextSize(0.04);
0462   int nplots = 0;
0463   TProfile *prof[nlaymax];
0464   for (int ii = 2; ii >= 0; ii--) {
0465     char hname[10], title[50];
0466     int idpl = istart + ihid[ii];
0467     sprintf(hname, "%i", idpl);
0468     gDirectory->GetObject(hname, prof[nplots]);
0469     prof[nplots]->GetXaxis()->SetTitle(xtit);
0470     prof[nplots]->GetYaxis()->SetTitle(ytit);
0471     prof[nplots]->GetXaxis()->SetLabelSize(0.05);
0472     prof[nplots]->GetYaxis()->SetLabelSize(0.05);
0473     prof[nplots]->GetXaxis()->SetTitleSize(0.05);
0474     prof[nplots]->GetYaxis()->SetTitleSize(0.05);
0475     prof[nplots]->GetYaxis()->SetTitleOffset(0.8);
0476     prof[nplots]->GetYaxis()->SetRangeUser(ymin, ymax);
0477     prof[nplots]->SetLineColor(colorLayer[ii]);
0478     prof[nplots]->SetFillColor(colorLayer[ii]);
0479     if (xh < 0.8)
0480       prof[nplots]->GetYaxis()->SetTitleOffset(1.7);
0481     if (ifEta && maxEta > 0)
0482       prof[nplots]->GetXaxis()->SetRangeUser(0.0, maxEta);
0483     if (ii >= 2) {
0484       sprintf(title, "With HO");
0485     } else if (ii == 1) {
0486       sprintf(title, "Without HO");
0487     } else {
0488       sprintf(title, "Before HCAL");
0489     }
0490     leg->AddEntry(prof[nplots], title, "lf");
0491     nplots++;
0492     if (ii == 1) {
0493       for (int kk = 0; kk < 7; kk++) {
0494         idpl--;
0495         sprintf(hname, "%i", idpl);
0496         gDirectory->GetObject(hname, prof[nplots]);
0497         prof[nplots]->GetXaxis()->SetTitle(xtit);
0498         prof[nplots]->GetYaxis()->SetTitle(ytit);
0499         prof[nplots]->GetYaxis()->SetRangeUser(ymin, ymax);
0500         prof[nplots]->SetLineColor(colorLayer[ii]);
0501         prof[nplots]->SetFillColor(colorLayer[ii]);
0502         if (ifEta && maxEta > 0)
0503           prof[nplots]->GetXaxis()->SetRangeUser(0.0, maxEta);
0504         nplots++;
0505       }
0506     }
0507   }
0508 
0509   TString cname = "c_HO" + plot + ztit + tag;
0510   new TCanvas(cname, cname, 700, 400);
0511 
0512   prof[0]->Draw("h");
0513   for (int i = 1; i < nplots; i++)
0514     prof[i]->Draw("h sames");
0515   if (drawLeg > 0)
0516     leg->Draw("sames");
0517 }
0518 
0519 void etaPhiPlotEC(TString fileName, TString plot, int drawLeg, bool ifEta, double maxEta, std::string tag) {
0520   TFile *hcalFile = new TFile(fileName);
0521   hcalFile->cd("g4SimHits");
0522   setStyle();
0523 
0524   int ihid[3] = {0, 1, 2};
0525   TString xtit = TString("#eta");
0526   TString ztit = TString("eta");
0527   TString ytit = "none";
0528   int ymin = 0, ymax = 7, istart = 200, ymax1 = 5;
0529   double xh = 0.90, xh1 = 0.90;
0530   if (plot.CompareTo("RadLen") == 0) {
0531     ytit = TString("Material Budget (X_{0})");
0532     ymin = 0;
0533     ymax = 70;
0534     istart = 100;
0535     ymax1 = 50;
0536   } else if (plot.CompareTo("StepLen") == 0) {
0537     ytit = TString("Material Budget (Step Length)");
0538     ymin = 0;
0539     ymax = 6000;
0540     istart = 300;
0541     xh = 0.35;
0542     ymax1 = 2500;
0543   } else {
0544     ytit = TString("Material Budget (#lambda)");
0545     ymin = 0;
0546     ymax = 7;
0547     istart = 200;
0548   }
0549   if (!ifEta) {
0550     istart += 400;
0551     xtit = TString("#phi");
0552     ztit = TString("phi");
0553   }
0554 
0555   TLegend *leg = new TLegend(xh - 0.25, 0.60, xh, 0.90);
0556   leg->SetBorderSize(1);
0557   leg->SetFillColor(10);
0558   leg->SetMargin(0.30);
0559   leg->SetTextSize(0.04);
0560   int nplots = 0;
0561   TProfile *prof[nlaymax];
0562   for (int ii = 2; ii >= 0; ii--) {
0563     char hname[10], title[50];
0564     int idpl = istart + ihid[ii];
0565     sprintf(hname, "%i", idpl);
0566     gDirectory->GetObject(hname, prof[nplots]);
0567     prof[nplots]->GetXaxis()->SetTitle(xtit);
0568     prof[nplots]->GetYaxis()->SetTitle(ytit);
0569     prof[nplots]->GetXaxis()->SetLabelSize(0.05);
0570     prof[nplots]->GetYaxis()->SetLabelSize(0.05);
0571     prof[nplots]->GetXaxis()->SetTitleSize(0.05);
0572     prof[nplots]->GetYaxis()->SetTitleSize(0.05);
0573     prof[nplots]->GetYaxis()->SetTitleOffset(0.8);
0574     prof[nplots]->GetYaxis()->SetRangeUser(ymin, ymax);
0575     prof[nplots]->SetLineColor(colorLayer[ii]);
0576     prof[nplots]->SetFillColor(colorLayer[ii]);
0577     if (xh < 0.8)
0578       prof[nplots]->GetYaxis()->SetTitleOffset(1.05);
0579     if (ifEta && maxEta > 0)
0580       prof[nplots]->GetXaxis()->SetRangeUser(0.0, maxEta);
0581     if (ii >= 2) {
0582       sprintf(title, "Front of HCAL");
0583     } else if (ii == 1) {
0584       sprintf(title, "After Crystals");
0585     } else {
0586       sprintf(title, "Before Crystals");
0587     }
0588     leg->AddEntry(prof[nplots], title, "lf");
0589     nplots++;
0590   }
0591 
0592   TString cname = "c_EC1" + plot + ztit + tag;
0593   new TCanvas(cname, cname, 700, 400);
0594 
0595   prof[0]->Draw("h");
0596   for (int i = 1; i < nplots; i++)
0597     prof[i]->Draw("h sames");
0598   if (drawLeg > 0)
0599     leg->Draw("sames");
0600 
0601   double xmin = prof[2]->GetXaxis()->GetXmin();
0602   double xmax = prof[2]->GetXaxis()->GetXmax();
0603   int nbins = prof[2]->GetNbinsX();
0604   TH1D *prof1 = new TH1D("Temp01", "Temp01", nbins, xmin, xmax);
0605   for (int ii = 0; ii < nbins; ii++) {
0606     double x1 = prof[0]->GetBinLowEdge(ii + 1);
0607     double x2 = prof[0]->GetBinWidth(ii + 1);
0608     double v0 = prof[0]->GetBinContent(ii + 1);
0609     double v1 = prof[1]->GetBinContent(ii + 1);
0610     double v2 = prof[2]->GetBinContent(ii + 1);
0611     double xx = x1 + 0.5 * x2;
0612     double cont = v0 - v1;
0613     prof1->Fill(xx, cont);
0614     std::cout << "Bin " << ii << " Eta/Phi " << std::setw(4) << xx << " Material " << std::setw(6) << cont << "  "
0615               << std::setw(6) << (v0 - v2) << "\n";
0616   }
0617   prof1->GetXaxis()->SetTitle(xtit);
0618   prof1->GetYaxis()->SetTitle(ytit);
0619   prof1->GetXaxis()->SetLabelSize(0.05);
0620   prof1->GetYaxis()->SetLabelSize(0.05);
0621   prof1->GetXaxis()->SetTitleSize(0.05);
0622   prof1->GetYaxis()->SetTitleSize(0.05);
0623   prof1->GetYaxis()->SetTitleOffset(0.8);
0624   prof1->GetYaxis()->SetRangeUser(ymin, ymax1);
0625   prof1->SetLineColor(colorLayer[3]);
0626   prof1->SetFillColor(colorLayer[3]);
0627   if (ifEta && maxEta > 0)
0628     prof1->GetXaxis()->SetRangeUser(0.0, maxEta);
0629   if (xh < 0.8)
0630     prof1->GetYaxis()->SetTitleOffset(1.05);
0631 
0632   TLegend *mlg = new TLegend(xh1 - 0.3, 0.80, xh1, 0.90);
0633   char title[100];
0634   sprintf(title, "End crystal to Layer 0");
0635   mlg->SetBorderSize(1);
0636   mlg->SetFillColor(10);
0637   mlg->SetMargin(0.30);
0638   mlg->SetTextSize(0.04);
0639   mlg->AddEntry(prof1, title, "lf");
0640 
0641   cname = "c_EC2" + plot + ztit + tag;
0642   new TCanvas(cname, cname, 700, 400);
0643   prof1->Draw();
0644   if (drawLeg > 0)
0645     mlg->Draw("sames");
0646 }
0647 
0648 void etaPhiPlotHC(TString fileName, TString plot, int drawLeg, bool ifEta, double maxEta, std::string tag) {
0649   TFile *hcalFile = new TFile(fileName);
0650   hcalFile->cd("g4SimHits");
0651   setStyle();
0652 
0653   int ihid[20] = {2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21};
0654   TString xtit = TString("#eta");
0655   TString ztit = TString("eta");
0656   TString ytit = "none";
0657   int ymin = 0, ymax = 20, istart = 200;
0658   double xh = 0.90;
0659   if (plot.CompareTo("RadLen") == 0) {
0660     ytit = TString("Material Budget (X_{0})");
0661     ymin = 0;
0662     ymax = 200;
0663     istart = 100;
0664   } else if (plot.CompareTo("StepLen") == 0) {
0665     ytit = TString("Material Budget (Step Length)");
0666     ymin = 0;
0667     ymax = 15000;
0668     istart = 300;
0669     xh = 0.35;
0670   } else {
0671     ytit = TString("Material Budget (#lambda)");
0672     ymin = 0;
0673     ymax = 20;
0674     istart = 200;
0675   }
0676   if (!ifEta) {
0677     istart += 400;
0678     xtit = TString("#phi");
0679     ztit = TString("phi");
0680   }
0681 
0682   TLegend *leg = new TLegend(xh - 0.25, 0.70, xh, 0.90);
0683   leg->SetBorderSize(1);
0684   leg->SetFillColor(10);
0685   leg->SetMargin(0.30);
0686   leg->SetTextSize(0.04);
0687   int nplots = 0;
0688   TProfile *prof[nlaymax];
0689   for (int ii = 19; ii >= 0; ii--) {
0690     char hname[10], title[50];
0691     int idpl = istart + ihid[ii];
0692     sprintf(hname, "%i", idpl);
0693     int icol = colorLayer[0];
0694     if (ii >= 18)
0695       icol = colorLayer[2];
0696     else if (ii > 0)
0697       icol = colorLayer[1];
0698     gDirectory->GetObject(hname, prof[nplots]);
0699     prof[nplots]->GetXaxis()->SetTitle(xtit);
0700     prof[nplots]->GetYaxis()->SetTitle(ytit);
0701     prof[nplots]->GetXaxis()->SetLabelSize(0.05);
0702     prof[nplots]->GetYaxis()->SetLabelSize(0.05);
0703     prof[nplots]->GetXaxis()->SetTitleSize(0.05);
0704     prof[nplots]->GetYaxis()->SetTitleSize(0.05);
0705     prof[nplots]->GetYaxis()->SetTitleOffset(0.8);
0706     prof[nplots]->GetYaxis()->SetRangeUser(ymin, ymax);
0707     prof[nplots]->SetLineColor(icol);
0708     prof[nplots]->SetFillColor(icol);
0709     if (xh < 0.8)
0710       prof[nplots]->GetYaxis()->SetTitleOffset(1.05);
0711     if (ifEta && maxEta > 0)
0712       prof[nplots]->GetXaxis()->SetRangeUser(0.0, maxEta);
0713     if (ii == 19) {
0714       sprintf(title, "End of HF");
0715       leg->AddEntry(prof[nplots], title, "lf");
0716     } else if (ii == 1) {
0717       sprintf(title, "After HB/HE/HO");
0718       leg->AddEntry(prof[nplots], title, "lf");
0719     } else if (ii == 0) {
0720       sprintf(title, "Before HCAL");
0721       leg->AddEntry(prof[nplots], title, "lf");
0722     }
0723     nplots++;
0724   }
0725 
0726   TString cname = "c_HC" + plot + ztit + tag;
0727   new TCanvas(cname, cname, 700, 400);
0728 
0729   prof[0]->Draw("h");
0730   for (int i = 1; i < nplots; i++)
0731     prof[i]->Draw("h sames");
0732   if (drawLeg > 0)
0733     leg->Draw("sames");
0734 }
0735 
0736 void etaPhi2DPlot(TString fileName, TString plot, int ifirst, int ilast, int drawLeg, std::string tag) {
0737   TFile *hcalFile = new TFile(fileName);
0738   hcalFile->cd("g4SimHits");
0739   setStyle();
0740 
0741   TString xtit = TString("#eta");
0742   TString ytit = TString("#phi");
0743   TString xytit = TString("etaphi");
0744   TString ztit = TString("HCal Material Budget (#lambda)");
0745   int ymin = 0, ymax = 20, istart = 1000;
0746   double xh = 0.95, yh = 0.95;
0747   if (plot.CompareTo("RadLen") == 0) {
0748     ztit = TString("HCal Material Budget (X_{0})");
0749     ymin = 0;
0750     ymax = 200;
0751     istart = 900;
0752   } else if (plot.CompareTo("StepLen") == 0) {
0753     ztit = TString("HCal Material Budget (Step Length)");
0754     ymin = 0;
0755     ymax = 15000;
0756     istart = 1100;
0757   }
0758 
0759   TLegend *leg = new TLegend(xh - 0.13, yh - 0.30, xh, yh);
0760   leg->SetBorderSize(1);
0761   leg->SetFillColor(10);
0762   leg->SetMargin(0.25);
0763   leg->SetTextSize(0.018);
0764 
0765   int nplots = 0;
0766   TProfile2D *prof[nlaymax];
0767   for (int ii = ilast; ii >= ifirst; ii--) {
0768     char hname[10], title[50];
0769     sprintf(hname, "%i", istart + ii);
0770     gDirectory->GetObject(hname, prof[nplots]);
0771     prof[nplots]->GetXaxis()->SetTitle(xtit);
0772     prof[nplots]->GetYaxis()->SetTitle(ytit);
0773     prof[nplots]->GetZaxis()->SetTitle(ztit);
0774     prof[nplots]->GetZaxis()->SetRangeUser(ymin, ymax);
0775     prof[nplots]->SetLineColor(colorLayer[ii]);
0776     prof[nplots]->SetFillColor(colorLayer[ii]);
0777     prof[nplots]->GetZaxis()->SetTitleOffset(1.4);
0778     int lay = ii - 1;
0779     if (lay > 0 && lay < 20) {
0780       sprintf(title, "Layer %d", lay);
0781     } else if (lay == 0) {
0782       sprintf(title, "After Crystal");
0783     } else if (lay >= 20) {
0784       sprintf(title, "After HF");
0785     } else {
0786       sprintf(title, "Before Crystal");
0787     }
0788     leg->AddEntry(prof[nplots], title, "lf");
0789     nplots++;
0790   }
0791 
0792   TString cname = "c_" + plot + xytit + tag;
0793   TCanvas *cc1 = new TCanvas(cname, cname, 700, 600);
0794   cc1->SetLeftMargin(0.15);
0795   cc1->SetRightMargin(0.05);
0796 
0797   prof[0]->Draw("lego fb bb");
0798   for (int i = 1; i < nplots; i++)
0799     prof[i]->Draw("lego fb bb sames");
0800   if (drawLeg > 0)
0801     leg->Draw("sames");
0802 }
0803 
0804 void etaPhi2DPlot(
0805     int nslice, int kslice, TString fileName, TString plot, int ifirst, int ilast, int drawLeg, std::string tag) {
0806   char hname[200], title[50];
0807 
0808   TFile *hcalFile = new TFile(fileName);
0809   hcalFile->cd("g4SimHits");
0810   setStyle();
0811 
0812   TString xtit = TString("#eta");
0813   TString ytit;
0814   int ymin, ymax, istart;
0815   double xh = 0.95, yh = 0.90;
0816   char type[10];
0817   if (plot.CompareTo("RadLen") == 0) {
0818     ytit = TString("HCal Material Budget (X_{0})");
0819     ymin = 0;
0820     ymax = 200;
0821     istart = 900;
0822     sprintf(type, "Radlen");
0823   } else if (plot.CompareTo("StepLen") == 0) {
0824     ytit = TString("HCal Material Budget (Step Length)");
0825     ymin = 0;
0826     ymax = 15000;
0827     istart = 1100;
0828     sprintf(type, "Steplen");
0829   } else {
0830     ytit = TString("HCal Material Budget (#lambda)");
0831     ymin = 0;
0832     ymax = 20;
0833     istart = 1000;
0834     sprintf(type, "Intlen");
0835   }
0836 
0837   int nplots = 0;
0838   TProfile2D *prof[nlaymax];
0839   for (int ii = ilast; ii >= ifirst; ii--) {
0840     sprintf(hname, "%i", istart + ii);
0841     gDirectory->GetObject(hname, prof[nplots]);
0842     nplots++;
0843   }
0844 
0845   double xmin = prof[0]->GetXaxis()->GetXmin();
0846   double xmax = prof[0]->GetXaxis()->GetXmax();
0847   int nbinX = prof[0]->GetNbinsX();
0848   double ymin1 = prof[0]->GetYaxis()->GetXmin();
0849   double ymax1 = prof[0]->GetYaxis()->GetXmax();
0850   int nbinY = prof[0]->GetNbinsY();
0851   int ngroup = nbinY / nslice;
0852   double dy = (ymax1 - ymin1) / nbinY;
0853   cout << "X " << nbinX << "/" << xmin << "/" << xmax << " Slice " << nbinY << "/" << nslice << "/" << ngroup << " "
0854        << nplots << " " << dy << "\n";
0855 
0856   istart = 0;
0857   TLegend *leg[360];
0858   TH1D *hist[nlaymax][360];
0859   for (int is = 0; is < nslice; is++) {
0860     leg[is] = new TLegend(xh - 0.13, yh - 0.43, xh, yh);
0861     leg[is]->SetBorderSize(1);
0862     leg[is]->SetFillColor(10);
0863     leg[is]->SetMargin(0.25);
0864     leg[is]->SetTextSize(0.023);
0865     double y1 = (ymin1 + istart * dy) * 180. / 3.1415926;
0866     double y2 = y1 + ngroup * dy * 180. / 3.1415926;
0867     if (y1 < 0) {
0868       y1 += 360.;
0869       y2 += 360.;
0870     }
0871     sprintf(title, "#phi = %6.1f :%6.1f", y1, y2);
0872     leg[is]->SetHeader(title);
0873 
0874     for (int ii = 0; ii < nplots; ii++) {
0875       sprintf(hname, "Hist%iSlice%i", ii, is);
0876       hist[ii][is] = new TH1D(hname, hname, nbinX, xmin, xmax);
0877       //      cout << "Hist " << ii;
0878       for (int ibx = 0; ibx < nbinX; ibx++) {
0879         double contb = 0;
0880         //  cout << " / " << ibx;
0881         for (int iby = 0; iby < ngroup; iby++) {
0882           int ibin = iby + istart;
0883           double cont = prof[ii]->GetBinContent(ibx + 1, ibin + 1);
0884           //      cout << "/" << ibin << " " << cont;
0885           contb += cont;
0886         }
0887         contb /= ngroup;
0888         //  cout << " " << contb;
0889         hist[ii][is]->SetBinContent(ibx + 1, contb);
0890       }
0891       //      cout << "\n";
0892       hist[ii][is]->GetXaxis()->SetTitle(xtit);
0893       hist[ii][is]->GetYaxis()->SetTitle(ytit);
0894       hist[ii][is]->GetYaxis()->SetRangeUser(ymin, ymax);
0895       hist[ii][is]->SetLineColor(colorLayer[ilast - ii]);
0896       hist[ii][is]->SetFillColor(colorLayer[ilast - ii]);
0897       hist[ii][is]->GetYaxis()->SetTitleOffset(0.8);
0898       int lay = ilast - ii - 1;
0899       if (lay > 0 && lay < 20) {
0900         sprintf(title, "Layer %d", lay);
0901       } else if (lay == 0) {
0902         sprintf(title, "After Crystal");
0903       } else if (lay >= 20) {
0904         sprintf(title, "After HF");
0905       } else {
0906         sprintf(title, "Before Crystal");
0907       }
0908       leg[is]->AddEntry(hist[ii][is], title, "lf");
0909     }
0910     istart += ngroup;
0911   }
0912 
0913   cout << "All histograms created now plot\n";
0914   TCanvas *cc1[360];
0915   int ismin = kslice, ismax = nslice;
0916   for (int is = ismin; is < ismax; is++) {
0917     sprintf(hname, "c_%s%i%s", type, is, tag.c_str());
0918     cc1[is] = new TCanvas(hname, hname, 700, 400);
0919     cc1[is]->SetLeftMargin(0.15);
0920     cc1[is]->SetRightMargin(0.05);
0921     hist[0][is]->Draw();
0922     for (int i = 1; i < nplots; i++)
0923       hist[i][is]->Draw("sames");
0924     if (drawLeg > 0)
0925       leg[is]->Draw("sames");
0926   }
0927 }
0928 
0929 void etaSlicePlot(
0930     TString fileName, TString plot, int ifirst, int ilast, int ietaRange, int drawLeg, std::string tag, bool debug) {
0931   TFile *hcalFile = new TFile(fileName);
0932   hcalFile->cd("g4SimHits");
0933   setStyle();
0934 
0935   TString xtit = TString("#phi (i#eta = 28)");
0936   TString ztit = TString("phi");
0937   TString ytit = "none";
0938   int ymin = 0, ymax = 20, istart = 1500;
0939   double xh = 0.90;
0940   if (plot.CompareTo("RadLen") == 0) {
0941     ytit = TString("HCal Material Budget (X_{0})");
0942     ymin = 0;
0943     ymax = 200;
0944     istart = 1600;
0945   } else if (plot.CompareTo("StepLen") == 0) {
0946     ytit = TString("HCal Material Budget (Step Length)");
0947     ymin = 0;
0948     ymax = 15000;
0949     istart = 1800;
0950     xh = 0.70;
0951   } else {
0952     ytit = TString("HCal Material Budget (#lambda)");
0953     ymin = 0;
0954     ymax = 18;
0955     istart = 1700;
0956   }
0957   if (ietaRange == 1) {
0958     istart += 300;
0959     xtit = TString("#phi (i#eta = 29)");
0960   } else if (ietaRange == 2) {
0961     istart += 600;
0962     xtit = TString("#phi (i#eta = 9)");
0963   }
0964 
0965   TLegend *leg = new TLegend(xh - 0.13, 0.60, xh, 0.90);
0966   leg->SetBorderSize(1);
0967   leg->SetFillColor(10);
0968   leg->SetMargin(0.25);
0969   leg->SetTextSize(0.018);
0970 
0971   int nplots = 0;
0972   TProfile *prof[nlaymax];
0973   for (int ii = ilast; ii >= ifirst; ii--) {
0974     char hname[10], title[50];
0975     sprintf(hname, "%i", istart + ii);
0976     gDirectory->GetObject(hname, prof[nplots]);
0977     prof[nplots]->GetXaxis()->SetTitle(xtit);
0978     prof[nplots]->GetYaxis()->SetTitle(ytit);
0979     prof[nplots]->GetYaxis()->SetRangeUser(ymin, ymax);
0980     prof[nplots]->SetLineColor(colorLayer[ii]);
0981     prof[nplots]->SetFillColor(colorLayer[ii]);
0982     if (xh < 0.8)
0983       prof[nplots]->GetYaxis()->SetTitleOffset(1.7);
0984     else
0985       prof[nplots]->GetYaxis()->SetTitleOffset(1.0);
0986     int lay = ii - 1;
0987     if (lay > 0 && lay < 20) {
0988       sprintf(title, "Layer %d", lay);
0989     } else if (lay == 0) {
0990       sprintf(title, "After Crystal");
0991     } else if (lay >= 20) {
0992       sprintf(title, "After HF");
0993     } else {
0994       sprintf(title, "Before Crystal");
0995     }
0996     leg->AddEntry(prof[nplots], title, "lf");
0997     nplots++;
0998     if (ii == ilast && debug) {
0999       int nbinX = prof[0]->GetNbinsX();
1000       double xmin = prof[0]->GetXaxis()->GetXmin();
1001       double xmax = prof[0]->GetXaxis()->GetXmax();
1002       double dx = (xmax - xmin) / nbinX;
1003       cout << "Hist " << ii;
1004       for (int ibx = 0; ibx < nbinX; ibx++) {
1005         double xx1 = xmin + ibx * dx;
1006         double cont = prof[0]->GetBinContent(ibx + 1);
1007         cout << " | " << ibx << "(" << xx1 << ":" << (xx1 + dx) << ") " << cont;
1008       }
1009       cout << "\n";
1010     }
1011   }
1012 
1013   TString cname = "c_" + plot + ztit + tag;
1014   TCanvas *cc1 = new TCanvas(cname, cname, 700, 600);
1015   if (xh < 0.8) {
1016     cc1->SetLeftMargin(0.15);
1017     cc1->SetRightMargin(0.05);
1018   }
1019 
1020   prof[0]->Draw("h");
1021   for (int i = 1; i < nplots; i++)
1022     prof[i]->Draw("h sames");
1023   if (drawLeg > 0)
1024     leg->Draw("sames");
1025 }
1026 
1027 void printTable(TString fileName, TString outputFileName, TString inputFileName) {
1028   double radl[nlaymax][nbinmax], intl[nlaymax][nbinmax];
1029   bool compare = false;
1030   if (inputFileName != "None") {
1031     ifstream inp(inputFileName, ios::in);
1032     cout << "Opens " << inputFileName << "\n";
1033     if (inp) {
1034       TString line;
1035       int tower;
1036       double eta;
1037       for (int i = 0; i < 23; i++)
1038         inp >> line;
1039       for (int itow = 0; itow < nbinmax; itow++) {
1040         inp >> tower >> eta;
1041         int laymax = 19;
1042         if (itow > 29)
1043           laymax = 2;
1044         else if (itow > 3)
1045           laymax = 18;
1046         for (int ilay = 0; ilay < laymax; ilay++)
1047           inp >> intl[ilay][tower];
1048       }
1049       for (int i = 0; i < 23; i++)
1050         inp >> line;
1051       for (int itow = 0; itow < nbinmax; itow++) {
1052         inp >> tower >> eta;
1053         int laymax = 19;
1054         if (itow > 29)
1055           laymax = 2;
1056         else if (itow > 3)
1057           laymax = 18;
1058         for (int ilay = 0; ilay < laymax; ilay++)
1059           inp >> radl[ilay][tower];
1060       }
1061       compare = true;
1062       inp.close();
1063     }
1064   }
1065   std::ofstream os;
1066   os.open(outputFileName);
1067 
1068   int nbadI = 0;
1069   getDiff(fileName, "IntLen");
1070   os << "Interaction Length\n"
1071      << "==================\n"
1072      << "Eta Tower/Layer   0      1       2       3       4       5     "
1073      << "  6       7       8       9     10      11      12      13     "
1074      << " 14      15      16      17\n";
1075   for (int itow = 0; itow < nbinmax; itow++) {
1076     os << setw(3) << itow << setw(7) << setprecision(3) << 0.5 * (towLow[itow] + towHigh[itow]);
1077     int laymax = 19, ioff = 1;
1078     if (itow > 29) {
1079       laymax = 2;
1080       ioff = 0;
1081     } else if (itow > 3) {
1082       laymax = 18;
1083       ioff = 1;
1084     }
1085     for (int ilay = 0; ilay < laymax; ilay++) {
1086       os << setw(8) << setprecision(4) << diff[ilay + ioff][itow];
1087       if (compare) {
1088         double num = (diff[ilay + ioff][itow] - intl[ilay][itow]);
1089         double den = (diff[ilay + ioff][itow] + intl[ilay][itow]);
1090         double dd = (den == 0. ? 0. : 2.0 * num / den);
1091         if (dd > 0.01) {
1092           nbadI++;
1093           cout << "Lambda::Tower " << setw(3) << itow << " Layer " << setw(3) << ilay << " Old" << setw(8)
1094                << setprecision(4) << intl[ilay][itow] << " New" << setw(8) << setprecision(4) << diff[ilay + ioff][itow]
1095                << " Diff" << setw(8) << setprecision(4) << dd << "\n";
1096         }
1097       }
1098     }
1099     os << "\n";
1100   }
1101 
1102   int nbadR = 0;
1103   getDiff(fileName, "RadLen");
1104   os << "\n\nRadiation Length\n"
1105      << "================\n"
1106      << "Eta Tower/Layer   0      1       2       3       4       5     "
1107      << "  6       7       8       9     10      11      12      13     "
1108      << " 14      15      16      17\n";
1109   for (int itow = 0; itow < nbinmax; itow++) {
1110     os << setw(3) << itow << setw(7) << setprecision(3) << 0.5 * (towLow[itow] + towHigh[itow]);
1111     int laymax = 19, ioff = 1;
1112     if (itow > 29) {
1113       laymax = 2;
1114       ioff = 0;
1115     } else if (itow > 3) {
1116       laymax = 18;
1117       ioff = 1;
1118     }
1119     for (int ilay = 0; ilay < laymax; ilay++) {
1120       os << setw(8) << setprecision(4) << diff[ilay + ioff][itow];
1121       if (compare) {
1122         double num = (diff[ilay + ioff][itow] - radl[ilay][itow]);
1123         double den = (diff[ilay + ioff][itow] + radl[ilay][itow]);
1124         double dd = (den == 0. ? 0. : 2.0 * num / den);
1125         if (dd > 0.01) {
1126           nbadR++;
1127           cout << "X0::Tower " << setw(3) << itow << " Layer " << setw(3) << ilay << " Old" << setw(8)
1128                << setprecision(4) << radl[ilay][itow] << " New" << setw(8) << setprecision(4) << diff[ilay + ioff][itow]
1129                << " Diff" << setw(8) << setprecision(4) << dd << "\n";
1130         }
1131       }
1132     }
1133     os << "\n";
1134   }
1135   os.close();
1136 
1137   cout << "Comparison Results " << nbadI << " discrepancies for Lambda and " << nbadR << " discrepancies for X0\n";
1138 }
1139 
1140 void plotDiff(TString fileName, TString plot, std::string tag) {
1141   setStyle();
1142   gStyle->SetTitleOffset(1.0, "Y");
1143   getDiff(fileName, plot);
1144   TString xtit = TString("Layer Number");
1145   TString ytit = TString("HCal Material Budget (#lambda)");
1146   if (plot.CompareTo("RadLen") == 0)
1147     ytit = TString("HCal Material Budget (X_{0})");
1148 
1149   TMultiGraph *mg = new TMultiGraph();
1150   TLegend *leg_mg = new TLegend(.5, .5, .75, .80);
1151   leg_mg->SetFillColor(10);
1152   leg_mg->SetBorderSize(1);
1153   leg_mg->SetMargin(0.25);
1154   leg_mg->SetTextSize(0.04);
1155   leg_mg->SetHeader(ytit);
1156 
1157   double diff_lay[19], idx[19];
1158   for (int ilay = 1; ilay < 20; ilay++) {
1159     diff_lay[ilay - 1] = diff[ilay][0];
1160     idx[ilay - 1] = ilay - 1;
1161   }
1162   TGraph *gr_eta1 = new TGraph(19, idx, diff_lay);
1163   gr_eta1->SetMarkerStyle(20);
1164   gr_eta1->SetMarkerColor(2);
1165   gr_eta1->SetLineColor(2);
1166   gr_eta1->GetXaxis()->SetTitle(xtit);
1167   gr_eta1->GetYaxis()->SetTitle(ytit);
1168   mg->Add(gr_eta1, "pc");
1169   leg_mg->AddEntry(gr_eta1, "HB #eta = 1");
1170 
1171   for (int ilay = 1; ilay < 20; ilay++)
1172     diff_lay[ilay - 1] = diff[ilay][6];
1173   TGraph *gr_eta7 = new TGraph(18, idx, diff_lay);
1174   gr_eta7->SetMarkerStyle(22);
1175   gr_eta7->SetMarkerColor(4);
1176   gr_eta7->SetLineColor(4);
1177   mg->Add(gr_eta7, "pc");
1178   gr_eta7->GetXaxis()->SetTitle(xtit);
1179   gr_eta7->GetYaxis()->SetTitle(ytit);
1180   leg_mg->AddEntry(gr_eta7, "HB #eta = 7");
1181 
1182   for (int ilay = 1; ilay < 20; ilay++)
1183     diff_lay[ilay - 1] = diff[ilay][12];
1184   TGraph *gr_eta13 = new TGraph(18, idx, diff_lay);
1185   gr_eta13->SetMarkerStyle(29);
1186   gr_eta13->SetMarkerColor(kGreen);
1187   gr_eta13->SetLineColor(kGreen);
1188   gr_eta13->GetXaxis()->SetTitle(xtit);
1189   gr_eta13->GetYaxis()->SetTitle(ytit);
1190   mg->Add(gr_eta13, "pc");
1191   leg_mg->AddEntry(gr_eta13, "HB #eta = 13");
1192 
1193   for (int ilay = 1; ilay < 20; ilay++)
1194     diff_lay[ilay - 1] = diff[ilay][19];
1195   TGraph *gr_eta19 = new TGraph(18, idx, diff_lay);
1196   gr_eta19->SetMarkerStyle(24);
1197   gr_eta19->SetMarkerColor(kCyan);
1198   gr_eta19->SetLineColor(kCyan);
1199   gr_eta19->GetXaxis()->SetTitle(xtit);
1200   gr_eta19->GetYaxis()->SetTitle(ytit);
1201   mg->Add(gr_eta19, "pc");
1202   leg_mg->AddEntry(gr_eta19, "HE #eta = 20");
1203 
1204   for (int ilay = 1; ilay < 20; ilay++)
1205     diff_lay[ilay - 1] = diff[ilay][25];
1206   TGraph *gr_eta25 = new TGraph(18, idx, diff_lay);
1207   gr_eta25->SetMarkerStyle(26);
1208   gr_eta25->SetMarkerColor(kCyan);
1209   gr_eta25->SetLineColor(kCyan);
1210   gr_eta25->GetXaxis()->SetTitle(xtit);
1211   gr_eta25->GetYaxis()->SetTitle(ytit);
1212   mg->Add(gr_eta25, "pc");
1213   leg_mg->AddEntry(gr_eta25, "HE #eta = 26");
1214 
1215   TString cname = "c_diff_" + plot + tag;
1216   new TCanvas(cname, cname, 700, 400);
1217   mg->Draw("a");
1218   mg->GetXaxis()->SetTitle(xtit);
1219   mg->GetYaxis()->SetTitle(ytit);
1220   leg_mg->Draw("same");
1221 }
1222 
1223 void getDiff(TString fileName, TString plot) {
1224   TFile *hcalFile = new TFile(fileName);
1225   hcalFile->cd("g4SimHits");
1226 
1227   int istart = 200;
1228   if (plot.CompareTo("RadLen") == 0) {
1229     istart = 100;
1230   } else if (plot.CompareTo("StepLen") == 0) {
1231     istart = 300;
1232   }
1233 
1234   for (int ilay = 0; ilay < 22; ilay++) {
1235     char hname[10];
1236     sprintf(hname, "%i", istart + ilay + 1);
1237     TProfile *prof;
1238     gDirectory->GetObject(hname, prof);
1239     int nbins = prof->GetNbinsX();
1240     for (int itow = 0; itow < nbinmax; itow++) {
1241       double ent = 0, value = 0;
1242       for (int ii = 0; ii < nbins; ii++) {
1243         double xl = prof->GetBinLowEdge(ii + 1);
1244         double xu = prof->GetBinWidth(ii + 1);
1245         if (xl >= 0) {
1246           xu += xl;
1247         } else {
1248           double tmp = xu;
1249           xu = -xl;
1250           xl = xu - tmp;
1251         }
1252         double cont = (prof->GetBinContent(ii + 1));
1253         double dx = 1;
1254         if (cont > 0) {
1255           if (xl >= towLow[itow] && xu <= towHigh[itow]) {
1256             ent += dx;
1257             value += cont;
1258           } else if (xl < towLow[itow] && xu > towLow[itow]) {
1259             dx = (xu - towLow[itow]) / (xu - xl);
1260             ent += dx;
1261             value += dx * cont;
1262           } else if (xu > towHigh[itow] && xl < towHigh[itow]) {
1263             dx = (towHigh[itow] - xl) / (xu - xl);
1264             ent += dx;
1265             value += dx * cont;
1266           }
1267         }
1268       }
1269       if (ent > 0)
1270         mean[ilay][itow] = value / ent;
1271       else
1272         mean[ilay][itow] = 0.;
1273     }
1274   }
1275 
1276   for (int itow = 30; itow < nbinmax; itow++) {
1277     mean[0][itow] = mean[19][itow];
1278     mean[1][itow] = mean[20][itow];
1279     mean[19][itow] = 0;
1280     mean[20][itow] = 0;
1281   }
1282 
1283   mean[0][15] = 0.5 * (mean[0][14] + mean[0][17]);
1284   mean[1][15] = 0.5 * (mean[1][14] + mean[1][17]);
1285   mean[2][15] = 0.5 * (mean[2][14] + mean[2][17]);
1286   /*
1287   for (int itow=0; itow<nbinmax; itow++) {
1288     std::cout << "Tower " << itow;
1289     for (int ilay=0; ilay<22; ilay++) 
1290       cout << " " << ilay << " " << mean[ilay][itow];
1291     std::cout << "\n";
1292   }
1293   */
1294   for (int itow = 0; itow < 30; itow++) {
1295     for (int ilay = 20; ilay > 0; ilay--) {
1296       if (mean[ilay - 1][itow] <= 0) {
1297         mean[ilay - 1][itow] = mean[ilay][itow];
1298       }
1299     }
1300   }
1301 
1302   for (int itow = 0; itow < nbinmax; itow++) {
1303     if (itow > 4 && itow < 26)
1304       mean[19][itow] = 0;
1305     diff[0][itow] = mean[0][itow];
1306   }
1307 
1308   for (int itow = 15; itow < 17; itow++) {
1309     for (int ilay = 1; ilay < 22; ilay++) {
1310       if (mean[ilay][itow] > mean[ilay + 1][itow]) {
1311         for (int jlay = ilay + 1; jlay < 22; jlay++)
1312           mean[jlay][itow] = 0;
1313         break;
1314       }
1315     }
1316   }
1317 
1318   for (int ilay = 1; ilay < 22; ilay++) {
1319     for (int itow = 0; itow < nbinmax; itow++) {
1320       diff[ilay][itow] = mean[ilay][itow] - mean[ilay - 1][itow];
1321       if (diff[ilay][itow] < 0)
1322         diff[ilay][itow] = 0;
1323     }
1324   }
1325   /*  
1326   for (int itow=0; itow<nbinmax; itow++) {
1327     std::cout << "Tower " << itow;
1328     for (int ilay=0; ilay<22; ilay++) {
1329       cout << " " << ilay << " " << mean[ilay][itow] << " " << diff[ilay][itow];
1330     }
1331     std::cout << "\n";
1332   }
1333   */
1334 }
1335 
1336 void plotHE(int flag, int logy, std::string tag, int save) {
1337   double angle[31] = {-2.5,  -2.25, -2.00, -1.75,  -1.50, -1.25, -1.00, -0.75, -0.50, -0.25, -0.20,
1338                       -0.15, -0.10, -0.05, -0.025, 0,     0.025, 0.05,  0.10,  0.15,  0.20,  0.25,
1339                       0.50,  0.75,  1.00,  1.25,   1.50,  1.75,  2.00,  2.25,  2.50};
1340   double lAir[31] = {11440, 11440, 11440, 11440, 11440, 11440, 11440, 11440, 11440, 11440, 11440,
1341                      11440, 11440, 11440, 11440, 11569, 11569, 11440, 11440, 11440, 11440, 11440,
1342                      11440, 11440, 11440, 11440, 11440, 11440, 11440, 11440, 11440};
1343   double xAir[31] = {0.037941, 0.037941, 0.037941, 0.037941, 0.037941, 0.037941, 0.037941, 0.037941,
1344                      0.037941, 0.037941, 0.037941, 0.037941, 0.037941, 0.037941, 0.037941, 0.038369,
1345                      0.038369, 0.037941, 0.037941, 0.037941, 0.037941, 0.037941, 0.037941, 0.037941,
1346                      0.037941, 0.037941, 0.037941, 0.037941, 0.037941, 0.037941, 0.037941};
1347   double iAir[31] = {0.0162481, 0.0162481, 0.0162481, 0.0162481, 0.0162481, 0.0162481, 0.0162481, 0.0162481,
1348                      0.0162481, 0.0162481, 0.0162481, 0.0162481, 0.0162481, 0.0162481, 0.0162481, 0.0164314,
1349                      0.0164314, 0.0162481, 0.0162481, 0.0162481, 0.0162481, 0.0162481, 0.0162481, 0.0162481,
1350                      0.0162481, 0.0162481, 0.0162481, 0.0162481, 0.0162481, 0.0162481, 0.0162481};
1351   double lPol[31] = {60.8381, 60.8381, 60.8381, 60.8381, 60.8381, 60.8381, 60.8381, 60.8381, 60.8381, 60.8381, 60.8381,
1352                      60.8381, 60.8381, 60.8381, 60.8381, 0.00000, 0.00000, 60.8381, 60.8381, 60.8381, 60.8381, 60.8381,
1353                      60.8381, 60.8381, 60.8381, 60.8381, 60.8381, 60.8381, 60.8381, 60.8381, 60.8381};
1354   double xPol[31] = {0.129083, 0.129083, 0.129083, 0.129083, 0.129083, 0.129083, 0.129083, 0.129083,
1355                      0.129083, 0.129083, 0.129083, 0.129083, 0.129083, 0.129083, 0.129083, 0.000000,
1356                      0.000000, 0.129083, 0.129083, 0.129083, 0.129083, 0.129083, 0.129083, 0.129083,
1357                      0.129083, 0.129083, 0.129083, 0.129083, 0.129083, 0.129083, 0.129083};
1358   double iPol[31] = {0.0854129, 0.0854129, 0.0854129, 0.0854129, 0.0854129, 0.0854129, 0.0854129, 0.0854129,
1359                      0.0854129, 0.0854129, 0.0854129, 0.0854129, 0.0854129, 0.0854129, 0.0854129, 0.0000000,
1360                      0.0000000, 0.0854129, 0.0854129, 0.0854129, 0.0854129, 0.0854129, 0.0854129, 0.0854129,
1361                      0.0854129, 0.0854129, 0.0854129, 0.0854129, 0.0854129, 0.0854129, 0.0854129};
1362   double lScn[31] = {68.2125, 68.2125, 68.2125, 68.2125, 68.2125, 68.2125, 68.2125, 68.2125, 68.2125, 68.2125, 68.2125,
1363                      68.2125, 68.2125, 68.2125, 68.2125, 0.00000, 0.00000, 68.2125, 68.2125, 68.2125, 68.2125, 68.2125,
1364                      68.2125, 68.2125, 68.2125, 68.2125, 68.2125, 68.2125, 68.2125, 68.2125, 68.2125};
1365   double xScn[31] = {0.160352, 0.160352, 0.160352, 0.160352, 0.160352, 0.160352, 0.160352, 0.160352,
1366                      0.160352, 0.160352, 0.160352, 0.160352, 0.160352, 0.160352, 0.160352, 0.000000,
1367                      0.000000, 0.160352, 0.160352, 0.160352, 0.160352, 0.160352, 0.160352, 0.160352,
1368                      0.160352, 0.160352, 0.160352, 0.160352, 0.160352, 0.160352, 0.160352};
1369   double iScn[31] = {0.0973967, 0.0973967, 0.0973967, 0.0973967, 0.0973967, 0.0973967, 0.0973967, 0.0973967,
1370                      0.0973967, 0.0973967, 0.0973967, 0.0973967, 0.0973967, 0.0973967, 0.0973967, 0.0000000,
1371                      0.0000000, 0.0973967, 0.0973967, 0.0973967, 0.0973967, 0.0973967, 0.0973967, 0.0973967,
1372                      0.0973967, 0.0973967, 0.0973967, 0.0973967, 0.0973967, 0.0973967, 0.0973967};
1373   double lBra[31] = {1444.41, 1444.41, 1444.41, 1444.41, 1444.41, 1444.41, 1444.41, 1444.41, 1444.41, 1444.41, 1444.41,
1374                      1444.41, 1444.41, 1444.41, 1444.41, 1444.41, 1444.41, 1444.41, 1444.41, 1444.41, 1444.41, 1444.41,
1375                      1444.41, 1444.41, 1444.41, 1444.41, 1444.41, 1444.41, 1444.41, 1444.41, 1444.41};
1376   double xBra[31] = {96.7848, 96.7848, 96.7848, 96.7848, 96.7848, 96.7848, 96.7848, 96.7848, 96.7848, 96.7848, 96.7848,
1377                      96.7848, 96.7848, 96.7848, 96.7848, 96.7848, 96.7848, 96.7848, 96.7848, 96.7848, 96.7848, 96.7848,
1378                      96.7848, 96.7848, 96.7848, 96.7848, 96.7848, 96.7848, 96.7848, 96.7848, 96.7848};
1379   double iBra[31] = {8.79637, 8.79637, 8.79637, 8.79637, 8.79637, 8.79637, 8.79637, 8.79637, 8.79637, 8.79637, 8.79637,
1380                      8.79637, 8.79637, 8.79637, 8.79637, 8.79637, 8.79637, 8.79637, 8.79637, 8.79637, 8.79637, 8.79637,
1381                      8.79637, 8.79637, 8.79637, 8.79637, 8.79637, 8.79637, 8.79637, 8.79637, 8.79637};
1382   std::string nameMat[4] = {"Air", "Polythene", "Scintillator", "Brass"};
1383   int colMat[4] = {1, 2, 6, 4};
1384   int symbMat[4] = {24, 29, 25, 27};
1385 
1386   setStyle();
1387   gStyle->SetTitleOffset(1.2, "Y");
1388   char name[30], title[60], gname[12];
1389   TGraph *gr[4];
1390   int kfirst = 3;
1391   double ymi = 0, ymx = 100;
1392   if (flag < 0) {
1393     sprintf(name, "Step Length");
1394     sprintf(title, "Step Length (mm)");
1395     sprintf(gname, "stepLength");
1396     gr[0] = new TGraph(31, angle, lAir);
1397     gr[1] = new TGraph(31, angle, lPol);
1398     gr[2] = new TGraph(31, angle, lScn);
1399     gr[3] = new TGraph(31, angle, lBra);
1400     kfirst = 0;
1401     if (logy == 0) {
1402       ymx = 12000;
1403     } else {
1404       ymi = 10.0;
1405       ymx = 20000;
1406     }
1407   } else if (flag > 0) {
1408     sprintf(name, "# Interaction Length");
1409     sprintf(title, "Number of Interaction Length");
1410     sprintf(gname, "intLength");
1411     gr[0] = new TGraph(31, angle, iAir);
1412     gr[1] = new TGraph(31, angle, iPol);
1413     gr[2] = new TGraph(31, angle, iScn);
1414     gr[3] = new TGraph(31, angle, iBra);
1415     if (logy == 0) {
1416       ymx = 10;
1417     } else {
1418       ymi = 0.01;
1419       ymx = 20;
1420     }
1421   } else {
1422     sprintf(name, "# Radiation Length");
1423     sprintf(title, "Number of Radiation Length");
1424     sprintf(gname, "radLength");
1425     gr[0] = new TGraph(31, angle, xAir);
1426     gr[1] = new TGraph(31, angle, xPol);
1427     gr[2] = new TGraph(31, angle, xScn);
1428     gr[3] = new TGraph(31, angle, xBra);
1429     if (logy == 0) {
1430       ymx = 100;
1431     } else {
1432       ymi = 0.01;
1433       ymx = 200;
1434     }
1435   }
1436   gr[kfirst]->GetXaxis()->SetTitle("#phi ( ^{o})");
1437   gr[kfirst]->GetYaxis()->SetTitle(title);
1438   gr[kfirst]->SetTitle("");
1439   for (int i = 0; i < 4; i++) {
1440     int icol = colMat[i];
1441     int type = symbMat[i];
1442     gr[i]->SetMarkerSize(1.2);
1443     gr[i]->SetLineColor(icol);
1444     gr[i]->SetLineStyle(i + 1);
1445     gr[i]->SetLineWidth(2);
1446     gr[i]->SetMarkerColor(icol);
1447     gr[i]->SetMarkerStyle(type);
1448     gr[i]->GetYaxis()->SetRangeUser(ymi, ymx);
1449     gr[i]->GetXaxis()->SetRangeUser(-3.0, +3.0);
1450   }
1451 
1452   TString cname = "c1" + tag;
1453   TCanvas *c1 = new TCanvas(cname, name, 800, 500);
1454   if (logy != 0)
1455     gPad->SetLogy(1);
1456   gr[kfirst]->Draw("alp");
1457   for (int i = 0; i < 4; i++) {
1458     if (i != kfirst)
1459       gr[i]->Draw("lp");
1460   }
1461 
1462   double ylow = 0.4;
1463   char list[20];
1464   TLegend *leg1 = new TLegend(0.60, ylow, 0.90, ylow + 0.2);
1465   for (int i = 0; i < 4; i++) {
1466     sprintf(list, "%s", nameMat[i].c_str());
1467     leg1->AddEntry(gr[i], list, "LP");
1468   }
1469   leg1->SetHeader(name);
1470   leg1->SetFillColor(0);
1471   leg1->SetTextSize(0.04);
1472   leg1->Draw();
1473 
1474   if (save != 0) {
1475     char fname[20];
1476     if (save > 0)
1477       sprintf(fname, "%s.eps", gname);
1478     else
1479       sprintf(fname, "%s.gif", gname);
1480     c1->SaveAs(fname);
1481   }
1482 }
1483 
1484 void etaPhiCastorPlot(
1485     TString fileName, TString plot, TString type, bool etaPlus, int drawLeg, bool ifEta, std::string tag, bool debug) {
1486   TFile *hcalFile = new TFile(fileName);
1487   hcalFile->cd("g4SimHits");
1488   setStyle();
1489   if (debug)
1490     std::cout << fileName << " is opened at " << hcalFile << std::endl;
1491 
1492   TString xtit = TString("#eta");
1493   TString ztit = TString("eta");
1494   char ytit[80], ytpart[10];
1495   int ymin = 0, ymax = 20, istart = 200, ifirst = 0;
1496   double xh = 0.90;
1497   if (!etaPlus)
1498     ifirst = 10;
1499   if (type.CompareTo("EC") == 0) {
1500     sprintf(ytpart, "(EC)");
1501     ifirst += 2;
1502   } else if (type.CompareTo("HC") == 0) {
1503     sprintf(ytpart, "(HC)");
1504     ifirst += 4;
1505   } else if (type.CompareTo("ED") == 0) {
1506     sprintf(ytpart, "(Dead EC)");
1507     ifirst += 6;
1508   } else if (type.CompareTo("HD") == 0) {
1509     sprintf(ytpart, "(Dead HC)");
1510     ifirst += 8;
1511   } else {
1512     sprintf(ytpart, "(All)");
1513   }
1514   if (debug)
1515     std::cout << type << " Gives " << ifirst << " Title " << ytpart << std::endl;
1516   if (plot.CompareTo("RadLen") == 0) {
1517     sprintf(ytit, "Castor %s Material Budget (X_{0})", ytpart);
1518     ymin = 0;
1519     ymax = 200;
1520     istart = 100;
1521   } else if (plot.CompareTo("StepLen") == 0) {
1522     sprintf(ytit, "Castor %s Material Budget (Step Length)", ytpart);
1523     ymin = 0;
1524     ymax = 15000;
1525     istart = 300;
1526     xh = 0.70;
1527   } else {
1528     sprintf(ytit, "Castor %s Material Budget (#lambda)", ytpart);
1529     ymin = 0;
1530     ymax = 20;
1531     istart = 200;
1532   }
1533   if (!ifEta) {
1534     istart += 400;
1535     xtit = TString("#phi");
1536     ztit = TString("phi");
1537   }
1538   if (debug)
1539     std::cout << "Title (x) " << xtit << " (y) " << ytit << " First " << ifirst << ":" << istart << std::endl;
1540 
1541   TLegend *leg = new TLegend(xh - 0.13, 0.80, xh, 0.90);
1542   leg->SetBorderSize(1);
1543   leg->SetFillColor(10);
1544   leg->SetMargin(0.25);
1545   leg->SetTextSize(0.025);
1546 
1547   int nplots = 0;
1548   TProfile *prof[2];
1549   for (int ii = ifirst + 1; ii >= ifirst; ii--) {
1550     char hname[10], title[50];
1551     sprintf(hname, "%i", istart + ii);
1552     if (debug)
1553       std::cout << "[" << nplots << "] " << ii << " " << hname << "\n";
1554     gDirectory->GetObject(hname, prof[nplots]);
1555     if (debug)
1556       std::cout << "Histogram[" << nplots << "] : " << hname << " at " << prof[nplots] << std::endl;
1557     prof[nplots]->GetXaxis()->SetTitle(xtit);
1558     prof[nplots]->GetYaxis()->SetTitle(ytit);
1559     prof[nplots]->GetYaxis()->SetRangeUser(ymin, ymax);
1560     prof[nplots]->SetLineColor(colorLayer[nplots]);
1561     prof[nplots]->SetFillColor(colorLayer[nplots]);
1562     if (xh < 0.8)
1563       prof[nplots]->GetYaxis()->SetTitleOffset(1.7);
1564     if (ii == ifirst) {
1565       sprintf(title, "Front");
1566     } else {
1567       sprintf(title, "End");
1568     }
1569     leg->AddEntry(prof[nplots], title, "lf");
1570     nplots++;
1571   }
1572 
1573   TString cname = "c_" + plot + ztit + tag;
1574   TCanvas *cc1 = new TCanvas(cname, cname, 700, 600);
1575   if (xh < 0.8) {
1576     cc1->SetLeftMargin(0.15);
1577     cc1->SetRightMargin(0.05);
1578   }
1579 
1580   prof[0]->Draw("h");
1581   for (int i = 1; i < nplots; i++)
1582     prof[i]->Draw("h sames");
1583   if (drawLeg > 0)
1584     leg->Draw("sames");
1585 }
1586 
1587 void efficiencyPlot(TString fileName, TString type, bool ifEtaPhi, double maxEta, std::string tag, bool debug) {
1588   TFile *hcalFile = new TFile(fileName);
1589   hcalFile->cd("g4SimHits");
1590   setStyle();
1591 
1592   int id0 = 1300, idpl1 = 8, idpl2 = 0;
1593   char hname[100], title[100];
1594   if (type.CompareTo("HB") == 0) {
1595     idpl1 = 1;
1596     idpl2 = 2;
1597     sprintf(title, "Efficiency for HB (%s)", tag.c_str());
1598   } else if (type.CompareTo("HE") == 0) {
1599     idpl1 = 3;
1600     idpl2 = 4;
1601     sprintf(title, "Efficiency for HE (%s)", tag.c_str());
1602   } else if (type.CompareTo("HO") == 0) {
1603     idpl1 = 5;
1604     idpl2 = 6;
1605     sprintf(title, "Efficiency for HO (%s)", tag.c_str());
1606   } else if (type.CompareTo("HF") == 0) {
1607     idpl1 = 7;
1608     sprintf(title, "Efficiency for HF (%s)", tag.c_str());
1609   } else {
1610     sprintf(title, "Efficiency for HCAL (%s)", tag.c_str());
1611   }
1612   TLegend *leg = new TLegend(0.70, 0.82, 0.90, 0.90);
1613   leg->SetBorderSize(1);
1614   leg->SetFillColor(10);
1615   leg->SetMargin(0.25);
1616   leg->SetTextSize(0.03);
1617 
1618   if (ifEtaPhi) {
1619     id0 = 1400;
1620     TH2F *hist0, *hist1, *hist2;
1621     sprintf(hname, "%i", id0);
1622     gDirectory->GetObject(hname, hist0);
1623     sprintf(hname, "%i", id0 + idpl1);
1624     gDirectory->GetObject(hname, hist1);
1625     if (idpl2 > 0) {
1626       sprintf(hname, "%i", id0 + idpl2);
1627       gDirectory->GetObject(hname, hist2);
1628     } else {
1629       hist2 = 0;
1630     }
1631     if (debug)
1632       std::cout << "Get Histos at " << hist0 << " and " << hist1 << "\n";
1633     if (hist0 && hist1) {
1634       double xmin = hist0->GetXaxis()->GetXmin();
1635       double xmax = hist0->GetXaxis()->GetXmax();
1636       int nbinX = hist0->GetNbinsX();
1637       double ymin = hist0->GetYaxis()->GetXmin();
1638       double ymax = hist0->GetYaxis()->GetXmax();
1639       int nbinY = hist0->GetNbinsY();
1640       if (debug)
1641         std::cout << "NbinX " << nbinX << " range " << std::setprecision(5) << xmin << ":" << std::setprecision(5)
1642                   << xmax << " "
1643                   << "NbinY " << nbinY << " range " << std::setprecision(5) << ymin << ":" << std::setprecision(5)
1644                   << ymax << "\n";
1645       TH2D *hist = new TH2D("hist", title, nbinX, xmin, xmax, nbinY, ymin, ymax);
1646       TH2D *histe = 0;
1647       if (hist2)
1648         histe = new TH2D("histe", title, nbinX, xmin, xmax, nbinY, ymin, ymax);
1649       for (int ibx = 0; ibx < nbinX; ++ibx) {
1650         for (int iby = 0; iby < nbinY; ++iby) {
1651           double contN = hist1->GetBinContent(ibx + 1, iby + 1);
1652           double contD = hist0->GetBinContent(ibx + 1, iby + 1);
1653           double cont = contN / std::max(contD, 1.0);
1654           hist->SetBinContent(ibx + 1, iby + 1, cont);
1655           if (hist2) {
1656             contN = hist2->GetBinContent(ibx + 1, iby + 1);
1657             cont = contN / std::max(contD, 1.0);
1658             histe->SetBinContent(ibx + 1, iby + 1, cont);
1659           }
1660         }
1661       }
1662       hist->GetXaxis()->SetTitle("#eta");
1663       hist->GetYaxis()->SetTitle("#phi");
1664       hist->GetZaxis()->SetTitle(title);
1665       hist->GetZaxis()->SetTitleOffset(.8);
1666       new TCanvas(title, title, 700, 400);
1667       hist->SetLineColor(2);
1668       hist->SetLineStyle(1);
1669       hist->SetLineWidth(1);
1670       if (maxEta > 0)
1671         hist->GetXaxis()->SetRangeUser(-maxEta, maxEta);
1672       hist->Draw("lego fb bb");
1673       leg->AddEntry(hist, "At least 1 layer", "l");
1674       if (histe) {
1675         histe->SetLineColor(4);
1676         histe->SetLineStyle(2);
1677         histe->SetLineWidth(1);
1678         if (maxEta > 0)
1679           histe->GetXaxis()->SetRangeUser(-maxEta, maxEta);
1680         histe->Draw("lego fb bb sames");
1681         leg->AddEntry(histe, "All layers", "l");
1682       }
1683       leg->Draw("sames");
1684     }
1685   } else {
1686     TH1F *hist0, *hist1, *hist2;
1687     sprintf(hname, "%i", id0);
1688     gDirectory->GetObject(hname, hist0);
1689     sprintf(hname, "%i", id0 + idpl1);
1690     gDirectory->GetObject(hname, hist1);
1691     if (idpl2 > 0) {
1692       sprintf(hname, "%i", id0 + idpl2);
1693       gDirectory->GetObject(hname, hist2);
1694     } else {
1695       hist2 = 0;
1696     }
1697     if (debug)
1698       std::cout << "Get Histos at " << hist0 << " and " << hist1 << "\n";
1699     if (hist0 && hist1) {
1700       double xmin = hist0->GetXaxis()->GetXmin();
1701       double xmax = hist0->GetXaxis()->GetXmax();
1702       int nbinX = hist0->GetNbinsX();
1703       if (debug)
1704         std::cout << "Nbin " << nbinX << " range " << std::setprecision(5) << xmin << ":" << std::setprecision(5)
1705                   << xmax << "\n";
1706       TH1D *hist = new TH1D("hist", title, nbinX, xmin, xmax);
1707       TH1D *histe = 0;
1708       if (hist2)
1709         histe = new TH1D("histe", title, nbinX, xmin, xmax);
1710       for (int ib = 0; ib < nbinX; ++ib) {
1711         double contN = hist1->GetBinContent(ib + 1);
1712         double contD = hist0->GetBinContent(ib + 1);
1713         double cont = contN / std::max(contD, 1.0);
1714         hist->SetBinContent(ib + 1, cont);
1715         /*
1716     double eror  = std::sqrt(contN)/std::max(contD,1.0);
1717     hist->SetBinError(ib+1, eror);
1718     */
1719         if (hist2) {
1720           contN = hist2->GetBinContent(ib + 1);
1721           cont = contN / std::max(contD, 1.0);
1722           histe->SetBinContent(ib + 1, cont);
1723         }
1724       }
1725       hist->GetXaxis()->SetTitle("#eta");
1726       hist->GetYaxis()->SetTitle(title);
1727       hist->GetYaxis()->SetTitleOffset(0.8);
1728       new TCanvas(title, title, 700, 400);
1729       hist->SetLineColor(2);
1730       hist->SetLineStyle(1);
1731       hist->SetLineWidth(1);
1732       if (maxEta > 0)
1733         hist->GetXaxis()->SetRangeUser(-maxEta, maxEta);
1734       hist->Draw();
1735       leg->AddEntry(hist, "At least 1 layer", "l");
1736       if (histe) {
1737         histe->SetLineColor(4);
1738         histe->SetLineStyle(2);
1739         histe->SetLineWidth(1);
1740         if (maxEta > 0)
1741           histe->GetXaxis()->SetRangeUser(-maxEta, maxEta);
1742         histe->Draw("sames");
1743         leg->AddEntry(histe, "All layers", "l");
1744       }
1745       leg->Draw("sames");
1746     }
1747   }
1748 }
1749 
1750 void etaPhiFwdPlot(TString fileName, TString plot, int first, int last, int drawLeg, std::string tag, bool debug) {
1751   TFile *hcalFile = new TFile(fileName);
1752   hcalFile->cd("g4SimHits");
1753   setStyle();
1754 
1755   TString xtit = TString("#eta");
1756   TString ztit = TString("eta");
1757   TString ytit = "none";
1758   int ymin = 0, ymax = 20, istart = 200;
1759   double xh = 0.90, xl = 0.1;
1760   if (plot.CompareTo("RadLen") == 0) {
1761     ytit = TString("Material Budget (X_{0})");
1762     ymin = 0;
1763     ymax = 400;
1764     istart = 100;
1765   } else if (plot.CompareTo("StepLen") == 0) {
1766     ytit = TString("Material Budget (Step Length)");
1767     ymin = 0;
1768     ymax = 20000;
1769     istart = 300;
1770     xh = 0.70, xl = 0.15;
1771   } else {
1772     ytit = TString("Material Budget (#lambda)");
1773     ymin = 0;
1774     ymax = 30;
1775     istart = 200;
1776   }
1777 
1778   int index[10] = {9, 0, 1, 2, 3, 8, 4, 7, 5, 6};
1779   std::string label[10] = {"Empty",
1780                            "Beam Pipe",
1781                            "Tracker",
1782                            "EM Calorimeter",
1783                            "Hadron Calorimeter",
1784                            "Muon System",
1785                            "Forward Hadron Calorimeter",
1786                            "Shielding",
1787                            "TOTEM",
1788                            "CASTOR"};
1789 
1790   TLegend *leg = new TLegend(xl, 0.75, xl + 0.3, 0.90);
1791   leg->SetBorderSize(1);
1792   leg->SetFillColor(10);
1793   leg->SetMargin(0.25);
1794   leg->SetTextSize(0.018);
1795 
1796   int nplots = 0;
1797   TProfile *prof[nlaymax];
1798   for (int ii = last; ii >= first; ii--) {
1799     char hname[10], title[50];
1800     sprintf(hname, "%i", istart + index[ii]);
1801     gDirectory->GetObject(hname, prof[nplots]);
1802     prof[nplots]->GetXaxis()->SetTitle(xtit);
1803     prof[nplots]->GetYaxis()->SetTitle(ytit);
1804     prof[nplots]->GetYaxis()->SetRangeUser(ymin, ymax);
1805     prof[nplots]->SetLineColor(colorLayer[ii]);
1806     prof[nplots]->SetFillColor(colorLayer[ii]);
1807     if (xh < 0.8)
1808       prof[nplots]->GetYaxis()->SetTitleOffset(1.7);
1809     sprintf(title, "%s", label[ii].c_str());
1810     leg->AddEntry(prof[nplots], title, "lf");
1811     if (debug) {
1812       int nbinX = prof[nplots]->GetNbinsX();
1813       double xmin = prof[nplots]->GetXaxis()->GetXmin();
1814       double xmax = prof[nplots]->GetXaxis()->GetXmax();
1815       double dx = (xmax - xmin) / nbinX;
1816       cout << "Hist " << ii;
1817       for (int ibx = 0; ibx < nbinX; ibx++) {
1818         double xx1 = xmin + ibx * dx;
1819         double cont = prof[nplots]->GetBinContent(ibx + 1);
1820         cout << " | " << ibx << "(" << xx1 << ":" << (xx1 + dx) << ") " << cont;
1821       }
1822       cout << "\n";
1823     }
1824     nplots++;
1825   }
1826 
1827   TString cname = "c_" + plot + ztit + tag;
1828   TCanvas *cc1 = new TCanvas(cname, cname, 700, 600);
1829   if (xh < 0.8) {
1830     cc1->SetLeftMargin(0.15);
1831     cc1->SetRightMargin(0.05);
1832   }
1833 
1834   prof[0]->Draw("h");
1835   for (int i = 1; i < nplots; i++)
1836     prof[i]->Draw("h sames");
1837   if (drawLeg > 0)
1838     leg->Draw("sames");
1839 }
1840 
1841 void setStyle() {
1842   gStyle->SetCanvasBorderMode(0);
1843   gStyle->SetCanvasColor(kWhite);
1844   gStyle->SetPadColor(kWhite);
1845   gStyle->SetFrameBorderMode(0);
1846   gStyle->SetFrameBorderSize(1);
1847   gStyle->SetFrameFillColor(0);
1848   gStyle->SetFrameFillStyle(0);
1849   gStyle->SetFrameLineColor(1);
1850   gStyle->SetFrameLineStyle(1);
1851   gStyle->SetFrameLineWidth(1);
1852   gStyle->SetOptStat(0);
1853   gStyle->SetLegendBorderSize(1);
1854   gStyle->SetOptTitle(0);
1855   gStyle->SetTitleOffset(2.5, "Y");
1856 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!