File indexing completed on 2024-04-06 12:29:48
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
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 }
