Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 10:03:55

0001 //////////////////////////////////////////////////////////////////////////////
0002 //
0003 // Usage:
0004 // .L MakeHitStudyPlots.C+g
0005 //
0006 //   To make plot of various quantities which are created by CaloSimHitStudy
0007 //   from one or two different settings
0008 //
0009 //     makeHitStudyPlots(file1, tag1, file2, tag2, todomin, todomax, gtitle,
0010 //                       ratio, save, dirnm)
0011 //
0012 //   To make plots of ratios of hits produced with identical geometry and
0013 //   generator level files created by HGCalHitCheck
0014 //
0015 //      makeDDDvsDD4hepPlots(dirnm, inType, geometry, layer, ratio, save)
0016 //
0017 //   where (for makeHitStudyPlots)
0018 //     file1    std::string   Name of the first ROOT file [old/analRun3.root]
0019 //     file2    std::string   Name of the second ROOT file [new/analRun3.root]
0020 //     tag1     std::string   Tag for the first file [Bug]
0021 //     tag2     std::string   Tag for the second file [Fix]
0022 //     gtitle   std::string   Overall Titile [""]
0023 //     todomin  int           Minimum type # of histograms to be plotted [0]
0024 //     todomax  int           Maximum type # of histograms to be plotted [0]
0025 //     dirnm    std::string   Name of the directory [CaloSimHitStudy]
0026 //
0027 //   where (for makeHitStudyPlots)
0028 //     dirnm    std::string   Directory name (EE/HEF/HEB)
0029 //     inType   std::string   Name of the input data (Muon/MinBias)
0030 //     geometry std::string   Tag for the geometry (D98/D99)
0031 //     layer    int           Layer number (if 0; all layers combined)
0032 //
0033 //   where (common to both macros)
0034 //     ratio    bool          if the ratio to be plotted [true]
0035 //                            (works when both files are active)
0036 //     save     bool          If the canvas is to be saved as jpg file [false]
0037 //
0038 //////////////////////////////////////////////////////////////////////////////
0039 
0040 #include <TCanvas.h>
0041 #include <TChain.h>
0042 #include <TFile.h>
0043 #include <TFitResult.h>
0044 #include <TFitResultPtr.h>
0045 #include <TH1D.h>
0046 #include <TH2D.h>
0047 #include <TLegend.h>
0048 #include <TPaveStats.h>
0049 #include <TPaveText.h>
0050 #include <TProfile.h>
0051 #include <TROOT.h>
0052 #include <TStyle.h>
0053 #include <vector>
0054 #include <string>
0055 #include <iomanip>
0056 #include <iostream>
0057 #include <fstream>
0058 
0059 void makeHitStudyPlots(std::string file1 = "uncorr/analRun3.root",
0060                        std::string file2 = "corr/analRun3.root",
0061                        std::string tag1 = "Bug",
0062                        std::string tag2 = "Fix",
0063                        std::string gtitle = "",
0064                        int todomin = 0,
0065                        int todomax = 0,
0066                        bool ratio = true,
0067                        bool save = false,
0068                        std::string dirnm = "CaloSimHitStudy") {
0069   const int plots = 20;
0070   std::string names[plots] = {"Etot",   "Hit",     "EtotG",  "Time",    "EdepTk", "Edep", "HitHigh",
0071                               "HitLow", "HitMu",   "HitTk",  "TimeAll", "TimeTk", "eta",  "phi",
0072                               "EdepEM", "EdepHad", "EneInc", "EtaInc",  "PhiInc", "PtInc"};
0073   int numb[plots] = {9, 9, 9, 9, 16, 9, 1, 1, 1, 16, 9, 16, 9, 9, 9, 9, 1, 1, 1, 1};
0074   int rebin[plots] = {10, 10, 10, 10, 1, 10, 10, 10, 10, 10, 10, 10, 2, 4, 10, 10, 1, 1, 1, 1};
0075   bool debug(false);
0076 
0077   gStyle->SetCanvasBorderMode(0);
0078   gStyle->SetCanvasColor(kWhite);
0079   gStyle->SetPadColor(kWhite);
0080   gStyle->SetFillColor(kWhite);
0081   if (ratio)
0082     gStyle->SetOptStat(0);
0083   else
0084     gStyle->SetOptStat(1110);
0085   TFile* file[2];
0086   int nfile(0);
0087   std::string tag(""), tags[2];
0088   if (file1 != "") {
0089     file[nfile] = new TFile(file1.c_str());
0090     if (file[nfile]) {
0091       tags[nfile] = tag1;
0092       ++nfile;
0093       tag += tag1;
0094     }
0095   }
0096   if (file2 != "") {
0097     file[nfile] = new TFile(file2.c_str());
0098     if (file[nfile]) {
0099       tags[nfile] = tag2;
0100       ++nfile;
0101       tag += tag2;
0102     }
0103   }
0104   if ((todomin < 0) || (todomin >= plots))
0105     todomin = 0;
0106   if (todomax < todomin) {
0107     todomax = todomin;
0108   } else if (todomax >= plots) {
0109     todomax = plots - 1;
0110   }
0111   std::cout << "Use " << nfile << " files from " << file1 << " and " << file2 << " and look for " << todomin << ":"
0112             << todomax << std::endl;
0113   for (int todo = todomin; todo <= todomax; ++todo) {
0114     for (int i1 = 0; i1 < numb[todo]; ++i1) {
0115       double y1(0.90), dy(0.12);
0116       double y2 = y1 - dy * nfile - 0.01;
0117       TLegend* leg = (ratio) ? (new TLegend(0.10, 0.86, 0.90, 0.90)) : (new TLegend(0.65, y2 - nfile * 0.04, 0.90, y2));
0118       leg->SetBorderSize(1);
0119       leg->SetFillColor(10);
0120       TH1D* hist0[nfile];
0121       char name[100], namec[100];
0122       if (numb[todo] == 1) {
0123         sprintf(name, "%s", names[todo].c_str());
0124         sprintf(namec, "c_%s%s%s", names[todo].c_str(), tag.c_str(), gtitle.c_str());
0125       } else {
0126         sprintf(name, "%s%d", names[todo].c_str(), i1);
0127         sprintf(namec, "c_%s%d%s%s", names[todo].c_str(), i1, tag.c_str(), gtitle.c_str());
0128       }
0129       TCanvas* pad = new TCanvas(namec, namec, 500, 500);
0130       for (int ifile = 0; ifile < nfile; ++ifile) {
0131         TDirectory* dir = (TDirectory*)file[ifile]->FindObjectAny(dirnm.c_str());
0132         hist0[ifile] = static_cast<TH1D*>(dir->FindObjectAny(name));
0133         if (debug)
0134           std::cout << name << " read out at " << hist0[ifile] << " for " << tags[ifile] << std::endl;
0135       }
0136       if (!ratio) {
0137         int first(0);
0138         for (int ifile = 0; ifile < nfile; ++ifile) {
0139           TH1D* hist(hist0[ifile]);
0140           if (hist != nullptr) {
0141             hist->SetLineColor(first + 1);
0142             hist->SetLineStyle(first + 1);
0143             hist->GetYaxis()->SetTitleOffset(1.4);
0144             if (rebin[todo] > 1)
0145               hist->Rebin(rebin[todo]);
0146             hist->SetTitle(gtitle.c_str());
0147             if (first == 0) {
0148               pad = new TCanvas(namec, namec, 500, 500);
0149               pad->SetRightMargin(0.10);
0150               pad->SetTopMargin(0.10);
0151               pad->SetLogy();
0152               hist->Draw();
0153             } else {
0154               hist->Draw("sames");
0155             }
0156             leg->AddEntry(hist, tags[ifile].c_str(), "lp");
0157             pad->Update();
0158             ++first;
0159             TPaveStats* st = ((TPaveStats*)hist->GetListOfFunctions()->FindObject("stats"));
0160             if (st != NULL) {
0161               st->SetLineColor(first);
0162               st->SetTextColor(first);
0163               st->SetY1NDC(y1 - dy);
0164               st->SetY2NDC(y1);
0165               st->SetX1NDC(0.65);
0166               st->SetX2NDC(0.90);
0167               y1 -= dy;
0168             }
0169             pad->Modified();
0170             pad->Update();
0171             leg->Draw("same");
0172             pad->Update();
0173             if (save) {
0174               sprintf(name, "%s.pdf", pad->GetName());
0175               pad->Print(name);
0176             }
0177           }
0178         }
0179       } else {
0180         if (nfile == 2) {
0181           int nbin = hist0[0]->GetNbinsX();
0182           int nbinR = nbin / rebin[todo];
0183           double xlow = hist0[0]->GetXaxis()->GetBinLowEdge(1);
0184           double xhigh = hist0[0]->GetXaxis()->GetBinLowEdge(nbin) + hist0[0]->GetXaxis()->GetBinWidth(nbin);
0185           if (numb[todo] == 1) {
0186             sprintf(name, "%sRatio", names[todo].c_str());
0187             sprintf(namec, "c_%sRatio%s%s", names[todo].c_str(), tag.c_str(), gtitle.c_str());
0188           } else {
0189             sprintf(name, "%s%dRatio", names[todo].c_str(), i1);
0190             sprintf(namec, "c_%s%dRatio%s%s", names[todo].c_str(), i1, tag.c_str(), gtitle.c_str());
0191           }
0192           pad = new TCanvas(namec, namec, 500, 500);
0193           TH1D* histr = new TH1D(name, hist0[0]->GetTitle(), nbinR, xlow, xhigh);
0194           sprintf(name, "Ratio (%s/%s)", tags[0].c_str(), tags[1].c_str());
0195           histr->SetTitle(gtitle.c_str());
0196           histr->GetXaxis()->SetTitle(hist0[0]->GetXaxis()->GetTitle());
0197           histr->GetYaxis()->SetTitle(name);
0198           histr->GetXaxis()->SetLabelOffset(0.005);
0199           histr->GetXaxis()->SetTitleOffset(1.00);
0200           histr->GetYaxis()->SetLabelOffset(0.005);
0201           histr->GetYaxis()->SetTitleOffset(1.20);
0202           histr->GetYaxis()->SetRangeUser(0.0, 2.0);
0203           double sumNum(0), sumDen(0), maxDev(0);
0204           for (int j = 0; j < nbinR; ++j) {
0205             double tnum(0), tden(0), rnum(0), rden(0);
0206             for (int i = 0; i < rebin[todo]; ++i) {
0207               int ib = j * rebin[todo] + i + 1;
0208               tnum += hist0[0]->GetBinContent(ib);
0209               tden += hist0[1]->GetBinContent(ib);
0210               rnum += ((hist0[0]->GetBinError(ib)) * (hist0[0]->GetBinError(ib)));
0211               rden += ((hist0[1]->GetBinError(ib)) * (hist0[1]->GetBinError(ib)));
0212             }
0213             if (tden > 0 && tnum > 0) {
0214               double rat = tnum / tden;
0215               double err = rat * sqrt((rnum / (tnum * tnum)) + (rden / (tden * tden)));
0216               histr->SetBinContent(j + 1, rat);
0217               histr->SetBinError(j + 1, err);
0218               double temp1 = (rat > 1.0) ? 1.0 / rat : rat;
0219               double temp2 = (rat > 1.0) ? err / (rat * rat) : err;
0220               sumNum += (fabs(1 - temp1) / (temp2 * temp2));
0221               sumDen += (1.0 / (temp2 * temp2));
0222               if (fabs(1 - temp1) > maxDev)
0223                 maxDev = fabs(1 - temp1);
0224             }
0225           }
0226           histr->Draw();
0227           sprintf(name, "%s vs %s", tag1.c_str(), tag2.c_str());
0228           leg->AddEntry(histr, name, "lp");
0229           leg->Draw("same");
0230           pad->Update();
0231           TLine* line = new TLine(xlow, 1.0, xhigh, 1.0);
0232           line->SetLineColor(2);
0233           line->SetLineWidth(2);
0234           line->SetLineStyle(2);
0235           line->Draw("same");
0236           pad->Modified();
0237           pad->Update();
0238           sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0;
0239           sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0;
0240           if (sumNum == 0)
0241             sumDen = 0;
0242           std::cout << tag1 << " vs " << tag2 << " " << hist0[0]->GetXaxis()->GetTitle() << " Mean deviation " << sumNum
0243                     << " +- " << sumDen << " maximum " << maxDev << std::endl;
0244           if (save) {
0245             sprintf(name, "%s.pdf", pad->GetName());
0246             pad->Print(name);
0247           }
0248         }
0249       }
0250     }
0251   }
0252 }
0253 
0254 void makeDDDvsDD4hepPlots(std::string dirnm = "EE",
0255                           std::string inType = "Muon",
0256                           std::string geometry = "D98",
0257                           int layer = 0,
0258                           bool ratio = true,
0259                           bool save = false) {
0260   const int plots = 3;
0261   std::string types[2] = {"DDD", "DD4hep"};
0262   std::string plotf[plots] = {"L", "F", "P"};
0263   std::string plotp[plots] = {"All", "Full|SiPM 2mm", "Partial|SiPM 4mm"};
0264   int rebins[2] = {4, 20};
0265   double xmaxs[2] = {1000, 5000};
0266   bool debug(false);
0267 
0268   gStyle->SetCanvasBorderMode(0);
0269   gStyle->SetCanvasColor(kWhite);
0270   gStyle->SetPadColor(kWhite);
0271   gStyle->SetFillColor(kWhite);
0272   if (ratio)
0273     gStyle->SetOptStat(0);
0274   else
0275     gStyle->SetOptStat(1110);
0276   TFile* file[2];
0277   int nfile(0);
0278   std::string tag(""), tags[2], filex[2];
0279   for (int i = 0; i < 2; ++i) {
0280     filex[i] = types[i] + geometry + inType + ".root";
0281     file[nfile] = new TFile(filex[i].c_str());
0282     if (file[nfile]) {
0283       tags[nfile] = types[i];
0284       ++nfile;
0285       tag += tags[nfile];
0286     }
0287   }
0288   char name[80], nameD[80], title[80];
0289   int rebin = (inType == "Muon") ? rebins[0] : rebins[1];
0290   double xmax = (inType == "Muon") ? xmaxs[0] : xmaxs[1];
0291   sprintf(nameD, "hgcalHitCheck%s", dirnm.c_str());
0292   sprintf(title, "%s vs %s for %s", types[0].c_str(), types[1].c_str(), inType.c_str());
0293   std::cout << "Use " << nfile << " files from " << filex[0] << " and " << filex[1] << " and look for " << plots
0294             << " plots in " << nameD << " with rebin " << rebin << " Max " << xmax << std::endl;
0295   for (int i = 0; i < plots; ++i) {
0296     if (layer == 0)
0297       sprintf(name, "Hits%s", plotf[i].c_str());
0298     else
0299       sprintf(name, "Hits%s%d", plotf[i].c_str(), layer);
0300     double y1(0.90), dy(0.12);
0301     double y2 = y1 - dy * nfile - 0.01;
0302     TLegend* leg = (ratio) ? (new TLegend(0.40, 0.86, 0.90, 0.90)) : (new TLegend(0.65, y2 - nfile * 0.04, 0.90, y2));
0303     leg->SetBorderSize(1);
0304     leg->SetFillColor(10);
0305     TH1D* hist0[nfile];
0306     for (int ifile = 0; ifile < nfile; ++ifile) {
0307       TDirectory* dir = (TDirectory*)file[ifile]->FindObjectAny(nameD);
0308       hist0[ifile] = static_cast<TH1D*>(dir->FindObjectAny(name));
0309       if (debug)
0310         std::cout << name << " read out at " << hist0[ifile] << " for " << tags[ifile] << std::endl;
0311     }
0312     char namec[160];
0313     if (!ratio) {
0314       sprintf(namec, "c_%s%s%s%s", geometry.c_str(), inType.c_str(), dirnm.c_str(), name);
0315       TCanvas* pad;
0316       int first(0);
0317       for (int ifile = 0; ifile < nfile; ++ifile) {
0318         TH1D* hist(hist0[ifile]);
0319         if (hist != nullptr) {
0320           if (rebin > 1)
0321             hist->Rebin(rebin);
0322           hist->SetTitle(title);
0323           hist->SetLineColor(first + 1);
0324           hist->SetLineStyle(first + 1);
0325           hist->GetYaxis()->SetTitleOffset(1.4);
0326           hist->GetXaxis()->SetRangeUser(0, xmax);
0327           hist->GetXaxis()->SetTitleSize(0.025);
0328           hist->GetXaxis()->SetTitleOffset(1.2);
0329           hist->SetMarkerStyle(first + 20);
0330           hist->SetMarkerColor(first + 1);
0331           hist->SetMarkerSize(0.7);
0332           if (first == 0) {
0333             pad = new TCanvas(namec, namec, 500, 500);
0334             pad->SetRightMargin(0.10);
0335             pad->SetTopMargin(0.10);
0336             pad->SetLogy();
0337             hist->Draw();
0338           } else {
0339             hist->Draw("sames");
0340           }
0341           leg->AddEntry(hist, tags[ifile].c_str(), "lp");
0342           pad->Update();
0343           ++first;
0344           TPaveStats* st = ((TPaveStats*)hist->GetListOfFunctions()->FindObject("stats"));
0345           if (st != NULL) {
0346             st->SetLineColor(first);
0347             st->SetTextColor(first);
0348             st->SetY1NDC(y1 - dy);
0349             st->SetY2NDC(y1);
0350             st->SetX1NDC(0.65);
0351             st->SetX2NDC(0.90);
0352             y1 -= dy;
0353           }
0354           pad->Modified();
0355           pad->Update();
0356           leg->Draw("same");
0357           pad->Update();
0358         }
0359         if (save) {
0360           sprintf(name, "%s.pdf", pad->GetName());
0361           pad->Print(name);
0362         }
0363       }
0364     } else if (nfile == 2) {
0365       sprintf(namec, "cR_%s%s%s%s", geometry.c_str(), inType.c_str(), dirnm.c_str(), name);
0366       TCanvas* pad = new TCanvas(namec, namec, 500, 500);
0367       int nbin = hist0[0]->GetNbinsX();
0368       int nbinR = nbin / rebin;
0369       double xlow = hist0[0]->GetXaxis()->GetBinLowEdge(1);
0370       double xhigh = xmax;
0371       TH1D* histr = new TH1D(name, hist0[0]->GetTitle(), nbinR, xlow, xhigh);
0372       histr->SetTitle(title);
0373       if (layer == 0)
0374         sprintf(name, "Number of hits (%s)", plotp[i].c_str());
0375       else
0376         sprintf(name, "Number of hits in Layer %d (%s)", layer, plotp[i].c_str());
0377       histr->GetXaxis()->SetTitle(name);
0378       sprintf(name, "Ratio (%s/%s)", tags[0].c_str(), tags[1].c_str());
0379       histr->GetYaxis()->SetTitle(name);
0380       histr->GetXaxis()->SetLabelOffset(0.005);
0381       histr->GetXaxis()->SetTitleOffset(1.30);
0382       histr->GetXaxis()->SetTitleSize(0.036);
0383       histr->GetYaxis()->SetLabelOffset(0.005);
0384       histr->GetYaxis()->SetTitleOffset(1.20);
0385       histr->GetYaxis()->SetTitleSize(0.036);
0386       histr->GetYaxis()->SetRangeUser(0.0, 5.0);
0387       histr->SetMarkerStyle(20);
0388       histr->SetMarkerColor(1);
0389       histr->SetMarkerSize(0.7);
0390       double sumNum(0), sumDen(0), maxDev(0);
0391       for (int j = 0; j < nbinR; ++j) {
0392         double tnum(0), tden(0), rnum(0), rden(0);
0393         for (int i = 0; i < rebin; ++i) {
0394           int ib = j * rebin + i + 1;
0395           tnum += hist0[0]->GetBinContent(ib);
0396           tden += hist0[1]->GetBinContent(ib);
0397           rnum += ((hist0[0]->GetBinError(ib)) * (hist0[0]->GetBinError(ib)));
0398           rden += ((hist0[1]->GetBinError(ib)) * (hist0[1]->GetBinError(ib)));
0399         }
0400         if (tden > 0 && tnum > 0) {
0401           double rat = tnum / tden;
0402           double err = rat * sqrt((rnum / (tnum * tnum)) + (rden / (tden * tden)));
0403           histr->SetBinContent(j + 1, rat);
0404           histr->SetBinError(j + 1, err);
0405           double temp1 = (rat > 1.0) ? 1.0 / rat : rat;
0406           double temp2 = (rat > 1.0) ? err / (rat * rat) : err;
0407           sumNum += (fabs(1 - temp1) / (temp2 * temp2));
0408           sumDen += (1.0 / (temp2 * temp2));
0409           if (fabs(1 - temp1) > maxDev)
0410             maxDev = fabs(1 - temp1);
0411         }
0412       }
0413       histr->Draw();
0414       if (layer == 0)
0415         sprintf(name, "%s %s", dirnm.c_str(), plotp[i].c_str());
0416       else
0417         sprintf(name, "%s (Layer %d) %s", dirnm.c_str(), layer, plotp[i].c_str());
0418       leg->AddEntry(histr, name, "lp");
0419       leg->Draw("same");
0420       pad->Update();
0421       TLine* line = new TLine(xlow, 1.0, xhigh, 1.0);
0422       line->SetLineColor(2);
0423       line->SetLineWidth(2);
0424       line->SetLineStyle(2);
0425       line->Draw("same");
0426       pad->Modified();
0427       pad->Update();
0428       sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0;
0429       sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0;
0430       if (sumNum == 0)
0431         sumDen = 0;
0432       std::cout << tags[0] << " vs " << tags[1] << " " << hist0[0]->GetXaxis()->GetTitle() << " Mean deviation "
0433                 << sumNum << " +- " << sumDen << " maximum " << maxDev << std::endl;
0434       if (save) {
0435         sprintf(name, "%s.pdf", pad->GetName());
0436         pad->Print(name);
0437       }
0438     }
0439   }
0440 }