File indexing completed on 2023-01-25 23:30:02
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
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138
0139
0140
0141
0142
0143
0144
0145
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164
0165
0166
0167
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189
0190
0191
0192
0193
0194
0195
0196
0197 #include <TCanvas.h>
0198 #include <TChain.h>
0199 #include <TProfile.h>
0200 #include <TF1.h>
0201 #include <TFile.h>
0202 #include <TFitResult.h>
0203 #include <TFitResultPtr.h>
0204 #include <TH1D.h>
0205 #include <TLegend.h>
0206 #include <TGraph.h>
0207 #include <TGraphErrors.h>
0208 #include <TGraphAsymmErrors.h>
0209 #include <TMath.h>
0210 #include <TPaveStats.h>
0211 #include <TPaveText.h>
0212 #include <TROOT.h>
0213 #include <TStyle.h>
0214 #include <cstdlib>
0215 #include <fstream>
0216 #include <iomanip>
0217 #include <iostream>
0218 #include <string>
0219 #include <vector>
0220
0221 #include "CalibCorr.C"
0222
0223 const double fitrangeFactor = 1.5;
0224
0225 struct cfactors {
0226 int ieta, depth;
0227 double corrf, dcorr;
0228 cfactors(int ie = 0, int dp = 0, double cor = 1, double dc = 0) : ieta(ie), depth(dp), corrf(cor), dcorr(dc){};
0229 };
0230
0231 struct results {
0232 double mean, errmean, width, errwidth;
0233 results(double v1 = 0, double er1 = 0, double v2 = 0, double er2 = 0)
0234 : mean(v1), errmean(er1), width(v2), errwidth(er2){};
0235 };
0236
0237 std::pair<double, double> GetMean(TH1D* hist, double xmin, double xmax, double& rms) {
0238 double mean(0), err(0), wt(0);
0239 rms = 0;
0240 for (int i = 1; i <= hist->GetNbinsX(); ++i) {
0241 double xlow = hist->GetBinLowEdge(i);
0242 double xhigh = xlow + hist->GetBinWidth(i);
0243 if ((xlow >= xmin) && (xhigh <= xmax)) {
0244 double cont = hist->GetBinContent(i);
0245 double valu = 0.5 * (xlow + xhigh);
0246 wt += cont;
0247 mean += (valu * cont);
0248 rms += (valu * valu * cont);
0249 }
0250 }
0251 if (wt > 0) {
0252 mean /= wt;
0253 rms /= wt;
0254 err = std::sqrt((rms - mean * mean) / wt);
0255 }
0256 return std::pair<double, double>(mean, err);
0257 }
0258
0259 std::pair<double, double> GetWidth(TH1D* hist, double xmin, double xmax) {
0260 double mean(0), mom2(0), rms(0), err(0), wt(0);
0261 for (int i = 1; i <= hist->GetNbinsX(); ++i) {
0262 double xlow = hist->GetBinLowEdge(i);
0263 double xhigh = xlow + hist->GetBinWidth(i);
0264 if ((xlow >= xmin) && (xhigh <= xmax)) {
0265 double cont = hist->GetBinContent(i);
0266 double valu = 0.5 * (xlow + xhigh);
0267 wt += cont;
0268 mean += (valu * cont);
0269 mom2 += (valu * valu * cont);
0270 }
0271 }
0272 if (wt > 0) {
0273 mean /= wt;
0274 mom2 /= wt;
0275 rms = std::sqrt(mom2 - mean * mean);
0276 err = rms / std::sqrt(2 * wt);
0277 }
0278 return std::pair<double, double>(rms, err);
0279 }
0280
0281 Double_t langaufun(Double_t* x, Double_t* par) {
0282
0283
0284
0285
0286
0287
0288
0289
0290
0291
0292
0293 Double_t invsq2pi = 0.3989422804014;
0294 Double_t mpshift = -0.22278298;
0295
0296
0297 Double_t np = 100.0;
0298 Double_t sc = 5.0;
0299
0300
0301 Double_t xx;
0302 Double_t mpc;
0303 Double_t fland;
0304 Double_t sum = 0.0;
0305 Double_t xlow, xupp;
0306 Double_t step;
0307 Double_t i;
0308 Double_t par0 = 0.2;
0309
0310
0311 mpc = par[0] - mpshift * par0 * par[0];
0312
0313
0314 xlow = x[0] - sc * par[2];
0315 xupp = x[0] + sc * par[2];
0316
0317 step = (xupp - xlow) / np;
0318
0319
0320 for (i = 1.0; i <= np / 2; i++) {
0321 xx = xlow + (i - .5) * step;
0322 fland = TMath::Landau(xx, mpc, par0 * par[0], kTRUE);
0323 sum += fland * TMath::Gaus(x[0], xx, par[2]);
0324
0325 xx = xupp - (i - .5) * step;
0326 fland = TMath::Landau(xx, mpc, par0 * par[0], kTRUE);
0327 sum += fland * TMath::Gaus(x[0], xx, par[2]);
0328 }
0329
0330 return (par[1] * step * sum * invsq2pi / par[2]);
0331 }
0332
0333 Double_t doubleGauss(Double_t* x, Double_t* par) {
0334 double x1 = x[0] - par[1];
0335 double sig1 = par[2];
0336 double x2 = x[0] - par[4];
0337 double sig2 = par[5];
0338 double yval =
0339 (par[0] * std::exp(-0.5 * (x1 / sig1) * (x1 / sig1)) + par[3] * std::exp(-0.5 * (x2 / sig2) * (x2 / sig2)));
0340 return yval;
0341 }
0342
0343 TFitResultPtr functionFit(TH1D* hist, double* fitrange, double* startvalues, double* parlimitslo, double* parlimitshi) {
0344 char FunName[100];
0345 sprintf(FunName, "Fitfcn_%s", hist->GetName());
0346 TF1* ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
0347 if (ffitold)
0348 delete ffitold;
0349
0350 int npar(6);
0351 TF1* ffit = new TF1(FunName, doubleGauss, fitrange[0], fitrange[1], npar);
0352 ffit->SetParameters(startvalues);
0353 ffit->SetLineColor(kBlue);
0354 ffit->SetParNames("Area1", "Mean1", "Width1", "Area2", "Mean2", "Width2");
0355 for (int i = 0; i < npar; i++)
0356 ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
0357 TFitResultPtr Fit = hist->Fit(FunName, "QRWLS");
0358 return Fit;
0359 }
0360
0361 std::pair<double, double> fitLanGau(TH1D* hist, bool debug) {
0362 double rms;
0363 std::pair<double, double> mrms = GetMean(hist, 0.005, 0.25, rms);
0364 double mean = mrms.first;
0365 double LowEdge = 0.005, HighEdge = mean + 3 * rms;
0366 TFitResultPtr Fit1 = hist->Fit("gaus", "+0wwqRS", "", LowEdge, HighEdge);
0367 if (debug)
0368 std::cout << hist->GetName() << " 0 " << Fit1->Value(0) << " 1 " << Fit1->Value(1) << " 2 " << Fit1->Value(2)
0369 << std::endl;
0370 double startvalues[3];
0371 startvalues[0] = Fit1->Value(1);
0372 startvalues[1] = Fit1->Value(0);
0373 startvalues[2] = Fit1->Value(2);
0374
0375 char FunName[100];
0376 sprintf(FunName, "Fitfcn_%s", hist->GetName());
0377 TF1* ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
0378 if (ffitold)
0379 delete ffitold;
0380
0381 TF1* ffit = new TF1(FunName, langaufun, LowEdge, HighEdge, 3);
0382 ffit->SetParameters(startvalues);
0383 ffit->SetParNames("MP", "Area", "GSigma");
0384 TFitResultPtr Fit2 = hist->Fit(FunName, "QRWLS");
0385 double value = Fit2->Value(0);
0386 double error = Fit2->FitResult::Error(0);
0387 if (debug)
0388 std::cout << hist->GetName() << " 0 " << Fit2->Value(0) << " 1 " << Fit2->Value(1) << " 2 " << Fit2->Value(2)
0389 << std::endl;
0390 return std::pair<double, double>(value, error);
0391 }
0392
0393 results fitTwoGauss(TH1D* hist, bool debug) {
0394 double rms;
0395 std::pair<double, double> mrms = GetMean(hist, 0.2, 2.0, rms);
0396 double mean = mrms.first;
0397 double LowEdge = mean - fitrangeFactor * rms;
0398 double HighEdge = mean + fitrangeFactor * rms;
0399 if (LowEdge < 0.15)
0400 LowEdge = 0.15;
0401 std::string option = (hist->GetEntries() > 100) ? "QRS" : "QRWLS";
0402 TF1* g1 = new TF1("g1", "gaus", LowEdge, HighEdge);
0403 g1->SetLineColor(kGreen);
0404 TFitResultPtr Fit = hist->Fit(g1, option.c_str(), "");
0405
0406 if (debug) {
0407 for (int k = 0; k < 3; ++k)
0408 std::cout << "Initial Parameter[" << k << "] = " << Fit->Value(k) << " +- " << Fit->FitResult::Error(k)
0409 << std::endl;
0410 }
0411 double startvalues[6], fitrange[2], lowValue[6], highValue[6];
0412 startvalues[0] = Fit->Value(0);
0413 lowValue[0] = 0.5 * startvalues[0];
0414 highValue[0] = 2. * startvalues[0];
0415 startvalues[1] = Fit->Value(1);
0416 lowValue[1] = 0.5 * startvalues[1];
0417 highValue[1] = 2. * startvalues[1];
0418 startvalues[2] = Fit->Value(2);
0419 lowValue[2] = 0.5 * startvalues[2];
0420 highValue[2] = 2. * startvalues[2];
0421 startvalues[3] = 0.1 * Fit->Value(0);
0422 lowValue[3] = 0.0;
0423 highValue[3] = 10. * startvalues[3];
0424 startvalues[4] = Fit->Value(1);
0425 lowValue[4] = 0.5 * startvalues[4];
0426 highValue[4] = 2. * startvalues[4];
0427 startvalues[5] = 2.0 * Fit->Value(2);
0428 lowValue[5] = 0.5 * startvalues[5];
0429 highValue[5] = 100. * startvalues[5];
0430
0431 fitrange[0] = Fit->Value(1) - fitrangeFactor * Fit->Value(2);
0432 fitrange[1] = Fit->Value(1) + fitrangeFactor * Fit->Value(2);
0433 TFitResultPtr Fitfun = functionFit(hist, fitrange, startvalues, lowValue, highValue);
0434 double wt1 = (Fitfun->Value(0)) * (Fitfun->Value(2));
0435 double value1 = Fitfun->Value(1);
0436 double error1 = Fitfun->FitResult::Error(1);
0437 double wval1 = Fitfun->Value(2);
0438 double werr1 = Fitfun->FitResult::Error(2);
0439 double wt2 = (Fitfun->Value(3)) * (Fitfun->Value(5));
0440 double value2 = Fitfun->Value(4);
0441 double error2 = Fitfun->FitResult::Error(4);
0442 double wval2 = Fitfun->Value(5);
0443 double werr2 = Fitfun->FitResult::Error(5);
0444 double value = (wt1 * value1 + wt2 * value2) / (wt1 + wt2);
0445 double wval = (wt1 * wval1 + wt2 * wval2) / (wt1 + wt2);
0446 double error = (sqrt((wt1 * error1) * (wt1 * error1) + (wt2 * error2) * (wt2 * error2)) / (wt1 + wt2));
0447 double werror = (sqrt((wt1 * werr1) * (wt1 * werr1) + (wt2 * werr2) * (wt2 * werr2)) / (wt1 + wt2));
0448 std::cout << hist->GetName() << " Fit " << value << "+-" << error << " width " << wval << " +- " << werror
0449 << " First " << value1 << "+-" << error1 << ":" << wval1 << "+-" << werr1 << ":" << wt1 << " Second "
0450 << value2 << "+-" << error2 << ":" << wval2 << "+-" << werr2 << ":" << wt2 << std::endl;
0451 if (debug) {
0452 for (int k = 0; k < 6; ++k)
0453 std::cout << hist->GetName() << ":Parameter[" << k << "] = " << Fitfun->Value(k) << " +- "
0454 << Fitfun->FitResult::Error(k) << std::endl;
0455 }
0456 return results(value, error, wval, werror);
0457 }
0458
0459 results fitOneGauss(TH1D* hist, bool fitTwice, bool debug) {
0460 double rms;
0461 std::pair<double, double> mrms = GetMean(hist, 0.2, 2.0, rms);
0462 double mean = mrms.first;
0463 double LowEdge = ((mean - fitrangeFactor * rms) < 0.5) ? 0.5 : (mean - fitrangeFactor * rms);
0464 double HighEdge = (mean + fitrangeFactor * rms);
0465 if (debug)
0466 std::cout << hist->GetName() << " Mean " << mean << " RMS " << rms << " Range " << LowEdge << ":" << HighEdge
0467 << "\n";
0468 std::string option = (hist->GetEntries() > 100) ? "QRS" : "QRWLS";
0469 TF1* g1 = new TF1("g1", "gaus", LowEdge, HighEdge);
0470 g1->SetLineColor(kGreen);
0471 TFitResultPtr Fit1 = hist->Fit(g1, option.c_str(), "");
0472 double value = Fit1->Value(1);
0473 double error = Fit1->FitResult::Error(1);
0474 double width = Fit1->Value(2);
0475 double werror = Fit1->FitResult::Error(2);
0476 if (fitTwice) {
0477 LowEdge = Fit1->Value(1) - fitrangeFactor * Fit1->Value(2);
0478 HighEdge = Fit1->Value(1) + fitrangeFactor * Fit1->Value(2);
0479 if (LowEdge < 0.5)
0480 LowEdge = 0.5;
0481 if (HighEdge > 5.0)
0482 HighEdge = 5.0;
0483 if (debug)
0484 std::cout << " Range for second Fit " << LowEdge << ":" << HighEdge << std::endl;
0485 TFitResultPtr Fit = hist->Fit("gaus", option.c_str(), "", LowEdge, HighEdge);
0486 value = Fit->Value(1);
0487 error = Fit->FitResult::Error(1);
0488 width = Fit->Value(2);
0489 werror = Fit->FitResult::Error(2);
0490 }
0491
0492 std::pair<double, double> meaner = GetMean(hist, 0.2, 2.0, rms);
0493 std::pair<double, double> rmserr = GetWidth(hist, 0.2, 2.0);
0494 if (debug) {
0495 std::cout << "Fit " << value << ":" << error << ":" << hist->GetMeanError() << " Mean " << meaner.first << ":"
0496 << meaner.second << " Width " << rmserr.first << ":" << rmserr.second;
0497 }
0498 double minvalue(0.30);
0499 if (value < minvalue || value > 2.0 || error > 0.5) {
0500 value = meaner.first;
0501 error = meaner.second;
0502 width = rmserr.first;
0503 werror = rmserr.second;
0504 }
0505 if (debug) {
0506 std::cout << " Final " << value << "+-" << error << ":" << width << "+-" << werror << std::endl;
0507 }
0508 return results(value, error, width, werror);
0509 }
0510
0511 void readCorrFactors(char* infile,
0512 double scale,
0513 std::map<int, cfactors>& cfacs,
0514 int& etamin,
0515 int& etamax,
0516 int& maxdepth,
0517 int iformat = 0,
0518 bool debug = false) {
0519 cfacs.clear();
0520 std::ifstream fInput(infile);
0521 if (!fInput.good()) {
0522 std::cout << "Cannot open file " << infile << std::endl;
0523 } else {
0524 char buffer[1024];
0525 unsigned int all(0), good(0);
0526 while (fInput.getline(buffer, 1024)) {
0527 ++all;
0528 if (buffer[0] == '#')
0529 continue;
0530 std::vector<std::string> items = splitString(std::string(buffer));
0531 if (items.size() != 5) {
0532 std::cout << "Ignore line: " << buffer << std::endl;
0533 } else {
0534 ++good;
0535 int ieta = (iformat == 1) ? std::atoi(items[0].c_str()) : std::atoi(items[1].c_str());
0536 int depth = (iformat == 1) ? std::atoi(items[1].c_str()) : std::atoi(items[2].c_str());
0537 float corrf = std::atof(items[3].c_str());
0538 float dcorr = (iformat == 1) ? (0.02 * corrf) : std::atof(items[4].c_str());
0539 cfactors cfac(ieta, depth, scale * corrf, scale * dcorr);
0540 int detId = (iformat == 1) ? repackId(items[2], ieta, depth) : repackId(ieta, depth);
0541 cfacs[detId] = cfactors(ieta, depth, corrf, dcorr);
0542 if (ieta > etamax)
0543 etamax = ieta;
0544 if (ieta < etamin)
0545 etamin = ieta;
0546 if (depth > maxdepth)
0547 maxdepth = depth;
0548 }
0549 }
0550 fInput.close();
0551 std::cout << "Reads total of " << all << " and " << good << " good records"
0552 << " from " << infile << std::endl;
0553 }
0554 if (debug) {
0555 unsigned k(0);
0556 std::cout << "Eta Range " << etamin << ":" << etamax << " Max Depth " << maxdepth << std::endl;
0557 for (std::map<int, cfactors>::const_iterator itr = cfacs.begin(); itr != cfacs.end(); ++itr, ++k)
0558 std::cout << "[" << k << "] " << std::hex << itr->first << std::dec << ": " << (itr->second).ieta << " "
0559 << (itr->second).depth << " " << (itr->second).corrf << " " << (itr->second).dcorr << std::endl;
0560 }
0561 }
0562
0563 void FitHistStandard(std::string infile,
0564 std::string outfile,
0565 std::string prefix,
0566 int mode = 11111,
0567 int type = 0,
0568 bool append = true,
0569 bool saveAll = false,
0570 bool debug = false) {
0571 int iname[6] = {0, 1, 2, 3, 4, 5};
0572 int checkmode[6] = {1000, 10, 10000, 1, 100000, 100};
0573 double xbin0[9] = {-21.0, -16.0, -12.0, -6.0, 0.0, 6.0, 12.0, 16.0, 21.0};
0574 double xbins[11] = {-25.0, -20.0, -15.0, -10.0, -5.0, 0.0, 5.0, 10.0, 15.0, 20.0, 25.0};
0575 double vbins[6] = {0.0, 7.0, 10.0, 13.0, 16.0, 50.0};
0576 double dlbins[9] = {0.0, 0.10, 0.20, 0.50, 1.0, 2.0, 2.5, 3.0, 10.0};
0577 std::string sname[4] = {"ratio", "etaR", "dl1R", "nvxR"};
0578 std::string lname[4] = {"Z", "E", "L", "V"};
0579 std::string wname[4] = {"W", "F", "M", "X"};
0580 std::string xname[4] = {"i#eta", "i#eta", "d_{L1}", "# Vertex"};
0581 int numb[4] = {10, 8, 8, 5};
0582
0583 if (type == 0) {
0584 numb[0] = 8;
0585 for (int i = 0; i < 9; ++i)
0586 xbins[i] = xbin0[i];
0587 }
0588 TFile* file = new TFile(infile.c_str());
0589 std::vector<TH1D*> hists;
0590 char name[100], namw[100];
0591 if (file != nullptr) {
0592 for (int m1 = 0; m1 < 4; ++m1) {
0593 for (int m2 = 0; m2 < 6; ++m2) {
0594 sprintf(name, "%s%s%d0", prefix.c_str(), sname[m1].c_str(), iname[m2]);
0595 TH1D* hist0 = (TH1D*)file->FindObjectAny(name);
0596 bool ok = ((hist0 != nullptr) && (hist0->GetEntries() > 25));
0597 if ((mode / checkmode[m2]) % 10 > 0 && ok) {
0598 TH1D *histo(0), *histw(0);
0599 sprintf(name, "%s%s%d", prefix.c_str(), lname[m1].c_str(), iname[m2]);
0600 sprintf(namw, "%s%s%d", prefix.c_str(), wname[m1].c_str(), iname[m2]);
0601 if (m1 <= 1) {
0602 histo = new TH1D(name, hist0->GetTitle(), numb[m1], xbins);
0603 histw = new TH1D(namw, hist0->GetTitle(), numb[m1], xbins);
0604 } else if (m1 == 2) {
0605 histo = new TH1D(name, hist0->GetTitle(), numb[m1], dlbins);
0606 histw = new TH1D(namw, hist0->GetTitle(), numb[m1], dlbins);
0607 } else {
0608 histo = new TH1D(name, hist0->GetTitle(), numb[m1], vbins);
0609 histw = new TH1D(namw, hist0->GetTitle(), numb[m1], vbins);
0610 }
0611 int jmin(numb[m1]), jmax(0);
0612 for (int j = 0; j <= numb[m1]; ++j) {
0613 sprintf(name, "%s%s%d%d", prefix.c_str(), sname[m1].c_str(), iname[m2], j);
0614 TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
0615 TH1D* hist = (TH1D*)hist1->Clone();
0616 double value(0), error(0), width(0), werror(0);
0617 if (hist->GetEntries() > 0) {
0618 value = hist->GetMean();
0619 error = hist->GetRMS();
0620 std::pair<double, double> rmserr = GetWidth(hist, 0.2, 2.0);
0621 width = rmserr.first;
0622 werror = rmserr.second;
0623 }
0624 if (hist->GetEntries() > 4) {
0625 bool flag = (j == 0) ? true : false;
0626 results meaner = fitOneGauss(hist, flag, debug);
0627 value = meaner.mean;
0628 error = meaner.errmean;
0629 width = meaner.width;
0630 werror = meaner.errwidth;
0631 if (j != 0) {
0632 if (j < jmin)
0633 jmin = j;
0634 if (j > jmax)
0635 jmax = j;
0636 }
0637 }
0638 if (j == 0) {
0639 hists.push_back(hist);
0640 } else {
0641 double wbyv = width / value;
0642 double wverr = wbyv * std::sqrt((werror * werror) / (width * width) + (error * error) / (value * value));
0643 if (saveAll)
0644 hists.push_back(hist);
0645 histo->SetBinContent(j, value);
0646 histo->SetBinError(j, error);
0647 histw->SetBinContent(j, wbyv);
0648 histw->SetBinError(j, wverr);
0649 }
0650 }
0651 if (histo->GetEntries() > 2) {
0652 double LowEdge = histo->GetBinLowEdge(jmin);
0653 double HighEdge = histo->GetBinLowEdge(jmax) + histo->GetBinWidth(jmax);
0654 TFitResultPtr Fit = histo->Fit("pol0", "+QRWLS", "", LowEdge, HighEdge);
0655 if (debug) {
0656 std::cout << "Fit to Pol0: " << Fit->Value(0) << " +- " << Fit->FitResult::Error(0) << std::endl;
0657 }
0658 histo->GetXaxis()->SetTitle(xname[m1].c_str());
0659 histo->GetYaxis()->SetTitle("MPV(E_{HCAL}/(p-E_{ECAL}))");
0660 histo->GetYaxis()->SetRangeUser(0.4, 1.6);
0661 histw->GetXaxis()->SetTitle(xname[m1].c_str());
0662 histw->GetYaxis()->SetTitle("Width/MPV of (E_{HCAL}/(p-E_{ECAL}))");
0663 histw->GetYaxis()->SetRangeUser(0.0, 0.5);
0664 }
0665 hists.push_back(histo);
0666 hists.push_back(histw);
0667 }
0668 }
0669 }
0670 TFile* theFile(0);
0671 if (append) {
0672 theFile = new TFile(outfile.c_str(), "UPDATE");
0673 } else {
0674 theFile = new TFile(outfile.c_str(), "RECREATE");
0675 }
0676 theFile->cd();
0677 for (unsigned int i = 0; i < hists.size(); ++i) {
0678 TH1D* hnew = (TH1D*)hists[i]->Clone();
0679 hnew->Write();
0680 }
0681 theFile->Close();
0682 file->Close();
0683 }
0684 }
0685
0686 void FitHistExtended(const char* infile,
0687 const char* outfile,
0688 std::string prefix,
0689 int numb = 50,
0690 int type = 3,
0691 bool append = true,
0692 bool fiteta = true,
0693 int iname = 3,
0694 bool debug = false) {
0695 std::string sname("ratio"), lname("Z"), wname("W"), ename("etaB");
0696 double xbins[99];
0697 double xbin[23] = {-23.0, -21.0, -19.0, -17.0, -15.0, -13.0, -11.0, -9.0, -7.0, -5.0, -3.0, 0.0,
0698 3.0, 5.0, 7.0, 9.0, 11.0, 13.0, 15.0, 17.0, 19.0, 21.0, 23.0};
0699 if (type == 2) {
0700 numb = 22;
0701 for (int k = 0; k <= numb; ++k)
0702 xbins[k] = xbin[k];
0703 } else if (type == 1) {
0704 numb = 1;
0705 xbins[0] = -25;
0706 xbins[1] = 25;
0707 } else {
0708 int neta = numb / 2;
0709 for (int k = 0; k < neta; ++k) {
0710 xbins[k] = (k - neta) - 0.5;
0711 xbins[numb - k] = (neta - k) + 0.5;
0712 }
0713 xbins[neta] = 0;
0714 }
0715 if (debug) {
0716 for (int k = 0; k <= numb; ++k)
0717 std::cout << " " << xbins[k];
0718 std::cout << std::endl;
0719 }
0720 TFile* file = new TFile(infile);
0721 std::vector<TH1D*> hists;
0722 char name[200];
0723 if (debug) {
0724 std::cout << infile << " " << file << std::endl;
0725 }
0726 if (file != nullptr) {
0727 sprintf(name, "%s%s%d0", prefix.c_str(), sname.c_str(), iname);
0728 TH1D* hist0 = (TH1D*)file->FindObjectAny(name);
0729 bool ok = (hist0 != nullptr);
0730 if (debug) {
0731 std::cout << name << " Pointer " << hist0 << " " << ok << std::endl;
0732 }
0733 if (ok) {
0734 TH1D *histo(0), *histw(0);
0735 if (numb > 0) {
0736 sprintf(name, "%s%s%d", prefix.c_str(), lname.c_str(), iname);
0737 histo = new TH1D(name, hist0->GetTitle(), numb, xbins);
0738 sprintf(name, "%s%s%d", prefix.c_str(), wname.c_str(), iname);
0739 histw = new TH1D(name, hist0->GetTitle(), numb, xbins);
0740 if (debug)
0741 std::cout << name << " " << histo->GetNbinsX() << std::endl;
0742 }
0743 if (hist0->GetEntries() > 10) {
0744 double rms;
0745 results meaner0 = fitOneGauss(hist0, true, debug);
0746 std::pair<double, double> meaner1 = GetMean(hist0, 0.2, 2.0, rms);
0747 std::pair<double, double> meaner2 = GetWidth(hist0, 0.2, 2.0);
0748 if (debug) {
0749 std::cout << "Fit " << meaner0.mean << ":" << meaner0.errmean << " Mean1 " << hist0->GetMean() << ":"
0750 << hist0->GetMeanError() << " Mean2 " << meaner1.first << ":" << meaner1.second << " Width "
0751 << meaner2.first << ":" << meaner2.second << std::endl;
0752 }
0753 }
0754 int nv1(100), nv2(0);
0755 int jmin(numb), jmax(0);
0756 for (int j = 0; j <= numb; ++j) {
0757 sprintf(name, "%s%s%d%d", prefix.c_str(), sname.c_str(), iname, j);
0758 TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
0759 if (debug) {
0760 std::cout << "Get Histogram for " << name << " at " << hist1 << std::endl;
0761 }
0762 double value(0), error(0), total(0), width(0), werror(0);
0763 if (hist1 == nullptr) {
0764 value = 1.0;
0765 } else {
0766 TH1D* hist = (TH1D*)hist1->Clone();
0767 if (debug)
0768 std::cout << "Histogram " << name << ":" << (hist->GetName()) << " with " << (hist->GetEntries())
0769 << " entries" << std::endl;
0770 if (hist->GetEntries() > 0) {
0771 value = hist->GetMean();
0772 error = hist->GetRMS();
0773 for (int i = 1; i <= hist->GetNbinsX(); ++i)
0774 total += hist->GetBinContent(i);
0775 std::pair<double, double> rmserr = GetWidth(hist, 0.2, 2.0);
0776 width = rmserr.first;
0777 werror = rmserr.second;
0778 }
0779 if (total > 4) {
0780 if (nv1 > j)
0781 nv1 = j;
0782 if (nv2 < j)
0783 nv2 = j;
0784 if (j == 0) {
0785 sprintf(name, "%sOne", hist1->GetName());
0786 TH1D* hist2 = (TH1D*)hist1->Clone(name);
0787 fitOneGauss(hist2, true, debug);
0788 hists.push_back(hist2);
0789 results meaner = fitOneGauss(hist, true, debug);
0790 value = meaner.mean;
0791 error = meaner.errmean;
0792 width = meaner.width;
0793 werror = meaner.errwidth;
0794 } else {
0795 results meaner = fitOneGauss(hist, true, debug);
0796 value = meaner.mean;
0797 error = meaner.errmean;
0798 width = meaner.width;
0799 werror = meaner.errwidth;
0800 }
0801 if (j != 0) {
0802 if (j < jmin)
0803 jmin = j;
0804 if (j > jmax)
0805 jmax = j;
0806 }
0807 }
0808 hists.push_back(hist);
0809 }
0810 if (debug) {
0811 std::cout << "Hist " << j << " Value " << value << " +- " << error << std::endl;
0812 }
0813 if (j != 0) {
0814 double wbyv = width / value;
0815 double wverr = wbyv * std::sqrt((werror * werror) / (width * width) + (error * error) / (value * value));
0816 histo->SetBinContent(j, value);
0817 histo->SetBinError(j, error);
0818 histw->SetBinContent(j, wbyv);
0819 histw->SetBinError(j, wverr);
0820 }
0821 }
0822 if (histo != nullptr) {
0823 if (histo->GetEntries() > 2 && fiteta) {
0824 if (debug) {
0825 std::cout << "Jmin/max " << jmin << ":" << jmax << ":" << histo->GetNbinsX() << std::endl;
0826 }
0827 double LowEdge = histo->GetBinLowEdge(jmin);
0828 double HighEdge = histo->GetBinLowEdge(jmax) + histo->GetBinWidth(jmax);
0829 TFitResultPtr Fit = histo->Fit("pol0", "+QRWLS", "", LowEdge, HighEdge);
0830 if (debug) {
0831 std::cout << "Fit to Pol0: " << Fit->Value(0) << " +- " << Fit->FitResult::Error(0) << " in range " << nv1
0832 << ":" << xbins[nv1] << ":" << nv2 << ":" << xbins[nv2] << std::endl;
0833 }
0834 histo->GetXaxis()->SetTitle("i#eta");
0835 histo->GetYaxis()->SetTitle("MPV(E_{HCAL}/(p-E_{ECAL}))");
0836 histo->GetYaxis()->SetRangeUser(0.4, 1.6);
0837 histw->GetXaxis()->SetTitle("i#eta");
0838 histw->GetYaxis()->SetTitle("MPV/Width(E_{HCAL}/(p-E_{ECAL}))");
0839 histw->GetYaxis()->SetRangeUser(0.0, 0.5);
0840 }
0841 hists.push_back(histo);
0842 hists.push_back(histw);
0843 } else {
0844 hists.push_back(hist0);
0845 }
0846
0847
0848 for (int j = 1; j <= 4; ++j) {
0849 sprintf(name, "%s%s%d%d", prefix.c_str(), ename.c_str(), iname, j);
0850 TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
0851 if (debug) {
0852 std::cout << "Get Histogram for " << name << " at " << hist1 << std::endl;
0853 }
0854 if (hist1 != nullptr) {
0855 TH1D* hist = (TH1D*)hist1->Clone();
0856 double value(0), error(0), total(0), width(0), werror(0);
0857 if (hist->GetEntries() > 0) {
0858 value = hist->GetMean();
0859 error = hist->GetRMS();
0860 for (int i = 1; i <= hist->GetNbinsX(); ++i)
0861 total += hist->GetBinContent(i);
0862 }
0863 if (total > 4) {
0864 sprintf(name, "%sOne", hist1->GetName());
0865 TH1D* hist2 = (TH1D*)hist1->Clone(name);
0866 results meanerr = fitOneGauss(hist2, true, debug);
0867 value = meanerr.mean;
0868 error = meanerr.errmean;
0869 width = meanerr.width;
0870 werror = meanerr.errwidth;
0871 double wbyv = width / value;
0872 double wverr = wbyv * std::sqrt((werror * werror) / (width * width) + (error * error) / (value * value));
0873 std::cout << hist2->GetName() << " MPV " << value << " +- " << error << " Width " << width << " +- "
0874 << werror << " W/M " << wbyv << " +- " << wverr << std::endl;
0875 hists.push_back(hist2);
0876 if (hist1->GetBinLowEdge(1) < 0.1) {
0877 sprintf(name, "%sTwo", hist1->GetName());
0878 TH1D* hist3 = (TH1D*)hist1->Clone(name);
0879 fitLanGau(hist3, debug);
0880 hists.push_back(hist3);
0881 }
0882
0883 results meaner0 = fitOneGauss(hist, true, debug);
0884 value = meaner0.mean;
0885 error = meaner0.errmean;
0886 double rms;
0887 std::pair<double, double> meaner = GetMean(hist, 0.2, 2.0, rms);
0888 if (debug) {
0889 std::cout << "Fit " << value << ":" << error << ":" << hist->GetMeanError() << " Mean " << meaner.first
0890 << ":" << meaner.second << std::endl;
0891 }
0892 }
0893 hists.push_back(hist);
0894 }
0895 }
0896 }
0897 TFile* theFile(0);
0898 if (append) {
0899 if (debug) {
0900 std::cout << "Open file " << outfile << " in append mode" << std::endl;
0901 }
0902 theFile = new TFile(outfile, "UPDATE");
0903 } else {
0904 if (debug) {
0905 std::cout << "Open file " << outfile << " in recreate mode" << std::endl;
0906 }
0907 theFile = new TFile(outfile, "RECREATE");
0908 }
0909
0910 theFile->cd();
0911 for (unsigned int i = 0; i < hists.size(); ++i) {
0912 TH1D* hnew = (TH1D*)hists[i]->Clone();
0913 if (debug) {
0914 std::cout << "Write Histogram " << hnew->GetTitle() << std::endl;
0915 }
0916 hnew->Write();
0917 }
0918 theFile->Close();
0919 file->Close();
0920 }
0921 }
0922
0923 void FitHistRBX(const char* infile, const char* outfile, std::string prefix, bool append = true, int iname = 3) {
0924 std::string sname("RBX"), lname("R");
0925 int numb(18);
0926 bool debug(false);
0927 char name[200];
0928
0929 TFile* file = new TFile(infile);
0930 std::vector<TH1D*> hists;
0931 sprintf(name, "%s%s%d", prefix.c_str(), lname.c_str(), iname);
0932 TH1D* histo = new TH1D(name, "", numb, 0, numb);
0933 if (debug) {
0934 std::cout << infile << " " << file << std::endl;
0935 }
0936 if (file != nullptr) {
0937 for (int j = 0; j < numb; ++j) {
0938 sprintf(name, "%s%s%d%d", prefix.c_str(), sname.c_str(), iname, j + 1);
0939 TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
0940 if (debug) {
0941 std::cout << "Get Histogram for " << name << " at " << hist1 << std::endl;
0942 }
0943 TH1D* hist = (TH1D*)hist1->Clone();
0944 double value(0), error(0), total(0);
0945 if (hist->GetEntries() > 0) {
0946 value = hist->GetMean();
0947 error = hist->GetRMS();
0948 for (int i = 1; i <= hist->GetNbinsX(); ++i)
0949 total += hist->GetBinContent(i);
0950 }
0951 if (total > 4) {
0952 results meaner = fitOneGauss(hist, false, debug);
0953 value = meaner.mean;
0954 error = meaner.errmean;
0955 }
0956 hists.push_back(hist);
0957 histo->SetBinContent(j + 1, value);
0958 histo->SetBinError(j + 1, error);
0959 }
0960 histo->GetXaxis()->SetTitle("RBX #");
0961 histo->GetYaxis()->SetTitle("MPV(E_{HCAL}/(p-E_{ECAL}))");
0962 histo->GetYaxis()->SetRangeUser(0.75, 1.20);
0963 hists.push_back(histo);
0964
0965 TFile* theFile(0);
0966 if (append) {
0967 if (debug) {
0968 std::cout << "Open file " << outfile << " in append mode" << std::endl;
0969 }
0970 theFile = new TFile(outfile, "UPDATE");
0971 } else {
0972 if (debug) {
0973 std::cout << "Open file " << outfile << " in recreate mode" << std::endl;
0974 }
0975 theFile = new TFile(outfile, "RECREATE");
0976 }
0977
0978 theFile->cd();
0979 for (unsigned int i = 0; i < hists.size(); ++i) {
0980 TH1D* hnew = (TH1D*)hists[i]->Clone();
0981 if (debug) {
0982 std::cout << "Write Histogram " << hnew->GetTitle() << std::endl;
0983 }
0984 hnew->Write();
0985 }
0986 theFile->Close();
0987 file->Close();
0988 }
0989 }
0990
0991 void PlotHist(const char* infile,
0992 std::string prefix,
0993 std::string text,
0994 int mode = 4,
0995 int kopt = 100,
0996 double lumi = 0,
0997 double ener = 13.0,
0998 bool dataMC = false,
0999 bool drawStatBox = true,
1000 int save = 0) {
1001 std::string name0[6] = {"ratio00", "ratio10", "ratio20", "ratio30", "ratio40", "ratio50"};
1002 std::string name1[5] = {"Z0", "Z1", "Z2", "Z3", "Z4"};
1003 std::string name2[5] = {"L0", "L1", "L2", "L3", "L4"};
1004 std::string name3[5] = {"V0", "V1", "V2", "V3", "V4"};
1005 std::string name4[20] = {"etaB41", "etaB42", "etaB43", "etaB44", "etaB31", "etaB32", "etaB33",
1006 "etaB34", "etaB21", "etaB22", "etaB23", "etaB24", "etaB11", "etaB12",
1007 "etaB13", "etaB14", "etaB01", "etaB02", "etaB03", "etaB04"};
1008 std::string name5[5] = {"W0", "W1", "W2", "W3", "W4"};
1009 std::string title[6] = {"Tracks with p = 10:20 GeV",
1010 "Tracks with p = 20:30 GeV",
1011 "Tracks with p = 30:40 GeV",
1012 "Tracks with p = 40:60 GeV",
1013 "Tracks with p = 60:100 GeV",
1014 "Tracks with p = 20:100 GeV"};
1015 std::string title1[20] = {"Tracks with p = 60:100 GeV (Barrel)", "Tracks with p = 60:100 GeV (Transition)",
1016 "Tracks with p = 60:100 GeV (Endcap)", "Tracks with p = 60:100 GeV",
1017 "Tracks with p = 40:60 GeV (Barrel)", "Tracks with p = 40:60 GeV (Transition)",
1018 "Tracks with p = 40:60 GeV (Endcap)", "Tracks with p = 40:60 GeV",
1019 "Tracks with p = 30:40 GeV (Barrel)", "Tracks with p = 30:40 GeV (Transition)",
1020 "Tracks with p = 30:40 GeV (Endcap)", "Tracks with p = 30:40 GeV",
1021 "Tracks with p = 20:30 GeV (Barrel)", "Tracks with p = 20:30 GeV (Transition)",
1022 "Tracks with p = 20:30 GeV (Endcap)", "Tracks with p = 20:30 GeV",
1023 "Tracks with p = 10:20 GeV (Barrel)", "Tracks with p = 10:20 GeV (Transition)",
1024 "Tracks with p = 10:20 GeV (Endcap)", "Tracks with p = 10:20 GeV"};
1025 std::string xtitl[5] = {"E_{HCAL}/(p-E_{ECAL})", "i#eta", "d_{L1}", "# Vertex", "E_{HCAL}/(p-E_{ECAL})"};
1026 std::string ytitl[5] = {
1027 "Tracks", "MPV(E_{HCAL}/(p-E_{ECAL}))", "MPV(E_{HCAL}/(p-E_{ECAL}))", "MPV(E_{HCAL}/(p-E_{ECAL}))", "Tracks"};
1028
1029 gStyle->SetCanvasBorderMode(0);
1030 gStyle->SetCanvasColor(kWhite);
1031 gStyle->SetPadColor(kWhite);
1032 gStyle->SetFillColor(kWhite);
1033 gStyle->SetOptTitle(0);
1034 if (mode < 0 || mode > 5)
1035 mode = 0;
1036 if (drawStatBox) {
1037 int iopt(1110);
1038 if (mode != 0)
1039 iopt = 10;
1040 gStyle->SetOptStat(iopt);
1041 gStyle->SetOptFit(1);
1042 } else {
1043 gStyle->SetOptStat(0);
1044 gStyle->SetOptFit(0);
1045 }
1046 TFile* file = new TFile(infile);
1047 TLine* line(0);
1048 char name[100], namep[100];
1049 int kmax = (mode == 4) ? 20 : (((mode < 1) && (mode > 5)) ? 6 : 5);
1050 for (int k = 0; k < kmax; ++k) {
1051 if (mode == 1) {
1052 sprintf(name, "%s%s", prefix.c_str(), name1[k].c_str());
1053 } else if (mode == 2) {
1054 sprintf(name, "%s%s", prefix.c_str(), name2[k].c_str());
1055 } else if (mode == 3) {
1056 sprintf(name, "%s%s", prefix.c_str(), name3[k].c_str());
1057 } else if (mode == 4) {
1058 if ((kopt / 100) % 10 == 0) {
1059 sprintf(name, "%s%s", prefix.c_str(), name4[k].c_str());
1060 } else if ((kopt / 100) % 10 == 2) {
1061 sprintf(name, "%s%sTwo", prefix.c_str(), name4[k].c_str());
1062 } else {
1063 sprintf(name, "%s%sOne", prefix.c_str(), name4[k].c_str());
1064 }
1065 } else if (mode == 5) {
1066 sprintf(name, "%s%s", prefix.c_str(), name5[k].c_str());
1067 } else {
1068 if ((kopt / 100) % 10 == 0) {
1069 sprintf(name, "%s%s", prefix.c_str(), name0[k].c_str());
1070 } else if ((kopt / 100) % 10 == 2) {
1071 sprintf(name, "%s%sTwo", prefix.c_str(), name0[k].c_str());
1072 } else {
1073 sprintf(name, "%s%sOne", prefix.c_str(), name0[k].c_str());
1074 }
1075 }
1076 TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1077 if (hist1 != nullptr) {
1078 TH1D* hist = (TH1D*)(hist1->Clone());
1079 double p0(1);
1080 double ymin(0.90);
1081 sprintf(namep, "c_%s", name);
1082 TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1083 pad->SetRightMargin(0.10);
1084 pad->SetTopMargin(0.10);
1085 if ((kopt / 10) % 10 > 0)
1086 gPad->SetGrid();
1087 hist->GetXaxis()->SetTitleSize(0.04);
1088 hist->GetXaxis()->SetTitle(xtitl[mode].c_str());
1089 hist->GetYaxis()->SetTitle(ytitl[mode].c_str());
1090 hist->GetYaxis()->SetLabelOffset(0.005);
1091 hist->GetYaxis()->SetTitleSize(0.04);
1092 hist->GetYaxis()->SetLabelSize(0.035);
1093 hist->GetYaxis()->SetTitleOffset(1.10);
1094 if (mode == 0 || mode == 4) {
1095 if ((kopt / 100) % 10 == 2)
1096 hist->GetXaxis()->SetRangeUser(0.0, 0.30);
1097 else
1098 hist->GetXaxis()->SetRangeUser(0.25, 2.25);
1099 } else {
1100 if (mode == 5)
1101 hist->GetYaxis()->SetRangeUser(0.1, 0.50);
1102 else if (dataMC)
1103 hist->GetYaxis()->SetRangeUser(0.5, 1.50);
1104 else
1105 hist->GetYaxis()->SetRangeUser(0.8, 1.20);
1106 if (kopt % 10 > 0) {
1107 int nbin = hist->GetNbinsX();
1108 double LowEdge = (kopt % 10 == 1) ? hist->GetBinLowEdge(1) : -20;
1109 double HighEdge = (kopt % 10 == 1) ? hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin) : 20;
1110 TFitResultPtr Fit = hist->Fit("pol0", "+QRWLS", "", LowEdge, HighEdge);
1111 p0 = Fit->Value(0);
1112 }
1113 }
1114 hist->SetMarkerStyle(20);
1115 hist->SetMarkerColor(2);
1116 hist->SetLineColor(2);
1117 hist->Draw();
1118 pad->Update();
1119 TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
1120 if (st1 != nullptr) {
1121 ymin = (mode == 0 || mode == 4) ? 0.70 : 0.80;
1122 st1->SetY1NDC(ymin);
1123 st1->SetY2NDC(0.90);
1124 st1->SetX1NDC(0.65);
1125 st1->SetX2NDC(0.90);
1126 }
1127 if (mode != 0 && mode != 4 && kopt % 10 > 0) {
1128 double xmin = hist->GetBinLowEdge(1);
1129 int nbin = hist->GetNbinsX();
1130 double xmax = hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin);
1131 double mean(0), rms(0), total(0);
1132 int kount(0);
1133 for (int k = 2; k < nbin; ++k) {
1134 double x = hist->GetBinContent(k);
1135 double w = hist->GetBinError(k);
1136 mean += (x * w);
1137 rms += (x * x * w);
1138 total += w;
1139 ++kount;
1140 }
1141 mean /= total;
1142 rms /= total;
1143 double error = sqrt(rms - mean * mean) / total;
1144 line = new TLine(xmin, p0, xmax, p0);
1145 std::cout << xmin << ":" << xmax << ":" << p0 << " Mean " << nbin << ":" << kount << ":" << total << ":" << mean
1146 << ":" << rms << ":" << error << std::endl;
1147 line->SetLineWidth(2);
1148 line->SetLineStyle(2);
1149 line->Draw("same");
1150 }
1151 pad->Update();
1152 double ymx(0.96), xmi(0.25), xmx(0.90);
1153 char txt[100];
1154 if (lumi > 0.1) {
1155 ymx = ymin - 0.005;
1156 xmi = 0.45;
1157 TPaveText* txt0 = new TPaveText(0.65, 0.91, 0.90, 0.96, "blNDC");
1158 txt0->SetFillColor(0);
1159 sprintf(txt, "%4.1f TeV %5.1f fb^{-1}", ener, lumi);
1160 txt0->AddText(txt);
1161 txt0->Draw("same");
1162 }
1163 double ymi = ymx - 0.05;
1164 TPaveText* txt1 = new TPaveText(xmi, ymi, xmx, ymx, "blNDC");
1165 txt1->SetFillColor(0);
1166 if (text == "") {
1167 if (mode == 4)
1168 sprintf(txt, "%s", title1[k].c_str());
1169 else
1170 sprintf(txt, "%s", title[k].c_str());
1171 } else {
1172 if (mode == 4)
1173 sprintf(txt, "%s (%s)", title1[k].c_str(), text.c_str());
1174 else
1175 sprintf(txt, "%s (%s)", title[k].c_str(), text.c_str());
1176 }
1177 txt1->AddText(txt);
1178 txt1->Draw("same");
1179 double xmax = (dataMC) ? 0.33 : 0.44;
1180 ymi = (lumi > 0.1) ? 0.91 : 0.84;
1181 ymx = ymi + 0.05;
1182 TPaveText* txt2 = new TPaveText(0.11, ymi, xmax, ymx, "blNDC");
1183 txt2->SetFillColor(0);
1184 if (dataMC)
1185 sprintf(txt, "CMS Preliminary");
1186 else
1187 sprintf(txt, "CMS Simulation Preliminary");
1188 txt2->AddText(txt);
1189 txt2->Draw("same");
1190 pad->Modified();
1191 pad->Update();
1192 if (save > 0) {
1193 sprintf(name, "%s.pdf", pad->GetName());
1194 pad->Print(name);
1195 } else if (save < 0) {
1196 sprintf(name, "%s.C", pad->GetName());
1197 pad->Print(name);
1198 }
1199 }
1200 }
1201 }
1202
1203 void PlotHistEta(const char* infile,
1204 std::string prefix,
1205 std::string text,
1206 int iene = 3,
1207 int numb = 50,
1208 int ieta = 0,
1209 double lumi = 0,
1210 double ener = 13.0,
1211 bool dataMC = false,
1212 bool drawStatBox = true,
1213 int save = 0) {
1214 std::string name0 = "ratio";
1215 std::string title[5] = {"10:20", "20:30", "30:40", "40:60", "60:100"};
1216 std::string xtitl = "E_{HCAL}/(p-E_{ECAL})";
1217 std::string ytitl = "Tracks";
1218
1219 gStyle->SetCanvasBorderMode(0);
1220 gStyle->SetCanvasColor(kWhite);
1221 gStyle->SetPadColor(kWhite);
1222 gStyle->SetFillColor(kWhite);
1223 gStyle->SetOptTitle(0);
1224 if (drawStatBox) {
1225 int iopt(1110);
1226 gStyle->SetOptStat(iopt);
1227 gStyle->SetOptFit(1);
1228 } else {
1229 gStyle->SetOptStat(0);
1230 gStyle->SetOptFit(0);
1231 }
1232 if (iene < 0 || iene >= 5)
1233 iene = 3;
1234 int numb2 = numb / 2;
1235 if (ieta < -numb2 || ieta > numb2)
1236 ieta = 0;
1237 int ietaMin = ((ieta == 0) ? 1 : ((ieta > 0) ? (numb2 + ieta) : (numb2 + ieta + 1)));
1238 int ietaMax = (ieta == 0) ? numb : ietaMin;
1239 TFile* file = new TFile(infile);
1240 char name[100], namep[100];
1241 for (int k = ietaMin; k <= ietaMax; ++k) {
1242 int eta = (k > numb2) ? (k - numb2) : (k - numb2 - 1);
1243 sprintf(name, "%s%s%d%d", prefix.c_str(), name0.c_str(), iene, k);
1244 TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1245 std::cout << name << " at " << hist1 << std::endl;
1246 if (hist1 != nullptr) {
1247 TH1D* hist = (TH1D*)(hist1->Clone());
1248 double ymin(0.90);
1249 sprintf(namep, "c_%s", name);
1250 TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1251 pad->SetRightMargin(0.10);
1252 pad->SetTopMargin(0.10);
1253 hist->GetXaxis()->SetTitleSize(0.04);
1254 hist->GetXaxis()->SetTitle(xtitl.c_str());
1255 hist->GetYaxis()->SetTitle(ytitl.c_str());
1256 hist->GetYaxis()->SetLabelOffset(0.005);
1257 hist->GetYaxis()->SetTitleSize(0.04);
1258 hist->GetYaxis()->SetLabelSize(0.035);
1259 hist->GetYaxis()->SetTitleOffset(1.10);
1260 hist->GetXaxis()->SetRangeUser(0.25, 2.25);
1261 hist->SetMarkerStyle(20);
1262 hist->SetMarkerColor(2);
1263 hist->SetLineColor(2);
1264 hist->Draw();
1265 pad->Update();
1266 TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
1267 if (st1 != nullptr) {
1268 ymin = 0.70;
1269 st1->SetY1NDC(ymin);
1270 st1->SetY2NDC(0.90);
1271 st1->SetX1NDC(0.65);
1272 st1->SetX2NDC(0.90);
1273 }
1274 double ymx(0.96), xmi(0.25), xmx(0.90);
1275 char txt[100];
1276 if (lumi > 0.1) {
1277 ymx = ymin - 0.005;
1278 xmi = 0.45;
1279 TPaveText* txt0 = new TPaveText(0.65, 0.91, 0.90, 0.96, "blNDC");
1280 txt0->SetFillColor(0);
1281 sprintf(txt, "%4.1f TeV %5.1f fb^{-1}", ener, lumi);
1282 txt0->AddText(txt);
1283 txt0->Draw("same");
1284 }
1285 double ymi = ymx - 0.05;
1286 TPaveText* txt1 = new TPaveText(xmi, ymi, xmx, ymx, "blNDC");
1287 txt1->SetFillColor(0);
1288 if (text == "") {
1289 sprintf(txt, "Tracks with p = %s GeV at i#eta = %d", title[iene].c_str(), eta);
1290 } else {
1291 sprintf(txt, "Tracks with p = %s GeV at i#eta = %d (%s)", title[iene].c_str(), eta, text.c_str());
1292 }
1293 txt1->AddText(txt);
1294 txt1->Draw("same");
1295 double xmax = (dataMC) ? 0.33 : 0.44;
1296 ymi = (lumi > 0.1) ? 0.91 : 0.84;
1297 ymx = ymi + 0.05;
1298 TPaveText* txt2 = new TPaveText(0.11, ymi, xmax, ymx, "blNDC");
1299 txt2->SetFillColor(0);
1300 if (dataMC)
1301 sprintf(txt, "CMS Preliminary");
1302 else
1303 sprintf(txt, "CMS Simulation Preliminary");
1304 txt2->AddText(txt);
1305 txt2->Draw("same");
1306 pad->Modified();
1307 pad->Update();
1308 if (save > 0) {
1309 sprintf(name, "%s.pdf", pad->GetName());
1310 pad->Print(name);
1311 } else if (save < 0) {
1312 sprintf(name, "%s.C", pad->GetName());
1313 pad->Print(name);
1314 }
1315 }
1316 }
1317 }
1318
1319 void PlotHists(std::string infile, std::string prefix, std::string text, bool drawStatBox = true, int save = 0) {
1320 int colors[6] = {1, 6, 4, 7, 2, 9};
1321 std::string types[6] = {"B", "C", "D", "E", "F", "G"};
1322 std::string names[3] = {"ratio20", "Z2", "W2"};
1323 std::string xtitl[3] = {"E_{HCAL}/(p-E_{ECAL})", "i#eta", "i#eta"};
1324 std::string ytitl[3] = {"Tracks", "MPV(E_{HCAL}/(p-E_{ECAL}))", "MPV/Width(E_{HCAL}/(p-E_{ECAL}))"};
1325
1326 gStyle->SetCanvasBorderMode(0);
1327 gStyle->SetCanvasColor(kWhite);
1328 gStyle->SetPadColor(kWhite);
1329 gStyle->SetFillColor(kWhite);
1330 gStyle->SetOptTitle(0);
1331 if (drawStatBox)
1332 gStyle->SetOptFit(10);
1333 else
1334 gStyle->SetOptFit(0);
1335
1336 char name[100], namep[100];
1337 TFile* file = new TFile(infile.c_str());
1338 for (int i = 0; i < 3; ++i) {
1339 std::vector<TH1D*> hists;
1340 std::vector<int> kks;
1341 double ymax(0.77);
1342 if (drawStatBox) {
1343 if (i == 0)
1344 gStyle->SetOptStat(1100);
1345 else
1346 gStyle->SetOptStat(10);
1347 } else {
1348 gStyle->SetOptStat(0);
1349 ymax = 0.89;
1350 }
1351 for (int k = 0; k < 6; ++k) {
1352 sprintf(name, "%s%s%s", prefix.c_str(), types[k].c_str(), names[i].c_str());
1353 TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1354 if (hist1 != nullptr) {
1355 hists.push_back((TH1D*)(hist1->Clone()));
1356 kks.push_back(k);
1357 }
1358 }
1359 if (hists.size() > 0) {
1360 sprintf(namep, "c_%s%s", prefix.c_str(), names[i].c_str());
1361 TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1362 TLegend* legend = new TLegend(0.44, 0.89 - 0.055 * hists.size(), 0.69, ymax);
1363 legend->SetFillColor(kWhite);
1364 pad->SetRightMargin(0.10);
1365 pad->SetTopMargin(0.10);
1366 double ymax(0.90);
1367 double dy = (i == 0) ? 0.13 : 0.08;
1368 for (unsigned int jk = 0; jk < hists.size(); ++jk) {
1369 int k = kks[jk];
1370 hists[jk]->GetXaxis()->SetTitle(xtitl[i].c_str());
1371 hists[jk]->GetYaxis()->SetTitle(ytitl[i].c_str());
1372 hists[jk]->GetYaxis()->SetLabelOffset(0.005);
1373 hists[jk]->GetYaxis()->SetLabelSize(0.035);
1374 hists[jk]->GetYaxis()->SetTitleOffset(1.15);
1375 if (i == 0) {
1376 hists[jk]->GetXaxis()->SetRangeUser(0.0, 2.5);
1377 } else if (i == 1) {
1378 hists[jk]->GetYaxis()->SetRangeUser(0.5, 2.0);
1379 } else {
1380 hists[jk]->GetYaxis()->SetRangeUser(0.0, 0.5);
1381 }
1382 hists[jk]->SetMarkerStyle(20);
1383 hists[jk]->SetMarkerColor(colors[k]);
1384 hists[jk]->SetLineColor(colors[k]);
1385 if (jk == 0)
1386 hists[jk]->Draw();
1387 else
1388 hists[jk]->Draw("sames");
1389 pad->Update();
1390 TPaveStats* st1 = (TPaveStats*)hists[jk]->GetListOfFunctions()->FindObject("stats");
1391 if (st1 != nullptr) {
1392 double ymin = ymax - dy;
1393 st1->SetLineColor(colors[k]);
1394 st1->SetTextColor(colors[k]);
1395 st1->SetY1NDC(ymin);
1396 st1->SetY2NDC(ymax);
1397 st1->SetX1NDC(0.70);
1398 st1->SetX2NDC(0.90);
1399 ymax = ymin;
1400 }
1401 sprintf(name, "%s%s", text.c_str(), types[k].c_str());
1402 legend->AddEntry(hists[jk], name, "lp");
1403 }
1404 legend->Draw("same");
1405 pad->Update();
1406 TPaveText* txt1 = new TPaveText(0.34, 0.825, 0.69, 0.895, "blNDC");
1407 txt1->SetFillColor(0);
1408 char txt[100];
1409 sprintf(txt, "Tracks with p = 40:60 GeV");
1410 txt1->AddText(txt);
1411 txt1->Draw("same");
1412 TPaveText* txt2 = new TPaveText(0.11, 0.825, 0.33, 0.895, "blNDC");
1413 txt2->SetFillColor(0);
1414 sprintf(txt, "CMS Preliminary");
1415 txt2->AddText(txt);
1416 txt2->Draw("same");
1417 if (!drawStatBox && i == 1) {
1418 double xmin = hists[0]->GetBinLowEdge(1);
1419 int nbin = hists[0]->GetNbinsX();
1420 double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
1421 TLine line = TLine(xmin, 1.0, xmax, 1.0);
1422 line.SetLineWidth(4);
1423 line.Draw("same");
1424 pad->Update();
1425 }
1426 pad->Modified();
1427 pad->Update();
1428 if (save > 0) {
1429 sprintf(name, "%s.pdf", pad->GetName());
1430 pad->Print(name);
1431 } else if (save < 0) {
1432 sprintf(name, "%s.C", pad->GetName());
1433 pad->Print(name);
1434 }
1435 }
1436 }
1437 }
1438
1439 void PlotTwoHists(std::string infile,
1440 std::string prefix1,
1441 std::string text1,
1442 std::string prefix2,
1443 std::string text2,
1444 std::string text0,
1445 int type = 0,
1446 int iname = 3,
1447 double lumi = 0,
1448 double ener = 13.0,
1449 int drawStatBox = 0,
1450 int save = 0) {
1451 int colors[2] = {2, 4};
1452 int numb[2] = {5, 1};
1453 std::string names0[5] = {"ratio00", "ratio00One", "etaB04One", "Z0", "W0"};
1454 std::string names1[5] = {"ratio10", "ratio10One", "etaB14One", "Z1", "W1"};
1455 std::string names2[5] = {"ratio30", "ratio30One", "etaB34One", "Z3", "W3"};
1456 std::string xtitl1[5] = {"E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "i#eta", "i#eta"};
1457 std::string ytitl1[5] = {
1458 "Tracks", "Tracks", "Tracks", "MPV(E_{HCAL}/(p-E_{ECAL}))", "MPV/Width(E_{HCAL}/(p-E_{ECAL}))"};
1459 std::string names3[1] = {"R"};
1460 std::string xtitl2[1] = {"RBX #"};
1461 std::string ytitl2[1] = {"MPV(E_{HCAL}/(p-E_{ECAL}))"};
1462
1463 gStyle->SetCanvasBorderMode(0);
1464 gStyle->SetCanvasColor(kWhite);
1465 gStyle->SetPadColor(kWhite);
1466 gStyle->SetFillColor(kWhite);
1467 gStyle->SetOptTitle(0);
1468 if ((drawStatBox / 10) % 10 > 0)
1469 gStyle->SetOptFit(10);
1470 else
1471 gStyle->SetOptFit(0);
1472
1473 if (type != 1)
1474 type = 0;
1475 char name[100], namep[100];
1476 TFile* file = new TFile(infile.c_str());
1477 for (int i = 0; i < numb[type]; ++i) {
1478 std::vector<TH1D*> hists;
1479 std::vector<int> kks;
1480 double ymax(0.77);
1481 if (drawStatBox % 10 > 0) {
1482 if (i != 2)
1483 gStyle->SetOptStat(1100);
1484 else
1485 gStyle->SetOptStat(10);
1486 } else {
1487 gStyle->SetOptStat(0);
1488 ymax = 0.82;
1489 }
1490 for (int k = 0; k < 2; ++k) {
1491 std::string prefix = (k == 0) ? prefix1 : prefix2;
1492 if (type == 0) {
1493 if (iname == 0)
1494 sprintf(name, "%s%s", prefix.c_str(), names0[i].c_str());
1495 else if (iname == 1)
1496 sprintf(name, "%s%s", prefix.c_str(), names1[i].c_str());
1497 else
1498 sprintf(name, "%s%s", prefix.c_str(), names2[i].c_str());
1499 } else {
1500 sprintf(name, "%s%s%d", prefix.c_str(), names3[i].c_str(), iname);
1501 }
1502 TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1503 if (hist1 != nullptr) {
1504 hists.push_back((TH1D*)(hist1->Clone()));
1505 kks.push_back(k);
1506 }
1507 }
1508 if (hists.size() == 2) {
1509 if (type == 0) {
1510 if (iname == 0)
1511 sprintf(namep, "c_%s%s%s", prefix1.c_str(), prefix2.c_str(), names0[i].c_str());
1512 else if (iname == 1)
1513 sprintf(namep, "c_%s%s%s", prefix1.c_str(), prefix2.c_str(), names1[i].c_str());
1514 else
1515 sprintf(namep, "c_%s%s%s", prefix1.c_str(), prefix2.c_str(), names2[i].c_str());
1516 } else {
1517 sprintf(namep, "c_%s%s%s%d", prefix1.c_str(), prefix2.c_str(), names3[i].c_str(), iname);
1518 }
1519 double ymax(0.90);
1520 double dy = (i == 0) ? 0.13 : 0.08;
1521 double ymx0 = (drawStatBox == 0) ? (ymax - .01) : (ymax - dy * hists.size() - .01);
1522 TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1523 TLegend* legend = new TLegend(0.64, ymx0 - 0.05 * hists.size(), 0.89, ymx0);
1524 legend->SetFillColor(kWhite);
1525 pad->SetRightMargin(0.10);
1526 pad->SetTopMargin(0.10);
1527 for (unsigned int jk = 0; jk < hists.size(); ++jk) {
1528 int k = kks[jk];
1529 hists[jk]->GetXaxis()->SetTitleSize(0.040);
1530 if (type == 0) {
1531 hists[jk]->GetXaxis()->SetTitle(xtitl1[i].c_str());
1532 hists[jk]->GetYaxis()->SetTitle(ytitl1[i].c_str());
1533 } else {
1534 hists[jk]->GetXaxis()->SetTitle(xtitl2[i].c_str());
1535 hists[jk]->GetYaxis()->SetTitle(ytitl2[i].c_str());
1536 }
1537 hists[jk]->GetYaxis()->SetLabelOffset(0.005);
1538 hists[jk]->GetYaxis()->SetLabelSize(0.035);
1539 hists[jk]->GetYaxis()->SetTitleSize(0.040);
1540 hists[jk]->GetYaxis()->SetTitleOffset(1.15);
1541 if ((type == 0) && (i != 3) && (i != 4))
1542 hists[jk]->GetXaxis()->SetRangeUser(0.0, 5.0);
1543 if (type == 0) {
1544 if (i == 3)
1545 hists[jk]->GetYaxis()->SetRangeUser(0.8, 1.2);
1546 else if (i == 4)
1547 hists[jk]->GetYaxis()->SetRangeUser(0.0, 0.5);
1548 }
1549 if (type != 0)
1550 hists[jk]->GetYaxis()->SetRangeUser(0.75, 1.2);
1551 hists[jk]->SetMarkerStyle(20);
1552 hists[jk]->SetMarkerColor(colors[k]);
1553 hists[jk]->SetLineColor(colors[k]);
1554 if (jk == 0)
1555 hists[jk]->Draw();
1556 else
1557 hists[jk]->Draw("sames");
1558 pad->Update();
1559 TPaveStats* st1 = (TPaveStats*)hists[jk]->GetListOfFunctions()->FindObject("stats");
1560 if (st1 != nullptr) {
1561 double ymin = ymax - dy;
1562 st1->SetLineColor(colors[k]);
1563 st1->SetTextColor(colors[k]);
1564 st1->SetY1NDC(ymin);
1565 st1->SetY2NDC(ymax);
1566 st1->SetX1NDC(0.70);
1567 st1->SetX2NDC(0.90);
1568 ymax = ymin;
1569 }
1570 if (k == 0)
1571 sprintf(name, "%s", text1.c_str());
1572 else
1573 sprintf(name, "%s", text2.c_str());
1574 legend->AddEntry(hists[jk], name, "lp");
1575 }
1576 legend->Draw("same");
1577 pad->Update();
1578 char txt[100];
1579 double xmi(0.10), xmx(0.895), ymx(0.95);
1580 if (lumi > 0.01) {
1581 xmx = 0.70;
1582 xmi = 0.30;
1583 TPaveText* txt0 = new TPaveText(0.705, 0.905, 0.90, 0.95, "blNDC");
1584 txt0->SetFillColor(0);
1585 sprintf(txt, "%4.1f TeV %5.1f fb^{-1}", ener, lumi);
1586 txt0->AddText(txt);
1587 txt0->Draw("same");
1588 }
1589 double ymi = ymx - 0.045;
1590 TPaveText* txt1 = new TPaveText(xmi, ymi, xmx, ymx, "blNDC");
1591 txt1->SetFillColor(0);
1592 if (iname == 0)
1593 sprintf(txt, "p = 20:30 GeV %s", text0.c_str());
1594 else
1595 sprintf(txt, "p = 40:60 GeV %s", text0.c_str());
1596 txt1->AddText(txt);
1597 txt1->Draw("same");
1598 ymi = (lumi > 0.1) ? 0.905 : 0.85;
1599 ymx = ymi + 0.045;
1600 TPaveText* txt2 = new TPaveText(0.105, ymi, 0.295, ymx, "blNDC");
1601 txt2->SetFillColor(0);
1602 sprintf(txt, "CMS Preliminary");
1603 txt2->AddText(txt);
1604 txt2->Draw("same");
1605 pad->Modified();
1606 if ((drawStatBox == 0) && (i == 3)) {
1607 double xmin = hists[0]->GetBinLowEdge(1);
1608 int nbin = hists[0]->GetNbinsX();
1609 double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
1610 TLine* line = new TLine(xmin, 1.0, xmax, 1.0);
1611 line->SetLineWidth(2);
1612 line->Draw("same");
1613 pad->Update();
1614 }
1615 pad->Update();
1616 if (save > 0) {
1617 sprintf(name, "%s.pdf", pad->GetName());
1618 pad->Print(name);
1619 } else if (save < 0) {
1620 sprintf(name, "%s.C", pad->GetName());
1621 pad->Print(name);
1622 }
1623 }
1624 }
1625 }
1626
1627 void PlotFiveHists(std::string infile,
1628 std::string text0,
1629 std::string prefix0,
1630 int type = 0,
1631 int iname = 3,
1632 int drawStatBox = 0,
1633 bool normalize = false,
1634 int save = 0,
1635 std::string prefix1 = "",
1636 std::string text1 = "",
1637 std::string prefix2 = "",
1638 std::string text2 = "",
1639 std::string prefix3 = "",
1640 std::string text3 = "",
1641 std::string prefix4 = "",
1642 std::string text4 = "",
1643 std::string prefix5 = "",
1644 std::string text5 = "") {
1645 int colors[5] = {2, 4, 6, 1, 7};
1646 int numb[3] = {5, 1, 4};
1647 std::string names0[5] = {"ratio00", "ratio00One", "etaB04", "Z0", "W0"};
1648 std::string names1[5] = {"ratio10", "ratio10One", "etaB14", "Z1", "W1"};
1649 std::string names2[5] = {"ratio30", "ratio30One", "etaB34", "Z3", "W3"};
1650 std::string xtitl1[5] = {"E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "i#eta", "i#eta"};
1651 std::string ytitl1[5] = {
1652 "Tracks", "Tracks", "Tracks", "MPV(E_{HCAL}/(p-E_{ECAL}))", "MPV/Width(E_{HCAL}/(p-E_{ECAL}))"};
1653 std::string names3[1] = {"R"};
1654 std::string xtitl2[1] = {"RBX #"};
1655 std::string ytitl2[1] = {"MPV(E_{HCAL}/(p-E_{ECAL}))"};
1656 std::string names4[4] = {"pp21", "pp22", "pp23", "pp24"};
1657 std::string xtitl3[4] = {"p (GeV)", "p (GeV)", "p (GeV)", "p (GeV)"};
1658 std::string ytitl3[4] = {"Tracks", "Tracks", "Tracks", "Tracks"};
1659 std::string title3[4] = {"Barrel", "Transition", "Endcap", "Combined"};
1660
1661 gStyle->SetCanvasBorderMode(0);
1662 gStyle->SetCanvasColor(kWhite);
1663 gStyle->SetPadColor(kWhite);
1664 gStyle->SetFillColor(kWhite);
1665 gStyle->SetOptTitle(0);
1666 if ((drawStatBox / 10) % 10 > 0)
1667 gStyle->SetOptFit(10);
1668 else
1669 gStyle->SetOptFit(0);
1670
1671 if (type != 1 && type != 2)
1672 type = 0;
1673 char name[100], namep[100];
1674 TFile* file = new TFile(infile.c_str());
1675 for (int i = 0; i < numb[type]; ++i) {
1676 std::vector<TH1D*> hists;
1677 std::vector<int> kks;
1678 std::vector<std::string> texts;
1679 double ymax(0.77);
1680 if (drawStatBox % 10 > 0) {
1681 if (type == 2)
1682 gStyle->SetOptStat(1110);
1683 else if (i != 3)
1684 gStyle->SetOptStat(1100);
1685 else
1686 gStyle->SetOptStat(10);
1687 } else {
1688 gStyle->SetOptStat(0);
1689 ymax = 0.82;
1690 }
1691 for (int k = 0; k < 5; ++k) {
1692 std::string prefix, text;
1693 if (k == 0) {
1694 prefix = prefix1;
1695 text = text1;
1696 } else if (k == 1) {
1697 prefix = prefix2;
1698 text = text2;
1699 } else if (k == 2) {
1700 prefix = prefix3;
1701 text = text3;
1702 } else if (k == 3) {
1703 prefix = prefix4;
1704 text = text4;
1705 } else {
1706 prefix = prefix5;
1707 text = text5;
1708 }
1709 if (prefix != "") {
1710 if (type == 0) {
1711 if (iname == 0)
1712 sprintf(name, "%s%s", prefix.c_str(), names0[i].c_str());
1713 else if (iname == 1)
1714 sprintf(name, "%s%s", prefix.c_str(), names1[i].c_str());
1715 else
1716 sprintf(name, "%s%s", prefix.c_str(), names2[i].c_str());
1717 } else if (type == 1) {
1718 sprintf(name, "%s%s%d", prefix.c_str(), names3[i].c_str(), iname);
1719 } else {
1720 sprintf(name, "%s%s", prefix.c_str(), names4[i].c_str());
1721 }
1722 TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
1723 if (hist1 != nullptr) {
1724 hists.push_back((TH1D*)(hist1->Clone()));
1725 kks.push_back(k);
1726 texts.push_back(text);
1727 }
1728 }
1729 }
1730 if (hists.size() > 0) {
1731 if (type == 0) {
1732 if (iname == 0)
1733 sprintf(namep, "c_%s%s", prefix0.c_str(), names0[i].c_str());
1734 else if (iname == 1)
1735 sprintf(namep, "c_%s%s", prefix0.c_str(), names1[i].c_str());
1736 else
1737 sprintf(namep, "c_%s%s", prefix0.c_str(), names2[i].c_str());
1738 } else if (type == 1) {
1739 sprintf(namep, "c_%s%s%d", prefix0.c_str(), names3[i].c_str(), iname);
1740 } else {
1741 sprintf(namep, "c_%s%s", prefix0.c_str(), names4[i].c_str());
1742 }
1743 double ymax(0.90);
1744 double dy = (i == 0 && type == 0) ? 0.13 : 0.08;
1745 double ymx0 = (drawStatBox == 0) ? (ymax - .01) : (ymax - dy * hists.size() - .01);
1746 TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1747 TLegend* legend = new TLegend(0.64, ymx0 - 0.05 * hists.size(), 0.89, ymx0);
1748 legend->SetFillColor(kWhite);
1749 pad->SetRightMargin(0.10);
1750 pad->SetTopMargin(0.10);
1751 for (unsigned int jk = 0; jk < hists.size(); ++jk) {
1752 int k = kks[jk];
1753 hists[jk]->GetXaxis()->SetTitleSize(0.040);
1754 if (type == 0) {
1755 hists[jk]->GetXaxis()->SetTitle(xtitl1[i].c_str());
1756 hists[jk]->GetYaxis()->SetTitle(ytitl1[i].c_str());
1757 } else if (type == 1) {
1758 hists[jk]->GetXaxis()->SetTitle(xtitl2[i].c_str());
1759 hists[jk]->GetYaxis()->SetTitle(ytitl2[i].c_str());
1760 } else {
1761 hists[jk]->GetXaxis()->SetTitle(xtitl3[i].c_str());
1762 hists[jk]->GetYaxis()->SetTitle(ytitl3[i].c_str());
1763 }
1764 hists[jk]->GetYaxis()->SetLabelOffset(0.005);
1765 hists[jk]->GetYaxis()->SetLabelSize(0.035);
1766 hists[jk]->GetYaxis()->SetTitleSize(0.040);
1767 hists[jk]->GetYaxis()->SetTitleOffset(1.15);
1768 if ((type == 0) && (i != 3) && (i != 4))
1769 hists[jk]->GetXaxis()->SetRangeUser(0.0, 2.5);
1770 if (type == 0) {
1771 if (i == 3)
1772 hists[jk]->GetYaxis()->SetRangeUser(0.8, 1.2);
1773 else if (i == 4)
1774 hists[jk]->GetYaxis()->SetRangeUser(0.0, 0.5);
1775 }
1776 if (type == 1)
1777 hists[jk]->GetYaxis()->SetRangeUser(0.75, 1.2);
1778 hists[jk]->SetMarkerStyle(20);
1779 hists[jk]->SetMarkerColor(colors[k]);
1780 hists[jk]->SetLineColor(colors[k]);
1781 if (normalize && ((type == 2) || ((type == 0) && (i != 3)))) {
1782 if (jk == 0)
1783 hists[jk]->DrawNormalized();
1784 else
1785 hists[jk]->DrawNormalized("sames");
1786 } else {
1787 if (jk == 0)
1788 hists[jk]->Draw();
1789 else
1790 hists[jk]->Draw("sames");
1791 }
1792 pad->Update();
1793 TPaveStats* st1 = (TPaveStats*)hists[jk]->GetListOfFunctions()->FindObject("stats");
1794 if (st1 != nullptr) {
1795 double ymin = ymax - dy;
1796 st1->SetLineColor(colors[k]);
1797 st1->SetTextColor(colors[k]);
1798 st1->SetY1NDC(ymin);
1799 st1->SetY2NDC(ymax);
1800 st1->SetX1NDC(0.70);
1801 st1->SetX2NDC(0.90);
1802 ymax = ymin;
1803 }
1804 sprintf(name, "%s", texts[jk].c_str());
1805 legend->AddEntry(hists[jk], name, "lp");
1806 }
1807 legend->Draw("same");
1808 pad->Update();
1809 TPaveText* txt1 = new TPaveText(0.10, 0.905, 0.80, 0.95, "blNDC");
1810 txt1->SetFillColor(0);
1811 char txt[100];
1812 if (type == 2) {
1813 sprintf(txt, "p = 40:60 GeV (%s)", title3[i].c_str());
1814 } else if (((type == 0) && (iname == 0))) {
1815 sprintf(txt, "p = 20:30 GeV %s", text0.c_str());
1816 } else {
1817 sprintf(txt, "p = 40:60 GeV %s", text0.c_str());
1818 }
1819 txt1->AddText(txt);
1820 txt1->Draw("same");
1821 TPaveText* txt2 = new TPaveText(0.11, 0.825, 0.33, 0.895, "blNDC");
1822 txt2->SetFillColor(0);
1823 sprintf(txt, "CMS Preliminary");
1824 txt2->AddText(txt);
1825 txt2->Draw("same");
1826 pad->Modified();
1827 if ((drawStatBox == 0) && (i == 3)) {
1828 double xmin = hists[0]->GetBinLowEdge(1);
1829 int nbin = hists[0]->GetNbinsX();
1830 double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
1831 TLine* line = new TLine(xmin, 1.0, xmax, 1.0);
1832 line->SetLineWidth(2);
1833 line->Draw("same");
1834 pad->Update();
1835 }
1836 pad->Update();
1837 if (save > 0) {
1838 sprintf(name, "%s.pdf", pad->GetName());
1839 pad->Print(name);
1840 } else if (save < 0) {
1841 sprintf(name, "%s.C", pad->GetName());
1842 pad->Print(name);
1843 }
1844 }
1845 }
1846 }
1847
1848 void PlotHistCorrResults(std::string infile, std::string text, std::string prefixF, int save = 0) {
1849 std::string name[5] = {"Eta1Bf", "Eta2Bf", "Eta1Af", "Eta2Af", "Cvg0"};
1850 std::string title[5] = {"Mean at the start of itertions",
1851 "Median at the start of itertions",
1852 "Mean at the end of itertions",
1853 "Median at the end of itertions",
1854 ""};
1855 int type[5] = {0, 0, 0, 0, 1};
1856
1857 gStyle->SetCanvasBorderMode(0);
1858 gStyle->SetCanvasColor(kWhite);
1859 gStyle->SetPadColor(kWhite);
1860 gStyle->SetFillColor(kWhite);
1861 gStyle->SetOptTitle(0);
1862 gStyle->SetOptStat(10);
1863 gStyle->SetOptFit(10);
1864 TFile* file = new TFile(infile.c_str());
1865 char namep[100];
1866 for (int k = 0; k < 5; ++k) {
1867 TH1D* hist1 = (TH1D*)file->FindObjectAny(name[k].c_str());
1868 if (hist1 != nullptr) {
1869 TH1D* hist = (TH1D*)(hist1->Clone());
1870 sprintf(namep, "c_%s%s", prefixF.c_str(), name[k].c_str());
1871 TCanvas* pad = new TCanvas(namep, namep, 700, 500);
1872 pad->SetRightMargin(0.10);
1873 pad->SetTopMargin(0.10);
1874 hist->GetYaxis()->SetLabelOffset(0.005);
1875 hist->GetYaxis()->SetTitleOffset(1.20);
1876 double xmin = hist->GetBinLowEdge(1);
1877 int nbin = hist->GetNbinsX();
1878 double xmax = hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin);
1879 std::cout << hist->GetTitle() << " Bins " << nbin << ":" << xmin << ":" << xmax << std::endl;
1880 double xlow(0.12), ylow(0.82);
1881 char txt[100], option[2];
1882 if (type[k] == 0) {
1883 sprintf(namep, "f_%s", name[k].c_str());
1884 TF1* func = new TF1(namep, "pol0", xmin, xmax);
1885 hist->Fit(func, "+QWL", "");
1886 if (text == "")
1887 sprintf(txt, "%s", title[k].c_str());
1888 else
1889 sprintf(txt, "%s (balancing the %s)", title[k].c_str(), text.c_str());
1890 sprintf(option, "%s", "");
1891 } else {
1892 hist->GetXaxis()->SetTitle("Iterations");
1893 hist->GetYaxis()->SetTitle("Convergence");
1894 hist->SetMarkerStyle(20);
1895 hist->SetMarkerColor(2);
1896 hist->SetMarkerSize(0.8);
1897 xlow = 0.50;
1898 ylow = 0.86;
1899 sprintf(txt, "(%s)", text.c_str());
1900 sprintf(option, "%s", "p");
1901 }
1902 hist->Draw(option);
1903 pad->Update();
1904 TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
1905 if (st1 != nullptr) {
1906 st1->SetY1NDC(ylow);
1907 st1->SetY2NDC(0.90);
1908 st1->SetX1NDC(0.70);
1909 st1->SetX2NDC(0.90);
1910 }
1911 TPaveText* txt1 = new TPaveText(xlow, 0.82, 0.68, 0.88, "blNDC");
1912 txt1->SetFillColor(0);
1913 txt1->AddText(txt);
1914 txt1->Draw("same");
1915 pad->Modified();
1916 pad->Update();
1917 if (save > 0) {
1918 sprintf(namep, "%s.pdf", pad->GetName());
1919 pad->Print(namep);
1920 } else if (save < 0) {
1921 sprintf(namep, "%s.C", pad->GetName());
1922 pad->Print(namep);
1923 }
1924 }
1925 }
1926 }
1927
1928 void PlotHistCorrFactor(char* infile,
1929 std::string text,
1930 std::string prefixF = "",
1931 double scale = 1.0,
1932 int nmin = 100,
1933 bool dataMC = false,
1934 bool drawStatBox = true,
1935 int iformat = 0,
1936 int save = 0) {
1937 std::map<int, cfactors> cfacs;
1938 int etamin(100), etamax(-100), maxdepth(0);
1939 readCorrFactors(infile, scale, cfacs, etamin, etamax, maxdepth, iformat);
1940
1941 gStyle->SetCanvasBorderMode(0);
1942 gStyle->SetCanvasColor(kWhite);
1943 gStyle->SetPadColor(kWhite);
1944 gStyle->SetFillColor(kWhite);
1945 gStyle->SetOptTitle(0);
1946 if (drawStatBox) {
1947 gStyle->SetOptStat(10);
1948 gStyle->SetOptFit(10);
1949 } else {
1950 gStyle->SetOptStat(0);
1951 gStyle->SetOptFit(0);
1952 }
1953 int colors[7] = {1, 6, 4, 7, 2, 9, 3};
1954 int mtype[7] = {20, 21, 22, 23, 24, 33, 25};
1955 int nbin = etamax - etamin + 1;
1956 std::vector<TH1D*> hists;
1957 std::vector<int> entries;
1958 char name[100];
1959 double dy(0);
1960 int fits(0);
1961 for (int j = 0; j < maxdepth; ++j) {
1962 sprintf(name, "hd%d", j + 1);
1963 TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
1964 int nent(0);
1965 for (std::map<int, cfactors>::const_iterator itr = cfacs.begin(); itr != cfacs.end(); ++itr) {
1966 if ((itr->second).depth == j + 1) {
1967 int ieta = (itr->second).ieta;
1968 int bin = ieta - etamin + 1;
1969 float val = (itr->second).corrf;
1970 float dvl = (itr->second).dcorr;
1971 h->SetBinContent(bin, val);
1972 h->SetBinError(bin, dvl);
1973 nent++;
1974 }
1975 }
1976 if (nent > nmin) {
1977 fits++;
1978 dy += 0.025;
1979 sprintf(name, "hdf%d", j + 1);
1980 TF1* func = new TF1(name, "pol0", etamin, etamax);
1981 h->Fit(func, "+QWLR", "");
1982 }
1983 h->SetLineColor(colors[j]);
1984 h->SetMarkerColor(colors[j]);
1985 h->SetMarkerStyle(mtype[j]);
1986 h->GetXaxis()->SetTitle("i#eta");
1987 h->GetYaxis()->SetTitle("Correction Factor");
1988 h->GetYaxis()->SetLabelOffset(0.005);
1989 h->GetYaxis()->SetTitleOffset(1.20);
1990 h->GetYaxis()->SetRangeUser(0.0, 2.0);
1991 hists.push_back(h);
1992 entries.push_back(nent);
1993 dy += 0.025;
1994 }
1995 sprintf(name, "c_%sCorrFactor", prefixF.c_str());
1996 TCanvas* pad = new TCanvas(name, name, 700, 500);
1997 pad->SetRightMargin(0.10);
1998 pad->SetTopMargin(0.10);
1999 double yh = 0.90;
2000
2001 double yl = 0.15;
2002 TLegend* legend = new TLegend(0.35, yl, 0.65, yl + 0.04 * hists.size());
2003 legend->SetFillColor(kWhite);
2004 for (unsigned int k = 0; k < hists.size(); ++k) {
2005 if (k == 0)
2006 hists[k]->Draw("");
2007 else
2008 hists[k]->Draw("sames");
2009 pad->Update();
2010 TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
2011 if (st1 != nullptr) {
2012 dy = (entries[k] > nmin) ? 0.05 : 0.025;
2013 st1->SetLineColor(colors[k]);
2014 st1->SetTextColor(colors[k]);
2015 st1->SetY1NDC(yh - dy);
2016 st1->SetY2NDC(yh);
2017 st1->SetX1NDC(0.70);
2018 st1->SetX2NDC(0.90);
2019 yh -= dy;
2020 }
2021 sprintf(name, "Depth %d (%s)", k + 1, text.c_str());
2022 legend->AddEntry(hists[k], name, "lp");
2023 }
2024 legend->Draw("same");
2025 pad->Update();
2026 if (fits < 1) {
2027 double xmin = hists[0]->GetBinLowEdge(1);
2028 int nbin = hists[0]->GetNbinsX();
2029 double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
2030 TLine* line = new TLine(xmin, 1.0, xmax, 1.0);
2031 line->SetLineColor(9);
2032 line->SetLineWidth(2);
2033 line->SetLineStyle(2);
2034 line->Draw("same");
2035 pad->Modified();
2036 pad->Update();
2037 }
2038 char txt1[30];
2039 double xmax = (dataMC) ? 0.33 : 0.44;
2040 TPaveText* txt2 = new TPaveText(0.11, 0.85, xmax, 0.89, "blNDC");
2041 txt2->SetFillColor(0);
2042 if (dataMC)
2043 sprintf(txt1, "CMS Preliminary");
2044 else
2045 sprintf(txt1, "CMS Simulation Preliminary");
2046 txt2->AddText(txt1);
2047 txt2->Draw("same");
2048 pad->Modified();
2049 pad->Update();
2050 if (save > 0) {
2051 sprintf(name, "%s.pdf", pad->GetName());
2052 pad->Print(name);
2053 } else if (save < 0) {
2054 sprintf(name, "%s.C", pad->GetName());
2055 pad->Print(name);
2056 }
2057 }
2058
2059 void PlotHistCorrAsymmetry(char* infile, std::string text, std::string prefixF = "", int iformat = 0, int save = 0) {
2060 std::map<int, cfactors> cfacs;
2061 int etamin(100), etamax(-100), maxdepth(0);
2062 double scale(1.0);
2063 readCorrFactors(infile, scale, cfacs, etamin, etamax, maxdepth, iformat);
2064
2065 gStyle->SetCanvasBorderMode(0);
2066 gStyle->SetCanvasColor(kWhite);
2067 gStyle->SetPadColor(kWhite);
2068 gStyle->SetFillColor(kWhite);
2069 gStyle->SetOptTitle(0);
2070 gStyle->SetOptStat(0);
2071 gStyle->SetOptFit(10);
2072 int colors[6] = {1, 6, 4, 7, 2, 9};
2073 int mtype[6] = {20, 21, 22, 23, 24, 33};
2074 int nbin = etamax + 1;
2075 std::vector<TH1D*> hists;
2076 std::vector<int> entries;
2077 char name[100];
2078 double dy(0);
2079 for (int j = 0; j < maxdepth; ++j) {
2080 sprintf(name, "hd%d", j + 1);
2081 TH1D* h = new TH1D(name, name, nbin, 0, etamax);
2082 int nent(0);
2083 for (std::map<int, cfactors>::const_iterator itr = cfacs.begin(); itr != cfacs.end(); ++itr) {
2084 if ((itr->second).depth == j + 1) {
2085 int ieta = (itr->second).ieta;
2086 float vl1 = (itr->second).corrf;
2087 float dv1 = (itr->second).dcorr;
2088 if (ieta > 0) {
2089 for (std::map<int, cfactors>::const_iterator ktr = cfacs.begin(); ktr != cfacs.end(); ++ktr) {
2090 if (((ktr->second).depth == j + 1) && ((ktr->second).ieta == -ieta)) {
2091 float vl2 = (ktr->second).corrf;
2092 float dv2 = (ktr->second).dcorr;
2093 float val = 2.0 * (vl1 - vl2) / (vl1 + vl2);
2094 float dvl = (4.0 * sqrt(vl1 * vl1 * dv2 * dv2 + vl2 * vl2 * dv1 * dv1) / ((vl1 + vl2) * (vl1 + vl2)));
2095 int bin = ieta;
2096 h->SetBinContent(bin, val);
2097 h->SetBinError(bin, dvl);
2098 nent++;
2099 }
2100 }
2101 }
2102 }
2103 }
2104 h->SetLineColor(colors[j]);
2105 h->SetMarkerColor(colors[j]);
2106 h->SetMarkerStyle(mtype[j]);
2107 h->GetXaxis()->SetTitle("i#eta");
2108 h->GetYaxis()->SetTitle("Asymmetry in Correction Factor");
2109 h->GetYaxis()->SetLabelOffset(0.005);
2110 h->GetYaxis()->SetTitleOffset(1.20);
2111 h->GetYaxis()->SetRangeUser(-0.25, 0.25);
2112 hists.push_back(h);
2113 entries.push_back(nent);
2114 dy += 0.025;
2115 }
2116 sprintf(name, "c_%sCorrAsymmetry", prefixF.c_str());
2117 TCanvas* pad = new TCanvas(name, name, 700, 500);
2118 pad->SetRightMargin(0.10);
2119 pad->SetTopMargin(0.10);
2120 double yh = 0.90;
2121 double yl = yh - 0.035 * hists.size() - dy - 0.01;
2122 TLegend* legend = new TLegend(0.60, yl, 0.90, yl + 0.035 * hists.size());
2123 legend->SetFillColor(kWhite);
2124 for (unsigned int k = 0; k < hists.size(); ++k) {
2125 if (k == 0)
2126 hists[k]->Draw("");
2127 else
2128 hists[k]->Draw("sames");
2129 pad->Update();
2130 sprintf(name, "Depth %d (%s)", k + 1, text.c_str());
2131 legend->AddEntry(hists[k], name, "lp");
2132 }
2133 legend->Draw("same");
2134 pad->Update();
2135
2136 TLine* line = new TLine(0.0, 0.0, etamax, 0.0);
2137 line->SetLineColor(9);
2138 line->SetLineWidth(2);
2139 line->SetLineStyle(2);
2140 line->Draw("same");
2141 pad->Update();
2142
2143 if (save > 0) {
2144 sprintf(name, "%s.pdf", pad->GetName());
2145 pad->Print(name);
2146 } else if (save < 0) {
2147 sprintf(name, "%s.C", pad->GetName());
2148 pad->Print(name);
2149 }
2150 }
2151
2152 void PlotHistCorrFactors(char* infile1,
2153 std::string text1,
2154 char* infile2,
2155 std::string text2,
2156 char* infile3,
2157 std::string text3,
2158 char* infile4,
2159 std::string text4,
2160 char* infile5,
2161 std::string text5,
2162 std::string prefixF,
2163 bool ratio = false,
2164 bool drawStatBox = true,
2165 int nmin = 100,
2166 bool dataMC = false,
2167 int year = 2018,
2168 int iformat = 0,
2169 int save = 0) {
2170 std::map<int, cfactors> cfacs[5];
2171 std::vector<std::string> texts;
2172 int nfile(0), etamin(100), etamax(-100), maxdepth(0);
2173 const char* blank("");
2174 if (infile1 != blank) {
2175 readCorrFactors(infile1, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2176 if (cfacs[nfile].size() > 0) {
2177 texts.push_back(text1);
2178 ++nfile;
2179 }
2180 }
2181 if (infile2 != blank) {
2182 readCorrFactors(infile2, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2183 if (cfacs[nfile].size() > 0) {
2184 texts.push_back(text2);
2185 ++nfile;
2186 }
2187 }
2188 if (infile3 != blank) {
2189 readCorrFactors(infile3, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2190 if (cfacs[nfile].size() > 0) {
2191 texts.push_back(text3);
2192 ++nfile;
2193 }
2194 }
2195 if (infile4 != blank) {
2196 readCorrFactors(infile4, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2197 if (cfacs[nfile].size() > 0) {
2198 texts.push_back(text4);
2199 ++nfile;
2200 }
2201 }
2202 if (infile5 != blank) {
2203 readCorrFactors(infile5, 1.0, cfacs[nfile], etamin, etamax, maxdepth, iformat);
2204 if (cfacs[nfile].size() > 0) {
2205 texts.push_back(text5);
2206 ++nfile;
2207 }
2208 }
2209
2210 if (nfile > 0) {
2211 gStyle->SetCanvasBorderMode(0);
2212 gStyle->SetCanvasColor(kWhite);
2213 gStyle->SetPadColor(kWhite);
2214 gStyle->SetFillColor(kWhite);
2215 gStyle->SetOptTitle(0);
2216 if ((!ratio) && drawStatBox) {
2217 gStyle->SetOptStat(10);
2218 gStyle->SetOptFit(10);
2219 } else {
2220 gStyle->SetOptStat(0);
2221 gStyle->SetOptFit(0);
2222 }
2223 int colors[6] = {1, 6, 4, 2, 7, 9};
2224 int mtype[6] = {20, 24, 22, 23, 21, 33};
2225 int nbin = etamax - etamin + 1;
2226 std::vector<TH1D*> hists;
2227 std::vector<int> entries, htype, depths;
2228 std::vector<double> fitr;
2229 char name[100];
2230 double dy(0);
2231 int fits(0);
2232 int nline(0);
2233 if (ratio) {
2234 for (int ih = 1; ih < nfile; ++ih) {
2235 for (int j = 0; j < maxdepth; ++j) {
2236 sprintf(name, "h%dd%d", ih, j + 1);
2237 TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2238 double sumNum(0), sumDen(0);
2239 std::map<int, cfactors>::const_iterator ktr = cfacs[ih].begin();
2240 for (std::map<int, cfactors>::const_iterator itr = cfacs[0].begin(); itr != cfacs[0].end(); ++itr, ++ktr) {
2241 int dep = (itr->second).depth;
2242 if (dep == j + 1) {
2243 int ieta = (itr->second).ieta;
2244 int bin = ieta - etamin + 1;
2245 float val = (itr->second).corrf / (ktr->second).corrf;
2246 float dvl =
2247 val *
2248 sqrt((((itr->second).dcorr * (itr->second).dcorr) / ((itr->second).corrf * (itr->second).corrf)) +
2249 (((ktr->second).dcorr * (ktr->second).dcorr) / ((ktr->second).corrf * (ktr->second).corrf)));
2250 h->SetBinContent(bin, val);
2251 h->SetBinError(bin, dvl);
2252 sumNum += (val / (dvl * dvl));
2253 sumDen += (1.0 / (dvl * dvl));
2254 }
2255 }
2256 double fit = (sumDen > 0) ? (sumNum / sumDen) : 1.0;
2257 std::cout << "Fit to Pol0: " << fit << std::endl;
2258 h->SetLineColor(colors[ih]);
2259 h->SetMarkerColor(colors[ih]);
2260 h->SetMarkerStyle(mtype[j]);
2261 h->SetMarkerSize(0.9);
2262 h->GetXaxis()->SetTitle("i#eta");
2263 sprintf(name, "CF_{%s}/CF_{Set}", texts[0].c_str());
2264 h->GetYaxis()->SetTitle(name);
2265 h->GetYaxis()->SetLabelOffset(0.005);
2266 h->GetYaxis()->SetTitleSize(0.036);
2267 h->GetYaxis()->SetTitleOffset(1.20);
2268 h->GetYaxis()->SetRangeUser(0.80, 1.20);
2269 hists.push_back(h);
2270 fitr.push_back(fit);
2271 htype.push_back(ih);
2272 depths.push_back(j + 1);
2273 }
2274 if (ih == 1)
2275 nline += hists.size();
2276 else
2277 ++nline;
2278 }
2279 } else {
2280 for (int k1 = 0; k1 < nfile; ++k1) {
2281 for (int j = 0; j < maxdepth; ++j) {
2282 sprintf(name, "h%dd%d", k1, j + 1);
2283 TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2284 int nent(0);
2285 for (std::map<int, cfactors>::const_iterator itr = cfacs[k1].begin(); itr != cfacs[k1].end(); ++itr) {
2286 int dep = (itr->second).depth;
2287 if (dep == j + 1) {
2288 int ieta = (itr->second).ieta;
2289 int bin = ieta - etamin + 1;
2290 float val = (itr->second).corrf;
2291 float dvl = (itr->second).dcorr;
2292 h->SetBinContent(bin, val);
2293 h->SetBinError(bin, dvl);
2294 nent++;
2295 }
2296 }
2297 if (nent > nmin) {
2298 fits++;
2299 if (drawStatBox)
2300 dy += 0.025;
2301 sprintf(name, "h%ddf%d", k1, j + 1);
2302 TF1* func = new TF1(name, "pol0", etamin, etamax);
2303 h->Fit(func, "+QWLR", "");
2304 }
2305 h->SetLineColor(colors[k1]);
2306 h->SetMarkerColor(colors[k1]);
2307 h->SetMarkerStyle(mtype[j]);
2308 h->SetMarkerSize(0.9);
2309 h->GetXaxis()->SetTitle("i#eta");
2310 h->GetYaxis()->SetTitle("Correction Factor");
2311 h->GetYaxis()->SetLabelOffset(0.005);
2312 h->GetYaxis()->SetTitleOffset(1.20);
2313 h->GetYaxis()->SetRangeUser(0.5, 1.5);
2314 hists.push_back(h);
2315 entries.push_back(nent);
2316 if (drawStatBox)
2317 dy += 0.025;
2318 htype.push_back(k1);
2319 depths.push_back(j + 1);
2320 }
2321 if (k1 == 0)
2322 nline += hists.size();
2323 else
2324 ++nline;
2325 }
2326 }
2327 if (ratio)
2328 sprintf(name, "c_Corr%sRatio", prefixF.c_str());
2329 else
2330 sprintf(name, "c_Corr%s", prefixF.c_str());
2331 TCanvas* pad = new TCanvas(name, name, 700, 500);
2332 pad->SetRightMargin(0.10);
2333 pad->SetTopMargin(0.10);
2334 double yh = 0.90;
2335 double yl = yh - 0.035 * hists.size() - dy - 0.01;
2336 TLegend* legend = new TLegend(0.55, yl, 0.90, yl + 0.035 * nline);
2337 legend->SetFillColor(kWhite);
2338 for (unsigned int k = 0; k < hists.size(); ++k) {
2339 if (k == 0)
2340 hists[k]->Draw("");
2341 else
2342 hists[k]->Draw("sames");
2343 pad->Update();
2344 int k1 = htype[k];
2345 if (!ratio) {
2346 TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
2347 if (st1 != nullptr) {
2348 dy = (entries[k] > nmin) ? 0.05 : 0.025;
2349 st1->SetLineColor(colors[k1]);
2350 st1->SetTextColor(colors[k1]);
2351 st1->SetY1NDC(yh - dy);
2352 st1->SetY2NDC(yh);
2353 st1->SetX1NDC(0.70);
2354 st1->SetX2NDC(0.90);
2355 yh -= dy;
2356 }
2357 sprintf(name, "Depth %d (%s)", depths[k], texts[k1].c_str());
2358 } else {
2359 sprintf(name, "Depth %d (%s Mean = %5.3f)", depths[k], texts[k1].c_str(), fitr[k]);
2360 }
2361 if ((depths[k] == 1) || (k1 == 0))
2362 legend->AddEntry(hists[k], name, "lp");
2363 }
2364 legend->Draw("same");
2365 TPaveText* txt0 = new TPaveText(0.12, 0.84, 0.49, 0.89, "blNDC");
2366 txt0->SetFillColor(0);
2367 char txt[40];
2368 if (dataMC)
2369 sprintf(txt, "CMS Preliminary (%d)", year);
2370 else
2371 sprintf(txt, "CMS Simulation Preliminary (%d)", year);
2372 txt0->AddText(txt);
2373 txt0->Draw("same");
2374 pad->Update();
2375 if (fits < 1) {
2376 double xmin = hists[0]->GetBinLowEdge(1);
2377 int nbin = hists[0]->GetNbinsX();
2378 double xmax = hists[0]->GetBinLowEdge(nbin) + hists[0]->GetBinWidth(nbin);
2379 TLine* line = new TLine(xmin, 1.0, xmax, 1.0);
2380 line->SetLineColor(9);
2381 line->SetLineWidth(2);
2382 line->SetLineStyle(2);
2383 line->Draw("same");
2384 pad->Update();
2385 }
2386 if (save > 0) {
2387 sprintf(name, "%s.pdf", pad->GetName());
2388 pad->Print(name);
2389 } else if (save < 0) {
2390 sprintf(name, "%s.C", pad->GetName());
2391 pad->Print(name);
2392 }
2393 }
2394 }
2395
2396 void PlotHistCorrSys(std::string infilec, int conds, std::string text, int save = 0) {
2397 char fname[100];
2398 int iformat(0);
2399 sprintf(fname, "%s_cond0.txt", infilec.c_str());
2400 int etamin(100), etamax(-100), maxdepth(0);
2401 std::map<int, cfactors> cfacs;
2402 readCorrFactors(fname, 1.0, cfacs, etamin, etamax, maxdepth, iformat);
2403
2404 if (cfacs.size() > 0) {
2405
2406 std::map<int, cfactors> errfacs;
2407 for (int i = 0; i < conds; ++i) {
2408 sprintf(fname, "%s_cond%d.txt", infilec.c_str(), i + 1);
2409 std::map<int, cfactors> cfacx;
2410 int etamin1(100), etamax1(-100), maxdepth1(0);
2411 readCorrFactors(fname, 1.0, cfacx, etamin1, etamax1, maxdepth1, iformat);
2412 for (std::map<int, cfactors>::const_iterator itr1 = cfacx.begin(); itr1 != cfacx.end(); ++itr1) {
2413 std::map<int, cfactors>::iterator itr2 = errfacs.find(itr1->first);
2414 float corrf = (itr1->second).corrf;
2415 if (itr2 == errfacs.end()) {
2416 errfacs[itr1->first] = cfactors(1, 0, corrf, corrf * corrf);
2417 } else {
2418 int nent = (itr2->second).ieta + 1;
2419 float c1 = (itr2->second).corrf + corrf;
2420 float c2 = (itr2->second).dcorr + (corrf * corrf);
2421 errfacs[itr1->first] = cfactors(nent, 0, c1, c2);
2422 }
2423 }
2424 }
2425
2426 for (std::map<int, cfactors>::iterator itr = errfacs.begin(); itr != errfacs.end(); ++itr) {
2427 int nent = (itr->second).ieta;
2428 float mean = (itr->second).corrf / nent;
2429 float rms2 = (itr->second).dcorr / nent - (mean * mean);
2430 float rms = rms2 > 0 ? sqrt(rms2) : 0;
2431 errfacs[itr->first] = cfactors(nent, 0, mean, rms);
2432 }
2433
2434 int k(0);
2435 gStyle->SetCanvasBorderMode(0);
2436 gStyle->SetCanvasColor(kWhite);
2437 gStyle->SetPadColor(kWhite);
2438 gStyle->SetFillColor(kWhite);
2439 gStyle->SetOptTitle(0);
2440 gStyle->SetOptStat(10);
2441 gStyle->SetOptFit(10);
2442 int colors[6] = {1, 6, 4, 7, 2, 9};
2443 int mtype[6] = {20, 21, 22, 23, 24, 33};
2444 std::vector<TH1D*> hists;
2445 char name[100];
2446 int nbin = etamax - etamin + 1;
2447 for (int j = 0; j < maxdepth; ++j) {
2448 sprintf(name, "hd%d", j + 1);
2449 TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2450 h->SetLineColor(colors[j]);
2451 h->SetMarkerColor(colors[j]);
2452 h->SetMarkerStyle(mtype[j]);
2453 h->GetXaxis()->SetTitle("i#eta");
2454 h->GetYaxis()->SetTitle("Correction Factor");
2455 h->GetYaxis()->SetLabelOffset(0.005);
2456 h->GetYaxis()->SetTitleOffset(1.20);
2457 h->GetYaxis()->SetRangeUser(0.0, 2.0);
2458 hists.push_back(h);
2459 }
2460 for (std::map<int, cfactors>::iterator itr = cfacs.begin(); itr != cfacs.end(); ++itr, ++k) {
2461 std::map<int, cfactors>::iterator itr2 = errfacs.find(itr->first);
2462 float mean2 = (itr2 == errfacs.end()) ? 0 : (itr2->second).corrf;
2463 float ersys = (itr2 == errfacs.end()) ? 0 : (itr2->second).dcorr;
2464 float erstt = (itr->second).dcorr;
2465 float ertot = sqrt(erstt * erstt + ersys * ersys);
2466 float mean = (itr->second).corrf;
2467 int ieta = (itr->second).ieta;
2468 int depth = (itr->second).depth;
2469 std::cout << "[" << k << "] " << ieta << " " << depth << " " << mean << ":" << mean2 << " " << erstt << ":"
2470 << ersys << ":" << ertot << std::endl;
2471 int bin = ieta - etamin + 1;
2472 hists[depth - 1]->SetBinContent(bin, mean);
2473 hists[depth - 1]->SetBinError(bin, ertot);
2474 }
2475 TCanvas* pad = new TCanvas("CorrFactorSys", "CorrFactorSys", 700, 500);
2476 pad->SetRightMargin(0.10);
2477 pad->SetTopMargin(0.10);
2478 double yh = 0.90;
2479 double yl = yh - 0.050 * hists.size() - 0.01;
2480 TLegend* legend = new TLegend(0.60, yl, 0.90, yl + 0.025 * hists.size());
2481 legend->SetFillColor(kWhite);
2482 for (unsigned int k = 0; k < hists.size(); ++k) {
2483 if (k == 0)
2484 hists[k]->Draw("");
2485 else
2486 hists[k]->Draw("sames");
2487 pad->Update();
2488 TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
2489 if (st1 != nullptr) {
2490 st1->SetLineColor(colors[k]);
2491 st1->SetTextColor(colors[k]);
2492 st1->SetY1NDC(yh - 0.025);
2493 st1->SetY2NDC(yh);
2494 st1->SetX1NDC(0.70);
2495 st1->SetX2NDC(0.90);
2496 yh -= 0.025;
2497 }
2498 sprintf(name, "Depth %d (%s)", k + 1, text.c_str());
2499 legend->AddEntry(hists[k], name, "lp");
2500 }
2501 legend->Draw("same");
2502 pad->Update();
2503 if (save > 0) {
2504 sprintf(name, "%s.pdf", pad->GetName());
2505 pad->Print(name);
2506 } else if (save < 0) {
2507 sprintf(name, "%s.C", pad->GetName());
2508 pad->Print(name);
2509 }
2510 }
2511 }
2512
2513 void PlotHistCorrLumis(std::string infilec, int conds, double lumi, int save = 0) {
2514 char fname[100];
2515 int iformat(0);
2516 sprintf(fname, "%s_0.txt", infilec.c_str());
2517 std::map<int, cfactors> cfacs;
2518 int etamin(100), etamax(-100), maxdepth(0);
2519 readCorrFactors(fname, 1.0, cfacs, etamin, etamax, maxdepth, iformat);
2520 int nbin = etamax - etamin + 1;
2521 std::cout << "Max Depth " << maxdepth << " and " << nbin << " eta bins for " << etamin << ":" << etamax << std::endl;
2522
2523
2524 int colors[8] = {4, 2, 6, 7, 1, 9, 3, 5};
2525 int mtype[8] = {20, 21, 22, 23, 24, 25, 26, 27};
2526 if (cfacs.size() > 0) {
2527
2528 std::vector<TH1D*> hists;
2529 char name[100];
2530 for (int i = 0; i < conds; ++i) {
2531 int ih = (int)(hists.size());
2532 for (int j = 0; j < maxdepth; ++j) {
2533 sprintf(name, "hd%d%d", j + 1, i);
2534 TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2535 h->SetLineColor(colors[j]);
2536 h->SetMarkerColor(colors[j]);
2537 h->SetMarkerStyle(mtype[i]);
2538 h->SetMarkerSize(0.9);
2539 h->GetXaxis()->SetTitle("i#eta");
2540 h->GetYaxis()->SetTitle("Fractional Error");
2541 h->GetYaxis()->SetLabelOffset(0.005);
2542 h->GetYaxis()->SetTitleOffset(1.20);
2543 h->GetYaxis()->SetRangeUser(0.0, 0.10);
2544 hists.push_back(h);
2545 }
2546 sprintf(fname, "%s_%d.txt", infilec.c_str(), i);
2547 int etamin1(100), etamax1(-100), maxdepth1(0);
2548 readCorrFactors(fname, 1.0, cfacs, etamin1, etamax1, maxdepth1, iformat);
2549 for (std::map<int, cfactors>::const_iterator itr = cfacs.begin(); itr != cfacs.end(); ++itr) {
2550 double value = (itr->second).dcorr / (itr->second).corrf;
2551 int bin = (itr->second).ieta - etamin + 1;
2552 hists[ih + (itr->second).depth - 1]->SetBinContent(bin, value);
2553 hists[ih + (itr->second).depth - 1]->SetBinError(bin, 0.0001);
2554 }
2555 }
2556
2557
2558 gStyle->SetCanvasBorderMode(0);
2559 gStyle->SetCanvasColor(kWhite);
2560 gStyle->SetPadColor(kWhite);
2561 gStyle->SetFillColor(kWhite);
2562 gStyle->SetOptTitle(0);
2563 gStyle->SetOptStat(0);
2564 gStyle->SetOptFit(0);
2565 TCanvas* pad = new TCanvas("CorrFactorErr", "CorrFactorErr", 700, 500);
2566 pad->SetRightMargin(0.10);
2567 pad->SetTopMargin(0.10);
2568 double yh(0.89);
2569 TLegend* legend = new TLegend(0.60, yh - 0.04 * conds, 0.89, yh);
2570 legend->SetFillColor(kWhite);
2571 double lumic(lumi);
2572 for (unsigned int k = 0; k < hists.size(); ++k) {
2573 if (k == 0)
2574 hists[k]->Draw("");
2575 else
2576 hists[k]->Draw("sames");
2577 pad->Update();
2578 if (k % maxdepth == 0) {
2579 sprintf(name, "L = %5.2f fb^{-1}", lumic);
2580 legend->AddEntry(hists[k], name, "lp");
2581 lumic *= 0.5;
2582 }
2583 }
2584 legend->Draw("same");
2585 pad->Update();
2586 if (save > 0) {
2587 sprintf(name, "%s.pdf", pad->GetName());
2588 pad->Print(name);
2589 } else if (save < 0) {
2590 sprintf(name, "%s.C", pad->GetName());
2591 pad->Print(name);
2592 }
2593 }
2594 }
2595
2596 void PlotHistCorrRel(char* infile1,
2597 char* infile2,
2598 std::string text1,
2599 std::string text2,
2600 int iformat1 = 0,
2601 int iformat2 = 0,
2602 int save = 0) {
2603 std::map<int, cfactors> cfacs1, cfacs2;
2604 int etamin(100), etamax(-100), maxdepth(0);
2605 readCorrFactors(infile1, 1.0, cfacs1, etamin, etamax, maxdepth, iformat1);
2606 readCorrFactors(infile2, 1.0, cfacs2, etamin, etamax, maxdepth, iformat2);
2607 std::map<int, std::pair<cfactors, cfactors> > cfacs;
2608 for (std::map<int, cfactors>::iterator itr = cfacs1.begin(); itr != cfacs1.end(); ++itr) {
2609 std::map<int, cfactors>::iterator ktr = cfacs2.find(itr->first);
2610 if (ktr == cfacs2.end()) {
2611 cfactors fac2(((itr->second).ieta), ((itr->second).depth), 0, -1);
2612 cfacs[itr->first] = std::pair<cfactors, cfactors>((itr->second), fac2);
2613 } else {
2614 cfactors fac2(ktr->second);
2615 cfacs[itr->first] = std::pair<cfactors, cfactors>((itr->second), fac2);
2616 }
2617 }
2618 for (std::map<int, cfactors>::iterator itr = cfacs2.begin(); itr != cfacs2.end(); ++itr) {
2619 std::map<int, cfactors>::const_iterator ktr = cfacs1.find(itr->first);
2620 if (ktr == cfacs1.end()) {
2621 cfactors fac1(((itr->second).ieta), ((itr->second).depth), 0, -1);
2622 cfacs[itr->first] = std::pair<cfactors, cfactors>(fac1, (itr->second));
2623 }
2624 }
2625
2626
2627 if ((cfacs1.size() > 0) && (cfacs2.size() > 0)) {
2628 int k(0);
2629 gStyle->SetCanvasBorderMode(0);
2630 gStyle->SetCanvasColor(kWhite);
2631 gStyle->SetPadColor(kWhite);
2632 gStyle->SetFillColor(kWhite);
2633 gStyle->SetOptTitle(0);
2634 gStyle->SetOptStat(10);
2635 gStyle->SetOptFit(10);
2636 int colors[6] = {1, 6, 4, 7, 2, 9};
2637 int mtype[6] = {20, 21, 22, 23, 24, 33};
2638 std::vector<TH1D*> hists;
2639 char name[100];
2640 int nbin = etamax - etamin + 1;
2641 for (int i = 0; i < 2; ++i) {
2642 for (int j = 0; j < maxdepth; ++j) {
2643 int j1 = (i == 0) ? j : maxdepth + j;
2644 sprintf(name, "hd%d%d", i, j + 1);
2645 TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2646 h->SetLineColor(colors[j1]);
2647 h->SetMarkerColor(colors[j1]);
2648 h->SetMarkerStyle(mtype[i]);
2649 h->GetXaxis()->SetTitle("i#eta");
2650 h->GetYaxis()->SetTitle("Correction Factor");
2651 h->GetYaxis()->SetLabelOffset(0.005);
2652 h->GetYaxis()->SetTitleOffset(1.20);
2653 h->GetYaxis()->SetRangeUser(0.0, 2.0);
2654 hists.push_back(h);
2655 }
2656 }
2657 for (std::map<int, std::pair<cfactors, cfactors> >::iterator it = cfacs.begin(); it != cfacs.end(); ++it, ++k) {
2658 float mean1 = (it->second).first.corrf;
2659 float error1 = (it->second).first.dcorr;
2660 float mean2 = (it->second).second.corrf;
2661 float error2 = (it->second).second.dcorr;
2662 int ieta = (it->second).first.ieta;
2663 int depth = (it->second).first.depth;
2664
2665
2666
2667
2668
2669 int bin = ieta - etamin + 1;
2670 if (error1 >= 0) {
2671 hists[depth - 1]->SetBinContent(bin, mean1);
2672 hists[depth - 1]->SetBinError(bin, error1);
2673 }
2674 if (error2 >= 0) {
2675 hists[maxdepth + depth - 1]->SetBinContent(bin, mean2);
2676 hists[maxdepth + depth - 1]->SetBinError(bin, error2);
2677 }
2678 }
2679 TCanvas* pad = new TCanvas("CorrFactors", "CorrFactors", 700, 500);
2680 pad->SetRightMargin(0.10);
2681 pad->SetTopMargin(0.10);
2682 double yh = 0.90;
2683 double yl = yh - 0.050 * hists.size() - 0.01;
2684 TLegend* legend = new TLegend(0.60, yl, 0.90, yl + 0.025 * hists.size());
2685 legend->SetFillColor(kWhite);
2686 for (unsigned int k = 0; k < hists.size(); ++k) {
2687 if (k == 0)
2688 hists[k]->Draw("");
2689 else
2690 hists[k]->Draw("sames");
2691 pad->Update();
2692 TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
2693 if (st1 != nullptr) {
2694 st1->SetLineColor(colors[k]);
2695 st1->SetTextColor(colors[k]);
2696 st1->SetY1NDC(yh - 0.025);
2697 st1->SetY2NDC(yh);
2698 st1->SetX1NDC(0.70);
2699 st1->SetX2NDC(0.90);
2700 yh -= 0.025;
2701 }
2702 if (k < (unsigned int)(maxdepth)) {
2703 sprintf(name, "Depth %d (%s)", k + 1, text1.c_str());
2704 } else {
2705 sprintf(name, "Depth %d (%s)", k - maxdepth + 1, text2.c_str());
2706 }
2707 legend->AddEntry(hists[k], name, "lp");
2708 }
2709 legend->Draw("same");
2710 pad->Update();
2711 if (save > 0) {
2712 sprintf(name, "%s.pdf", pad->GetName());
2713 pad->Print(name);
2714 } else if (save < 0) {
2715 sprintf(name, "%s.C", pad->GetName());
2716 pad->Print(name);
2717 }
2718 }
2719 }
2720
2721 void PlotHistCorrDepth(char* infile1,
2722 char* infile2,
2723 std::string text1,
2724 std::string text2,
2725 int depth,
2726 int ietamax,
2727 int iformat1 = 0,
2728 int iformat2 = 0,
2729 int save = 0,
2730 int debug = 1) {
2731 std::map<int, cfactors> cfacs1, cfacs2;
2732 int etamin(100), etamax(-100), maxdepth(0);
2733 readCorrFactors(infile1, 1.0, cfacs1, etamin, etamax, maxdepth, iformat1, (debug > 1));
2734 readCorrFactors(infile2, 1.0, cfacs2, etamin, etamax, maxdepth, iformat2, (debug > 1));
2735
2736 double sumNum(0), sumDen(0), ratMax(0);
2737 int npt0(0), npt1(0);
2738 for (std::map<int, cfactors>::iterator itr = cfacs1.begin(); itr != cfacs1.end(); ++itr) {
2739 std::map<int, cfactors>::iterator ktr = cfacs2.find(itr->first);
2740 if ((ktr != cfacs2.end()) && ((itr->second).depth == depth)) {
2741 double er1 = (itr->second).dcorr / (itr->second).corrf;
2742 double er2 = (ktr->second).dcorr / (ktr->second).corrf;
2743 double tmp = (ktr->second).corrf / (itr->second).corrf;
2744 double dtmp = tmp * sqrt(er1 * er1 + er2 * er2);
2745 double rat = (tmp > 1.0) ? 1.0 / tmp : tmp;
2746 double drt = (tmp > 1.0) ? dtmp / (tmp * tmp) : dtmp;
2747 rat = fabs(1.0 - rat);
2748 ratMax = std::max(ratMax, rat);
2749 ++npt0;
2750 if (debug > 0)
2751 std::cout << std::hex << (itr->first) << std::dec << " ieta:depth" << (itr->second).ieta << ":"
2752 << (itr->second).depth << " Corr " << (itr->second).corrf << ":" << (ktr->second).corrf << " Ratio "
2753 << rat << ":" << drt << std::endl;
2754 if (std::abs((itr->second).ieta) <= ietamax) {
2755 sumNum += (rat / (drt * drt));
2756 sumDen += (1.0 / (drt * drt));
2757 ++npt1;
2758 }
2759 }
2760 }
2761 sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0;
2762 sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0;
2763 std::cout << "Get Ratio of mean for " << npt0 << ":" << npt1 << " points for depth " << depth << " Mean " << sumNum
2764 << " +- " << sumDen << " (Maximum " << ratMax << ")" << std::endl;
2765
2766 gStyle->SetCanvasBorderMode(0);
2767 gStyle->SetCanvasColor(kWhite);
2768 gStyle->SetPadColor(kWhite);
2769 gStyle->SetFillColor(kWhite);
2770 gStyle->SetOptTitle(0);
2771 gStyle->SetOptStat(0);
2772 gStyle->SetOptFit(0);
2773 int colors[2] = {1, 2};
2774 int mtype[2] = {20, 24};
2775 int nbin = etamax - etamin + 1;
2776 std::vector<TH1D*> hists;
2777 char name[100];
2778 for (int j = 0; j < 2; ++j) {
2779 sprintf(name, "hd%d", (j + 1));
2780 TH1D* h = new TH1D(name, name, nbin, etamin, etamax);
2781 if (j == 0) {
2782 for (std::map<int, cfactors>::const_iterator itr = cfacs1.begin(); itr != cfacs1.end(); ++itr) {
2783 if ((itr->second).depth == depth) {
2784 int ieta = (itr->second).ieta;
2785 int bin = ieta - etamin + 1;
2786 float val = (itr->second).corrf;
2787 float dvl = (itr->second).dcorr;
2788 h->SetBinContent(bin, val);
2789 h->SetBinError(bin, dvl);
2790 }
2791 }
2792 } else {
2793 for (std::map<int, cfactors>::const_iterator itr = cfacs2.begin(); itr != cfacs2.end(); ++itr) {
2794 if ((itr->second).depth == depth) {
2795 int ieta = (itr->second).ieta;
2796 int bin = ieta - etamin + 1;
2797 float val = (itr->second).corrf;
2798 float dvl = (itr->second).dcorr;
2799 h->SetBinContent(bin, val);
2800 h->SetBinError(bin, dvl);
2801 }
2802 }
2803 }
2804 h->SetLineColor(colors[j]);
2805 h->SetMarkerColor(colors[j]);
2806 h->SetMarkerStyle(mtype[j]);
2807 h->GetXaxis()->SetTitle("i#eta");
2808 h->GetYaxis()->SetTitle("Correction Factor");
2809 h->GetYaxis()->SetLabelOffset(0.005);
2810 h->GetYaxis()->SetTitleOffset(1.20);
2811 h->GetYaxis()->SetRangeUser(0.0, 2.0);
2812 hists.push_back(h);
2813 }
2814 sprintf(name, "c_CorrFactorDepth%d", depth);
2815 TCanvas* pad = new TCanvas(name, name, 700, 500);
2816 pad->SetRightMargin(0.10);
2817 pad->SetTopMargin(0.10);
2818 double yl = 0.25;
2819 TLegend* legend = new TLegend(0.25, yl, 0.85, yl + 0.04 * hists.size());
2820 legend->SetFillColor(kWhite);
2821 for (unsigned int k = 0; k < hists.size(); ++k) {
2822 if (k == 0) {
2823 hists[k]->Draw("");
2824 sprintf(name, "Depth %d (%s)", depth, text1.c_str());
2825 } else {
2826 hists[k]->Draw("sames");
2827 sprintf(name, "Depth %d (%s)", depth, text2.c_str());
2828 }
2829 pad->Update();
2830 legend->AddEntry(hists[k], name, "lp");
2831 }
2832 legend->Draw("same");
2833 pad->Update();
2834 if (save > 0) {
2835 sprintf(name, "%s.pdf", pad->GetName());
2836 pad->Print(name);
2837 } else if (save < 0) {
2838 sprintf(name, "%s.C", pad->GetName());
2839 pad->Print(name);
2840 }
2841 }
2842
2843 void PlotFourHists(std::string infile,
2844 std::string prefix0,
2845 int type = 0,
2846 int drawStatBox = 0,
2847 bool normalize = false,
2848 int save = 0,
2849 std::string prefix1 = "",
2850 std::string text1 = "",
2851 std::string prefix2 = "",
2852 std::string text2 = "",
2853 std::string prefix3 = "",
2854 std::string text3 = "",
2855 std::string prefix4 = "",
2856 std::string text4 = "") {
2857 int colors[4] = {2, 4, 6, 1};
2858 std::string names[5] = {"eta03", "eta13", "eta23", "eta33", "eta43"};
2859 std::string xtitle[5] = {"i#eta", "i#eta", "i#eta", "i#eta", "i#eta"};
2860 std::string ytitle[5] = {"Tracks", "Tracks", "Tracks", "Tracks", "Tracks"};
2861 std::string title[5] = {"All Tracks (p = 40:60 GeV)",
2862 "Good Quality Tracks (p = 40:60 GeV)",
2863 "Selected Tracks (p = 40:60 GeV)",
2864 "Isolated Tracks (p = 40:60 GeV)",
2865 "Isolated MIP Tracks (p = 40:60 GeV)"};
2866
2867 gStyle->SetCanvasBorderMode(0);
2868 gStyle->SetCanvasColor(kWhite);
2869 gStyle->SetPadColor(kWhite);
2870 gStyle->SetFillColor(kWhite);
2871 gStyle->SetOptTitle(0);
2872 gStyle->SetOptFit(0);
2873 if (drawStatBox == 0)
2874 gStyle->SetOptStat(0);
2875 else
2876 gStyle->SetOptStat(1110);
2877
2878 if (type < 0 || type > 4)
2879 type = 0;
2880 char name[100], namep[100];
2881 TFile* file = new TFile(infile.c_str());
2882
2883 std::vector<TH1D*> hists;
2884 std::vector<int> kks;
2885 std::vector<std::string> texts;
2886 for (int k = 0; k < 4; ++k) {
2887 std::string prefix, text;
2888 if (k == 0) {
2889 prefix = prefix1;
2890 text = text1;
2891 } else if (k == 1) {
2892 prefix = prefix2;
2893 text = text2;
2894 } else if (k == 2) {
2895 prefix = prefix3;
2896 text = text3;
2897 } else {
2898 prefix = prefix4;
2899 text = text4;
2900 }
2901 if (prefix != "") {
2902 sprintf(name, "%s%s", prefix.c_str(), names[type].c_str());
2903 TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
2904 if (hist1 != nullptr) {
2905 hists.push_back((TH1D*)(hist1->Clone()));
2906 kks.push_back(k);
2907 texts.push_back(text);
2908 }
2909 }
2910 }
2911 if (hists.size() > 0) {
2912 sprintf(namep, "c_%s%s", prefix0.c_str(), names[type].c_str());
2913 double ymax(0.90), dy(0.13);
2914 double ymx0 = (drawStatBox == 0) ? (ymax - .01) : (ymax - dy * hists.size() - .01);
2915 TCanvas* pad = new TCanvas(namep, namep, 700, 500);
2916 TLegend* legend = new TLegend(0.64, ymx0 - 0.05 * hists.size(), 0.89, ymx0);
2917 legend->SetFillColor(kWhite);
2918 pad->SetRightMargin(0.10);
2919 pad->SetTopMargin(0.10);
2920 for (unsigned int jk = 0; jk < hists.size(); ++jk) {
2921 int k = kks[jk];
2922 hists[jk]->GetXaxis()->SetTitleSize(0.040);
2923 hists[jk]->GetXaxis()->SetTitle(xtitle[type].c_str());
2924 hists[jk]->GetYaxis()->SetTitle(ytitle[type].c_str());
2925 hists[jk]->GetYaxis()->SetLabelOffset(0.005);
2926 hists[jk]->GetYaxis()->SetLabelSize(0.035);
2927 hists[jk]->GetYaxis()->SetTitleSize(0.040);
2928 hists[jk]->GetYaxis()->SetTitleOffset(1.15);
2929 hists[jk]->SetMarkerStyle(20);
2930 hists[jk]->SetMarkerColor(colors[k]);
2931 hists[jk]->SetLineColor(colors[k]);
2932 if (normalize) {
2933 if (jk == 0)
2934 hists[jk]->DrawNormalized();
2935 else
2936 hists[jk]->DrawNormalized("sames");
2937 } else {
2938 if (jk == 0)
2939 hists[jk]->Draw();
2940 else
2941 hists[jk]->Draw("sames");
2942 }
2943 pad->Update();
2944 TPaveStats* st1 = (TPaveStats*)hists[jk]->GetListOfFunctions()->FindObject("stats");
2945 if (st1 != nullptr) {
2946 double ymin = ymax - dy;
2947 st1->SetLineColor(colors[k]);
2948 st1->SetTextColor(colors[k]);
2949 st1->SetY1NDC(ymin);
2950 st1->SetY2NDC(ymax);
2951 st1->SetX1NDC(0.70);
2952 st1->SetX2NDC(0.90);
2953 ymax = ymin;
2954 }
2955 sprintf(name, "%s", texts[jk].c_str());
2956 legend->AddEntry(hists[jk], name, "lp");
2957 }
2958 legend->Draw("same");
2959 pad->Update();
2960 TPaveText* txt1 = new TPaveText(0.10, 0.905, 0.80, 0.95, "blNDC");
2961 txt1->SetFillColor(0);
2962 char txt[100];
2963 sprintf(txt, "%s", title[type].c_str());
2964 txt1->AddText(txt);
2965 txt1->Draw("same");
2966
2967
2968
2969
2970
2971
2972
2973 pad->Modified();
2974 pad->Update();
2975 if (save > 0) {
2976 sprintf(name, "%s.pdf", pad->GetName());
2977 pad->Print(name);
2978 } else if (save < 0) {
2979 sprintf(name, "%s.C", pad->GetName());
2980 pad->Print(name);
2981 }
2982 }
2983 }
2984
2985 void PlotPUCorrHists(std::string infile = "corrfac.root",
2986 std::string prefix = "",
2987 int drawStatBox = 0,
2988 bool approve = true,
2989 int save = 0) {
2990 std::string name1[4] = {"W0", "W1", "W2", "P"};
2991 std::string name2[4] = {"All", "Barrel", "Endcap", ""};
2992 std::string name3[2] = {"", "p = 40:60 GeV"};
2993 std::string name4[2] = {"Loose Isolation", "Tight Isolation"};
2994 std::string xtitle[4] = {"Correction Factor", "Correction Factor", "Correction Factor", "i#eta"};
2995 std::string ytitle[4] = {"Tracks", "Tracks", "Tracks", "Correction Factor"};
2996
2997 gStyle->SetCanvasBorderMode(0);
2998 gStyle->SetCanvasColor(kWhite);
2999 gStyle->SetPadColor(kWhite);
3000 gStyle->SetFillColor(kWhite);
3001 gStyle->SetOptTitle(0);
3002 gStyle->SetOptFit(0);
3003 if (drawStatBox == 0)
3004 gStyle->SetOptStat(0);
3005 else
3006 gStyle->SetOptStat(1110);
3007
3008 char name[100], namep[100], title[100];
3009 TFile* file = new TFile(infile.c_str());
3010
3011 if (file != nullptr) {
3012 for (int i1 = 0; i1 < 4; ++i1) {
3013 for (int i2 = 0; i2 < 2; ++i2) {
3014 for (int i3 = 0; i3 < 2; ++i3) {
3015 sprintf(name, "%s%d%d", name1[i1].c_str(), i2, i3);
3016 if (i2 == 0)
3017 sprintf(title, "%s Tracks Selected with %s", name2[i1].c_str(), name4[i3].c_str());
3018 else
3019 sprintf(title, "%s Tracks Selected with %s (%s)", name2[i1].c_str(), name4[i3].c_str(), name3[i2].c_str());
3020 TH1D* hist1(nullptr);
3021 TProfile* hist2(nullptr);
3022 if (i1 != 3) {
3023 TH1D* hist = (TH1D*)file->FindObjectAny(name);
3024 if (hist != nullptr) {
3025 hist1 = (TH1D*)(hist->Clone());
3026 hist1->GetXaxis()->SetTitleSize(0.040);
3027 hist1->GetXaxis()->SetTitle(xtitle[i1].c_str());
3028 hist1->GetYaxis()->SetTitle(ytitle[i1].c_str());
3029 hist1->GetYaxis()->SetLabelOffset(0.005);
3030 hist1->GetYaxis()->SetLabelSize(0.035);
3031 hist1->GetYaxis()->SetTitleSize(0.040);
3032 hist1->GetYaxis()->SetTitleOffset(1.15);
3033 }
3034 } else {
3035 TProfile* hist = (TProfile*)file->FindObjectAny(name);
3036 if (hist != nullptr) {
3037 hist2 = (TProfile*)(hist->Clone());
3038 hist2->GetXaxis()->SetTitleSize(0.040);
3039 hist2->GetXaxis()->SetTitle(xtitle[i1].c_str());
3040 hist2->GetYaxis()->SetTitle(ytitle[i1].c_str());
3041 hist2->GetYaxis()->SetLabelOffset(0.005);
3042 hist2->GetYaxis()->SetLabelSize(0.035);
3043 hist2->GetYaxis()->SetTitleSize(0.040);
3044 hist2->GetYaxis()->SetTitleOffset(1.15);
3045
3046 hist2->SetMarkerStyle(20);
3047 }
3048 }
3049 if ((hist1 != nullptr) || (hist2 != nullptr)) {
3050 sprintf(namep, "c_%s%s", name, prefix.c_str());
3051 TCanvas* pad = new TCanvas(namep, namep, 700, 500);
3052 pad->SetRightMargin(0.10);
3053 pad->SetTopMargin(0.10);
3054 if (hist1 != nullptr) {
3055 pad->SetLogy();
3056 hist1->Draw();
3057 pad->Update();
3058 TPaveStats* st1 = (TPaveStats*)hist1->GetListOfFunctions()->FindObject("stats");
3059 if (st1 != nullptr) {
3060 st1->SetY1NDC(0.77);
3061 st1->SetY2NDC(0.90);
3062 st1->SetX1NDC(0.70);
3063 st1->SetX2NDC(0.90);
3064 }
3065 } else {
3066 hist2->Draw();
3067 pad->Update();
3068 }
3069 TPaveText* txt1 = new TPaveText(0.10, 0.905, 0.80, 0.95, "blNDC");
3070 txt1->SetFillColor(0);
3071 char txt[100];
3072 sprintf(txt, "%s", title);
3073 txt1->AddText(txt);
3074 txt1->Draw("same");
3075 if (approve) {
3076 double xoff = (i1 == 3) ? 0.11 : 0.22;
3077 TPaveText* txt2 = new TPaveText(xoff, 0.825, xoff + 0.22, 0.895, "blNDC");
3078 txt2->SetFillColor(0);
3079 sprintf(txt, "CMS Preliminary");
3080 txt2->AddText(txt);
3081 txt2->Draw("same");
3082 }
3083 pad->Modified();
3084 pad->Update();
3085 if (save > 0) {
3086 sprintf(name, "%s.pdf", pad->GetName());
3087 pad->Print(name);
3088 } else if (save < 0) {
3089 sprintf(name, "%s.C", pad->GetName());
3090 pad->Print(name);
3091 }
3092 }
3093 }
3094 }
3095 }
3096 }
3097 }
3098
3099 void PlotHistCorr(const char* infile,
3100 std::string prefix,
3101 std::string text0,
3102 int eta = 0,
3103 int mode = 1,
3104 bool drawStatBox = true,
3105 int save = 0) {
3106 gStyle->SetCanvasBorderMode(0);
3107 gStyle->SetCanvasColor(kWhite);
3108 gStyle->SetPadColor(kWhite);
3109 gStyle->SetFillColor(kWhite);
3110 gStyle->SetOptTitle(0);
3111 if (drawStatBox)
3112 gStyle->SetOptStat(1100);
3113 else
3114 gStyle->SetOptStat(0);
3115
3116 std::string tags[3] = {"UnNoPU", "UnPU", "Cor"};
3117 std::string text[3] = {"Uncorrected no PU", "Uncorrected PU", "Corrected PU"};
3118 int colors[3] = {1, 4, 2};
3119 int styles[3] = {1, 3, 2};
3120 TFile* file = new TFile(infile);
3121 if (mode < 0 || mode > 2)
3122 mode = 1;
3123 int etamin = (eta == 0) ? -27 : eta;
3124 int etamax = (eta == 0) ? 27 : eta;
3125 for (int ieta = etamin; ieta <= etamax; ++ieta) {
3126 char name[20];
3127 double yh(0.90), dy(0.09);
3128 double yh1 = drawStatBox ? (yh - 3 * dy - 0.01) : (yh - 0.01);
3129 TLegend* legend = new TLegend(0.55, yh1 - 0.15, 0.89, yh1);
3130 legend->SetFillColor(kWhite);
3131 sprintf(name, "c_%sEovp%d", prefix.c_str(), ieta);
3132 TCanvas* pad = new TCanvas(name, name, 700, 500);
3133 pad->SetRightMargin(0.10);
3134 pad->SetTopMargin(0.10);
3135 TH1D* hist[3];
3136 double ymax(0);
3137 for (int k = 0; k < 3; ++k) {
3138 if (k < 2)
3139 sprintf(name, "EovP_ieta%d%s", ieta, tags[k].c_str());
3140 else
3141 sprintf(name, "EovP_ieta%dCor%dPU", ieta, mode);
3142 TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
3143 if (hist1 != nullptr) {
3144 hist[k] = (TH1D*)(hist1->Clone());
3145 ymax = std::max(ymax, (hist1->GetMaximum()));
3146 }
3147 }
3148 int imax = 10 * (2 + int(0.1 * ymax));
3149 for (int k = 0; k < 3; ++k) {
3150 hist[k]->GetYaxis()->SetLabelOffset(0.005);
3151 hist[k]->GetYaxis()->SetTitleOffset(1.20);
3152 hist[k]->GetXaxis()->SetTitle("E/p");
3153 hist[k]->GetYaxis()->SetTitle("Tracks");
3154 hist[k]->SetLineColor(colors[k]);
3155 hist[k]->SetLineStyle(styles[k]);
3156 hist[k]->GetYaxis()->SetRangeUser(0.0, imax);
3157 if (k == 0)
3158 hist[k]->Draw();
3159 else
3160 hist[k]->Draw("sames");
3161 legend->AddEntry(hist[k], text[k].c_str(), "lp");
3162 pad->Update();
3163 if (drawStatBox) {
3164 TPaveStats* st1 = (TPaveStats*)hist[k]->GetListOfFunctions()->FindObject("stats");
3165 if (st1 != nullptr) {
3166 st1->SetLineColor(colors[k]);
3167 st1->SetTextColor(colors[k]);
3168 st1->SetY1NDC(yh - dy);
3169 st1->SetY2NDC(yh);
3170 st1->SetX1NDC(0.70);
3171 st1->SetX2NDC(0.90);
3172 yh -= dy;
3173 }
3174 }
3175 }
3176 pad->Update();
3177 legend->Draw("same");
3178 pad->Update();
3179 TPaveText* txt1 = new TPaveText(0.10, 0.905, 0.80, 0.95, "blNDC");
3180 txt1->SetFillColor(0);
3181 char title[100];
3182 sprintf(title, "%s for i#eta = %d", text0.c_str(), ieta);
3183 txt1->AddText(title);
3184 txt1->Draw("same");
3185 pad->Modified();
3186 pad->Update();
3187 if (save > 0) {
3188 sprintf(name, "%s.pdf", pad->GetName());
3189 pad->Print(name);
3190 } else if (save < 0) {
3191 sprintf(name, "%s.C", pad->GetName());
3192 pad->Print(name);
3193 }
3194 }
3195 }
3196
3197 void PlotPropertyHist(const char* infile,
3198 std::string prefix,
3199 std::string text,
3200 int etaMax = 25,
3201 double lumi = 0,
3202 double ener = 13.0,
3203 bool dataMC = false,
3204 bool drawStatBox = true,
3205 int save = 0) {
3206 std::string name0[3] = {"energyE2", "energyH2", "energyP2"};
3207 std::string title0[3] = {"Energy in ECAL", "Energy in HCAL", "Track Momentum"};
3208 std::string xtitl0[3] = {"Energy (GeV)", "Energy (GeV)", "p (GeV)"};
3209 std::string name1[5] = {"eta02", "eta12", "eta22", "eta32", "eta42"};
3210 std::string name10[5] = {"eta0", "eta1", "eta2", "eta3", "eta4"};
3211 std::string xtitl1 = "i#eta";
3212 std::string name2[5] = {"p0", "p1", "p2", "p3", "p4"};
3213 std::string xtitl2 = "p (GeV)";
3214 std::string title1[5] = {"Tracks with p=40:60 GeV",
3215 "Good Quality Tracks with p=40:60 GeV",
3216 "Selected Tracks with p=40:60 GeV",
3217 "Isolated Tracks with p=40:60 GeV",
3218 "Isolated Tracks with p=40:60 GeV and MIPS in ECAL"};
3219 std::string title2[5] = {
3220 "All Tracks", "Good Quality Tracks", "Selected Tracks", "Isolated Tracks", "Isolated Tracks with MIPS in ECAL"};
3221 std::string ytitle = "Tracks";
3222
3223 gStyle->SetCanvasBorderMode(0);
3224 gStyle->SetCanvasColor(kWhite);
3225 gStyle->SetPadColor(kWhite);
3226 gStyle->SetFillColor(kWhite);
3227 gStyle->SetOptTitle(0);
3228 if (drawStatBox)
3229 gStyle->SetOptStat(1110);
3230 else
3231 gStyle->SetOptStat(0);
3232 gStyle->SetOptFit(0);
3233
3234 TFile* file = new TFile(infile);
3235 char name[100], namep[100];
3236 for (int k = 1; k <= etaMax; ++k) {
3237 for (int j = 0; j < 3; ++j) {
3238 sprintf(name, "%s%s%d", prefix.c_str(), name0[j].c_str(), k);
3239 TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
3240 if (hist1 != nullptr) {
3241 TH1D* hist = (TH1D*)(hist1->Clone());
3242 double ymin(0.90);
3243 sprintf(namep, "c_%s", name);
3244 TCanvas* pad = new TCanvas(namep, namep, 700, 500);
3245 pad->SetLogy();
3246 pad->SetRightMargin(0.10);
3247 pad->SetTopMargin(0.10);
3248 hist->GetXaxis()->SetTitleSize(0.04);
3249 hist->GetXaxis()->SetTitle(xtitl0[j].c_str());
3250 hist->GetYaxis()->SetTitle(ytitle.c_str());
3251 hist->GetYaxis()->SetLabelOffset(0.005);
3252 hist->GetYaxis()->SetTitleSize(0.04);
3253 hist->GetYaxis()->SetLabelSize(0.035);
3254 hist->GetYaxis()->SetTitleOffset(1.10);
3255 hist->SetMarkerStyle(20);
3256 hist->Draw();
3257 pad->Update();
3258 TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
3259 if (st1 != nullptr) {
3260 st1->SetY1NDC(0.78);
3261 st1->SetY2NDC(0.90);
3262 st1->SetX1NDC(0.65);
3263 st1->SetX2NDC(0.90);
3264 }
3265 pad->Update();
3266 double ymx(0.96), xmi(0.35), xmx(0.95);
3267 char txt[100];
3268 if (lumi > 0.1) {
3269 ymx = ymin - 0.005;
3270 xmi = 0.45;
3271 TPaveText* txt0 = new TPaveText(0.65, 0.91, 0.90, 0.96, "blNDC");
3272 txt0->SetFillColor(0);
3273 sprintf(txt, "%4.1f TeV %5.1f fb^{-1}", ener, lumi);
3274 txt0->AddText(txt);
3275 txt0->Draw("same");
3276 }
3277 double ymi = ymx - 0.05;
3278 TPaveText* txt1 = new TPaveText(xmi, ymi, xmx, ymx, "blNDC");
3279 txt1->SetFillColor(0);
3280 if (text == "") {
3281 sprintf(txt, "%s", title0[j].c_str());
3282 } else {
3283 sprintf(txt, "%s (%s)", title0[j].c_str(), text.c_str());
3284 }
3285 txt1->AddText(txt);
3286 txt1->Draw("same");
3287 double xmax = (dataMC) ? 0.24 : 0.35;
3288 ymi = 0.91;
3289 ymx = ymi + 0.05;
3290 TPaveText* txt2 = new TPaveText(0.02, ymi, xmax, ymx, "blNDC");
3291 txt2->SetFillColor(0);
3292 if (dataMC)
3293 sprintf(txt, "CMS Preliminary");
3294 else
3295 sprintf(txt, "CMS Simulation Preliminary");
3296 txt2->AddText(txt);
3297 txt2->Draw("same");
3298 pad->Modified();
3299 pad->Update();
3300 if (save > 0) {
3301 sprintf(name, "%s.pdf", pad->GetName());
3302 pad->Print(name);
3303 } else if (save < 0) {
3304 sprintf(name, "%s.C", pad->GetName());
3305 pad->Print(name);
3306 }
3307 }
3308 }
3309 }
3310
3311 for (int k = 0; k < 3; ++k) {
3312 for (int j = 0; j < 5; ++j) {
3313 if (k == 0)
3314 sprintf(name, "%s%s", prefix.c_str(), name1[j].c_str());
3315 else if (k == 1)
3316 sprintf(name, "%s%s", prefix.c_str(), name10[j].c_str());
3317 else
3318 sprintf(name, "%s%s", prefix.c_str(), name2[j].c_str());
3319 TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
3320 if (hist1 != nullptr) {
3321 TH1D* hist = (TH1D*)(hist1->Clone());
3322 double ymin(0.90);
3323 sprintf(namep, "c_%s", name);
3324 TCanvas* pad = new TCanvas(namep, namep, 700, 500);
3325 pad->SetLogy();
3326 pad->SetRightMargin(0.10);
3327 pad->SetTopMargin(0.10);
3328 hist->GetXaxis()->SetTitleSize(0.04);
3329 if (k <= 1)
3330 hist->GetXaxis()->SetTitle(xtitl1.c_str());
3331 else
3332 hist->GetXaxis()->SetTitle(xtitl2.c_str());
3333 hist->GetYaxis()->SetTitle(ytitle.c_str());
3334 hist->GetYaxis()->SetLabelOffset(0.005);
3335 hist->GetYaxis()->SetTitleSize(0.04);
3336 hist->GetYaxis()->SetLabelSize(0.035);
3337 hist->GetYaxis()->SetTitleOffset(1.10);
3338 hist->SetMarkerStyle(20);
3339 hist->Draw();
3340 pad->Update();
3341 TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
3342 if (st1 != nullptr) {
3343 st1->SetY1NDC(0.78);
3344 st1->SetY2NDC(0.90);
3345 st1->SetX1NDC(0.65);
3346 st1->SetX2NDC(0.90);
3347 }
3348 pad->Update();
3349 double ymx(0.96), xmi(0.35), xmx(0.95);
3350 char txt[100];
3351 if (lumi > 0.1) {
3352 ymx = ymin - 0.005;
3353 xmi = 0.45;
3354 TPaveText* txt0 = new TPaveText(0.65, 0.91, 0.90, 0.96, "blNDC");
3355 txt0->SetFillColor(0);
3356 sprintf(txt, "%4.1f TeV %5.1f fb^{-1}", ener, lumi);
3357 txt0->AddText(txt);
3358 txt0->Draw("same");
3359 }
3360 double ymi = ymx - 0.05;
3361 TPaveText* txt1 = new TPaveText(xmi, ymi, xmx, ymx, "blNDC");
3362 txt1->SetFillColor(0);
3363 if (text == "") {
3364 if (k == 0)
3365 sprintf(txt, "%s", title1[j].c_str());
3366 else
3367 sprintf(txt, "%s", title2[j].c_str());
3368 } else {
3369 if (k == 0)
3370 sprintf(txt, "%s (%s)", title1[j].c_str(), text.c_str());
3371 else
3372 sprintf(txt, "%s (%s)", title2[j].c_str(), text.c_str());
3373 }
3374 txt1->AddText(txt);
3375 txt1->Draw("same");
3376 double xmax = (dataMC) ? 0.24 : 0.35;
3377 ymi = 0.91;
3378 ymx = ymi + 0.05;
3379 TPaveText* txt2 = new TPaveText(0.02, ymi, xmax, ymx, "blNDC");
3380 txt2->SetFillColor(0);
3381 if (dataMC)
3382 sprintf(txt, "CMS Preliminary");
3383 else
3384 sprintf(txt, "CMS Simulation Preliminary");
3385 txt2->AddText(txt);
3386 txt2->Draw("same");
3387 pad->Modified();
3388 pad->Update();
3389 if (save > 0) {
3390 sprintf(name, "%s.pdf", pad->GetName());
3391 pad->Print(name);
3392 } else if (save < 0) {
3393 sprintf(name, "%s.C", pad->GetName());
3394 pad->Print(name);
3395 }
3396 }
3397 }
3398 }
3399 }
3400
3401 void PlotMeanError(const std::string infilest, int reg = 3, bool resol = false, int save = 0, bool debug = false) {
3402 bool ok(false);
3403 const int ntypmx = 3;
3404 const int nregmx = 4;
3405 if (reg < 0 || reg >= nregmx)
3406 reg = nregmx - 1;
3407 int nEner(0), nType(0), nPts(0);
3408 std::vector<double> energy[ntypmx], denergy[ntypmx], value[ntypmx], dvalue[ntypmx];
3409
3410 std::ifstream fInput(infilest.c_str());
3411 if (!fInput.good()) {
3412 std::cout << "Cannot open file " << infilest << std::endl;
3413 } else {
3414 ok = true;
3415 fInput >> nEner >> nType >> nPts;
3416 int nmax = nEner * nType;
3417 int type, elow, ehigh;
3418 double v1[4], e1[4], v2[4], e2[4];
3419 for (int n = 0; n < nmax; ++n) {
3420 fInput >> type >> elow >> ehigh;
3421 fInput >> v1[0] >> e1[0] >> v1[1] >> e1[1] >> v1[2] >> e1[2] >> v1[3] >> e1[3];
3422 fInput >> v2[0] >> e2[0] >> v2[1] >> e2[1] >> v2[2] >> e2[2] >> v2[3] >> e2[3];
3423 double ener = 0.5 * (ehigh + elow);
3424 double dene = 0.5 * (ehigh - elow);
3425 energy[type].push_back(ener);
3426 denergy[type].push_back(dene);
3427 if (resol) {
3428 value[type].push_back(v2[reg]);
3429 dvalue[type].push_back(e2[reg]);
3430 } else {
3431 value[type].push_back(v1[reg]);
3432 dvalue[type].push_back(e1[reg]);
3433 }
3434 }
3435 fInput.close();
3436 std::cout << "Reads " << (nmax + 1) << " cards from " << infilest << " with measurements for " << nEner
3437 << " energies and " << nType << " types" << std::endl;
3438 if (debug) {
3439 for (int n = 0; n < nType; ++n) {
3440 std::cout << "Type " << n << " with " << energy[n].size() << " points\n";
3441 for (unsigned int k = 0; k < energy[n].size(); ++k)
3442 std::cout << " [" << k << "] " << energy[n][k] << " +- " << denergy[n][k] << " Value " << value[n][k]
3443 << " +- " << dvalue[n][k] << std::endl;
3444 }
3445 }
3446 }
3447
3448
3449 if (ok) {
3450 int mvsres = (resol) ? 1 : 0;
3451 std::string names[2] = {"Mean", "Resol"};
3452 std::string namet[nregmx] = {"Barrel", "Transition", "Endcap", "Combined"};
3453 char cname[100];
3454 sprintf(cname, "c_%s%s", names[mvsres].c_str(), namet[reg].c_str());
3455 int color[ntypmx] = {2, 4, 1};
3456 int mtype[ntypmx] = {20, 21, 22};
3457 double ymin[2] = {0.65, 0.10};
3458 double ymax[2] = {1.30, 0.50};
3459 gStyle->SetCanvasBorderMode(0);
3460 gStyle->SetCanvasColor(kWhite);
3461 gStyle->SetPadColor(kWhite);
3462 gStyle->SetFillColor(kWhite);
3463 gStyle->SetOptTitle(kFALSE);
3464 gStyle->SetPadBorderMode(0);
3465 gStyle->SetCanvasBorderMode(0);
3466 gStyle->SetOptStat(0);
3467 TCanvas* canvas = new TCanvas(cname, cname, 500, 500);
3468 canvas->SetTopMargin(0.05);
3469 canvas->SetBottomMargin(0.14);
3470 canvas->SetLeftMargin(0.15);
3471 canvas->SetRightMargin(0.10);
3472 TH1F* vFrame = canvas->DrawFrame(0.0, ymin[mvsres], 120.0, ymax[mvsres]);
3473 vFrame->GetXaxis()->SetRangeUser(0.0, 120.0);
3474 vFrame->GetYaxis()->SetRangeUser(ymin[mvsres], ymax[mvsres]);
3475 vFrame->GetXaxis()->SetLabelSize(0.04);
3476 vFrame->GetYaxis()->SetLabelSize(0.04);
3477 vFrame->GetXaxis()->SetTitleSize(0.04);
3478 vFrame->GetYaxis()->SetTitleSize(0.04);
3479 vFrame->GetXaxis()->SetTitleOffset(1.2);
3480 vFrame->GetYaxis()->SetTitleOffset(1.6);
3481 vFrame->GetXaxis()->SetLabelOffset(0.0);
3482 vFrame->GetXaxis()->SetTitle("p_{Beam} (GeV/c)");
3483 if (resol) {
3484 vFrame->GetYaxis()->SetTitle("Width(E_{HCAL}/(p-E_{ECAL}))");
3485 } else {
3486 vFrame->GetYaxis()->SetTitle("MPV(E_{HCAL}/(p-E_{ECAL}))");
3487 }
3488 TLegend* legend = new TLegend(0.70, 0.80, 0.90, 0.94);
3489 legend->SetFillColor(kWhite);
3490 std::string nameg[ntypmx] = {"MAHI", "M0", "M2"};
3491 for (int n = 0; n < nType; ++n) {
3492 unsigned int nmax0 = energy[n].size();
3493 double mom[nmax0], dmom[nmax0], mean[nmax0], dmean[nmax0];
3494 for (unsigned int k = 0; k < nmax0; ++k) {
3495 mom[k] = energy[n][k];
3496 dmom[k] = denergy[n][k];
3497 mean[k] = value[n][k];
3498 dmean[k] = dvalue[n][k];
3499 }
3500 TGraphErrors* graph = new TGraphErrors(nmax0, mom, mean, dmom, dmean);
3501 graph->SetMarkerStyle(mtype[n]);
3502 graph->SetMarkerColor(color[n]);
3503 graph->SetMarkerSize(1.4);
3504 graph->SetLineColor(color[n]);
3505 graph->SetLineWidth(2);
3506 sprintf(cname, "%s", nameg[n].c_str());
3507 legend->AddEntry(graph, cname, "lp");
3508 graph->Draw("P");
3509 }
3510 legend->Draw("same");
3511 std::string regions[nregmx] = {"20118B Barrel", "2018B Transition", "2018B Endcap", "2018B"};
3512 sprintf(cname, "%s", regions[reg].c_str());
3513 TPaveText* txt0 = new TPaveText(0.16, 0.90, 0.40, 0.94, "blNDC");
3514 txt0->SetFillColor(0);
3515 txt0->AddText(cname);
3516 txt0->Draw("same");
3517 TPaveText* txt1 = new TPaveText(0.15, 0.95, 0.40, 0.99, "blNDC");
3518 txt1->SetFillColor(0);
3519 sprintf(cname, "CMS Preliminary");
3520 txt1->AddText(cname);
3521 txt1->Draw("same");
3522 canvas->Update();
3523 if (save > 0) {
3524 sprintf(cname, "%s.pdf", canvas->GetName());
3525 canvas->Print(cname);
3526 } else if (save < 0) {
3527 sprintf(cname, "%s.C", canvas->GetName());
3528 canvas->Print(cname);
3529 }
3530 }
3531 }