File indexing completed on 2021-10-29 22:50:49
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 #include <TROOT.h>
0104 #include <TChain.h>
0105 #include <TFile.h>
0106 #include <TF1.h>
0107 #include <TH1D.h>
0108 #include <TH2F.h>
0109 #include <TProfile.h>
0110 #include <TFitResult.h>
0111 #include <TFitResultPtr.h>
0112 #include <TStyle.h>
0113 #include <TCanvas.h>
0114 #include <TLegend.h>
0115 #include <TPaveStats.h>
0116 #include <TPaveText.h>
0117
0118 #include <algorithm>
0119 #include <iomanip>
0120 #include <iostream>
0121 #include <fstream>
0122 #include <map>
0123 #include <vector>
0124 #include <string>
0125
0126 #include "CalibCorr.C"
0127
0128 namespace CalibPlots {
0129 static const int npbin = 4;
0130 static const int netabin = 4;
0131 static const int ndepth = 7;
0132 static const int ntitles = 5;
0133 static const int npbin0 = (npbin + 1);
0134 int getP(int k) {
0135 const int ipbin[npbin0] = {20, 30, 40, 60, 100};
0136 return ((k >= 0 && k < npbin0) ? ipbin[k] : 0);
0137 }
0138 double getMomentum(int k) { return (double)(getP(k)); }
0139 int getEta(int k) {
0140 int ietas[netabin] = {0, 13, 18, 23};
0141 return ((k >= 0 && k < netabin) ? ietas[k] : -1);
0142 }
0143 std::string getTitle(int k) {
0144 std::string titl[ntitles] = {
0145 "all tracks", "good quality tracks", "selected tracks", "isolated good tracks", "tracks having MIP in ECAL"};
0146 return ((k >= 0 && k < ntitles) ? titl[k] : "");
0147 }
0148 }
0149
0150 class CalibPlotProperties {
0151 public:
0152 TChain *fChain;
0153 Int_t fCurrent;
0154
0155
0156 Int_t t_Run;
0157 Int_t t_Event;
0158 Int_t t_DataType;
0159 Int_t t_ieta;
0160 Int_t t_iphi;
0161 Double_t t_EventWeight;
0162 Int_t t_nVtx;
0163 Int_t t_nTrk;
0164 Int_t t_goodPV;
0165 Double_t t_l1pt;
0166 Double_t t_l1eta;
0167 Double_t t_l1phi;
0168 Double_t t_l3pt;
0169 Double_t t_l3eta;
0170 Double_t t_l3phi;
0171 Double_t t_p;
0172 Double_t t_pt;
0173 Double_t t_phi;
0174 Double_t t_mindR1;
0175 Double_t t_mindR2;
0176 Double_t t_eMipDR;
0177 Double_t t_eHcal;
0178 Double_t t_eHcal10;
0179 Double_t t_eHcal30;
0180 Double_t t_hmaxNearP;
0181 Double_t t_rhoh;
0182 Bool_t t_selectTk;
0183 Bool_t t_qltyFlag;
0184 Bool_t t_qltyMissFlag;
0185 Bool_t t_qltyPVFlag;
0186 Double_t t_gentrackP;
0187 std::vector<unsigned int> *t_DetIds;
0188 std::vector<double> *t_HitEnergies;
0189 std::vector<bool> *t_trgbits;
0190 std::vector<unsigned int> *t_DetIds1;
0191 std::vector<unsigned int> *t_DetIds3;
0192 std::vector<double> *t_HitEnergies1;
0193 std::vector<double> *t_HitEnergies3;
0194
0195
0196 TBranch *b_t_Run;
0197 TBranch *b_t_Event;
0198 TBranch *b_t_DataType;
0199 TBranch *b_t_ieta;
0200 TBranch *b_t_iphi;
0201 TBranch *b_t_EventWeight;
0202 TBranch *b_t_nVtx;
0203 TBranch *b_t_nTrk;
0204 TBranch *b_t_goodPV;
0205 TBranch *b_t_l1pt;
0206 TBranch *b_t_l1eta;
0207 TBranch *b_t_l1phi;
0208 TBranch *b_t_l3pt;
0209 TBranch *b_t_l3eta;
0210 TBranch *b_t_l3phi;
0211 TBranch *b_t_p;
0212 TBranch *b_t_pt;
0213 TBranch *b_t_phi;
0214 TBranch *b_t_mindR1;
0215 TBranch *b_t_mindR2;
0216 TBranch *b_t_eMipDR;
0217 TBranch *b_t_eHcal;
0218 TBranch *b_t_eHcal10;
0219 TBranch *b_t_eHcal30;
0220 TBranch *b_t_hmaxNearP;
0221 TBranch *b_t_rhoh;
0222 TBranch *b_t_selectTk;
0223 TBranch *b_t_qltyFlag;
0224 TBranch *b_t_qltyMissFlag;
0225 TBranch *b_t_qltyPVFlag;
0226 TBranch *b_t_gentrackP;
0227 TBranch *b_t_DetIds;
0228 TBranch *b_t_HitEnergies;
0229 TBranch *b_t_trgbits;
0230 TBranch *b_t_DetIds1;
0231 TBranch *b_t_DetIds3;
0232 TBranch *b_t_HitEnergies1;
0233 TBranch *b_t_HitEnergies3;
0234
0235 CalibPlotProperties(const char *fname,
0236 const std::string &dirname,
0237 const char *dupFileName,
0238 const std::string &prefix = "",
0239 const char *corrFileName = "",
0240 const char *rcorFileName = "",
0241 int puCorr = -1,
0242 int flag = 101111,
0243 bool dataMC = true,
0244 int truncateFlag = 0,
0245 bool useGen = false,
0246 double scale = 1.0,
0247 int useScale = 0,
0248 int etalo = 0,
0249 int etahi = 30,
0250 int runlo = 0,
0251 int runhi = 99999999,
0252 int phimin = 1,
0253 int phimax = 72,
0254 int zside = 1,
0255 int nvxlo = 0,
0256 int nvxhi = 1000,
0257 int rbx = 0,
0258 bool exclude = false,
0259 bool etamax = false);
0260 virtual ~CalibPlotProperties();
0261 virtual Int_t Cut(Long64_t entry);
0262 virtual Int_t GetEntry(Long64_t entry);
0263 virtual Long64_t LoadTree(Long64_t entry);
0264 virtual void Init(TChain *, const char *);
0265 virtual void Loop(Long64_t nentries = -1);
0266 virtual Bool_t Notify();
0267 virtual void Show(Long64_t entry = -1);
0268 bool goodTrack(double &eHcal, double &cut, bool debug);
0269 bool selectPhi(bool debug);
0270 void savePlot(const std::string &theName, bool append, bool all = false, bool debug = false);
0271 void correctEnergy(double &ener);
0272
0273 private:
0274 static const int kp50 = 2;
0275 CalibCorrFactor *corrFactor_;
0276 CalibCorr *cFactor_;
0277 CalibSelectRBX *cSelect_;
0278 const std::string fname_, dirnm_, prefix_, outFileName_;
0279 const int corrPU_, flag_;
0280 const bool dataMC_, useGen_;
0281 const int truncateFlag_;
0282 const int etalo_, etahi_;
0283 int runlo_, runhi_;
0284 const int phimin_, phimax_, zside_, nvxlo_, nvxhi_, rbx_;
0285 bool exclude_, corrE_, cutL1T_;
0286 bool includeRun_, getHist_;
0287 int flexibleSelect_;
0288 bool plotBasic_, plotEnergy_, plotHists_;
0289 double log2by18_;
0290 std::ofstream fileout_;
0291 std::vector<Long64_t> entries_;
0292 std::vector<std::pair<int, int> > events_;
0293 TH1D *h_p[CalibPlots::ntitles];
0294 TH1D *h_eta[CalibPlots::ntitles], *h_nvtx;
0295 std::vector<TH1D *> h_eta0, h_eta1, h_eta2, h_eta3, h_eta4;
0296 std::vector<TH1D *> h_dL1, h_vtx;
0297 std::vector<TH1D *> h_etaEH[CalibPlots::npbin0];
0298 std::vector<TH1D *> h_etaEp[CalibPlots::npbin0];
0299 std::vector<TH1D *> h_etaEE[CalibPlots::npbin0];
0300 std::vector<TH1D *> h_etaEE0[CalibPlots::npbin0];
0301 std::vector<TH1D *> h_mom, h_eEcal, h_eHcal;
0302 std::vector<TH1F *> h_bvlist, h_bvlist2, h_evlist, h_evlist2;
0303 std::vector<TH1F *> h_bvlist3, h_evlist3;
0304 TH2F *h_etaE;
0305 };
0306
0307 CalibPlotProperties::CalibPlotProperties(const char *fname,
0308 const std::string &dirnm,
0309 const char *dupFileName,
0310 const std::string &prefix,
0311 const char *corrFileName,
0312 const char *rcorFileName,
0313 int puCorr,
0314 int flag,
0315 bool dataMC,
0316 int truncate,
0317 bool useGen,
0318 double scl,
0319 int useScale,
0320 int etalo,
0321 int etahi,
0322 int runlo,
0323 int runhi,
0324 int phimin,
0325 int phimax,
0326 int zside,
0327 int nvxlo,
0328 int nvxhi,
0329 int rbx,
0330 bool exc,
0331 bool etam)
0332 : corrFactor_(nullptr),
0333 cFactor_(nullptr),
0334 cSelect_(nullptr),
0335 fname_(fname),
0336 dirnm_(dirnm),
0337 prefix_(prefix),
0338 corrPU_(puCorr),
0339 flag_(flag),
0340 dataMC_(dataMC),
0341 useGen_(useGen),
0342 truncateFlag_(truncate),
0343 etalo_(etalo),
0344 etahi_(etahi),
0345 runlo_(runlo),
0346 runhi_(runhi),
0347 phimin_(phimin),
0348 phimax_(phimax),
0349 zside_(zside),
0350 nvxlo_(nvxlo),
0351 nvxhi_(nvxhi),
0352 rbx_(rbx),
0353 exclude_(exc),
0354 includeRun_(true) {
0355
0356
0357
0358 flexibleSelect_ = ((flag_ / 1) % 10);
0359 plotBasic_ = (((flag_ / 10) % 10) > 0);
0360 plotEnergy_ = (((flag_ / 10) % 10) > 0);
0361 int oneplace = ((flag_ / 1000) % 10);
0362 cutL1T_ = (oneplace % 2);
0363 bool marina = ((oneplace / 2) % 2);
0364 bool ifDepth = (((flag_ / 10000) % 10) > 0);
0365 plotHists_ = (((flag_ / 100000) % 10) > 0);
0366 log2by18_ = std::log(2.5) / 18.0;
0367 if (runlo_ < 0 || runhi_ < 0) {
0368 runlo_ = std::abs(runlo_);
0369 runhi_ = std::abs(runhi_);
0370 includeRun_ = false;
0371 }
0372 char treeName[400];
0373 sprintf(treeName, "%s/CalibTree", dirnm.c_str());
0374 TChain *chain = new TChain(treeName);
0375 std::cout << "Create a chain for " << treeName << " from " << fname << " flags " << flexibleSelect_ << "|"
0376 << plotBasic_ << "|"
0377 << "|" << plotEnergy_ << "|" << plotHists_ << "|" << corrPU_ << " cons " << log2by18_ << " eta range "
0378 << etalo_ << ":" << etahi_ << " run range " << runlo_ << ":" << runhi_ << " (inclusion flag " << includeRun_
0379 << ") Vertex Range " << nvxlo_ << ":" << nvxhi_ << std::endl;
0380 corrFactor_ = new CalibCorrFactor(corrFileName, useScale, scl, etam, marina, false);
0381 if (!fillChain(chain, fname)) {
0382 std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
0383 } else {
0384 std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
0385 Init(chain, dupFileName);
0386 if (std::string(rcorFileName) != "")
0387 cFactor_ = new CalibCorr(rcorFileName, ifDepth, false);
0388 if (rbx != 0)
0389 cSelect_ = new CalibSelectRBX(rbx, false);
0390 }
0391 }
0392
0393 CalibPlotProperties::~CalibPlotProperties() {
0394 delete corrFactor_;
0395 delete cFactor_;
0396 delete cSelect_;
0397 if (!fChain)
0398 return;
0399 delete fChain->GetCurrentFile();
0400 }
0401
0402 Int_t CalibPlotProperties::GetEntry(Long64_t entry) {
0403
0404 if (!fChain)
0405 return 0;
0406 return fChain->GetEntry(entry);
0407 }
0408
0409 Long64_t CalibPlotProperties::LoadTree(Long64_t entry) {
0410
0411 if (!fChain)
0412 return -5;
0413 Long64_t centry = fChain->LoadTree(entry);
0414 if (centry < 0)
0415 return centry;
0416 if (!fChain->InheritsFrom(TChain::Class()))
0417 return centry;
0418 TChain *chain = (TChain *)fChain;
0419 if (chain->GetTreeNumber() != fCurrent) {
0420 fCurrent = chain->GetTreeNumber();
0421 Notify();
0422 }
0423 return centry;
0424 }
0425
0426 void CalibPlotProperties::Init(TChain *tree, const char *dupFileName) {
0427
0428
0429
0430
0431
0432
0433
0434
0435
0436 t_DetIds = 0;
0437 t_DetIds1 = 0;
0438 t_DetIds3 = 0;
0439 t_HitEnergies = 0;
0440 t_HitEnergies1 = 0;
0441 t_HitEnergies3 = 0;
0442 t_trgbits = 0;
0443
0444 fChain = tree;
0445 fCurrent = -1;
0446 if (!tree)
0447 return;
0448 fChain->SetMakeClass(1);
0449
0450 fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
0451 fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
0452 fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
0453 fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
0454 fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
0455 fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
0456 fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
0457 fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
0458 fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
0459 fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
0460 fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
0461 fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
0462 fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
0463 fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
0464 fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
0465 fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
0466 fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
0467 fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
0468 fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
0469 fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
0470 fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
0471 fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
0472 fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
0473 fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
0474 fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
0475 fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
0476 fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
0477 fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
0478 fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
0479 fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
0480 fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
0481 fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
0482 fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
0483 fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
0484 fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
0485 fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
0486 fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
0487 fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
0488 Notify();
0489
0490 if (std::string(dupFileName) != "") {
0491 ifstream infil1(dupFileName);
0492 if (!infil1.is_open()) {
0493 std::cout << "Cannot open duplicate file " << dupFileName << std::endl;
0494 } else {
0495 while (1) {
0496 Long64_t jentry;
0497 infil1 >> jentry;
0498 if (!infil1.good())
0499 break;
0500 entries_.push_back(jentry);
0501 }
0502 infil1.close();
0503 std::cout << "Reads a list of " << entries_.size() << " events from " << dupFileName << std::endl;
0504 }
0505 }
0506
0507 char name[20], title[200];
0508 unsigned int kk(0);
0509
0510 if (plotBasic_) {
0511 std::cout << "Book Basic Histos" << std::endl;
0512 for (int k = 0; k < CalibPlots::ntitles; ++k) {
0513 sprintf(name, "%sp%d", prefix_.c_str(), k);
0514 sprintf(title, "Momentum for %s", CalibPlots::getTitle(k).c_str());
0515 h_p[k] = new TH1D(name, title, 150, 0.0, 150.0);
0516 sprintf(name, "%seta%d", prefix_.c_str(), k);
0517 sprintf(title, "#eta for %s", CalibPlots::getTitle(k).c_str());
0518 h_eta[k] = new TH1D(name, title, 60, -30.0, 30.0);
0519 }
0520 for (int k = 0; k < CalibPlots::npbin; ++k) {
0521 sprintf(name, "%seta0%d", prefix_.c_str(), k);
0522 sprintf(title,
0523 "#eta for %s (p = %d:%d GeV)",
0524 CalibPlots::getTitle(0).c_str(),
0525 CalibPlots::getP(k),
0526 CalibPlots::getP(k + 1));
0527 h_eta0.push_back(new TH1D(name, title, 60, -30.0, 30.0));
0528 kk = h_eta0.size() - 1;
0529 h_eta0[kk]->Sumw2();
0530 sprintf(name, "%seta1%d", prefix_.c_str(), k);
0531 sprintf(title,
0532 "#eta for %s (p = %d:%d GeV)",
0533 CalibPlots::getTitle(1).c_str(),
0534 CalibPlots::getP(k),
0535 CalibPlots::getP(k + 1));
0536 h_eta1.push_back(new TH1D(name, title, 60, -30.0, 30.0));
0537 kk = h_eta1.size() - 1;
0538 h_eta1[kk]->Sumw2();
0539 sprintf(name, "%seta2%d", prefix_.c_str(), k);
0540 sprintf(title,
0541 "#eta for %s (p = %d:%d GeV)",
0542 CalibPlots::getTitle(2).c_str(),
0543 CalibPlots::getP(k),
0544 CalibPlots::getP(k + 1));
0545 h_eta2.push_back(new TH1D(name, title, 60, -30.0, 30.0));
0546 kk = h_eta2.size() - 1;
0547 h_eta2[kk]->Sumw2();
0548 sprintf(name, "%seta3%d", prefix_.c_str(), k);
0549 sprintf(title,
0550 "#eta for %s (p = %d:%d GeV)",
0551 CalibPlots::getTitle(3).c_str(),
0552 CalibPlots::getP(k),
0553 CalibPlots::getP(k + 1));
0554 h_eta3.push_back(new TH1D(name, title, 60, -30.0, 30.0));
0555 kk = h_eta3.size() - 1;
0556 h_eta3[kk]->Sumw2();
0557 sprintf(name, "%seta4%d", prefix_.c_str(), k);
0558 sprintf(title,
0559 "#eta for %s (p = %d:%d GeV)",
0560 CalibPlots::getTitle(4).c_str(),
0561 CalibPlots::getP(k),
0562 CalibPlots::getP(k + 1));
0563 h_eta4.push_back(new TH1D(name, title, 60, -30.0, 30.0));
0564 kk = h_eta4.size() - 1;
0565 h_eta4[kk]->Sumw2();
0566 sprintf(name, "%sdl1%d", prefix_.c_str(), k);
0567 sprintf(title, "Distance from L1 (p = %d:%d GeV)", CalibPlots::getP(k), CalibPlots::getP(k + 1));
0568 h_dL1.push_back(new TH1D(name, title, 160, 0.0, 8.0));
0569 kk = h_dL1.size() - 1;
0570 h_dL1[kk]->Sumw2();
0571 sprintf(name, "%svtx%d", prefix_.c_str(), k);
0572 sprintf(title, "N_{Vertex} (p = %d:%d GeV)", CalibPlots::getP(k), CalibPlots::getP(k + 1));
0573 h_vtx.push_back(new TH1D(name, title, 100, 0.0, 100.0));
0574 kk = h_vtx.size() - 1;
0575 h_vtx[kk]->Sumw2();
0576 }
0577 }
0578
0579 if (plotEnergy_) {
0580 std::cout << "Make plots for good tracks" << std::endl;
0581 for (int k = 0; k < CalibPlots::npbin; ++k) {
0582 for (int j = etalo_; j <= etahi_ + 1; ++j) {
0583 sprintf(name, "%senergyH%d%d", prefix_.c_str(), k, j);
0584 if (j > etahi_)
0585 sprintf(title,
0586 "HCAL energy for %s (p = %d:%d GeV |#eta| = %d:%d)",
0587 CalibPlots::getTitle(3).c_str(),
0588 CalibPlots::getP(k),
0589 CalibPlots::getP(k + 1),
0590 etalo_,
0591 etahi_);
0592 else
0593 sprintf(title,
0594 "HCAL energy for %s (p = %d:%d GeV |#eta| = %d)",
0595 CalibPlots::getTitle(3).c_str(),
0596 CalibPlots::getP(k),
0597 CalibPlots::getP(k + 1),
0598 j);
0599 h_etaEH[k].push_back(new TH1D(name, title, 200, 0, 200));
0600 kk = h_etaEH[k].size() - 1;
0601 h_etaEH[k][kk]->Sumw2();
0602 sprintf(name, "%senergyP%d%d", prefix_.c_str(), k, j);
0603 if (j > etahi_)
0604 sprintf(title,
0605 "momentum for %s (p = %d:%d GeV |#eta| = %d:%d)",
0606 CalibPlots::getTitle(3).c_str(),
0607 CalibPlots::getP(k),
0608 CalibPlots::getP(k + 1),
0609 etalo_,
0610 etahi_);
0611 else
0612 sprintf(title,
0613 "momentum for %s (p = %d:%d GeV |#eta| = %d)",
0614 CalibPlots::getTitle(3).c_str(),
0615 CalibPlots::getP(k),
0616 CalibPlots::getP(k + 1),
0617 j);
0618 h_etaEp[k].push_back(new TH1D(name, title, 100, 0, 100));
0619 kk = h_etaEp[k].size() - 1;
0620 h_etaEp[k][kk]->Sumw2();
0621 sprintf(name, "%senergyE%d%d", prefix_.c_str(), k, j);
0622 if (j > etahi_)
0623 sprintf(title,
0624 "ECAL energy for %s (p = %d:%d GeV |#eta| = %d:%d)",
0625 CalibPlots::getTitle(3).c_str(),
0626 CalibPlots::getP(k),
0627 CalibPlots::getP(k + 1),
0628 etalo_,
0629 etahi_);
0630 else
0631 sprintf(title,
0632 "ECAL energy for %s (p = %d:%d GeV |#eta| = %d)",
0633 CalibPlots::getTitle(3).c_str(),
0634 CalibPlots::getP(k),
0635 CalibPlots::getP(k + 1),
0636 j);
0637 h_etaEE[k].push_back(new TH1D(name, title, 100, 0, 10));
0638 kk = h_etaEE[k].size() - 1;
0639 h_etaEE[k][kk]->Sumw2();
0640 sprintf(name, "%senergyER%d%d", prefix_.c_str(), k, j);
0641 h_etaEE0[k].push_back(new TH1D(name, title, 100, 0, 1));
0642 kk = h_etaEE0[k].size() - 1;
0643 h_etaEE0[k][kk]->Sumw2();
0644 }
0645 }
0646
0647 for (int j = 0; j < CalibPlots::netabin; ++j) {
0648 sprintf(name, "%senergyH%d", prefix_.c_str(), j);
0649 if (j == 0)
0650 sprintf(title, "HCAL energy for %s (All)", CalibPlots::getTitle(3).c_str());
0651 else
0652 sprintf(title,
0653 "HCAL energy for %s (|#eta| = %d:%d)",
0654 CalibPlots::getTitle(3).c_str(),
0655 CalibPlots::getEta(j - 1),
0656 CalibPlots::getEta(j));
0657 h_eHcal.push_back(new TH1D(name, title, 200, 0, 200));
0658 kk = h_eHcal.size() - 1;
0659 h_eHcal[kk]->Sumw2();
0660 sprintf(name, "%senergyP%d", prefix_.c_str(), j);
0661 if (j == 0)
0662 sprintf(title, "Track momentum for %s (All)", CalibPlots::getTitle(3).c_str());
0663 else
0664 sprintf(title,
0665 "Track momentum for %s (|#eta| = %d:%d)",
0666 CalibPlots::getTitle(3).c_str(),
0667 CalibPlots::getEta(j - 1),
0668 CalibPlots::getEta(j));
0669 h_mom.push_back(new TH1D(name, title, 100, 0, 100));
0670 kk = h_mom.size() - 1;
0671 h_mom[kk]->Sumw2();
0672 sprintf(name, "%senergyE%d", prefix_.c_str(), j);
0673 if (j == 0)
0674 sprintf(title, "ECAL energy for %s (All)", CalibPlots::getTitle(3).c_str());
0675 else
0676 sprintf(title,
0677 "ECAL energy for %s (|#eta| = %d:%d)",
0678 CalibPlots::getTitle(3).c_str(),
0679 CalibPlots::getEta(j - 1),
0680 CalibPlots::getEta(j));
0681 h_eEcal.push_back(new TH1D(name, title, 100, 0, 10));
0682 kk = h_eEcal.size() - 1;
0683 h_eEcal[kk]->Sumw2();
0684 }
0685 }
0686
0687 if (plotHists_) {
0688 h_nvtx = new TH1D("hnvtx", "Number of vertices", 10, 0, 100);
0689 h_nvtx->Sumw2();
0690 for (int i = 0; i < CalibPlots::ndepth; i++) {
0691 sprintf(name, "b_edepth%d", i);
0692 sprintf(title, "Total RecHit energy in depth %d (Barrel)", i + 1);
0693 h_bvlist.push_back(new TH1F(name, title, 1000, 0, 100));
0694 h_bvlist[i]->Sumw2();
0695 sprintf(name, "b_recedepth%d", i);
0696 sprintf(title, "RecHit energy in depth %d (Barrel)", i + 1);
0697 h_bvlist2.push_back(new TH1F(name, title, 1000, 0, 100));
0698 h_bvlist2[i]->Sumw2();
0699 sprintf(name, "b_nrecdepth%d", i);
0700 sprintf(title, "#RecHits in depth %d (Barrel)", i + 1);
0701 h_bvlist3.push_back(new TH1F(name, title, 1000, 0, 100));
0702 h_bvlist3[i]->Sumw2();
0703 sprintf(name, "e_edepth%d", i);
0704 sprintf(title, "Total RecHit energy in depth %d (Endcap)", i + 1);
0705 h_evlist.push_back(new TH1F(name, title, 1000, 0, 100));
0706 h_evlist[i]->Sumw2();
0707 sprintf(name, "e_recedepth%d", i);
0708 sprintf(title, "RecHit energy in depth %d (Endcap)", i + 1);
0709 h_evlist2.push_back(new TH1F(name, title, 1000, 0, 100));
0710 h_evlist2[i]->Sumw2();
0711 sprintf(name, "e_nrecdepth%d", i);
0712 sprintf(title, "#RecHits in depth %d (Endcap)", i + 1);
0713 h_evlist3.push_back(new TH1F(name, title, 1000, 0, 100));
0714 h_evlist3[i]->Sumw2();
0715 }
0716 h_etaE = new TH2F("heta", "", 60, -30, 30, 100, 0, 100);
0717 }
0718 }
0719
0720 Bool_t CalibPlotProperties::Notify() {
0721
0722
0723
0724
0725
0726
0727 return kTRUE;
0728 }
0729
0730 void CalibPlotProperties::Show(Long64_t entry) {
0731
0732
0733 if (!fChain)
0734 return;
0735 fChain->Show(entry);
0736 }
0737
0738 Int_t CalibPlotProperties::Cut(Long64_t) {
0739
0740
0741
0742 return 1;
0743 }
0744
0745 void CalibPlotProperties::Loop(Long64_t nentries) {
0746
0747
0748
0749
0750
0751
0752
0753
0754
0755
0756
0757
0758
0759
0760
0761
0762
0763
0764
0765
0766
0767
0768
0769 const bool debug(false);
0770
0771 if (fChain == 0)
0772 return;
0773
0774
0775 if (nentries < 0)
0776 nentries = fChain->GetEntriesFast();
0777 std::cout << "Total entries " << nentries << ":" << fChain->GetEntriesFast() << std::endl;
0778 Long64_t nbytes(0), nb(0);
0779 unsigned int duplicate(0), good(0), kount(0);
0780 double sel(0), selHB(0), selHE(0);
0781 for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0782 Long64_t ientry = LoadTree(jentry);
0783 if (ientry < 0)
0784 break;
0785 nb = fChain->GetEntry(jentry);
0786 nbytes += nb;
0787 if (jentry % 1000000 == 0)
0788 std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
0789 bool select = (std::find(entries_.begin(), entries_.end(), jentry) == entries_.end());
0790 if (!select) {
0791 ++duplicate;
0792 if (debug)
0793 std::cout << "Duplicate event " << t_Run << " " << t_Event << " " << t_p << std::endl;
0794 continue;
0795 }
0796 bool selRun = (includeRun_ ? ((t_Run >= runlo_) && (t_Run <= runhi_)) : ((t_Run < runlo_) || (t_Run > runhi_)));
0797 select =
0798 (selRun && (fabs(t_ieta) >= etalo_) && (fabs(t_ieta) <= etahi_) && (t_nVtx >= nvxlo_) && (t_nVtx <= nvxhi_));
0799 if (!select) {
0800 if (debug)
0801 std::cout << "Run # " << t_Run << " out of range of " << runlo_ << ":" << runhi_ << " or ieta " << t_ieta
0802 << " (" << etalo_ << ":" << etahi_ << ") or nvtx " << t_nVtx << " (" << nvxlo_ << ":" << nvxhi_
0803 << ") out of range" << std::endl;
0804 continue;
0805 }
0806 if (cSelect_ != nullptr) {
0807 if (exclude_) {
0808 if (cSelect_->isItRBX(t_DetIds))
0809 continue;
0810 } else {
0811 if (!(cSelect_->isItRBX(t_ieta, t_iphi)))
0812 continue;
0813 }
0814 }
0815 select = (!cutL1T_ || (t_mindR1 >= 0.5));
0816 if (!select) {
0817 if (debug)
0818 std::cout << "Reject Run # " << t_Run << " Event # " << t_Event << " too close to L1 trigger " << t_mindR1
0819 << std::endl;
0820 continue;
0821 }
0822 select = ((events_.size() == 0) ||
0823 (std::find(events_.begin(), events_.end(), std::pair<int, int>(t_Run, t_Event)) != events_.end()));
0824 if (!select) {
0825 if (debug)
0826 std::cout << "Reject Run # " << t_Run << " Event # " << t_Event << " not in the selection list" << std::endl;
0827 continue;
0828 }
0829
0830
0831 int kp = CalibPlots::npbin0;
0832 double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
0833 for (int k = 1; k < CalibPlots::npbin0; ++k) {
0834 if (pmom >= CalibPlots::getMomentum(k - 1) && pmom < CalibPlots::getMomentum(k)) {
0835 kp = k - 1;
0836 break;
0837 }
0838 }
0839 int jp1 = (((std::abs(t_ieta) >= etalo_) && (std::abs(t_ieta) <= etahi_)) ? (std::abs(t_ieta) - etalo_) : -1);
0840 int jp2 = (etahi_ - etalo_ + 1);
0841 int je1 = CalibPlots::netabin;
0842 for (int j = 1; j < CalibPlots::netabin; ++j) {
0843 if (std::abs(t_ieta) > CalibPlots::getEta(j - 1) && std::abs(t_ieta) <= CalibPlots::getEta(j)) {
0844 je1 = j;
0845 break;
0846 }
0847 }
0848 int je2 = (je1 != CalibPlots::netabin) ? 0 : -1;
0849 if (debug)
0850 std::cout << "Bin " << kp << ":" << je1 << ":" << je2 << ":" << jp1 << ":" << jp2 << std::endl;
0851 double cut = (pmom > 20) ? ((flexibleSelect_ == 0) ? 2.0 : 10.0) : 0.0;
0852 double rcut(-1000.0);
0853
0854
0855 double eHcal(t_eHcal);
0856 if (corrFactor_->doCorr()) {
0857 eHcal = 0;
0858 for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
0859
0860 unsigned int id = truncateId((*t_DetIds)[k], truncateFlag_, false);
0861 double cfac = corrFactor_->getCorr(id);
0862 if (cFactor_ != 0)
0863 cfac *= cFactor_->getCorr(t_Run, (*t_DetIds)[k]);
0864 eHcal += (cfac * ((*t_HitEnergies)[k]));
0865 if (debug) {
0866 int subdet, zside, ieta, iphi, depth;
0867 unpackDetId(id, subdet, zside, ieta, iphi, depth);
0868 std::cout << zside << ":" << ieta << ":" << depth << " Corr " << cfac << " " << (*t_HitEnergies)[k] << " Out "
0869 << eHcal << std::endl;
0870 }
0871 }
0872 }
0873 bool goodTk = goodTrack(eHcal, cut, debug);
0874 bool selPhi = selectPhi(debug);
0875 double rat = (pmom > 0) ? (eHcal / (pmom - t_eMipDR)) : 1.0;
0876 if (debug)
0877 std::cout << "Entry " << jentry << " p|eHcal|ratio " << pmom << "|" << t_eHcal << "|" << eHcal << "|" << rat
0878 << "|" << kp << " Cuts " << t_qltyFlag << "|" << t_selectTk << "|" << (t_hmaxNearP < cut) << "|"
0879 << (t_eMipDR < 1.0) << "|" << goodTk << "|" << (rat > rcut) << " Select Phi " << selPhi << std::endl;
0880 if (plotBasic_) {
0881 h_p[0]->Fill(pmom, t_EventWeight);
0882 h_eta[0]->Fill(t_ieta, t_EventWeight);
0883 if (kp < CalibPlots::npbin0)
0884 h_eta0[kp]->Fill(t_ieta, t_EventWeight);
0885 if (t_qltyFlag) {
0886 h_p[1]->Fill(pmom, t_EventWeight);
0887 h_eta[1]->Fill(t_ieta, t_EventWeight);
0888 if (kp < CalibPlots::npbin0)
0889 h_eta1[kp]->Fill(t_ieta, t_EventWeight);
0890 if (t_selectTk) {
0891 h_p[2]->Fill(pmom, t_EventWeight);
0892 h_eta[2]->Fill(t_ieta, t_EventWeight);
0893 if (kp < CalibPlots::npbin0)
0894 h_eta2[kp]->Fill(t_ieta, t_EventWeight);
0895 if (t_hmaxNearP < cut) {
0896 h_p[3]->Fill(pmom, t_EventWeight);
0897 h_eta[3]->Fill(t_ieta, t_EventWeight);
0898 if (kp < CalibPlots::npbin0)
0899 h_eta3[kp]->Fill(t_ieta, t_EventWeight);
0900 if (t_eMipDR < 1.0) {
0901 h_p[4]->Fill(pmom, t_EventWeight);
0902 h_eta[4]->Fill(t_ieta, t_EventWeight);
0903 if (kp < CalibPlots::npbin0) {
0904 h_eta4[kp]->Fill(t_ieta, t_EventWeight);
0905 h_dL1[kp]->Fill(t_mindR1, t_EventWeight);
0906 h_vtx[kp]->Fill(t_goodPV, t_EventWeight);
0907 }
0908 }
0909 }
0910 }
0911 }
0912 }
0913
0914 if (goodTk && kp != CalibPlots::npbin0 && selPhi) {
0915 if (rat > rcut) {
0916 if (plotEnergy_) {
0917 if (jp1 >= 0) {
0918 h_etaEH[kp][jp1]->Fill(eHcal, t_EventWeight);
0919 h_etaEH[kp][jp2]->Fill(eHcal, t_EventWeight);
0920 h_etaEp[kp][jp1]->Fill(pmom, t_EventWeight);
0921 h_etaEp[kp][jp2]->Fill(pmom, t_EventWeight);
0922 h_etaEE[kp][jp1]->Fill(t_eMipDR, t_EventWeight);
0923 h_etaEE[kp][jp2]->Fill(t_eMipDR, t_EventWeight);
0924 h_etaEE0[kp][jp1]->Fill(t_eMipDR, t_EventWeight);
0925 h_etaEE0[kp][jp2]->Fill(t_eMipDR, t_EventWeight);
0926 }
0927 if (kp == kp50) {
0928 if (je1 != CalibPlots::netabin) {
0929 h_eHcal[je1]->Fill(eHcal, t_EventWeight);
0930 h_eHcal[je2]->Fill(eHcal, t_EventWeight);
0931 h_mom[je1]->Fill(pmom, t_EventWeight);
0932 h_mom[je2]->Fill(pmom, t_EventWeight);
0933 h_eEcal[je1]->Fill(t_eMipDR, t_EventWeight);
0934 h_eEcal[je2]->Fill(t_eMipDR, t_EventWeight);
0935 }
0936 }
0937 }
0938
0939 if (plotHists_) {
0940 h_nvtx->Fill(t_nVtx);
0941 if ((std::fabs(rat - 1) < 0.15) && (kp == kp50) && ((std::abs(t_ieta) < 15) || (std::abs(t_ieta) > 17))) {
0942 float weight = (dataMC_ ? t_EventWeight : t_EventWeight * puweight(t_nVtx));
0943 h_etaE->Fill(t_ieta, eHcal, weight);
0944 sel += weight;
0945 std::vector<float> bv(7, 0.0f), ev(7, 0.0f);
0946 std::vector<int> bnrec(7, 0), enrec(7, 0);
0947 double eb(0), ee(0);
0948 for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
0949 unsigned int id = truncateId((*t_DetIds)[k], truncateFlag_, false);
0950 double cfac = corrFactor_->getCorr(id);
0951 if (cFactor_ != 0)
0952 cfac *= cFactor_->getCorr(t_Run, (*t_DetIds)[k]);
0953 double ener = cfac * (*t_HitEnergies)[k];
0954 if (corrPU_)
0955 correctEnergy(ener);
0956 unsigned int idx = (unsigned int)((*t_DetIds)[k]);
0957 int subdet, zside, ieta, iphi, depth;
0958 unpackDetId(idx, subdet, zside, ieta, iphi, depth);
0959 if (depth > 0 && depth <= CalibPlots::ndepth) {
0960 if (subdet == 1) {
0961 eb += ener;
0962 bv[depth - 1] += ener;
0963 h_bvlist2[depth - 1]->Fill(ener, weight);
0964 ++bnrec[depth - 1];
0965 } else if (subdet == 2) {
0966 ee += ener;
0967 ev[depth - 1] += ener;
0968 h_evlist2[depth - 1]->Fill(ener, weight);
0969 ++enrec[depth - 1];
0970 }
0971 }
0972 }
0973 bool barrel = (eb > ee);
0974 if (barrel)
0975 selHB += weight;
0976 else
0977 selHE += weight;
0978 for (int i = 0; i < CalibPlots::ndepth; i++) {
0979 if (barrel) {
0980 h_bvlist[i]->Fill(bv[i], weight);
0981 h_bvlist3[i]->Fill((bnrec[i] + 0.001), weight);
0982 } else {
0983 h_evlist[i]->Fill(ev[i], weight);
0984 h_evlist3[i]->Fill((enrec[i] + 0.001), weight);
0985 }
0986 }
0987 }
0988 }
0989 }
0990 good++;
0991 }
0992 ++kount;
0993 }
0994 std::cout << "Finds " << duplicate << " Duplicate events out of " << kount << " events in this file and " << good
0995 << " selected events" << std::endl;
0996 if (plotHists_)
0997 std::cout << "Number of weighted selected events " << sel << " HB " << selHB << " HE " << selHE << std::endl;
0998 }
0999
1000 bool CalibPlotProperties::goodTrack(double &eHcal, double &cuti, bool debug) {
1001 bool select(true);
1002 double cut(cuti);
1003 if (debug) {
1004 std::cout << "goodTrack input " << eHcal << ":" << cut;
1005 }
1006 if (flexibleSelect_ > 1) {
1007 double eta = (t_ieta > 0) ? t_ieta : -t_ieta;
1008 cut = 8.0 * exp(eta * log2by18_);
1009 }
1010 correctEnergy(eHcal);
1011 select = ((t_qltyFlag) && (t_selectTk) && (t_hmaxNearP < cut) && (t_eMipDR < 100.0));
1012 if (debug) {
1013 std::cout << " output " << eHcal << ":" << cut << ":" << select << std::endl;
1014 }
1015 return select;
1016 }
1017
1018 bool CalibPlotProperties::selectPhi(bool debug) {
1019 bool select(true);
1020 if (phimin_ > 1 || phimax_ < 72) {
1021 double eTotal(0), eSelec(0);
1022
1023 for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
1024 int iphi = ((*t_DetIds)[k]) & (0x3FF);
1025 int zside = ((*t_DetIds)[k] & 0x80000) ? (1) : (-1);
1026 eTotal += ((*t_HitEnergies)[k]);
1027 if (iphi >= phimin_ && iphi <= phimax_ && zside == zside_)
1028 eSelec += ((*t_HitEnergies)[k]);
1029 }
1030 if (eSelec < 0.9 * eTotal)
1031 select = false;
1032 if (debug) {
1033 std::cout << "Etotal " << eTotal << " and ESelec " << eSelec << " (phi " << phimin_ << ":" << phimax_ << " z "
1034 << zside_ << ") Selection " << select << std::endl;
1035 }
1036 }
1037 return select;
1038 }
1039
1040 void CalibPlotProperties::savePlot(const std::string &theName, bool append, bool all, bool debug) {
1041 TFile *theFile(0);
1042 if (append) {
1043 theFile = new TFile(theName.c_str(), "UPDATE");
1044 } else {
1045 theFile = new TFile(theName.c_str(), "RECREATE");
1046 }
1047
1048 theFile->cd();
1049
1050 if (plotBasic_) {
1051 for (int k = 0; k < CalibPlots::ntitles; ++k) {
1052 if (debug)
1053 std::cout << "[" << k << "] p " << h_p[k] << " eta " << h_eta[k] << std::endl;
1054 h_p[k]->Write();
1055 h_eta[k]->Write();
1056 }
1057 for (int k = 0; k < CalibPlots::npbin; ++k) {
1058 if (debug)
1059 std::cout << "[" << k << "] eta " << h_eta0[k] << ":" << h_eta1[k] << ":" << h_eta2[k] << ":" << h_eta3[k]
1060 << ":" << h_eta4[k] << " dl " << h_dL1[k] << " vtx " << h_vtx[k] << std::endl;
1061 if (h_eta0[k] != 0) {
1062 TH1D *h1 = (TH1D *)h_eta0[k]->Clone();
1063 h1->Write();
1064 }
1065 if (h_eta1[k] != 0) {
1066 TH1D *h2 = (TH1D *)h_eta1[k]->Clone();
1067 h2->Write();
1068 }
1069 if (h_eta2[k] != 0) {
1070 TH1D *h3 = (TH1D *)h_eta2[k]->Clone();
1071 h3->Write();
1072 }
1073 if (h_eta3[k] != 0) {
1074 TH1D *h4 = (TH1D *)h_eta3[k]->Clone();
1075 h4->Write();
1076 }
1077 if (h_eta4[k] != 0) {
1078 TH1D *h5 = (TH1D *)h_eta4[k]->Clone();
1079 h5->Write();
1080 }
1081 if (h_dL1[k] != 0) {
1082 TH1D *h6 = (TH1D *)h_dL1[k]->Clone();
1083 h6->Write();
1084 }
1085 if (h_vtx[k] != 0) {
1086 TH1D *h7 = (TH1D *)h_vtx[k]->Clone();
1087 h7->Write();
1088 }
1089 }
1090 }
1091
1092 if (plotEnergy_) {
1093 for (int k = 0; k < CalibPlots::npbin; ++k) {
1094 if (debug)
1095 std::cout << "Energy[" << k << "] "
1096 << " eta " << etalo_ << ":" << etahi_ << ":" << CalibPlots::netabin << " etaEH " << h_etaEH[k].size()
1097 << " etaEp " << h_etaEp[k].size() << " etaEE " << h_etaEE[k].size() << " eHcal " << h_eHcal.size()
1098 << " mom " << h_mom.size() << " eEcal " << h_eEcal.size() << std::endl;
1099 for (unsigned int j = 0; j <= (unsigned int)(etahi_ - etalo_); ++j) {
1100 if (h_etaEH[k].size() > j && h_etaEH[k][j] != nullptr && (all || (k == kp50))) {
1101 TH1D *hist = (TH1D *)h_etaEH[k][j]->Clone();
1102 hist->Write();
1103 }
1104 if (h_etaEp[k].size() > j && h_etaEp[k][j] != nullptr && (all || (k == kp50))) {
1105 TH1D *hist = (TH1D *)h_etaEp[k][j]->Clone();
1106 hist->Write();
1107 }
1108 if (h_etaEE[k].size() > j && h_etaEE[k][j] != nullptr && (all || (k == kp50))) {
1109 TH1D *hist = (TH1D *)h_etaEE[k][j]->Clone();
1110 hist->Write();
1111 }
1112 if (h_etaEE0[k].size() > j && h_etaEE0[k][j] != nullptr && (all || (k == kp50))) {
1113 TH1D *hist = (TH1D *)h_etaEE0[k][j]->Clone();
1114 hist->Write();
1115 }
1116 }
1117 }
1118
1119 for (int j = 0; j < CalibPlots::netabin; ++j) {
1120 if (h_eHcal.size() > (unsigned int)(j) && (h_eHcal[j] != nullptr)) {
1121 TH1D *hist = (TH1D *)h_eHcal[j]->Clone();
1122 hist->Write();
1123 }
1124 if (h_mom.size() > (unsigned int)(j) && (h_mom[j] != nullptr)) {
1125 TH1D *hist = (TH1D *)h_mom[j]->Clone();
1126 hist->Write();
1127 }
1128 if (h_eEcal.size() > (unsigned int)(j) && (h_eEcal[j] != nullptr)) {
1129 TH1D *hist = (TH1D *)h_eEcal[j]->Clone();
1130 hist->Write();
1131 }
1132 }
1133 }
1134
1135 if (plotHists_) {
1136 if (debug)
1137 std::cout << "nvtx " << h_nvtx << " etaE " << h_etaE << std::endl;
1138 h_nvtx->Write();
1139 h_etaE->Write();
1140 for (int i = 0; i < CalibPlots::ndepth; ++i) {
1141 if (debug)
1142 std::cout << "Depth[" << i << "] bv " << h_bvlist[i] << ":" << h_bvlist2[i] << ":" << h_bvlist3[i] << " ev "
1143 << h_evlist[i] << ":" << h_evlist2[i] << ":" << h_evlist3[i] << std::endl;
1144 h_bvlist[i]->Write();
1145 h_bvlist2[i]->Write();
1146 h_bvlist3[i]->Write();
1147 h_evlist[i]->Write();
1148 h_evlist2[i]->Write();
1149 h_evlist3[i]->Write();
1150 }
1151 }
1152 std::cout << "All done" << std::endl;
1153 theFile->Close();
1154 }
1155
1156 void CalibPlotProperties::correctEnergy(double &eHcal) {
1157 double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
1158 if ((corrPU_ < 0) && (pmom > 0)) {
1159 double ediff = (t_eHcal30 - t_eHcal10);
1160 if (t_DetIds1 != 0 && t_DetIds3 != 0) {
1161 double Etot1(0), Etot3(0);
1162
1163 for (unsigned int idet = 0; idet < (*t_DetIds1).size(); idet++) {
1164 unsigned int id = truncateId((*t_DetIds1)[idet], truncateFlag_, false);
1165 double cfac = corrFactor_->getCorr(id);
1166 if (cFactor_ != 0)
1167 cfac *= cFactor_->getCorr(t_Run, (*t_DetIds1)[idet]);
1168 double hitEn = cfac * (*t_HitEnergies1)[idet];
1169 Etot1 += hitEn;
1170 }
1171 for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
1172 unsigned int id = truncateId((*t_DetIds3)[idet], truncateFlag_, false);
1173 double cfac = corrFactor_->getCorr(id);
1174 if (cFactor_ != 0)
1175 cfac *= cFactor_->getCorr(t_Run, (*t_DetIds3)[idet]);
1176 double hitEn = cfac * (*t_HitEnergies3)[idet];
1177 Etot3 += hitEn;
1178 }
1179 ediff = (Etot3 - Etot1);
1180 }
1181 double fac = puFactor(-corrPU_, t_ieta, pmom, eHcal, ediff);
1182 eHcal *= fac;
1183 } else if (corrPU_ > 1) {
1184 eHcal = puFactorRho(corrPU_, t_ieta, t_rhoh, eHcal);
1185 }
1186 }
1187
1188 void PlotThisHist(TH1D *hist, const std::string &text, int save) {
1189 char namep[120];
1190 sprintf(namep, "c_%s", hist->GetName());
1191 TCanvas *pad = new TCanvas(namep, namep, 700, 500);
1192 pad->SetRightMargin(0.10);
1193 pad->SetTopMargin(0.10);
1194 hist->GetXaxis()->SetTitleSize(0.04);
1195 hist->GetYaxis()->SetTitle("Tracks");
1196 hist->GetYaxis()->SetLabelOffset(0.005);
1197 hist->GetYaxis()->SetTitleSize(0.04);
1198 hist->GetYaxis()->SetLabelSize(0.035);
1199 hist->GetYaxis()->SetTitleOffset(1.10);
1200 hist->SetMarkerStyle(20);
1201 hist->SetMarkerColor(2);
1202 hist->SetLineColor(2);
1203 hist->Draw("Hist");
1204 pad->Modified();
1205 pad->Update();
1206 TPaveStats *st1 = (TPaveStats *)hist->GetListOfFunctions()->FindObject("stats");
1207 TPaveText *txt0 = new TPaveText(0.12, 0.91, 0.49, 0.96, "blNDC");
1208 txt0->SetFillColor(0);
1209 char txt[100];
1210 sprintf(txt, "CMS Simulation Preliminary");
1211 txt0->AddText(txt);
1212 txt0->Draw("same");
1213 TPaveText *txt1 = new TPaveText(0.51, 0.91, 0.90, 0.96, "blNDC");
1214 txt1->SetFillColor(0);
1215 sprintf(txt, "%s", text.c_str());
1216 txt1->AddText(txt);
1217 txt1->Draw("same");
1218 if (st1 != nullptr) {
1219 st1->SetY1NDC(0.70);
1220 st1->SetY2NDC(0.90);
1221 st1->SetX1NDC(0.65);
1222 st1->SetX2NDC(0.90);
1223 }
1224 pad->Update();
1225 if (save != 0) {
1226 if (save > 0)
1227 sprintf(namep, "%s.pdf", pad->GetName());
1228 else
1229 sprintf(namep, "%s.jpg", pad->GetName());
1230 pad->Print(namep);
1231 }
1232 }
1233
1234 void PlotHist(const char *hisFileName,
1235 const std::string &prefix = "",
1236 const std::string &text = "",
1237 int flagC = 111,
1238 int etalo = 0,
1239 int etahi = 30,
1240 int save = 0) {
1241 gStyle->SetCanvasBorderMode(0);
1242 gStyle->SetCanvasColor(kWhite);
1243 gStyle->SetPadColor(kWhite);
1244 gStyle->SetFillColor(kWhite);
1245 gStyle->SetOptTitle(0);
1246 gStyle->SetOptStat(1110);
1247
1248 bool plotBasic = (((flagC / 1) % 10) > 0);
1249 bool plotEnergy = (((flagC / 10) % 10) > 0);
1250 bool plotHists = (((flagC / 100) % 10) > 0);
1251
1252 TFile *file = new TFile(hisFileName);
1253 char name[100], title[200];
1254 TH1D *hist;
1255 if ((file != nullptr) && plotBasic) {
1256 std::cout << "Plot Basic Histos" << std::endl;
1257 for (int k = 0; k < CalibPlots::ntitles; ++k) {
1258 sprintf(name, "%sp%d", prefix.c_str(), k);
1259 hist = (TH1D *)(file->FindObjectAny(name));
1260 if (hist != nullptr) {
1261 sprintf(title, "Momentum for %s (GeV)", CalibPlots::getTitle(k).c_str());
1262 hist->GetXaxis()->SetTitle(title);
1263 PlotThisHist(hist, text, save);
1264 }
1265 sprintf(name, "%seta%d", prefix.c_str(), k);
1266 hist = (TH1D *)(file->FindObjectAny(name));
1267 if (hist != nullptr) {
1268 sprintf(title, "#eta for %s", CalibPlots::getTitle(k).c_str());
1269 hist->GetXaxis()->SetTitle(title);
1270 PlotThisHist(hist, text, save);
1271 }
1272 }
1273 for (int k = 0; k < CalibPlots::npbin; ++k) {
1274 sprintf(name, "%seta0%d", prefix.c_str(), k);
1275 hist = (TH1D *)(file->FindObjectAny(name));
1276 if (hist != nullptr) {
1277 sprintf(title,
1278 "#eta for %s (p = %d:%d GeV)",
1279 CalibPlots::getTitle(0).c_str(),
1280 CalibPlots::getP(k),
1281 CalibPlots::getP(k + 1));
1282 hist->GetXaxis()->SetTitle(title);
1283 PlotThisHist(hist, text, save);
1284 }
1285 sprintf(name, "%seta1%d", prefix.c_str(), k);
1286 hist = (TH1D *)(file->FindObjectAny(name));
1287 if (hist != nullptr) {
1288 sprintf(title,
1289 "#eta for %s (p = %d:%d GeV)",
1290 CalibPlots::getTitle(1).c_str(),
1291 CalibPlots::getP(k),
1292 CalibPlots::getP(k + 1));
1293 hist->GetXaxis()->SetTitle(title);
1294 PlotThisHist(hist, text, save);
1295 }
1296 sprintf(name, "%seta2%d", prefix.c_str(), k);
1297 hist = (TH1D *)(file->FindObjectAny(name));
1298 if (hist != nullptr) {
1299 sprintf(title,
1300 "#eta for %s (p = %d:%d GeV)",
1301 CalibPlots::getTitle(2).c_str(),
1302 CalibPlots::getP(k),
1303 CalibPlots::getP(k + 1));
1304 hist->GetXaxis()->SetTitle(title);
1305 PlotThisHist(hist, text, save);
1306 }
1307 sprintf(name, "%seta3%d", prefix.c_str(), k);
1308 hist = (TH1D *)(file->FindObjectAny(name));
1309 if (hist != nullptr) {
1310 sprintf(title,
1311 "#eta for %s (p = %d:%d GeV)",
1312 CalibPlots::getTitle(3).c_str(),
1313 CalibPlots::getP(k),
1314 CalibPlots::getP(k + 1));
1315 hist->GetXaxis()->SetTitle(title);
1316 PlotThisHist(hist, text, save);
1317 }
1318 sprintf(name, "%seta4%d", prefix.c_str(), k);
1319 hist = (TH1D *)(file->FindObjectAny(name));
1320 if (hist != nullptr) {
1321 sprintf(title,
1322 "#eta for %s (p = %d:%d GeV)",
1323 CalibPlots::getTitle(4).c_str(),
1324 CalibPlots::getP(k),
1325 CalibPlots::getP(k + 1));
1326 hist->GetXaxis()->SetTitle(title);
1327 PlotThisHist(hist, text, save);
1328 }
1329 sprintf(name, "%sdl1%d", prefix.c_str(), k);
1330 hist = (TH1D *)(file->FindObjectAny(name));
1331 if (hist != nullptr) {
1332 sprintf(title, "Distance from L1 (p = %d:%d GeV)", CalibPlots::getP(k), CalibPlots::getP(k + 1));
1333 hist->GetXaxis()->SetTitle(title);
1334 PlotThisHist(hist, text, save);
1335 }
1336 sprintf(name, "%svtx%d", prefix.c_str(), k);
1337 hist = (TH1D *)(file->FindObjectAny(name));
1338 if (hist != nullptr) {
1339 sprintf(title, "N_{Vertex} (p = %d:%d GeV)", CalibPlots::getP(k), CalibPlots::getP(k + 1));
1340 hist->GetXaxis()->SetTitle(title);
1341 PlotThisHist(hist, text, save);
1342 }
1343 }
1344 }
1345
1346 if ((file != nullptr) && plotEnergy) {
1347 std::cout << "Make plots for good tracks" << std::endl;
1348 for (int k = 0; k < CalibPlots::npbin; ++k) {
1349 for (int j = etalo; j <= etahi + 1; ++j) {
1350 sprintf(name, "%senergyH%d%d", prefix.c_str(), k, j);
1351 hist = (TH1D *)(file->FindObjectAny(name));
1352 if (hist != nullptr) {
1353 if (j > etahi)
1354 sprintf(title,
1355 "HCAL energy for %s (p = %d:%d GeV |#eta| = %d:%d)",
1356 CalibPlots::getTitle(3).c_str(),
1357 CalibPlots::getP(k),
1358 CalibPlots::getP(k + 1),
1359 etalo,
1360 etahi);
1361 else
1362 sprintf(title,
1363 "HCAL energy for %s (p = %d:%d GeV |#eta| = %d)",
1364 CalibPlots::getTitle(3).c_str(),
1365 CalibPlots::getP(k),
1366 CalibPlots::getP(k + 1),
1367 j);
1368 hist->GetXaxis()->SetTitle(title);
1369 PlotThisHist(hist, text, save);
1370 }
1371 sprintf(name, "%senergyP%d%d", prefix.c_str(), k, j);
1372 hist = (TH1D *)(file->FindObjectAny(name));
1373 if (hist != nullptr) {
1374 if (j > etahi)
1375 sprintf(title,
1376 "momentum for %s (p = %d:%d GeV |#eta| = %d:%d)",
1377 CalibPlots::getTitle(3).c_str(),
1378 CalibPlots::getP(k),
1379 CalibPlots::getP(k + 1),
1380 etalo,
1381 etahi);
1382 else
1383 sprintf(title,
1384 "momentum for %s (p = %d:%d GeV |#eta| = %d)",
1385 CalibPlots::getTitle(3).c_str(),
1386 CalibPlots::getP(k),
1387 CalibPlots::getP(k + 1),
1388 j);
1389 hist->GetXaxis()->SetTitle(title);
1390 PlotThisHist(hist, text, save);
1391 }
1392 sprintf(name, "%senergyE%d%d", prefix.c_str(), k, j);
1393 hist = (TH1D *)(file->FindObjectAny(name));
1394 if (hist != nullptr) {
1395 if (j > etahi)
1396 sprintf(title,
1397 "ECAL energy for %s (p = %d:%d GeV |#eta| = %d:%d)",
1398 CalibPlots::getTitle(3).c_str(),
1399 CalibPlots::getP(k),
1400 CalibPlots::getP(k + 1),
1401 etalo,
1402 etahi);
1403 else
1404 sprintf(title,
1405 "ECAL energy for %s (p = %d:%d GeV |#eta| = %d)",
1406 CalibPlots::getTitle(3).c_str(),
1407 CalibPlots::getP(k),
1408 CalibPlots::getP(k + 1),
1409 j);
1410 hist->GetXaxis()->SetTitle(title);
1411 PlotThisHist(hist, text, save);
1412 }
1413 sprintf(name, "%senergyER%d%d", prefix.c_str(), k, j);
1414 hist = (TH1D *)(file->FindObjectAny(name));
1415 if (hist != nullptr) {
1416 std::cout << name << " Mean " << hist->GetMean() << " +- " << hist->GetMeanError() << " Entries "
1417 << hist->GetEntries() << " RMS " << hist->GetRMS() << std::endl;
1418 }
1419 }
1420 }
1421
1422 for (int j = 0; j < CalibPlots::netabin; ++j) {
1423 sprintf(name, "%senergyH%d", prefix.c_str(), j);
1424 hist = (TH1D *)(file->FindObjectAny(name));
1425 if (hist != nullptr) {
1426 if (j == 0)
1427 sprintf(title, "HCAL energy for %s (All)", CalibPlots::getTitle(3).c_str());
1428 else
1429 sprintf(title,
1430 "HCAL energy for %s (|#eta| = %d:%d)",
1431 CalibPlots::getTitle(3).c_str(),
1432 CalibPlots::getEta(j - 1),
1433 CalibPlots::getEta(j));
1434 hist->GetXaxis()->SetTitle(title);
1435 PlotThisHist(hist, text, save);
1436 }
1437 sprintf(name, "%senergyP%d", prefix.c_str(), j);
1438 hist = (TH1D *)(file->FindObjectAny(name));
1439 if (hist != nullptr) {
1440 if (j == 0)
1441 sprintf(title, "Track momentum for %s (All)", CalibPlots::getTitle(3).c_str());
1442 else
1443 sprintf(title,
1444 "Track momentum for %s (|#eta| = %d:%d)",
1445 CalibPlots::getTitle(3).c_str(),
1446 CalibPlots::getEta(j - 1),
1447 CalibPlots::getEta(j));
1448 hist->GetXaxis()->SetTitle(title);
1449 PlotThisHist(hist, text, save);
1450 }
1451 sprintf(name, "%senergyE%d", prefix.c_str(), j);
1452 hist = (TH1D *)(file->FindObjectAny(name));
1453 if (hist != nullptr) {
1454 if (j == 0)
1455 sprintf(title, "ECAL energy for %s (All)", CalibPlots::getTitle(3).c_str());
1456 else
1457 sprintf(title,
1458 "ECAL energy for %s (|#eta| = %d:%d)",
1459 CalibPlots::getTitle(3).c_str(),
1460 CalibPlots::getEta(j - 1),
1461 CalibPlots::getEta(j));
1462 hist->GetXaxis()->SetTitle(title);
1463 PlotThisHist(hist, text, save);
1464 }
1465 }
1466 }
1467
1468 if (plotHists) {
1469 hist = (TH1D *)(file->FindObjectAny("hnvtx"));
1470 if (hist != nullptr) {
1471 hist->GetXaxis()->SetTitle("Number of vertices");
1472 PlotThisHist(hist, text, save);
1473 }
1474 for (int i = 0; i < CalibPlots::ndepth; i++) {
1475 sprintf(name, "b_edepth%d", i);
1476 hist = (TH1D *)(file->FindObjectAny(name));
1477 if (hist != nullptr) {
1478 sprintf(title, "Total RecHit energy in depth %d (Barrel)", i + 1);
1479 hist->GetXaxis()->SetTitle(title);
1480 PlotThisHist(hist, text, save);
1481 }
1482 sprintf(name, "b_recedepth%d", i);
1483 hist = (TH1D *)(file->FindObjectAny(name));
1484 if (hist != nullptr) {
1485 sprintf(title, "RecHit energy in depth %d (Barrel)", i + 1);
1486 hist->GetXaxis()->SetTitle(title);
1487 PlotThisHist(hist, text, save);
1488 }
1489 sprintf(name, "b_nrecdepth%d", i);
1490 hist = (TH1D *)(file->FindObjectAny(name));
1491 if (hist != nullptr) {
1492 sprintf(title, "#RecHits in depth %d (Barrel)", i + 1);
1493 hist->GetXaxis()->SetTitle(title);
1494 PlotThisHist(hist, text, save);
1495 }
1496 sprintf(name, "e_edepth%d", i);
1497 hist = (TH1D *)(file->FindObjectAny(name));
1498 if (hist != nullptr) {
1499 sprintf(title, "Total RecHit energy in depth %d (Endcap)", i + 1);
1500 hist->GetXaxis()->SetTitle(title);
1501 PlotThisHist(hist, text, save);
1502 }
1503 sprintf(name, "e_recedepth%d", i);
1504 hist = (TH1D *)(file->FindObjectAny(name));
1505 if (hist != nullptr) {
1506 sprintf(title, "RecHit energy in depth %d (Endcap)", i + 1);
1507 hist->GetXaxis()->SetTitle(title);
1508 PlotThisHist(hist, text, save);
1509 }
1510 sprintf(name, "e_nrecdepth%d", i);
1511 hist = (TH1D *)(file->FindObjectAny(name));
1512 if (hist != nullptr) {
1513 sprintf(title, "#RecHits in depth %d (Endcap)", i + 1);
1514 hist->GetXaxis()->SetTitle(title);
1515 PlotThisHist(hist, text, save);
1516 }
1517 }
1518 TH2F *h_etaE = (TH2F *)(file->FindObjectAny("heta"));
1519 if (h_etaE != nullptr) {
1520 h_etaE->GetXaxis()->SetTitle("i#eta");
1521 h_etaE->GetYaxis()->SetTitle("Energy (GeV)");
1522 char namep[120];
1523 sprintf(namep, "c_%s", h_etaE->GetName());
1524 TCanvas *pad = new TCanvas(namep, namep, 700, 700);
1525 pad->SetRightMargin(0.10);
1526 pad->SetTopMargin(0.10);
1527 h_etaE->GetXaxis()->SetTitleSize(0.04);
1528 h_etaE->GetYaxis()->SetLabelOffset(0.005);
1529 h_etaE->GetYaxis()->SetTitleSize(0.04);
1530 h_etaE->GetYaxis()->SetLabelSize(0.035);
1531 h_etaE->GetYaxis()->SetTitleOffset(1.10);
1532 h_etaE->SetMarkerStyle(20);
1533 h_etaE->SetMarkerColor(2);
1534 h_etaE->SetLineColor(2);
1535 h_etaE->Draw();
1536 pad->Update();
1537 TPaveStats *st1 = (TPaveStats *)h_etaE->GetListOfFunctions()->FindObject("stats");
1538 if (st1 != nullptr) {
1539 st1->SetY1NDC(0.70);
1540 st1->SetY2NDC(0.90);
1541 st1->SetX1NDC(0.65);
1542 st1->SetX2NDC(0.90);
1543 }
1544 if (save != 0) {
1545 if (save > 0)
1546 sprintf(namep, "%s.pdf", pad->GetName());
1547 else
1548 sprintf(namep, "%s.jpg", pad->GetName());
1549 pad->Print(namep);
1550 }
1551 }
1552 }
1553 }