File indexing completed on 2023-03-17 11:27:23
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
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
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
0878 for (int ibx = 0; ibx < nbinX; ibx++) {
0879 double contb = 0;
0880
0881 for (int iby = 0; iby < ngroup; iby++) {
0882 int ibin = iby + istart;
0883 double cont = prof[ii]->GetBinContent(ibx + 1, ibin + 1);
0884
0885 contb += cont;
0886 }
0887 contb /= ngroup;
0888
0889 hist[ii][is]->SetBinContent(ibx + 1, contb);
0890 }
0891
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
1288
1289
1290
1291
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
1327
1328
1329
1330
1331
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
1717
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 }