Calo/macros/MakeGVPlots.C

0001 //////////////////////////////////////////////////////////////////////////////
0002 //
0003 // Usage:
0004 // .L MakeGVPlots.C+g
0005 //
0006 //   To make plot from a file created using CaloSimHitAnalysis
0007 //     makeGVPlots(fname, ifG4, todomin, todomax, tag, text, save, dirnm);
0008 //
0009 //   To plot on the same Canvas plots from Geant4 and GeantV done using
0010 //   CaloSimHitAnalysis
0011 //     makeGV2Plots(fnmG4, fnmGV, todomin, todomax, normalize, tag, text,
0012 //                  save, modet, mode, dirnm)
0013 //
0014 //   To plot from passive hit collection produced using CaloSimHitAnalysis
0015 //      makeGVSPlots(fnmG4, fnmGV, todomin, todomax, tag, text, save, dirnm)
0016 //
0017 //   where
0018 //     fname   std::string   Name of the ROOT file (analG4.root)
0019 //     fnmG4   std::string   Name of the ROOT file created using Geant4
0020 //                           (analG4.root)
0021 //     fnmGV   std::string   Name of the ROOT file created using GeantV
0022 //                           (analGV.root)
0023 //     ifGV    bool          If the file created using Geant4 (true)
0024 //     todomin int           The first one in the series to be created (0)
0025 //     todomax int           The last one in the series to be created
0026 //                           (3:2:7 for GVPlots:GV2Plots:GVSPlots)
0027 //     tag     std::string   To be added to the name of the canvas ("")
0028 //     text    std::string   To be added to the title of the histogram ("")
0029 //     save    bool          If the canvas is to be saved as jpg file (false)
0030 //     modet   int           Packed word (xhto) where non-zero value of a
0031 //                           digit indicate text message on the plot
0032 //                           (o for header "CMS Simulation"; t input type;
0033 //                           h for input specification; x for B-field/MT etc)
0034 //     mode    int           Flag if one file is G4 and other GV (0), or
0035 //                           both files are GV (1) or G4 (2) (defualt 0)
0036 //     dirnm   std::string   Name of the directory ("caloSimHitAnalysis")
0037 //
0038 //   The histogram series have different meanings for each function
0039 //
0040 //   GVPlots (17 in total):
0041 //     "Energy deposit per Hit", "Hit time", "Total energy deposit",
0042 //     "Energy deposit after 15 ns",  "R vs Z", "R vs Z for hits after 15 ns",
0043 //     "#phi vs #eta", "Energy deposit", "Time of Hit",
0044 //     "Energy deposit per Hit after 15 ns", "Total energy deposit in 100 ns",
0045 //     "Energy deposit for EM particles", "Energy deposit for non-EM particles",
0046 //     "R", "Z", "#eta", "phi"
0047 //
0048 //   GV2Plots (14 in total):
0049 //     "Energy deposit per Hit", "Hit time", "Total energy deposit",
0050 //     "Energy deposit after 15 ns", "Energy deposit", "Time of Hit",
0051 //     "Energy deposit per Hit after 15 ns", "Total energy deposit in 100 ns",
0052 //     "Energy deposit for EM particles", "Energy deposit for non-EM particles",
0053 //     "R", "Z", "#eta", "phi"
0054 //
0055 //   GVSPlots (8 in total):
0056 //     "Total Hits", "Tracks", "Step Length',
0057 //     "Energy Deposit in all components", "Time of all hits",
0058 //     "Energy Deposit in tracker subdetectors",
0059 //     "Time of hits in tracker subdetectors"
0060 //
0061 //   All plots in GVPlots or GV2Plots happen for EB, EE, HB and HE
0062 //   There are 6 subdetectors for tracker:
0063 //      Pixel Barrel, Pixel Forward, TIB, TID, TOB, TEC
0064 //
0065 //////////////////////////////////////////////////////////////////////////////
0066
0067 #include <TCanvas.h>
0068 #include <TChain.h>
0069 #include <TFile.h>
0070 #include <TFitResult.h>
0071 #include <TFitResultPtr.h>
0072 #include <TGraphAsymmErrors.h>
0073 #include <TH1D.h>
0074 #include <TH2D.h>
0075 #include <TLegend.h>
0076 #include <TPaveStats.h>
0077 #include <TPaveText.h>
0078 #include <TProfile.h>
0079 #include <TROOT.h>
0080 #include <TStyle.h>
0081
0082 #include <fstream>
0083 #include <iomanip>
0084 #include <iostream>
0085 #include <string>
0086 #include <vector>
0087
0088 void setTDRStyle();
0089
0090 void makeGVPlots(std::string fname = "analG4.root",
0091                  bool ifG4 = true,
0092                  int todomin = 0,
0093                  int todomax = 3,
0094                  std::string tag = "",
0095                  std::string text = "",
0096                  bool save = false,
0097                  std::string dirnm = "caloSimHitAnalysis") {
0098   std::string names[17] = {"EdepT",
0099                            "TimeT",
0100                            "Etot",
0101                            "Edep15",
0102                            "rz",
0103                            "rz2",
0104                            "etaphi",
0105                            "Edep",
0106                            "Time",
0107                            "EdepT15",
0108                            "EtotG",
0109                            "EdepEM",
0110                            "EdepHad",
0111                            "rr",
0112                            "zz",
0113                            "eta",
0114                            "phi"};
0115   std::string namex[17] = {"Edep",
0116                            "Time",
0117                            "Etot",
0118                            "Edep15",
0119                            "rz",
0120                            "rz2",
0121                            "etaphi",
0122                            "EdepX",
0123                            "TimeX",
0124                            "EdepT15",
0125                            "EtotG",
0126                            "EdepEM",
0127                            "EdepHad",
0128                            "rr",
0129                            "zz",
0130                            "eta",
0131                            "phi"};
0132   int types[17] = {1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
0133   int dets[4] = {0, 1, 2, 3};
0134   std::string detName[4] = {"EB", "EE", "HB", "HE"};
0135   std::string xtitle[17] = {"Energy deposit per Hit (GeV)",
0136                             "Hit time (ns)",
0137                             "Total energy deposit (GeV)",
0138                             "Energy deposit (GeV) after 15 ns",
0139                             "z (cm)",
0140                             "z (cm)",
0141                             "#eta",
0142                             "Energy deposit (GeV)",
0143                             "Time of Hit (ns)",
0144                             "Energy deposit (GeV) per Hit after 15 ns",
0145                             "Total energy deposit (GeV) in 100 ns",
0146                             "Energy deposit (GeV) for EM particles",
0147                             "Energy deposit (GeV) for non-EM particles",
0148                             "R (cm)",
0149                             "z (cm)",
0150                             "#eta",
0151                             "phi"};
0152   std::string ytitle[17] = {"Hits",
0153                             "Hits",
0154                             "Events",
0155                             "Hits",
0156                             "R (cm)",
0157                             "R (cm)",
0158                             "#phi",
0159                             "Hits",
0160                             "Hits",
0161                             "Hits",
0162                             "Events",
0163                             "Hits",
0164                             "Hits",
0165                             "Hits",
0166                             "Hits",
0167                             "Hits",
0168                             "Hits"};
0169
0170   gStyle->SetCanvasBorderMode(0);
0171   gStyle->SetCanvasColor(kWhite);
0172   gStyle->SetPadColor(kWhite);
0173   gStyle->SetFillColor(kWhite);
0174   gStyle->SetOptStat(111110);
0175   TFile *file = new TFile(fname.c_str());
0176   if (file) {
0177     TDirectory *dir = (TDirectory *)file->FindObjectAny(dirnm.c_str());
0178     char cname[100], name[100], name1[100], title[100];
0179     for (int i1 = 0; i1 < 4; ++i1) {
0180       for (int i2 = todomin; i2 <= todomax; ++i2) {
0181         if (types[i2] == 1) {
0182           sprintf(name, "%s%d", names[i2].c_str(), dets[i1]);
0183           sprintf(name1, "%s%d", namex[i2].c_str(), dets[i1]);
0184           if (ifG4)
0185             sprintf(title, "%s %s (Geant4 Simulation)", text.c_str(), detName[i1].c_str());
0186           else
0187             sprintf(title, "%s %s (GeantV Simulation)", text.c_str(), detName[i1].c_str());
0188         } else if (i1 == 0) {
0189           sprintf(name, "%s", names[i2].c_str());
0190           sprintf(name1, "%s", namex[i2].c_str());
0191           if (ifG4)
0192             sprintf(title, "%s (Geant4 Simulation)", text.c_str());
0193           else
0194             sprintf(title, "%s (GeantV Simulation)", text.c_str());
0195         } else {
0196           continue;
0197         }
0198         TH1D *hist1(nullptr);
0199         TH2D *hist2(nullptr);
0200         if (types[i2] == 1)
0201           hist1 = (TH1D *)dir->FindObjectAny(name);
0202         else
0203           hist2 = (TH2D *)dir->FindObjectAny(name);
0204         //      std::cout << name << " read out at " << hist1 << ":" << hist2 << std::endl;
0205         if ((hist1 != nullptr) || (hist2 != nullptr)) {
0206           if (ifG4)
0207             sprintf(cname, "%sG4%s", name1, tag.c_str());
0208           else
0209             sprintf(cname, "%sGV%s", name1, tag.c_str());
0210           TCanvas *pad = new TCanvas(cname, cname, 500, 500);
0211           pad->SetRightMargin(0.10);
0212           pad->SetTopMargin(0.10);
0213           if (types[i2] == 1) {
0214             hist1->GetYaxis()->SetTitleOffset(1.2);
0215             hist1->GetYaxis()->SetTitle(ytitle[i2].c_str());
0216             hist1->GetXaxis()->SetTitle(xtitle[i2].c_str());
0217             hist1->SetTitle(title);
0218             pad->SetLogy();
0219             hist1->Draw();
0220           } else {
0221             hist2->GetYaxis()->SetTitleOffset(1.2);
0222             hist2->GetYaxis()->SetTitle(ytitle[i2].c_str());
0223             hist2->GetXaxis()->SetTitle(xtitle[i2].c_str());
0224             hist2->SetTitle(title);
0225             hist2->Draw();
0226           }
0227           pad->Update();
0228           TPaveStats *st1 = ((hist1 != nullptr) ? ((TPaveStats *)hist1->GetListOfFunctions()->FindObject("stats"))
0229                                                 : ((TPaveStats *)hist2->GetListOfFunctions()->FindObject("stats")));
0230           if (st1 != NULL) {
0231             st1->SetY1NDC(0.70);
0232             st1->SetY2NDC(0.90);
0233             st1->SetX1NDC(0.65);
0234             st1->SetX2NDC(0.90);
0235           }
0236           pad->Modified();
0237           pad->Update();
0238           if (save) {
0239             sprintf(name, "c_%s.jpg", pad->GetName());
0240             pad->Print(name);
0241           }
0242         }
0243       }
0244     }
0245   }
0246 }
0247
0248 void makeGV2Plots(std::string fnmG4 = "analG4.root",
0249                   std::string fnmGV = "analGV.root",
0250                   int todomin = 0,
0251                   int todomax = 2,
0252                   bool normalize = true,
0253                   std::string tag = "",
0254                   std::string text = "",
0255                   int save = 0,
0256                   int modet = 0,
0257                   int mode = 0,
0258                   std::string dirnm = "caloSimHitAnalysis") {
0259   std::string names[14] = {"EdepT",
0260                            "TimeT",
0261                            "Etot",
0262                            "Edep15",
0263                            "Edep",
0264                            "Time",
0265                            "EdepT15",
0266                            "EtotG",
0267                            "EdepEM",
0268                            "EdepHad",
0269                            "rr",
0270                            "zz",
0271                            "eta",
0272                            "phi"};
0273   std::string namex[14] = {"Edep",
0274                            "Time",
0275                            "Etot",
0276                            "Edep15",
0277                            "EdepX",
0278                            "TimeX",
0279                            "EdepT15",
0280                            "EtotG",
0281                            "EdepEM",
0282                            "EdepHad",
0283                            "rr",
0284                            "zz",
0285                            "eta",
0286                            "phi"};
0287   int dets[4] = {0, 1, 2, 3};
0288   std::string detName[4] = {"EB", "EE", "HB", "HE"};
0289   std::string xtitle[14] = {"Energy deposit per Hit (GeV)",
0290                             "Hit time (ns)",
0291                             "Total energy deposit (GeV)",
0292                             "Energy deposit (GeV) after 15 ns",
0293                             "Energy deposit (GeV)",
0294                             "Time of Hit (ns)",
0295                             "Energy deposit (GeV) per Hit after 15 ns",
0296                             "Total energy deposit (GeV) in 100 ns",
0297                             "Energy deposit (GeV) for EM particles",
0298                             "Energy deposit (GeV) for non-EM particles",
0299                             "R (cm)",
0300                             "z (cm)",
0301                             "#eta",
0302                             "phi"};
0303   std::string ytitle[14] = {"Hits",
0304                             "Hits",
0305                             "Events",
0306                             "Hits",
0307                             "Hits",
0308                             "Hits",
0309                             "Hits",
0310                             "Events",
0311                             "Hits",
0312                             "Hits",
0313                             "Hits",
0314                             "Hits",
0315                             "Hits",
0316                             "Hits"};
0317   std::string title1[4] = {"2 GeV", "10 GeV", "50 GeV", "100 GeV"};
0318   std::string title2[2] = {"#eta = 1.0; #phi = 1.0", "|#eta| < 3.0; 0.0 < #phi < 2#pi"};
0319   std::string title3[4] = {
0320       "B = 0.0 Tesla", "B = 3.8 Tesla", "Multithreaded (B = 0.0 Tesla)", "Multithreaded (B = 3.8 Tesla)"};
0321
0322   if (modet >= 10) {
0323     setTDRStyle();
0324   } else {
0325     gStyle->SetCanvasBorderMode(0);
0326     gStyle->SetCanvasColor(kWhite);
0327     gStyle->SetPadColor(kWhite);
0328     gStyle->SetFillColor(kWhite);
0329   }
0330
0331   if (normalize)
0332     gStyle->SetOptStat(0);
0333   else
0334     gStyle->SetOptStat(111110);
0335   TFile *file1 = new TFile(fnmG4.c_str());
0336   TFile *file2 = new TFile(fnmGV.c_str());
0337   if (file1 && file2) {
0338     TDirectory *dir1 = (TDirectory *)file1->FindObjectAny(dirnm.c_str());
0339     TDirectory *dir2 = (TDirectory *)file2->FindObjectAny(dirnm.c_str());
0340     char name[100], cname[100], title[100];
0341     for (int i1 = 0; i1 < 4; ++i1) {
0342       for (int i2 = todomin; i2 <= todomax; ++i2) {
0343         sprintf(name, "%s%d", names[i2].c_str(), dets[i1]);
0344         sprintf(cname, "%s%s%s", namex[i2].c_str(), detName[i1].c_str(), tag.c_str());
0345         if (modet >= 10)
0346           sprintf(title, "");
0347         else if (mode == 1)
0348           sprintf(title, "%s  %s (2 runs of GeantV)", text.c_str(), detName[i1].c_str());
0349         else if (mode == 2)
0350           sprintf(title, "%s  %s (2 runs of Geant4)", text.c_str(), detName[i1].c_str());
0351         else
0352           sprintf(title, "%s  %s (Geant4 vs GeantV)", text.c_str(), detName[i1].c_str());
0353         TH1D *hist[2];
0354         hist[0] = (TH1D *)dir1->FindObjectAny(name);
0355         hist[1] = (TH1D *)dir2->FindObjectAny(name);
0356         if ((hist[0] != nullptr) && (hist[1] != nullptr)) {
0357           // Plot superimposed histograms
0358           double ymx(0.78);
0359           TCanvas *pad = new TCanvas(cname, cname, 500, 500);
0360           TLegend *legend = new TLegend(0.44, ymx, 0.64, 0.89);
0361           pad->SetRightMargin(0.10);
0362           pad->SetTopMargin(0.10);
0363           pad->SetLogy();
0364           pad->SetFillColor(kWhite);
0365           legend->SetFillColor(kWhite);
0366           int icol[2] = {1, 2};
0367           int isty[2] = {1, 2};
0368           int imty[2] = {20, 24};
0369           std::string type[2] = {"Geant4", "GeantV"};
0370           double ymax(0.90);
0371           double total[2] = {0, 0};
0372           for (int i = 0; i < 2; ++i) {
0373             hist[i]->GetYaxis()->SetTitleOffset(1.2);
0374             hist[i]->GetYaxis()->SetTitle(ytitle[i2].c_str());
0375             hist[i]->GetXaxis()->SetTitle(xtitle[i2].c_str());
0376             hist[i]->SetTitle(title);
0377             hist[i]->SetMarkerStyle(imty[i]);
0378             hist[i]->SetMarkerColor(icol[i]);
0379             hist[i]->SetLineColor(icol[i]);
0380             hist[i]->SetLineStyle(isty[i]);
0381             hist[i]->SetNdivisions(505, "X");
0382             total[i] = hist[i]->GetEntries();
0383             if (mode == 1)
0384               sprintf(name, "%s (run %d)", type[1].c_str(), i);
0385             else if (mode == 2)
0386               sprintf(name, "%s (run %d)", type[0].c_str(), i);
0387             else
0388               sprintf(name, "%s", type[i].c_str());
0389             legend->AddEntry(hist[i], name, "lp");
0390             if (i == 0) {
0391               if (normalize)
0392                 hist[i]->DrawNormalized("hist");
0393               else
0394                 hist[i]->Draw();
0395             } else {
0396               if (normalize)
0397                 hist[i]->DrawNormalized("sames hist");
0398               else
0399                 hist[i]->Draw("sames");
0400             }
0401             pad->Update();
0402           }
0403           legend->Draw("same");
0404           pad->Modified();
0405           pad->Update();
0406           for (int i = 0; i < 2; ++i) {
0407             TPaveStats *st = (TPaveStats *)hist[i]->GetListOfFunctions()->FindObject("stats");
0408             if (st != NULL) {
0409               st->SetLineColor(icol[i]);
0410               st->SetTextColor(icol[i]);
0411               st->SetY1NDC(ymax - 0.15);
0412               st->SetY2NDC(ymax);
0413               st->SetX1NDC(0.65);
0414               st->SetX2NDC(0.90);
0415               ymax -= 0.15;
0416             }
0417           }
0418           pad->Modified();
0419           pad->Update();
0420           if (ymx > ymax)
0421             ymx = ymax;
0422           ymx -= 0.005;
0423           int indx = (modet / 10) % 10;
0424           double ymi = ymx - 0.05;
0425           if (indx > 0 && indx <= 4) {
0426             sprintf(title, "%s Electrons", title1[indx - 1].c_str());
0427             TPaveText *txt0 = new TPaveText(0.70, ymi, 0.895, ymx, "blNDC");
0428             txt0->SetFillColor(0);
0429             txt0->AddText(title);
0430             txt0->Draw("same");
0431             ymx -= 0.05;
0432           }
0433           indx = (modet / 100) % 10;
0434           ymi = ymx - 0.05;
0435           if (indx > 0 && indx <= 2) {
0436             sprintf(title, "%s", title2[indx - 1].c_str());
0437             TPaveText *txt0 = new TPaveText(0.65, ymi, 0.895, ymx, "blNDC");
0438             txt0->SetFillColor(0);
0439             txt0->AddText(title);
0440             txt0->Draw("same");
0441             ymx -= 0.05;
0442           }
0443           indx = (modet / 1000) % 10;
0444           ymi = ymx - 0.05;
0445           if (indx > 0 && indx <= 4) {
0446             sprintf(title, "%s", title3[indx - 1].c_str());
0447             TPaveText *txt0 = new TPaveText(0.55, ymi, 0.895, ymx, "blNDC");
0448             txt0->SetFillColor(0);
0449             txt0->AddText(title);
0450             txt0->Draw("same");
0451             ymx -= 0.05;
0452           }
0453           indx = modet % 10;
0454           if (indx > 0) {
0455             sprintf(title, "CMS Simulation");
0456             TPaveText *txt0 = new TPaveText(0.10, 0.91, 0.35, 0.96, "blNDC");
0457             txt0->SetFillColor(0);
0458             txt0->AddText(title);
0459             txt0->Draw("same");
0460           }
0461           if (save != 0) {
0462             if (save > 0)
0463               sprintf(name, "c_%s.pdf", pad->GetName());
0464             else
0465               sprintf(name, "c_%s.jpg", pad->GetName());
0466             pad->Print(name);
0467           }
0468
0469           // Ratio plot
0470           if (normalize) {
0471             int nbin = hist[0]->GetNbinsX();
0472             double xmin = hist[0]->GetBinLowEdge(1);
0473             double dx = hist[0]->GetBinWidth(1);
0474             double xmax = xmin + nbin * dx;
0475             double fac = total[1] / total[0];
0476             int npt = 0;
0477             double sumNum(0), sumDen(0), xpt[200], dxp[200], ypt[200], dyp[200];
0478             for (int i = 0; i < nbin; ++i) {
0479               if (hist[0]->GetBinContent(i + 1) > 0 && hist[1]->GetBinContent(i + 1) > 0) {
0480                 ypt[npt] = (fac * hist[0]->GetBinContent(i + 1) / hist[1]->GetBinContent(i + 1));
0481                 double er1 = hist[0]->GetBinError(i + 1) / hist[0]->GetBinContent(i + 1);
0482                 double er2 = hist[1]->GetBinError(i + 1) / hist[1]->GetBinContent(i + 1);
0483                 dyp[npt] = ypt[npt] * sqrt(er1 * er1 + er2 * er2);
0484                 xpt[npt] = xmin + (i - 0.5) * dx;
0485                 dxp[npt] = 0;
0486                 double temp1 = (ypt[npt] > 1.0) ? 1.0 / ypt[npt] : ypt[npt];
0487                 double temp2 = (ypt[npt] > 1.0) ? dyp[npt] / (ypt[npt] * ypt[npt]) : dyp[npt];
0488                 sumNum += (fabs(1 - temp1) / (temp2 * temp2));
0489                 sumDen += (1.0 / (temp2 * temp2));
0490                 ++npt;
0491               }
0492             }
0493             sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0;
0494             sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0;
0495             TGraphAsymmErrors *graph = new TGraphAsymmErrors(npt, xpt, ypt, dxp, dxp, dyp, dyp);
0496             graph->SetMarkerStyle(24);
0497             graph->SetMarkerColor(1);
0498             graph->SetMarkerSize(0.8);
0499             graph->SetLineColor(1);
0500             graph->SetLineWidth(2);
0501             sprintf(name, "%sRatio", pad->GetName());
0502             TCanvas *canvas = new TCanvas(name, name, 500, 400);
0503             gStyle->SetOptStat(0);
0504             gPad->SetTopMargin(0.05);
0505             gPad->SetLeftMargin(0.15);
0506             gPad->SetRightMargin(0.025);
0507             gPad->SetBottomMargin(0.20);
0508             TH1F *vFrame = canvas->DrawFrame(xmin, 0.01, xmax, 0.5);
0509             vFrame->GetYaxis()->SetRangeUser(0.4, 1.5);
0510             vFrame->GetXaxis()->SetLabelSize(0.035);
0511             vFrame->GetYaxis()->SetLabelSize(0.04);
0512             vFrame->GetXaxis()->SetTitleSize(0.045);
0513             vFrame->GetYaxis()->SetTitleSize(0.045);
0514             vFrame->GetYaxis()->SetTitleOffset(1.2);
0515             vFrame->GetXaxis()->SetRangeUser(xmin, xmax);
0516             vFrame->GetYaxis()->SetTitle("Geant4/GeantV");
0517             sprintf(name, "%s in %s", xtitle[i2].c_str(), detName[i1].c_str());
0518             vFrame->GetXaxis()->SetTitle(name);
0519             graph->Draw("P");
0520             TLine *line = new TLine(xmin, 1.0, xmax, 1.0);
0521             line->SetLineStyle(2);
0522             line->SetLineWidth(2);
0523             line->SetLineColor(kRed);
0524             line->Draw();
0525             sprintf(title, "Mean Deviation = %5.3f #pm %5.3f", sumNum, sumDen);
0526             TPaveText *text = new TPaveText(0.16, 0.85, 0.60, 0.90, "brNDC");
0527             text->AddText(title);
0528             text->Draw("same");
0529             canvas->Modified();
0530             canvas->Update();
0531             if (save != 0) {
0532               if (save > 0)
0533                 sprintf(name, "c_%s.pdf", canvas->GetName());
0534               else
0535                 sprintf(name, "c_%s.jpg", canvas->GetName());
0536               canvas->Print(name);
0537             }
0538           }
0539         }
0540       }
0541     }
0542   }
0543 }
0544
0545 void makeGVSPlots(std::string fnmG4 = "analG4.root",
0546                   std::string fnmGV = "analGV.root",
0547                   int todomin = 0,
0548                   int todomax = 7,
0549                   std::string tag = "",
0550                   std::string text = "",
0551                   bool save = false,
0552                   std::string dirnm = "caloSimHitAnalysis") {
0553   std::string names[8] = {"hitp", "trackp", "stepp", "edepp", "timep", "volp", "edept", "timet"};
0554   std::string xtitle[8] = {"Hits",
0555                            "Tracks",
0556                            "Step Length (cm)",
0557                            "Energy Deposit (MeV)",
0558                            "Time (ns)",
0559                            "Volume elements",
0560                            "Energy Deposit (MeV)",
0561                            "Time (ns)"};
0562   std::string ytitle[8] = {"Events", "Events", "Hits", "Hits", "Hits", "Events", "Hits", "Hits"};
0563   std::string detName[6] = {"Barrel Pixel", "Forward Pixel", "TIB", "TID", "TOB", "TEC"};
0564   std::string particles[3] = {"Electrons", "Positrons", "Photons"};
0565   int boxp[8] = {0, 1, 0, 0, 0, 0, 0, 0};
0566   int nhis[8] = {-4, -4, -4, 1, 1, 1, 6, 6};
0567
0568   gStyle->SetCanvasBorderMode(0);
0569   gStyle->SetCanvasColor(kWhite);
0570   gStyle->SetPadColor(kWhite);
0571   gStyle->SetFillColor(kWhite);
0572   gStyle->SetOptStat(111110);
0573   TFile *file1 = new TFile(fnmG4.c_str());
0574   TFile *file2 = new TFile(fnmGV.c_str());
0575   if (file1 && file2) {
0576     TDirectory *dir1 = (TDirectory *)file1->FindObjectAny(dirnm.c_str());
0577     TDirectory *dir2 = (TDirectory *)file2->FindObjectAny(dirnm.c_str());
0578     char name[100], cname[100], title[100];
0579     for (int i2 = todomin; i2 <= todomax; ++i2) {
0580       int nhism = (nhis[i2] >= 0) ? nhis[i2] : -nhis[i2];
0581       for (int i1 = 0; i1 < nhism; ++i1) {
0582         if (nhis[i2] <= 1 && i1 == 0) {
0583           sprintf(name, "%s", names[i2].c_str());
0584           sprintf(cname, "%s%s", names[i2].c_str(), tag.c_str());
0585           sprintf(title, "%s (Geant4 vs GeantV)", text.c_str());
0586         } else if (nhis[i2] < 0) {
0587           sprintf(name, "%s%d", names[i2].c_str(), i1 - 1);
0588           sprintf(cname, "%s%d%s", names[i2].c_str(), i1 - 1, tag.c_str());
0589           sprintf(title, "%s in %s (Geant4 vs GeantV)", particles[i1 - 1].c_str(), text.c_str());
0590         } else {
0591           sprintf(name, "%s%d", names[i2].c_str(), i1);
0592           sprintf(cname, "%s%d%s", names[i2].c_str(), i1, tag.c_str());
0593           sprintf(title, "%s in %s (Geant4 vs GeantV)", text.c_str(), detName[i1].c_str());
0594         }
0595         TH1D *hist[2];
0596         hist[0] = (TH1D *)dir1->FindObjectAny(name);
0597         hist[1] = (TH1D *)dir2->FindObjectAny(name);
0598         if ((hist[0] != nullptr) && (hist[1] != nullptr)) {
0599           // Plot superimposed histograms
0600           TCanvas *pad = new TCanvas(cname, cname, 500, 500);
0601           TLegend *legend = new TLegend(0.44, 0.78, 0.64, 0.89);
0602           pad->SetRightMargin(0.10);
0603           pad->SetTopMargin(0.10);
0604           pad->SetLogy();
0605           pad->SetFillColor(kWhite);
0606           legend->SetFillColor(kWhite);
0607           int icol[2] = {1, 2};
0608           int isty[2] = {1, 2};
0609           int imty[2] = {20, 24};
0610           std::string type[2] = {"Geant4", "GeantV"};
0611           double ymax(0.90);
0612           double total[2] = {0, 0};
0613           double ymaxv[2] = {0, 0};
0614           for (int i = 0; i < 2; ++i) {
0615             hist[i]->GetYaxis()->SetTitleOffset(1.2);
0616             hist[i]->GetYaxis()->SetTitle(ytitle[i2].c_str());
0617             hist[i]->GetXaxis()->SetTitle(xtitle[i2].c_str());
0618             hist[i]->SetTitle(title);
0619             hist[i]->SetMarkerStyle(imty[i]);
0620             hist[i]->SetMarkerColor(icol[i]);
0621             hist[i]->SetLineColor(icol[i]);
0622             hist[i]->SetLineStyle(isty[i]);
0623             hist[i]->SetNdivisions(505, "X");
0624             total[i] = hist[i]->GetEntries();
0625             legend->AddEntry(hist[i], type[i].c_str(), "lp");
0626             ymaxv[i] = hist[i]->GetMaximum();
0627           }
0628           int first = (ymaxv[0] > ymaxv[1]) ? 0 : 1;
0629           int next = 1 - first;
0630           hist[first]->Draw();
0631           hist[next]->Draw("sames");
0632           pad->Update();
0633           legend->Draw("same");
0634           pad->Modified();
0635           pad->Update();
0636           for (int i = 0; i < 2; ++i) {
0637             TPaveStats *st = (TPaveStats *)hist[i]->GetListOfFunctions()->FindObject("stats");
0638             if (st != NULL) {
0639               double xl = (boxp[i2] == 0) ? 0.65 : 0.10;
0640               st->SetLineColor(icol[i]);
0641               st->SetTextColor(icol[i]);
0642               st->SetY1NDC(ymax - 0.15);
0643               st->SetY2NDC(ymax);
0644               st->SetX1NDC(xl);
0645               st->SetX2NDC(xl + 0.25);
0646               ymax -= 0.15;
0647             }
0648           }
0649           pad->Modified();
0650           pad->Update();
0651           if (save) {
0652             sprintf(name, "c_%s.jpg", pad->GetName());
0653             pad->Print(name);
0654           }
0655         }
0656       }
0657     }
0658   }
0659 }
0660
0661 void setTDRStyle() {
0662   TStyle *tdrStyle = new TStyle("tdrStyle", "Style for P-TDR");
0663
0664   // For the canvas:
0665   tdrStyle->SetCanvasBorderMode(0);
0666   tdrStyle->SetCanvasColor(kWhite);
0667   tdrStyle->SetCanvasDefH(600);  //Height of canvas
0668   tdrStyle->SetCanvasDefW(600);  //Width of canvas
0669   tdrStyle->SetCanvasDefX(0);    //POsition on screen
0670   tdrStyle->SetCanvasDefY(0);
0671
0672   // For the Pad:
0673   tdrStyle->SetPadBorderMode(0);
0674   tdrStyle->SetPadColor(kWhite);
0675   tdrStyle->SetPadGridX(false);
0676   tdrStyle->SetPadGridY(false);
0677   tdrStyle->SetGridColor(0);
0678   tdrStyle->SetGridStyle(3);
0679   tdrStyle->SetGridWidth(1);
0680
0681   // For the frame:
0682   tdrStyle->SetFrameBorderMode(0);
0683   tdrStyle->SetFrameBorderSize(1);
0684   tdrStyle->SetFrameFillColor(0);
0685   tdrStyle->SetFrameFillStyle(0);
0686   tdrStyle->SetFrameLineColor(1);
0687   tdrStyle->SetFrameLineStyle(1);
0688   tdrStyle->SetFrameLineWidth(1);
0689
0690   // For the histo:
0691   tdrStyle->SetHistLineColor(1);
0692   tdrStyle->SetHistLineStyle(0);
0693   tdrStyle->SetHistLineWidth(1);
0694   tdrStyle->SetEndErrorSize(2);
0695
0696   tdrStyle->SetMarkerStyle(20);
0697
0698   //For the fit/function:
0699   tdrStyle->SetOptFit(1);
0700   tdrStyle->SetFitFormat("5.4g");
0701   tdrStyle->SetFuncColor(2);
0702   tdrStyle->SetFuncStyle(1);
0703   tdrStyle->SetFuncWidth(1);
0704
0705   //For the date:
0706   tdrStyle->SetOptDate(0);
0707
0708   // For the statistics box:
0709   tdrStyle->SetOptFile(0);
0710   tdrStyle->SetOptStat(0);  // To display the mean and RMS:   SetOptStat("mr");
0711   tdrStyle->SetStatColor(kWhite);
0712   tdrStyle->SetStatFont(42);
0713   tdrStyle->SetStatFontSize(0.025);
0714   tdrStyle->SetStatTextColor(1);
0715   tdrStyle->SetStatFormat("6.4g");
0716   tdrStyle->SetStatBorderSize(1);
0717   tdrStyle->SetStatH(0.1);
0718   tdrStyle->SetStatW(0.15);
0719
0720   // Margins:
0721   tdrStyle->SetPadTopMargin(0.05);
0722   tdrStyle->SetPadBottomMargin(0.13);
0723   tdrStyle->SetPadLeftMargin(0.16);
0724   tdrStyle->SetPadRightMargin(0.02);
0725
0726   // For the Global title:
0727
0728   tdrStyle->SetOptTitle(0);
0729   tdrStyle->SetTitleFont(42);
0730   tdrStyle->SetTitleColor(1);
0731   tdrStyle->SetTitleTextColor(1);
0732   tdrStyle->SetTitleFillColor(10);
0733   tdrStyle->SetTitleFontSize(0.05);
0734
0735   // For the axis titles:
0736
0737   tdrStyle->SetTitleColor(1, "XYZ");
0738   tdrStyle->SetTitleFont(42, "XYZ");
0739   tdrStyle->SetTitleSize(0.06, "XYZ");
0740   tdrStyle->SetTitleXOffset(0.9);
0741   tdrStyle->SetTitleYOffset(1.25);
0742
0743   // For the axis labels:
0744
0745   tdrStyle->SetLabelColor(1, "XYZ");
0746   tdrStyle->SetLabelFont(42, "XYZ");
0747   tdrStyle->SetLabelOffset(0.007, "XYZ");
0748   tdrStyle->SetLabelSize(0.05, "XYZ");
0749
0750   // For the axis:
0751
0752   tdrStyle->SetAxisColor(1, "XYZ");
0753   tdrStyle->SetStripDecimals(kTRUE);
0754   tdrStyle->SetTickLength(0.03, "XYZ");
0755   tdrStyle->SetNdivisions(510, "XYZ");
0756   tdrStyle->SetPadTickX(1);  // To get tick marks on the opposite side of the frame
0757   tdrStyle->SetPadTickY(1);
0758
0759   // Change for log plots:
0760   tdrStyle->SetOptLogx(0);
0761   tdrStyle->SetOptLogy(0);
0762   tdrStyle->SetOptLogz(0);
0763
0764   // Postscript options:
0765   tdrStyle->SetPaperSize(20., 20.);
0766   tdrStyle->SetOptLogy(0);
0767   tdrStyle->SetOptLogz(0);
0768
0769   // Postscript options:
0770   tdrStyle->SetPaperSize(20., 20.);
0771
0772   tdrStyle->SetHatchesLineWidth(5);
0773   tdrStyle->SetHatchesSpacing(0.05);
0774
0775   tdrStyle->cd();
0776 }