Calo/macros/MakeHitStudyPlots.C

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, toDo, ratio, save, dirnm)
0010 //
0011 //   where
0012 //     file1   std::string   Name of the first ROOT file [analRun3Old.root]
0013 //     tag1    std::string   Tag for the first file [Old]
0014 //     file2   std::string   Name of the second ROOT file [analRun3New.root]
0015 //     tag2    std::string   Tag for the second file [New]
0016 //     todo    int           The plot type to be made [0]
0017 //                           if -1, 6 different types are plotted
0018 //                           (3, 5, 8, 9, 10, 11)
0019 //     ratio   bool          if the ratio to be plotted [true]
0020 //                           (works when both files are active)
0021 //     save    bool          If the canvas is to be saved as jpg file [false]
0022 //     dirnm   std::string   Name of the directory [CaloSimHitStudy]
0023 //
0024 //////////////////////////////////////////////////////////////////////////////
0025
0026 #include <TCanvas.h>
0027 #include <TChain.h>
0028 #include <TFile.h>
0029 #include <TFitResult.h>
0030 #include <TFitResultPtr.h>
0031 #include <TH1D.h>
0032 #include <TH2D.h>
0033 #include <TLegend.h>
0034 #include <TPaveStats.h>
0035 #include <TPaveText.h>
0036 #include <TProfile.h>
0037 #include <TROOT.h>
0038 #include <TStyle.h>
0039 #include <vector>
0040 #include <string>
0041 #include <iomanip>
0042 #include <iostream>
0043 #include <fstream>
0044
0045 void makeHitStudyPlots(std::string file1 = "analRun3Old.root",
0046                        std::string tag1 = "Old",
0047                        std::string file2 = "analRun3New.root",
0048                        std::string tag2 = "New",
0049                        int toDo = 0,
0050                        bool ratio = true,
0051                        bool save = false,
0052                        std::string dirnm = "CaloSimHitStudy") {
0053   std::string names[18] = {"Edep",
0054                            "EdepEM",
0055                            "EdepHad",
0056                            "EdepTk",
0057                            "Etot",
0058                            "EtotG",
0059                            "Hit",
0060                            "HitHigh",
0061                            "HitLow",
0062                            "HitMu",
0063                            "HitTk",
0064                            "Time",
0065                            "TimeAll",
0066                            "TimeTk",
0067                            "EneInc",
0068                            "EtaInc",
0069                            "PhiInc",
0070                            "PtInc"};
0071   int numb[18] = {9, 9, 9, 16, 9, 9, 9, 1, 1, 1, 16, 9, 9, 16, 1, 1, 1, 1};
0072   int rebin[18] = {10, 10, 10, 1, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 1, 1, 1, 1};
0073   int todos[6] = {3, 5, 8, 9, 10, 11};
0074   bool debug(false);
0075
0076   gStyle->SetCanvasBorderMode(0);
0077   gStyle->SetCanvasColor(kWhite);
0078   gStyle->SetPadColor(kWhite);
0079   gStyle->SetFillColor(kWhite);
0080   if (ratio)
0081     gStyle->SetOptStat(0);
0082   else
0083     gStyle->SetOptStat(1110);
0084   TFile* file[2];
0085   int nfile(0);
0086   std::string tag(""), tags[2];
0087   if (file1 != "") {
0088     file[nfile] = new TFile(file1.c_str());
0089     if (file[nfile]) {
0090       tags[nfile] = tag1;
0091       ++nfile;
0092       tag += tag1;
0093     }
0094   }
0095   if (file2 != "") {
0096     file[nfile] = new TFile(file2.c_str());
0097     if (file[nfile]) {
0098       tags[nfile] = tag2;
0099       ++nfile;
0100       tag += tag2;
0101     }
0102   }
0103   int todoMin = (toDo >= 0) ? 0 : 0;
0104   int todoMax = (toDo >= 0) ? 0 : 5;
0105   std::cout << "Use " << nfile << " files from " << file1 << " and " << file2 << " and look for " << todoMin << ":"
0106             << todoMax << std::endl;
0107   for (int i0 = todoMin; i0 <= todoMax; ++i0) {
0108     int todo = (todoMax == 0) ? toDo : todos[i0];
0109     for (int i1 = 0; i1 < numb[todo]; ++i1) {
0110       double y1(0.90), dy(0.12);
0111       double y2 = y1 - dy * nfile - 0.01;
0112       TLegend* leg = (ratio) ? (new TLegend(0.10, 0.86, 0.90, 0.90)) : (new TLegend(0.65, y2 - nfile * 0.04, 0.90, y2));
0113       leg->SetBorderSize(1);
0114       leg->SetFillColor(10);
0115       TCanvas* pad;
0116       TH1D* hist0[nfile];
0117       char name[100], namec[100];
0118       for (int ifile = 0; ifile < nfile; ++ifile) {
0119         TDirectory* dir = (TDirectory*)file[ifile]->FindObjectAny(dirnm.c_str());
0120         if (numb[todo] == 1) {
0121           sprintf(name, "%s", names[todo].c_str());
0122           sprintf(namec, "%s%s", names[todo].c_str(), tag.c_str());
0123         } else {
0124           sprintf(name, "%s%d", names[todo].c_str(), i1);
0125           sprintf(namec, "%s%d%s", names[todo].c_str(), i1, tag.c_str());
0126         }
0127         hist0[ifile] = static_cast<TH1D*>(dir->FindObjectAny(name));
0128         if (debug)
0129           std::cout << name << " read out at " << hist0[ifile] << " for " << tags[ifile] << std::endl;
0130       }
0131       if (!ratio) {
0132         int first(0);
0133         for (int ifile = 0; ifile < nfile; ++ifile) {
0134           TH1D* hist(hist0[ifile]);
0135           if (hist != nullptr) {
0136             hist->SetLineColor(first + 1);
0137             hist->SetLineStyle(first + 1);
0138             hist->GetYaxis()->SetTitleOffset(1.4);
0139             if (rebin[todo] > 1)
0140               hist->Rebin(rebin[todo]);
0141             if (first == 0) {
0142               pad = new TCanvas(namec, namec, 500, 500);
0143               pad->SetRightMargin(0.10);
0144               pad->SetTopMargin(0.10);
0145               pad->SetLogy();
0146               hist->Draw();
0147             } else {
0148               hist->Draw("sames");
0149             }
0150             leg->AddEntry(hist, tags[ifile].c_str(), "lp");
0151             pad->Update();
0152             ++first;
0153             TPaveStats* st = ((TPaveStats*)hist->GetListOfFunctions()->FindObject("stats"));
0154             if (st != NULL) {
0155               st->SetLineColor(first);
0156               st->SetTextColor(first);
0157               st->SetY1NDC(y1 - dy);
0158               st->SetY2NDC(y1);
0159               st->SetX1NDC(0.65);
0160               st->SetX2NDC(0.90);
0161               y1 -= dy;
0162             }
0163             pad->Modified();
0164             pad->Update();
0165             leg->Draw("same");
0166             pad->Update();
0167             if (save) {
0168               sprintf(name, "c_%s.pdf", pad->GetName());
0169               pad->Print(name);
0170             }
0171           }
0172         }
0173       } else {
0174         if (nfile == 2) {
0175           int nbin = hist0[0]->GetNbinsX();
0176           int nbinR = nbin / rebin[todo];
0177           double xlow = hist0[0]->GetXaxis()->GetBinLowEdge(1);
0178           double xhigh = hist0[0]->GetXaxis()->GetBinLowEdge(nbin) + hist0[0]->GetXaxis()->GetBinWidth(nbin);
0179           ;
0180           if (numb[todo] == 1) {
0181             sprintf(name, "%sRatio", names[todo].c_str());
0182             sprintf(namec, "%sRatio%s", names[todo].c_str(), tag.c_str());
0183           } else {
0184             sprintf(name, "%s%dRatio", names[todo].c_str(), i1);
0185             sprintf(namec, "%s%dRatio%s", names[todo].c_str(), i1, tag.c_str());
0186           }
0187           pad = new TCanvas(namec, namec, 500, 500);
0188           TH1D* histr = new TH1D(name, hist0[0]->GetTitle(), nbinR, xlow, xhigh);
0189           sprintf(name, "Ratio (%s/%s)", tags[0].c_str(), tags[1].c_str());
0190           histr->GetXaxis()->SetTitle(hist0[0]->GetXaxis()->GetTitle());
0191           histr->GetYaxis()->SetTitle(name);
0192           histr->GetXaxis()->SetLabelOffset(0.005);
0193           histr->GetXaxis()->SetTitleOffset(1.00);
0194           histr->GetYaxis()->SetLabelOffset(0.005);
0195           histr->GetYaxis()->SetTitleOffset(1.20);
0196           histr->GetYaxis()->SetRangeUser(0.0, 2.0);
0197           double sumNum(0), sumDen(0), maxDev(0);
0198           for (int j = 0; j < nbinR; ++j) {
0199             double tnum(0), tden(0), rnum(0), rden(0);
0200             for (int i = 0; i < rebin[todo]; ++i) {
0201               int ib = j * rebin[todo] + i + 1;
0202               tnum += hist0[0]->GetBinContent(ib);
0203               tden += hist0[1]->GetBinContent(ib);
0204               rnum += ((hist0[0]->GetBinError(ib)) * (hist0[0]->GetBinError(ib)));
0205               rden += ((hist0[1]->GetBinError(ib)) * (hist0[1]->GetBinError(ib)));
0206             }
0207             if (tden > 0 && tnum > 0) {
0208               double rat = tnum / tden;
0209               double err = rat * sqrt((rnum / (tnum * tnum)) + (rden / (tden * tden)));
0210               histr->SetBinContent(j + 1, rat);
0211               histr->SetBinError(j + 1, err);
0212               double temp1 = (rat > 1.0) ? 1.0 / rat : rat;
0213               double temp2 = (rat > 1.0) ? err / (rat * rat) : err;
0214               sumNum += (fabs(1 - temp1) / (temp2 * temp2));
0215               sumDen += (1.0 / (temp2 * temp2));
0216               if (fabs(1 - temp1) > maxDev)
0217                 maxDev = fabs(1 - temp1);
0218             }
0219           }
0220           histr->Draw();
0221           sprintf(name, "%s vs %s", tag1.c_str(), tag2.c_str());
0222           leg->AddEntry(histr, name, "lp");
0223           leg->Draw("same");
0224           pad->Update();
0225           TLine* line = new TLine(xlow, 1.0, xhigh, 1.0);
0226           line->SetLineColor(2);
0227           line->SetLineWidth(2);
0228           line->SetLineStyle(2);
0229           line->Draw("same");
0230           pad->Modified();
0231           pad->Update();
0232           sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0;
0233           sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0;
0234           if (sumNum == 0)
0235             sumDen = 0;
0236           std::cout << tag1 << " vs " << tag2 << " " << hist0[0]->GetXaxis()->GetTitle() << " Mean deviation " << sumNum
0237                     << " +- " << sumDen << " maximum " << maxDev << std::endl;
0238         }
0239       }
0240     }
0241   }
0242 }