File indexing completed on 2022-04-11 04:18:28
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 #include <TROOT.h>
0143 #include <TChain.h>
0144 #include <TFile.h>
0145 #include <TF1.h>
0146 #include <TH1D.h>
0147 #include <TH2F.h>
0148 #include <TProfile.h>
0149 #include <TFitResult.h>
0150 #include <TFitResultPtr.h>
0151 #include <TStyle.h>
0152 #include <TCanvas.h>
0153 #include <TLegend.h>
0154 #include <TPaveStats.h>
0155 #include <TPaveText.h>
0156
0157 #include <algorithm>
0158 #include <iomanip>
0159 #include <iostream>
0160 #include <fstream>
0161 #include <map>
0162 #include <vector>
0163 #include <string>
0164
0165 #include "CalibCorr.C"
0166
0167 class CalibMonitor {
0168 public:
0169 TChain *fChain;
0170 Int_t fCurrent;
0171
0172
0173 Int_t t_Run;
0174 Int_t t_Event;
0175 Int_t t_DataType;
0176 Int_t t_ieta;
0177 Int_t t_iphi;
0178 Double_t t_EventWeight;
0179 Int_t t_nVtx;
0180 Int_t t_nTrk;
0181 Int_t t_goodPV;
0182 Double_t t_l1pt;
0183 Double_t t_l1eta;
0184 Double_t t_l1phi;
0185 Double_t t_l3pt;
0186 Double_t t_l3eta;
0187 Double_t t_l3phi;
0188 Double_t t_p;
0189 Double_t t_pt;
0190 Double_t t_phi;
0191 Double_t t_mindR1;
0192 Double_t t_mindR2;
0193 Double_t t_eMipDR;
0194 Double_t t_eHcal;
0195 Double_t t_eHcal10;
0196 Double_t t_eHcal30;
0197 Double_t t_hmaxNearP;
0198 Double_t t_rhoh;
0199 Bool_t t_selectTk;
0200 Bool_t t_qltyFlag;
0201 Bool_t t_qltyMissFlag;
0202 Bool_t t_qltyPVFlag;
0203 Double_t t_gentrackP;
0204 std::vector<unsigned int> *t_DetIds;
0205 std::vector<double> *t_HitEnergies;
0206 std::vector<bool> *t_trgbits;
0207 std::vector<unsigned int> *t_DetIds1;
0208 std::vector<unsigned int> *t_DetIds3;
0209 std::vector<double> *t_HitEnergies1;
0210 std::vector<double> *t_HitEnergies3;
0211
0212
0213 TBranch *b_t_Run;
0214 TBranch *b_t_Event;
0215 TBranch *b_t_DataType;
0216 TBranch *b_t_ieta;
0217 TBranch *b_t_iphi;
0218 TBranch *b_t_EventWeight;
0219 TBranch *b_t_nVtx;
0220 TBranch *b_t_nTrk;
0221 TBranch *b_t_goodPV;
0222 TBranch *b_t_l1pt;
0223 TBranch *b_t_l1eta;
0224 TBranch *b_t_l1phi;
0225 TBranch *b_t_l3pt;
0226 TBranch *b_t_l3eta;
0227 TBranch *b_t_l3phi;
0228 TBranch *b_t_p;
0229 TBranch *b_t_pt;
0230 TBranch *b_t_phi;
0231 TBranch *b_t_mindR1;
0232 TBranch *b_t_mindR2;
0233 TBranch *b_t_eMipDR;
0234 TBranch *b_t_eHcal;
0235 TBranch *b_t_eHcal10;
0236 TBranch *b_t_eHcal30;
0237 TBranch *b_t_hmaxNearP;
0238 TBranch *b_t_rhoh;
0239 TBranch *b_t_selectTk;
0240 TBranch *b_t_qltyFlag;
0241 TBranch *b_t_qltyMissFlag;
0242 TBranch *b_t_qltyPVFlag;
0243 TBranch *b_t_gentrackP;
0244 TBranch *b_t_DetIds;
0245 TBranch *b_t_HitEnergies;
0246 TBranch *b_t_trgbits;
0247 TBranch *b_t_DetIds1;
0248 TBranch *b_t_DetIds3;
0249 TBranch *b_t_HitEnergies1;
0250 TBranch *b_t_HitEnergies3;
0251
0252 struct counter {
0253 static const int npsize = 4;
0254 counter() {
0255 total = 0;
0256 for (int k = 0; k < npsize; ++k)
0257 count[k] = 0;
0258 };
0259 unsigned int total, count[npsize];
0260 };
0261
0262 CalibMonitor(const char *fname,
0263 const std::string &dirname,
0264 const char *dupFileName,
0265 const char *comFileName,
0266 const char *outFileName,
0267 const std::string &prefix = "",
0268 const char *corrFileName = "",
0269 const char *rcorFileName = "",
0270 int puCorr = -2,
0271 int flag = 1031,
0272 int numb = 50,
0273 bool datMC = true,
0274 int truncateFlag = 0,
0275 bool useGen = false,
0276 double scale = 1.0,
0277 int useScale = 0,
0278 int etalo = 0,
0279 int etahi = 30,
0280 int runlo = 0,
0281 int runhi = 99999999,
0282 int phimin = 1,
0283 int phimax = 72,
0284 int zside = 1,
0285 int nvxlo = 0,
0286 int nvxhi = 1000,
0287 int rbx = 0,
0288 bool exclude = false,
0289 bool etamax = false);
0290 virtual ~CalibMonitor();
0291 virtual Int_t Cut(Long64_t entry);
0292 virtual Int_t GetEntry(Long64_t entry);
0293 virtual Long64_t LoadTree(Long64_t entry);
0294 virtual void Init(TChain *, const char *, const char *, const char *);
0295 virtual void Loop();
0296 virtual Bool_t Notify();
0297 virtual void Show(Long64_t entry = -1);
0298 bool goodTrack(double &eHcal, double &cut, const Long64_t &entry, bool debug);
0299 bool selectPhi(bool debug);
0300 void plotHist(int type, int num, bool save = false);
0301 template <class Hist>
0302 void drawHist(Hist *, TCanvas *);
0303 void savePlot(const std::string &theName, bool append = true, bool all = false);
0304 void correctEnergy(double &ener, const Long64_t &entry);
0305
0306 private:
0307 static const unsigned int npbin = 6, kp50 = 3;
0308 CalibCorrFactor *corrFactor_;
0309 CalibCorr *cFactor_;
0310 CalibSelectRBX *cSelect_;
0311 const std::string fname_, dirnm_, prefix_, outFileName_;
0312 const int corrPU_, flag_, numb_;
0313 const bool dataMC_, useGen_;
0314 const int truncateFlag_;
0315 const int etalo_, etahi_;
0316 int runlo_, runhi_;
0317 const int phimin_, phimax_, zside_, nvxlo_, nvxhi_, rbx_;
0318 bool exclude_, corrE_, cutL1T_, selRBX_;
0319 bool includeRun_;
0320 int coarseBin_, etamp_, etamn_, plotType_;
0321 int flexibleSelect_, ifDepth_;
0322 double log2by18_;
0323 std::ofstream fileout_;
0324 std::vector<Long64_t> entries_;
0325 std::vector<std::pair<int, int> > events_;
0326 std::vector<double> etas_, ps_, dl1_;
0327 std::vector<int> nvx_, ietasL_, ietasH_;
0328 std::vector<TH1D *> h_rbx, h_etaF[npbin], h_etaB[npbin];
0329 std::vector<TProfile *> h_etaX[npbin];
0330 std::vector<TH1D *> h_etaR[npbin], h_nvxR[npbin], h_dL1R[npbin];
0331 std::vector<TH1D *> h_pp[npbin];
0332 };
0333
0334 CalibMonitor::CalibMonitor(const char *fname,
0335 const std::string &dirnm,
0336 const char *dupFileName,
0337 const char *comFileName,
0338 const char *outFName,
0339 const std::string &prefix,
0340 const char *corrFileName,
0341 const char *rcorFileName,
0342 int puCorr,
0343 int flag,
0344 int numb,
0345 bool dataMC,
0346 int truncate,
0347 bool useGen,
0348 double scale,
0349 int useScale,
0350 int etalo,
0351 int etahi,
0352 int runlo,
0353 int runhi,
0354 int phimin,
0355 int phimax,
0356 int zside,
0357 int nvxlo,
0358 int nvxhi,
0359 int rbx,
0360 bool exc,
0361 bool etam)
0362 : corrFactor_(nullptr),
0363 cFactor_(nullptr),
0364 cSelect_(nullptr),
0365 fname_(fname),
0366 dirnm_(dirnm),
0367 prefix_(prefix),
0368 outFileName_(std::string(outFName)),
0369 corrPU_(puCorr),
0370 flag_(flag),
0371 numb_(numb),
0372 dataMC_(dataMC),
0373 useGen_(useGen),
0374 truncateFlag_(truncate),
0375 etalo_(etalo),
0376 etahi_(etahi),
0377 runlo_(runlo),
0378 runhi_(runhi),
0379 phimin_(phimin),
0380 phimax_(phimax),
0381 zside_(zside),
0382 nvxlo_(nvxlo),
0383 nvxhi_(nvxhi),
0384 rbx_(rbx),
0385 exclude_(exc),
0386 includeRun_(true) {
0387
0388
0389
0390 plotType_ = ((flag_ / 10) % 10);
0391 if (plotType_ < 0 || plotType_ > 3)
0392 plotType_ = 3;
0393 flexibleSelect_ = (((flag_ / 1) % 10));
0394 int oneplace = ((flag_ / 1000) % 10);
0395 cutL1T_ = (oneplace % 2);
0396 bool marina = ((oneplace / 2) % 2);
0397 ifDepth_ = ((flag_ / 10000) % 10);
0398 selRBX_ = (((flag_ / 100000) % 10) > 0);
0399 coarseBin_ = ((flag_ / 1000000) % 10);
0400 log2by18_ = std::log(2.5) / 18.0;
0401 if (runlo_ < 0 || runhi_ < 0) {
0402 runlo_ = std::abs(runlo_);
0403 runhi_ = std::abs(runhi_);
0404 includeRun_ = false;
0405 }
0406 char treeName[400];
0407 sprintf(treeName, "%s/CalibTree", dirnm.c_str());
0408 TChain *chain = new TChain(treeName);
0409 std::cout << "Create a chain for " << treeName << " from " << fname << " flags " << flexibleSelect_ << "|"
0410 << plotType_ << "|" << corrPU_ << "\n cons " << log2by18_ << " eta range " << etalo_ << ":" << etahi_
0411 << " run range " << runlo_ << ":" << runhi_ << " (inclusion flag " << includeRun_ << ")\n Selection of RBX "
0412 << selRBX_ << " Vertex Range " << nvxlo_ << ":" << nvxhi_ << "\n corrFileName: " << corrFileName
0413 << " useScale " << useScale << ":" << scale << ":" << etam << "\n rcorFileName: " << rcorFileName
0414 << " flag " << ifDepth_ << ":" << cutL1T_ << ":" << marina << std::endl;
0415 if (!fillChain(chain, fname)) {
0416 std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
0417 } else {
0418 std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
0419 corrFactor_ = new CalibCorrFactor(corrFileName, useScale, scale, etam, marina, false);
0420 Init(chain, dupFileName, comFileName, outFName);
0421 if (std::string(rcorFileName) != "") {
0422 cFactor_ = new CalibCorr(rcorFileName, ifDepth_, false);
0423 } else {
0424 ifDepth_ = 0;
0425 }
0426 if (rbx != 0)
0427 cSelect_ = new CalibSelectRBX(rbx, false);
0428 }
0429 }
0430
0431 CalibMonitor::~CalibMonitor() {
0432 delete corrFactor_;
0433 delete cFactor_;
0434 delete cSelect_;
0435 if (!fChain)
0436 return;
0437 delete fChain->GetCurrentFile();
0438 }
0439
0440 Int_t CalibMonitor::GetEntry(Long64_t entry) {
0441
0442 if (!fChain)
0443 return 0;
0444 return fChain->GetEntry(entry);
0445 }
0446
0447 Long64_t CalibMonitor::LoadTree(Long64_t entry) {
0448
0449 if (!fChain)
0450 return -5;
0451 Long64_t centry = fChain->LoadTree(entry);
0452 if (centry < 0)
0453 return centry;
0454 if (!fChain->InheritsFrom(TChain::Class()))
0455 return centry;
0456 TChain *chain = (TChain *)fChain;
0457 if (chain->GetTreeNumber() != fCurrent) {
0458 fCurrent = chain->GetTreeNumber();
0459 Notify();
0460 }
0461 return centry;
0462 }
0463
0464 void CalibMonitor::Init(TChain *tree, const char *dupFileName, const char *comFileName, const char *outFileName) {
0465
0466
0467
0468
0469
0470
0471
0472
0473
0474 t_DetIds = 0;
0475 t_DetIds1 = 0;
0476 t_DetIds3 = 0;
0477 t_HitEnergies = 0;
0478 t_HitEnergies1 = 0;
0479 t_HitEnergies3 = 0;
0480 t_trgbits = 0;
0481
0482 fChain = tree;
0483 fCurrent = -1;
0484 if (!tree)
0485 return;
0486 fChain->SetMakeClass(1);
0487
0488 fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
0489 fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
0490 fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
0491 fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
0492 fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
0493 fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
0494 fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
0495 fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
0496 fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
0497 fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
0498 fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
0499 fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
0500 fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
0501 fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
0502 fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
0503 fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
0504 fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
0505 fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
0506 fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
0507 fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
0508 fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
0509 fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
0510 fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
0511 fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
0512 fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
0513 fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
0514 fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
0515 fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
0516 fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
0517 fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
0518 fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
0519 fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
0520 fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
0521 fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
0522 fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
0523 fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
0524 fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
0525 fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
0526 Notify();
0527
0528 if (strcmp(dupFileName, "") != 0) {
0529 ifstream infil1(dupFileName);
0530 if (!infil1.is_open()) {
0531 std::cout << "Cannot open duplicate file " << dupFileName << std::endl;
0532 } else {
0533 while (1) {
0534 Long64_t jentry;
0535 infil1 >> jentry;
0536 if (!infil1.good())
0537 break;
0538 entries_.push_back(jentry);
0539 }
0540 infil1.close();
0541 std::cout << "Reads a list of " << entries_.size() << " events from " << dupFileName << std::endl;
0542 }
0543 } else {
0544 std::cout << "No duplicate events in the input file" << std::endl;
0545 }
0546
0547 if (strcmp(comFileName, "") != 0) {
0548 ifstream infil2(comFileName);
0549 if (!infil2.is_open()) {
0550 std::cout << "Cannot open selection file " << comFileName << std::endl;
0551 } else {
0552 while (1) {
0553 int irun, ievt;
0554 infil2 >> irun >> ievt;
0555 if (!infil2.good())
0556 break;
0557 std::pair<int, int> key(irun, ievt);
0558 events_.push_back(key);
0559 }
0560 infil2.close();
0561 std::cout << "Reads a list of " << events_.size() << " events from " << comFileName << std::endl;
0562 }
0563 } else {
0564 std::cout << "No event list provided for selection" << std::endl;
0565 }
0566
0567 if (((flag_ / 100) % 10) > 0) {
0568 if (((flag_ / 100) % 10) == 2) {
0569 fileout_.open(outFileName, std::ofstream::out);
0570 std::cout << "Opens " << outFileName << " in output mode" << std::endl;
0571 } else {
0572 fileout_.open(outFileName, std::ofstream::app);
0573 std::cout << "Opens " << outFileName << " in append mode" << std::endl;
0574 }
0575 fileout_ << "Input file: " << fname_ << " Directory: " << dirnm_ << " Prefix: " << prefix_ << std::endl;
0576 }
0577
0578 double xbins[99];
0579 int nbins(-1);
0580 if (plotType_ == 0) {
0581 double xbin[9] = {-21.0, -16.0, -12.0, -6.0, 0.0, 6.0, 12.0, 16.0, 21.0};
0582 for (int i = 0; i < 9; ++i) {
0583 etas_.push_back(xbin[i]);
0584 xbins[i] = xbin[i];
0585 }
0586 nbins = 8;
0587 } else if (plotType_ == 1) {
0588 double xbin[11] = {-25.0, -20.0, -15.0, -10.0, -5.0, 0.0, 5.0, 10.0, 15.0, 20.0, 25.0};
0589 for (int i = 0; i < 11; ++i) {
0590 etas_.push_back(xbin[i]);
0591 xbins[i] = xbin[i];
0592 }
0593 nbins = 10;
0594 } else if (plotType_ == 2) {
0595 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,
0596 3.0, 5.0, 7.0, 9.0, 11.0, 13.0, 15.0, 17.0, 19.0, 21.0, 23.0};
0597 for (int i = 0; i < 23; ++i) {
0598 etas_.push_back(xbin[i]);
0599 xbins[i] = xbin[i];
0600 }
0601 nbins = 22;
0602 } else {
0603 double xbina[99];
0604 int neta = numb_ / 2;
0605 for (int k = 0; k < neta; ++k) {
0606 xbina[k] = (k - neta) - 0.5;
0607 xbina[numb_ - k] = (neta - k) + 0.5;
0608 }
0609 xbina[neta] = 0;
0610 for (int i = 0; i < numb_ + 1; ++i) {
0611 etas_.push_back(xbina[i]);
0612 xbins[i] = xbina[i];
0613 ++nbins;
0614 }
0615 }
0616 int ipbin[npbin] = {10, 20, 30, 40, 60, 100};
0617 for (unsigned int i = 0; i < npbin; ++i)
0618 ps_.push_back((double)(ipbin[i]));
0619 int npvtx[6] = {0, 7, 10, 13, 16, 100};
0620 for (int i = 0; i < 6; ++i)
0621 nvx_.push_back(npvtx[i]);
0622 double dl1s[9] = {0, 0.10, 0.20, 0.50, 1.0, 2.0, 2.5, 3.0, 10.0};
0623 int ietasL[4] = {0, 13, 17, 0};
0624 int ietasH[4] = {14, 18, 25, 25};
0625 for (int i = 0; i < 4; ++i) {
0626 ietasL_.push_back(ietasL[i]);
0627 ietasH_.push_back(ietasH[i]);
0628 }
0629 int nxbin(100);
0630 double xlow(0.0), xhigh(5.0);
0631 if (coarseBin_ == 1) {
0632 xlow = 0.25;
0633 xhigh = 5.25;
0634 nxbin = 50;
0635 } else if (coarseBin_ > 1) {
0636 if (coarseBin_ == 2)
0637 nxbin = 500;
0638 else
0639 nxbin = 1000;
0640 }
0641
0642 char name[20], title[200];
0643 std::string titl[5] = {
0644 "All tracks", "Good quality tracks", "Selected tracks", "Tracks with charge isolation", "Tracks MIP in ECAL"};
0645 for (int i = 0; i < 9; ++i)
0646 dl1_.push_back(dl1s[i]);
0647 unsigned int kp = (ps_.size() - 1);
0648 for (unsigned int k = 0; k < kp; ++k) {
0649 for (unsigned int j = 0; j <= ietasL_.size(); ++j) {
0650 sprintf(name, "%spp%d%d", prefix_.c_str(), k, j);
0651 if (j == 0)
0652 sprintf(title, "E/p for %s (p = %d:%d GeV All)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0653 else
0654 sprintf(title,
0655 "E/p for %s (p = %d:%d GeV |#eta| %d:%d)",
0656 titl[4].c_str(),
0657 ipbin[k],
0658 ipbin[k + 1],
0659 (ietasL_[j - 1] + 1),
0660 ietasH_[j - 1]);
0661 h_pp[k].push_back(new TH1D(name, title, 100, 10.0, 110.0));
0662 int kk = h_pp[k].size() - 1;
0663 h_pp[k][kk]->Sumw2();
0664 }
0665 }
0666 if (plotType_ <= 1) {
0667 std::cout << "Book Histos for Standard" << std::endl;
0668 for (unsigned int k = 0; k < kp; ++k) {
0669 for (unsigned int i = 0; i < nvx_.size(); ++i) {
0670 sprintf(name, "%setaX%d%d", prefix_.c_str(), k, i);
0671 if (i == 0) {
0672 sprintf(title, "%s (p = %d:%d GeV all vertices)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0673 } else {
0674 sprintf(
0675 title, "%s (p = %d:%d GeV # Vtx %d:%d)", titl[4].c_str(), ipbin[k], ipbin[k + 1], nvx_[i - 1], nvx_[i]);
0676 }
0677 h_etaX[k].push_back(new TProfile(name, title, nbins, xbins));
0678 unsigned int kk = h_etaX[k].size() - 1;
0679 h_etaX[k][kk]->Sumw2();
0680 sprintf(name, "%snvxR%d%d", prefix_.c_str(), k, i);
0681 if (i == 0) {
0682 sprintf(title, "E/p for %s (p = %d:%d GeV all vertices)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0683 } else {
0684 sprintf(title,
0685 "E/p for %s (p = %d:%d GeV # Vtx %d:%d)",
0686 titl[4].c_str(),
0687 ipbin[k],
0688 ipbin[k + 1],
0689 nvx_[i - 1],
0690 nvx_[i]);
0691 }
0692 h_nvxR[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0693 kk = h_nvxR[k].size() - 1;
0694 h_nvxR[k][kk]->Sumw2();
0695 }
0696 for (unsigned int j = 0; j < etas_.size(); ++j) {
0697 sprintf(name, "%sratio%d%d", prefix_.c_str(), k, j);
0698 if (j == 0) {
0699 sprintf(title, "E/p for %s (p = %d:%d GeV)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0700 } else {
0701 sprintf(title,
0702 "E/p for %s (p = %d:%d GeV #eta %4.1f:%4.1f)",
0703 titl[4].c_str(),
0704 ipbin[k],
0705 ipbin[k + 1],
0706 etas_[j - 1],
0707 etas_[j]);
0708 }
0709 h_etaF[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0710 unsigned int kk = h_etaF[k].size() - 1;
0711 h_etaF[k][kk]->Sumw2();
0712 sprintf(name, "%setaR%d%d", prefix_.c_str(), k, j);
0713 h_etaR[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0714 kk = h_etaR[k].size() - 1;
0715 h_etaR[k][kk]->Sumw2();
0716 }
0717 for (unsigned int j = 1; j <= ietasL_.size(); ++j) {
0718 sprintf(name, "%setaB%d%d", prefix_.c_str(), k, j);
0719 sprintf(title,
0720 "E/p for %s (p = %d:%d GeV |#eta| %d:%d)",
0721 titl[4].c_str(),
0722 ipbin[k],
0723 ipbin[k + 1],
0724 (ietasL_[j - 1] + 1),
0725 ietasH_[j - 1]);
0726 h_etaB[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0727 unsigned int kk = h_etaB[k].size() - 1;
0728 h_etaB[k][kk]->Sumw2();
0729 }
0730 for (unsigned int j = 0; j < dl1_.size(); ++j) {
0731 sprintf(name, "%sdl1R%d%d", prefix_.c_str(), k, j);
0732 if (j == 0) {
0733 sprintf(title, "E/p for %s (p = %d:%d GeV All d_{L1})", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0734 } else {
0735 sprintf(title,
0736 "E/p for %s (p = %d:%d GeV d_{L1} %4.2f:%4.2f)",
0737 titl[4].c_str(),
0738 ipbin[k],
0739 ipbin[k + 1],
0740 dl1_[j - 1],
0741 dl1_[j]);
0742 }
0743 h_dL1R[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0744 unsigned int kk = h_dL1R[k].size() - 1;
0745 h_dL1R[k][kk]->Sumw2();
0746 }
0747 }
0748 for (unsigned int i = 0; i < nvx_.size(); ++i) {
0749 sprintf(name, "%setaX%d%d", prefix_.c_str(), kp, i);
0750 if (i == 0) {
0751 sprintf(title, "%s (All Momentum all vertices)", titl[4].c_str());
0752 } else {
0753 sprintf(title, "%s (All Momentum # Vtx %d:%d)", titl[4].c_str(), nvx_[i - 1], nvx_[i]);
0754 }
0755 h_etaX[npbin - 1].push_back(new TProfile(name, title, nbins, xbins));
0756 unsigned int kk = h_etaX[npbin - 1].size() - 1;
0757 h_etaX[npbin - 1][kk]->Sumw2();
0758 sprintf(name, "%snvxR%d%d", prefix_.c_str(), kp, i);
0759 if (i == 0) {
0760 sprintf(title, "E/p for %s (All Momentum all vertices)", titl[4].c_str());
0761 } else {
0762 sprintf(title, "E/p for %s (All Momentum # Vtx %d:%d)", titl[4].c_str(), nvx_[i - 1], nvx_[i]);
0763 }
0764 h_nvxR[npbin - 1].push_back(new TH1D(name, title, 200, 0., 10.));
0765 kk = h_nvxR[npbin - 1].size() - 1;
0766 h_nvxR[npbin - 1][kk]->Sumw2();
0767 }
0768 for (unsigned int j = 0; j < etas_.size(); ++j) {
0769 sprintf(name, "%sratio%d%d", prefix_.c_str(), kp, j);
0770 if (j == 0) {
0771 sprintf(title, "E/p for %s (All momentum)", titl[4].c_str());
0772 } else {
0773 sprintf(title, "E/p for %s (All momentum #eta %4.1f:%4.1f)", titl[4].c_str(), etas_[j - 1], etas_[j]);
0774 }
0775 h_etaF[npbin - 1].push_back(new TH1D(name, title, 200, 0., 10.));
0776 unsigned int kk = h_etaF[npbin - 1].size() - 1;
0777 h_etaF[npbin - 1][kk]->Sumw2();
0778 sprintf(name, "%setaR%d%d", prefix_.c_str(), kp, j);
0779 h_etaR[npbin - 1].push_back(new TH1D(name, title, 200, 0., 10.));
0780 kk = h_etaR[npbin - 1].size() - 1;
0781 h_etaR[npbin - 1][kk]->Sumw2();
0782 }
0783 for (unsigned int j = 1; j <= ietasL_.size(); ++j) {
0784 sprintf(name, "%setaB%d%d", prefix_.c_str(), kp, j);
0785 sprintf(title, "E/p for %s (All momentum |#eta| %d:%d)", titl[4].c_str(), (ietasL_[j - 1] + 1), ietasH_[j - 1]);
0786 h_etaB[npbin - 1].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0787 unsigned int kk = h_etaB[npbin - 1].size() - 1;
0788 h_etaB[npbin - 1][kk]->Sumw2();
0789 }
0790 for (unsigned int j = 0; j < dl1_.size(); ++j) {
0791 sprintf(name, "%sdl1R%d%d", prefix_.c_str(), kp, j);
0792 if (j == 0) {
0793 sprintf(title, "E/p for %s (All momentum)", titl[4].c_str());
0794 } else {
0795 sprintf(title, "E/p for %s (All momentum d_{L1} %4.2f:%4.2f)", titl[4].c_str(), dl1_[j - 1], dl1_[j]);
0796 }
0797 h_dL1R[npbin - 1].push_back(new TH1D(name, title, 200, 0., 10.));
0798 unsigned int kk = h_dL1R[npbin - 1].size() - 1;
0799 h_dL1R[npbin - 1][kk]->Sumw2();
0800 }
0801 } else {
0802 std::cout << "Book Histos for Non-Standard " << etas_.size() << ":" << kp50 << std::endl;
0803 unsigned int kp = (ps_.size() - 1);
0804 for (unsigned int k = 0; k < kp; ++k) {
0805 for (unsigned int j = 0; j < etas_.size(); ++j) {
0806 sprintf(name, "%sratio%d%d", prefix_.c_str(), k, j);
0807 if (j == 0) {
0808 sprintf(title, "E/p for %s (p = %d:%d GeV)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0809 } else {
0810 sprintf(title,
0811 "E/p for %s (p = %d:%d GeV #eta %4.1f:%4.1f)",
0812 titl[4].c_str(),
0813 ipbin[k],
0814 ipbin[k + 1],
0815 etas_[j - 1],
0816 etas_[j]);
0817 }
0818 h_etaF[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0819 unsigned int kk = h_etaF[k].size() - 1;
0820 h_etaF[k][kk]->Sumw2();
0821 }
0822 for (unsigned int j = 1; j <= ietasL_.size(); ++j) {
0823 sprintf(name, "%setaB%d%d", prefix_.c_str(), k, j);
0824 sprintf(title,
0825 "E/p for %s (p = %d:%d GeV |#eta| %d:%d)",
0826 titl[4].c_str(),
0827 ipbin[k],
0828 ipbin[k + 1],
0829 (ietasL_[j - 1] + 1),
0830 ietasH_[j - 1]);
0831 h_etaB[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0832 unsigned int kk = h_etaB[k].size() - 1;
0833 h_etaB[k][kk]->Sumw2();
0834 }
0835 }
0836 }
0837 if (selRBX_) {
0838 for (unsigned int j = 1; j <= 18; ++j) {
0839 sprintf(name, "%sRBX%d%d", prefix_.c_str(), kp50, j);
0840 sprintf(title, "E/p for RBX%d (p = %d:%d GeV |#eta| %d:%d)", j, ipbin[kp50], ipbin[kp50 + 1], etalo_, etahi_);
0841 h_rbx.push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0842 h_rbx[j - 1]->Sumw2();
0843 }
0844 }
0845 }
0846
0847 Bool_t CalibMonitor::Notify() {
0848
0849
0850
0851
0852
0853
0854 return kTRUE;
0855 }
0856
0857 void CalibMonitor::Show(Long64_t entry) {
0858
0859
0860 if (!fChain)
0861 return;
0862 fChain->Show(entry);
0863 }
0864
0865 Int_t CalibMonitor::Cut(Long64_t) {
0866
0867
0868
0869 return 1;
0870 }
0871
0872 void CalibMonitor::Loop() {
0873
0874
0875
0876
0877
0878
0879
0880
0881
0882
0883
0884
0885
0886
0887
0888
0889
0890
0891
0892
0893
0894
0895
0896 const bool debug(false);
0897
0898 if (fChain == 0)
0899 return;
0900 std::map<int, counter> runSum, runEn1, runEn2;
0901
0902
0903 Long64_t nentries = fChain->GetEntriesFast();
0904 std::cout << "Total entries " << nentries << std::endl;
0905 Long64_t nbytes(0), nb(0);
0906 unsigned int duplicate(0), good(0), kount(0);
0907 unsigned int kp1 = ps_.size() - 1;
0908 unsigned int kv1 = 0;
0909 std::vector<int> kounts(kp1, 0);
0910 std::vector<int> kount50(20, 0);
0911 std::vector<int> kount0(20, 0);
0912 std::vector<int> kount1(20, 0);
0913 std::vector<int> kount2(20, 0);
0914 std::vector<int> kount3(20, 0);
0915 std::vector<int> kount4(20, 0);
0916 std::vector<int> kount5(20, 0);
0917 for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0918
0919 Long64_t ientry = LoadTree(jentry);
0920 if (ientry < 0)
0921 break;
0922 nb = fChain->GetEntry(jentry);
0923 nbytes += nb;
0924 if (jentry % 1000000 == 0)
0925 std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
0926 double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
0927 int kp(-1);
0928 for (unsigned int k = 1; k < ps_.size(); ++k) {
0929 if (pmom >= ps_[k - 1] && pmom < ps_[k]) {
0930 kp = k - 1;
0931 break;
0932 }
0933 }
0934 bool p4060 = ((pmom >= 40.0) && (pmom <= 60.0));
0935 if (p4060)
0936 ++kount50[0];
0937 if (kp == 0) {
0938 ++kount0[0];
0939 } else if (kp == 1) {
0940 ++kount1[0];
0941 } else if (kp == 2) {
0942 ++kount2[0];
0943 } else if (kp == 3) {
0944 ++kount3[0];
0945 } else if (kp == 4) {
0946 ++kount4[0];
0947 } else if (kp == 5) {
0948 ++kount5[0];
0949 }
0950 bool select = (std::find(entries_.begin(), entries_.end(), jentry) == entries_.end());
0951 if (!select) {
0952 ++duplicate;
0953 if (debug)
0954 std::cout << "Duplicate event " << t_Run << " " << t_Event << " " << t_p << std::endl;
0955 continue;
0956 }
0957 if (p4060)
0958 ++kount50[1];
0959 if (kp == 0) {
0960 ++kount0[1];
0961 } else if (kp == 1) {
0962 ++kount1[1];
0963 } else if (kp == 2) {
0964 ++kount2[1];
0965 } else if (kp == 3) {
0966 ++kount3[1];
0967 } else if (kp == 4) {
0968 ++kount4[1];
0969 } else if (kp == 5) {
0970 ++kount5[1];
0971 }
0972 bool selRun = (includeRun_ ? ((t_Run >= runlo_) && (t_Run <= runhi_)) : ((t_Run < runlo_) || (t_Run > runhi_)));
0973 if (select) {
0974 if (p4060)
0975 ++kount50[2];
0976 if (kp == 0) {
0977 ++kount0[2];
0978 } else if (kp == 1) {
0979 ++kount1[2];
0980 } else if (kp == 2) {
0981 ++kount2[2];
0982 } else if (kp == 3) {
0983 ++kount3[2];
0984 } else if (kp == 4) {
0985 ++kount4[2];
0986 } else if (kp == 5) {
0987 ++kount5[2];
0988 }
0989 }
0990 select =
0991 (selRun && (fabs(t_ieta) >= etalo_) && (fabs(t_ieta) <= etahi_) && (t_nVtx >= nvxlo_) && (t_nVtx <= nvxhi_));
0992 if (select) {
0993 if (p4060)
0994 ++kount50[3];
0995 if (kp == 0) {
0996 ++kount0[3];
0997 } else if (kp == 1) {
0998 ++kount1[3];
0999 } else if (kp == 2) {
1000 ++kount2[3];
1001 } else if (kp == 3) {
1002 ++kount3[3];
1003 } else if (kp == 4) {
1004 ++kount4[3];
1005 } else if (kp == 5) {
1006 ++kount5[3];
1007 }
1008 }
1009 if (!select) {
1010 if (debug)
1011 std::cout << "Run # " << t_Run << " out of range of " << runlo_ << ":" << runhi_ << " or ieta " << t_ieta
1012 << " (" << etalo_ << ":" << etahi_ << ") or nvtx " << t_nVtx << " (" << nvxlo_ << ":" << nvxhi_
1013 << ") out of range" << std::endl;
1014 continue;
1015 }
1016 if (cSelect_ != nullptr) {
1017 if (exclude_) {
1018 if (cSelect_->isItRBX(t_DetIds))
1019 continue;
1020 } else {
1021 if (!(cSelect_->isItRBX(t_ieta, t_iphi)))
1022 continue;
1023 }
1024 }
1025 if (p4060)
1026 ++kount50[4];
1027 if (kp == 0) {
1028 ++kount0[4];
1029 } else if (kp == 1) {
1030 ++kount1[4];
1031 } else if (kp == 2) {
1032 ++kount2[4];
1033 } else if (kp == 3) {
1034 ++kount3[4];
1035 } else if (kp == 4) {
1036 ++kount4[4];
1037 } else if (kp == 5) {
1038 ++kount5[4];
1039 }
1040 select = (!cutL1T_ || (t_mindR1 >= 0.5));
1041 if (!select) {
1042 if (debug)
1043 std::cout << "Reject Run # " << t_Run << " Event # " << t_Event << " too close to L1 trigger " << t_mindR1
1044 << std::endl;
1045 continue;
1046 }
1047 if (p4060)
1048 ++kount50[5];
1049 if (kp == 0) {
1050 ++kount0[5];
1051 } else if (kp == 1) {
1052 ++kount1[5];
1053 } else if (kp == 2) {
1054 ++kount2[5];
1055 } else if (kp == 3) {
1056 ++kount3[5];
1057 } else if (kp == 4) {
1058 ++kount4[5];
1059 } else if (kp == 5) {
1060 ++kount5[5];
1061 }
1062 select = ((events_.size() == 0) ||
1063 (std::find(events_.begin(), events_.end(), std::pair<int, int>(t_Run, t_Event)) != events_.end()));
1064 if (!select) {
1065 if (debug)
1066 std::cout << "Reject Run # " << t_Run << " Event # " << t_Event << " not in the selection list" << std::endl;
1067 continue;
1068 }
1069 if (p4060)
1070 ++kount50[6];
1071 if (kp == 0) {
1072 ++kount0[6];
1073 } else if (kp == 1) {
1074 ++kount1[6];
1075 } else if (kp == 2) {
1076 ++kount2[6];
1077 } else if (kp == 3) {
1078 ++kount3[6];
1079 } else if (kp == 4) {
1080 ++kount4[6];
1081 } else if (kp == 5) {
1082 ++kount5[6];
1083 }
1084 if ((ifDepth_ == 3) && (cFactor_ != nullptr) && (cFactor_->absent(ientry)))
1085 continue;
1086
1087
1088 int jp(-1), jp1(-1), jp2(-1);
1089 unsigned int kv = nvx_.size() - 1;
1090 for (unsigned int k = 1; k < nvx_.size(); ++k) {
1091 if (t_goodPV >= nvx_[k - 1] && t_goodPV < nvx_[k]) {
1092 kv = k;
1093 break;
1094 }
1095 }
1096 unsigned int kd1 = 0;
1097 unsigned int kd = dl1_.size() - 1;
1098 for (unsigned int k = 1; k < dl1_.size(); ++k) {
1099 if (t_mindR1 >= dl1_[k - 1] && t_mindR1 < dl1_[k]) {
1100 kd = k;
1101 break;
1102 }
1103 }
1104 double eta = (t_ieta > 0) ? ((double)(t_ieta)-0.001) : ((double)(t_ieta) + 0.001);
1105 for (unsigned int j = 1; j < etas_.size(); ++j) {
1106 if (eta > etas_[j - 1] && eta < etas_[j]) {
1107 jp = j;
1108 break;
1109 }
1110 }
1111 for (unsigned int j = 0; j < (ietasL_.size() - 1); ++j) {
1112 if (std::abs(t_ieta) > ietasL_[j] && std::abs(t_ieta) <= ietasH_[j]) {
1113 if (jp1 >= 0)
1114 jp2 = j;
1115 if (jp1 < 0)
1116 jp1 = j;
1117 }
1118 }
1119 int jp3 = (jp1 >= 0) ? (int)(ietasL_.size() - 1) : jp1;
1120 if (debug)
1121 std::cout << "Bin " << kp << ":" << kp1 << ":" << kv << ":" << kv1 << ":" << kd << ":" << kd1 << ":" << jp << ":"
1122 << jp1 << ":" << jp2 << ":" << jp3 << " pmom:ieta:pv:mindR " << pmom << ":" << std::abs(t_ieta) << ":"
1123 << t_goodPV << ":" << t_mindR1 << std::endl;
1124 double cut = (pmom > 20) ? ((flexibleSelect_ == 0) ? 2.0 : 10.0) : 0.0;
1125
1126 double rcut(-1000.0);
1127
1128
1129 double rat(1.0), eHcal(t_eHcal);
1130 if (corrFactor_->doCorr() || (cFactor_ != nullptr)) {
1131 eHcal = 0;
1132 for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
1133
1134 double cfac(1.0);
1135 if (corrFactor_->doCorr()) {
1136 unsigned int id = truncateId((*t_DetIds)[k], truncateFlag_, false);
1137 cfac = corrFactor_->getCorr(id);
1138 }
1139 if ((cFactor_ != nullptr) && (ifDepth_ != 3))
1140 cfac *= cFactor_->getCorr(t_Run, (*t_DetIds)[k]);
1141 eHcal += (cfac * ((*t_HitEnergies)[k]));
1142 if (debug) {
1143 int subdet, zside, ieta, iphi, depth;
1144 unpackDetId((*t_DetIds)[k], subdet, zside, ieta, iphi, depth);
1145 std::cout << zside << ":" << ieta << ":" << depth << " Corr " << cfac << " " << (*t_HitEnergies)[k] << " Out "
1146 << eHcal << std::endl;
1147 }
1148 }
1149 }
1150 bool goodTk = goodTrack(eHcal, cut, jentry, debug);
1151 bool selPhi = selectPhi(debug);
1152 if (t_qltyFlag) {
1153 if (p4060)
1154 ++kount50[7];
1155 if (kp == 0) {
1156 ++kount0[7];
1157 } else if (kp == 1) {
1158 ++kount1[7];
1159 } else if (kp == 2) {
1160 ++kount2[7];
1161 } else if (kp == 3) {
1162 ++kount3[7];
1163 } else if (kp == 4) {
1164 ++kount4[7];
1165 } else if (kp == 5) {
1166 ++kount5[7];
1167 }
1168 if (t_selectTk) {
1169 if (p4060)
1170 ++kount50[8];
1171 if (kp == 0) {
1172 ++kount0[8];
1173 } else if (kp == 1) {
1174 ++kount1[8];
1175 } else if (kp == 2) {
1176 ++kount2[8];
1177 } else if (kp == 3) {
1178 ++kount3[8];
1179 } else if (kp == 4) {
1180 ++kount4[8];
1181 } else if (kp == 5) {
1182 ++kount5[8];
1183 }
1184 if (t_hmaxNearP < cut) {
1185 if (p4060)
1186 ++kount50[9];
1187 if (kp == 0) {
1188 ++kount0[9];
1189 } else if (kp == 1) {
1190 ++kount1[9];
1191 } else if (kp == 2) {
1192 ++kount2[9];
1193 } else if (kp == 3) {
1194 ++kount3[9];
1195 } else if (kp == 4) {
1196 ++kount4[9];
1197 } else if (kp == 5) {
1198 ++kount5[9];
1199 }
1200 if (t_eMipDR < 1.0) {
1201 if (p4060)
1202 ++kount50[10];
1203 if (kp == 0) {
1204 ++kount0[10];
1205 } else if (kp == 1) {
1206 ++kount1[10];
1207 } else if (kp == 2) {
1208 ++kount2[10];
1209 } else if (kp == 3) {
1210 ++kount3[10];
1211 } else if (kp == 4) {
1212 ++kount4[10];
1213 } else if (kp == 5) {
1214 ++kount5[10];
1215 }
1216 if (eHcal > 0.001) {
1217 if (p4060)
1218 ++kount50[11];
1219 if (kp == 0) {
1220 ++kount0[11];
1221 } else if (kp == 1) {
1222 ++kount1[11];
1223 } else if (kp == 2) {
1224 ++kount2[11];
1225 } else if (kp == 3) {
1226 ++kount3[11];
1227 } else if (kp == 4) {
1228 ++kount4[11];
1229 } else if (kp == 5) {
1230 ++kount5[11];
1231 }
1232 if (selPhi) {
1233 if (p4060)
1234 ++kount50[12];
1235 if (kp == 0) {
1236 ++kount0[12];
1237 } else if (kp == 1) {
1238 ++kount1[12];
1239 } else if (kp == 2) {
1240 ++kount2[12];
1241 } else if (kp == 3) {
1242 ++kount3[12];
1243 } else if (kp == 4) {
1244 ++kount4[12];
1245 } else if (kp == 5) {
1246 ++kount5[12];
1247 }
1248 }
1249 }
1250 }
1251 }
1252 }
1253 }
1254 if (pmom > 0)
1255 rat = (eHcal / (pmom - t_eMipDR));
1256 if (debug) {
1257 std::cout << "Entry " << jentry << " p|eHcal|ratio " << pmom << "|" << t_eHcal << "|" << eHcal << "|" << rat
1258 << "|" << kp << "|" << kv << "|" << jp << " Cuts " << t_qltyFlag << "|" << t_selectTk << "|"
1259 << (t_hmaxNearP < cut) << "|" << (t_eMipDR < 1.0) << "|" << goodTk << "|" << (rat > rcut)
1260 << " Select Phi " << selPhi << std::endl;
1261 std::cout << "D1 : " << kp << ":" << kp1 << ":" << kv << ":" << kv1 << ":" << kd << ":" << kd1 << ":" << jp
1262 << std::endl;
1263 }
1264 if (goodTk && (kp >= 0) && selPhi) {
1265 if (p4060)
1266 ++kount50[13];
1267 if (kp == 0) {
1268 ++kount0[13];
1269 } else if (kp == 1) {
1270 ++kount1[13];
1271 } else if (kp == 2) {
1272 ++kount2[13];
1273 } else if (kp == 3) {
1274 ++kount3[13];
1275 } else if (kp == 4) {
1276 ++kount4[13];
1277 } else if (kp == 5) {
1278 ++kount5[13];
1279 }
1280 if (t_eHcal < 0.01) {
1281 std::map<int, counter>::const_iterator itr = runEn1.find(t_Run);
1282 if (itr == runEn1.end()) {
1283 counter knt;
1284 if ((kp >= 0) && (kp < counter::npsize))
1285 knt.count[kp] = 1;
1286 knt.total = 1;
1287 runEn1[t_Run] = knt;
1288 } else {
1289 counter knt = runEn1[t_Run];
1290 if ((kp >= 0) && (kp < counter::npsize))
1291 ++knt.count[kp];
1292 ++knt.total;
1293 runEn1[t_Run] = knt;
1294 }
1295 }
1296 if (t_eMipDR < 0.01 && t_eHcal < 0.01) {
1297 if (p4060)
1298 ++kount50[14];
1299 if (kp == 0) {
1300 ++kount0[14];
1301 } else if (kp == 1) {
1302 ++kount1[14];
1303 } else if (kp == 2) {
1304 ++kount2[14];
1305 } else if (kp == 3) {
1306 ++kount3[14];
1307 } else if (kp == 4) {
1308 ++kount4[14];
1309 } else if (kp == 5) {
1310 ++kount5[14];
1311 }
1312 std::map<int, counter>::const_iterator itr = runEn2.find(t_Run);
1313 if (itr == runEn2.end()) {
1314 counter knt;
1315 if ((kp >= 0) && (kp < counter::npsize))
1316 knt.count[kp] = 1;
1317 knt.total = 1;
1318 runEn2[t_Run] = knt;
1319 } else {
1320 counter knt = runEn2[t_Run];
1321 if ((kp >= 0) && (kp < counter::npsize))
1322 ++knt.count[kp];
1323 ++knt.total;
1324 runEn2[t_Run] = knt;
1325 }
1326 }
1327 if (rat > rcut) {
1328 if (p4060)
1329 ++kount50[15];
1330 if (kp == 0) {
1331 ++kount0[15];
1332 } else if (kp == 1) {
1333 ++kount1[15];
1334 } else if (kp == 2) {
1335 ++kount2[15];
1336 } else if (kp == 3) {
1337 ++kount3[15];
1338 } else if (kp == 4) {
1339 ++kount4[15];
1340 } else if (kp == 5) {
1341 ++kount5[15];
1342 }
1343 if (plotType_ <= 1) {
1344 h_etaX[kp][kv]->Fill(eta, rat, t_EventWeight);
1345 h_etaX[kp][kv1]->Fill(eta, rat, t_EventWeight);
1346 h_nvxR[kp][kv]->Fill(rat, t_EventWeight);
1347 h_nvxR[kp][kv1]->Fill(rat, t_EventWeight);
1348 h_dL1R[kp][kd]->Fill(rat, t_EventWeight);
1349 h_dL1R[kp][kd1]->Fill(rat, t_EventWeight);
1350 if (jp > 0)
1351 h_etaR[kp][jp]->Fill(rat, t_EventWeight);
1352 h_etaR[kp][0]->Fill(rat, t_EventWeight);
1353 }
1354 h_pp[kp][0]->Fill(pmom, t_EventWeight);
1355 if (jp1 >= 0) {
1356 h_pp[kp][jp1 + 1]->Fill(pmom, t_EventWeight);
1357 h_pp[kp][jp3 + 1]->Fill(pmom, t_EventWeight);
1358 }
1359 if (jp2 >= 0)
1360 h_pp[kp][jp2 + 1]->Fill(pmom, t_EventWeight);
1361 if (kp == (int)(kp50)) {
1362 std::map<int, counter>::const_iterator itr = runSum.find(t_Run);
1363 if (itr == runSum.end()) {
1364 counter knt;
1365 if ((kp >= 0) && (kp < counter::npsize))
1366 knt.count[kp] = 1;
1367 knt.total = 1;
1368 runSum[t_Run] = knt;
1369 } else {
1370 counter knt = runSum[t_Run];
1371 if ((kp >= 0) && (kp < counter::npsize))
1372 ++knt.count[kp];
1373 ++knt.total;
1374 runSum[t_Run] = knt;
1375 }
1376 }
1377 if ((!dataMC_) || (t_mindR1 > 0.5) || (t_DataType == 1)) {
1378 if (p4060)
1379 ++kount50[16];
1380 if (kp == 0) {
1381 ++kount0[16];
1382 } else if (kp == 1) {
1383 ++kount1[16];
1384 } else if (kp == 2) {
1385 ++kount2[16];
1386 } else if (kp == 3) {
1387 ++kount3[16];
1388 } else if (kp == 4) {
1389 ++kount4[16];
1390 } else if (kp == 5) {
1391 ++kount5[16];
1392 }
1393 ++kounts[kp];
1394 if (plotType_ <= 1) {
1395 if (jp > 0)
1396 h_etaF[kp][jp]->Fill(rat, t_EventWeight);
1397 h_etaF[kp][0]->Fill(rat, t_EventWeight);
1398 } else {
1399 if (debug) {
1400 std::cout << "kp " << kp << h_etaF[kp].size() << std::endl;
1401 }
1402 if (jp > 0)
1403 h_etaF[kp][jp]->Fill(rat, t_EventWeight);
1404 h_etaF[kp][0]->Fill(rat, t_EventWeight);
1405 if (jp1 >= 0) {
1406 h_etaB[kp][jp1]->Fill(rat, t_EventWeight);
1407 h_etaB[kp][jp3]->Fill(rat, t_EventWeight);
1408 }
1409 if (jp2 >= 0)
1410 h_etaB[kp][jp2]->Fill(rat, t_EventWeight);
1411 }
1412 if (selRBX_ && (kp == (int)(kp50)) && ((t_ieta * zside_) > 0)) {
1413 int rbx = (t_iphi > 70) ? 0 : (t_iphi + 1) / 4;
1414 h_rbx[rbx]->Fill(rat, t_EventWeight);
1415 }
1416 }
1417 if (pmom > 10.0) {
1418 if (plotType_ <= 1) {
1419 h_etaX[kp1][kv]->Fill(eta, rat, t_EventWeight);
1420 h_etaX[kp1][kv1]->Fill(eta, rat, t_EventWeight);
1421 h_nvxR[kp1][kv]->Fill(rat, t_EventWeight);
1422 h_nvxR[kp1][kv1]->Fill(rat, t_EventWeight);
1423 h_dL1R[kp1][kd]->Fill(rat, t_EventWeight);
1424 h_dL1R[kp1][kd1]->Fill(rat, t_EventWeight);
1425 if (jp > 0)
1426 h_etaR[kp1][jp]->Fill(rat, t_EventWeight);
1427 h_etaR[kp1][0]->Fill(rat, t_EventWeight);
1428 if (jp1 >= 0) {
1429 h_etaB[kp][jp1]->Fill(rat, t_EventWeight);
1430 h_etaB[kp][jp3]->Fill(rat, t_EventWeight);
1431 }
1432 if (jp2 >= 0)
1433 h_etaB[kp][jp2]->Fill(rat, t_EventWeight);
1434 }
1435 if (p4060)
1436 ++kount50[17];
1437 if (kp == 0) {
1438 ++kount0[17];
1439 } else if (kp == 1) {
1440 ++kount1[17];
1441 } else if (kp == 2) {
1442 ++kount2[17];
1443 } else if (kp == 3) {
1444 ++kount3[17];
1445 } else if (kp == 4) {
1446 ++kount4[17];
1447 } else if (kp == 5) {
1448 ++kount5[17];
1449 }
1450 }
1451 }
1452 }
1453 if (pmom > 10.0) {
1454 kount++;
1455 if (((flag_ / 100) % 10) != 0) {
1456 good++;
1457 fileout_ << good << " " << jentry << " " << t_Run << " " << t_Event << " " << t_ieta << " " << pmom
1458 << std::endl;
1459 }
1460 }
1461 }
1462 unsigned int k(0);
1463 std::cout << "\nSummary of entries with " << runSum.size() << " runs\n";
1464 for (std::map<int, counter>::iterator itr = runSum.begin(); itr != runSum.end(); ++itr, ++k)
1465 std::cout << "[" << k << "] Run " << itr->first << " Total " << (itr->second).total << " in p-bins "
1466 << (itr->second).count[0] << ":" << (itr->second).count[1] << ":" << (itr->second).count[2] << ":"
1467 << (itr->second).count[3] << std::endl;
1468 k = 0;
1469 std::cout << runEn1.size() << " runs with 0 energy in HCAL\n";
1470 for (std::map<int, counter>::iterator itr = runEn1.begin(); itr != runEn1.end(); ++itr, ++k)
1471 std::cout << "[" << k << "] Run " << itr->first << " Total " << (itr->second).total << " in p-bins "
1472 << (itr->second).count[0] << ":" << (itr->second).count[1] << ":" << (itr->second).count[2] << ":"
1473 << (itr->second).count[3] << std::endl;
1474 k = 0;
1475 std::cout << runEn2.size() << " runs with 0 energy in ECAL and HCAL\n";
1476 for (std::map<int, counter>::iterator itr = runEn2.begin(); itr != runEn2.end(); ++itr, ++k)
1477 std::cout << "[" << k << "] Run " << itr->first << " Total " << (itr->second).total << " in p-bins "
1478 << (itr->second).count[0] << ":" << (itr->second).count[1] << ":" << (itr->second).count[2] << ":"
1479 << (itr->second).count[3] << std::endl;
1480 if (((flag_ / 100) % 10) > 0) {
1481 fileout_.close();
1482 std::cout << "Writes " << good << " events in the file " << outFileName_ << std::endl;
1483 }
1484 std::cout << "Finds " << duplicate << " Duplicate events out of " << kount << " events in this file with p>10 Gev"
1485 << std::endl;
1486 std::cout << "Number of selected events:" << std::endl;
1487 for (unsigned int k = 1; k < ps_.size(); ++k)
1488 std::cout << ps_[k - 1] << ":" << ps_[k] << " " << kounts[k - 1] << std::endl;
1489 std::cout << "Number in each step for tracks of momentum 40-60 GeV: ";
1490 for (unsigned int k = 0; k < 18; ++k)
1491 std::cout << " [" << k << "] " << kount50[k];
1492 std::cout << std::endl;
1493 for (unsigned int k = 1; k < ps_.size(); ++k) {
1494 std::cout << "Number in each step for tracks of momentum " << ps_[k - 1] << "-" << ps_[k] << " Gev: ";
1495 for (unsigned int k1 = 0; k1 < 18; ++k1) {
1496 if (k == 1) {
1497 std::cout << " [" << k1 << "] " << kount0[k1];
1498 } else if (k == 2) {
1499 std::cout << " [" << k1 << "] " << kount1[k1];
1500 } else if (k == 3) {
1501 std::cout << " [" << k1 << "] " << kount2[k1];
1502 } else if (k == 4) {
1503 std::cout << " [" << k1 << "] " << kount3[k1];
1504 } else if (k == 5) {
1505 std::cout << " [" << k1 << "] " << kount4[k1];
1506 } else if (k == 6) {
1507 std::cout << " [" << k1 << "] " << kount5[k1];
1508 }
1509 }
1510 std::cout << std::endl;
1511 }
1512 }
1513
1514 bool CalibMonitor::goodTrack(double &eHcal, double &cuti, const Long64_t &entry, bool debug) {
1515 bool select(true);
1516 double cut(cuti);
1517 if (debug) {
1518 std::cout << "goodTrack input " << eHcal << ":" << cut;
1519 }
1520 if (flexibleSelect_ > 1) {
1521 double eta = (t_ieta > 0) ? t_ieta : -t_ieta;
1522 cut = 8.0 * exp(eta * log2by18_);
1523 }
1524 correctEnergy(eHcal, entry);
1525 select = ((t_qltyFlag) && (t_selectTk) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (eHcal > 0.001));
1526 if (debug) {
1527 std::cout << " output " << select << " Based on " << t_qltyFlag << ":" << t_selectTk << ":" << t_hmaxNearP << ":"
1528 << cut << ":" << t_eMipDR << ":" << eHcal << std::endl;
1529 }
1530 return select;
1531 }
1532
1533 bool CalibMonitor::selectPhi(bool debug) {
1534 bool select(true);
1535 if (phimin_ > 1 || phimax_ < 72) {
1536 double eTotal(0), eSelec(0);
1537
1538 for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
1539 int iphi = ((*t_DetIds)[k]) & (0x3FF);
1540 int zside = ((*t_DetIds)[k] & 0x80000) ? (1) : (-1);
1541 eTotal += ((*t_HitEnergies)[k]);
1542 if (iphi >= phimin_ && iphi <= phimax_ && zside == zside_)
1543 eSelec += ((*t_HitEnergies)[k]);
1544 }
1545 if (eSelec < 0.9 * eTotal)
1546 select = false;
1547 if (debug) {
1548 std::cout << "Etotal " << eTotal << " and ESelec " << eSelec << " (phi " << phimin_ << ":" << phimax_ << " z "
1549 << zside_ << ") Selection " << select << std::endl;
1550 }
1551 }
1552 return select;
1553 }
1554
1555 void CalibMonitor::plotHist(int itype, int inum, bool save) {
1556 gStyle->SetCanvasBorderMode(0);
1557 gStyle->SetCanvasColor(kWhite);
1558 gStyle->SetPadColor(kWhite);
1559 gStyle->SetFillColor(kWhite);
1560 gStyle->SetOptStat(111110);
1561 gStyle->SetOptFit(1);
1562 char name[100];
1563 int itmin = (itype >= 0 && itype < 4) ? itype : 0;
1564 int itmax = (itype >= 0 && itype < 4) ? itype : 3;
1565 std::string types[4] = {
1566 "E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})"};
1567 int nmax[4] = {npbin, npbin, npbin, npbin};
1568 for (int type = itmin; type <= itmax; ++type) {
1569 int inmin = (inum >= 0 && inum < nmax[type]) ? inum : 0;
1570 int inmax = (inum >= 0 && inum < nmax[type]) ? inum : nmax[type] - 1;
1571 int kmax = 1;
1572 if (type == 0)
1573 kmax = (int)(etas_.size());
1574 else if (type == 1)
1575 kmax = (int)(etas_.size());
1576 else if (type == 2)
1577 kmax = (int)(nvx_.size());
1578 else
1579 kmax = (int)(dl1_.size());
1580 for (int num = inmin; num <= inmax; ++num) {
1581 for (int k = 0; k < kmax; ++k) {
1582 sprintf(name, "c_%d%d%d", type, num, k);
1583 TCanvas *pad = new TCanvas(name, name, 700, 500);
1584 pad->SetRightMargin(0.10);
1585 pad->SetTopMargin(0.10);
1586 sprintf(name, "%s", types[type].c_str());
1587 if (type != 7) {
1588 TH1D *hist(0);
1589 if (type == 0)
1590 hist = (TH1D *)(h_etaR[num][k]->Clone());
1591 else if (type == 1)
1592 hist = (TH1D *)(h_etaF[num][k]->Clone());
1593 else if (type == 2)
1594 hist = (TH1D *)(h_nvxR[num][k]->Clone());
1595 else
1596 hist = (TH1D *)(h_dL1R[num][k]->Clone());
1597 hist->GetXaxis()->SetTitle(name);
1598 hist->GetYaxis()->SetTitle("Tracks");
1599 drawHist(hist, pad);
1600 if (save) {
1601 sprintf(name, "c_%s%d%d%d.gif", prefix_.c_str(), type, num, k);
1602 pad->Print(name);
1603 }
1604 } else {
1605 TProfile *hist = (TProfile *)(h_etaX[num][k]->Clone());
1606 hist->GetXaxis()->SetTitle(name);
1607 hist->GetYaxis()->SetTitle("<E_{HCAL}/(p-E_{ECAL})>");
1608 hist->GetYaxis()->SetRangeUser(0.4, 1.6);
1609 hist->Fit("pol0", "q");
1610 drawHist(hist, pad);
1611 if (save) {
1612 sprintf(name, "c_%s%d%d%d.gif", prefix_.c_str(), type, num, k);
1613 pad->Print(name);
1614 }
1615 }
1616 }
1617 }
1618 }
1619 }
1620
1621 template <class Hist>
1622 void CalibMonitor::drawHist(Hist *hist, TCanvas *pad) {
1623 hist->GetYaxis()->SetLabelOffset(0.005);
1624 hist->GetYaxis()->SetTitleOffset(1.20);
1625 hist->Draw();
1626 pad->Update();
1627 TPaveStats *st1 = (TPaveStats *)hist->GetListOfFunctions()->FindObject("stats");
1628 if (st1 != NULL) {
1629 st1->SetY1NDC(0.70);
1630 st1->SetY2NDC(0.90);
1631 st1->SetX1NDC(0.55);
1632 st1->SetX2NDC(0.90);
1633 }
1634 pad->Modified();
1635 pad->Update();
1636 }
1637
1638 void CalibMonitor::savePlot(const std::string &theName, bool append, bool all) {
1639 TFile *theFile(0);
1640 if (append) {
1641 theFile = new TFile(theName.c_str(), "UPDATE");
1642 } else {
1643 theFile = new TFile(theName.c_str(), "RECREATE");
1644 }
1645
1646 theFile->cd();
1647 for (unsigned int k = 0; k < ps_.size(); ++k) {
1648 for (unsigned int j = 0; j <= ietasL_.size(); ++j) {
1649 if ((all || k == kp50) && h_pp[k].size() > j && h_pp[k][j] != 0) {
1650 TH1D *hist = (TH1D *)h_pp[k][j]->Clone();
1651 hist->Write();
1652 }
1653 }
1654 if (plotType_ <= 1) {
1655 for (unsigned int i = 0; i < nvx_.size(); ++i) {
1656 if (h_etaX[k][i] != 0) {
1657 TProfile *hnew = (TProfile *)h_etaX[k][i]->Clone();
1658 hnew->Write();
1659 }
1660 if (h_nvxR[k].size() > i && h_nvxR[k][i] != 0) {
1661 TH1D *hist = (TH1D *)h_nvxR[k][i]->Clone();
1662 hist->Write();
1663 }
1664 }
1665 }
1666 for (unsigned int j = 0; j < etas_.size(); ++j) {
1667 if ((plotType_ <= 1) && (h_etaR[k][j] != 0)) {
1668 TH1D *hist = (TH1D *)h_etaR[k][j]->Clone();
1669 hist->Write();
1670 }
1671 if (h_etaF[k].size() > j && h_etaF[k][j] != 0 && (all || (k == kp50))) {
1672 TH1D *hist = (TH1D *)h_etaF[k][j]->Clone();
1673 hist->Write();
1674 }
1675 }
1676 if (plotType_ <= 1) {
1677 for (unsigned int j = 0; j < dl1_.size(); ++j) {
1678 if (h_dL1R[k][j] != 0) {
1679 TH1D *hist = (TH1D *)h_dL1R[k][j]->Clone();
1680 hist->Write();
1681 }
1682 }
1683 }
1684 for (unsigned int j = 0; j < ietasL_.size(); ++j) {
1685 if (h_etaB[k].size() > j && h_etaB[k][j] != 0 && (all || (k == kp50))) {
1686 TH1D *hist = (TH1D *)h_etaB[k][j]->Clone();
1687 hist->Write();
1688 }
1689 }
1690 }
1691 if (selRBX_) {
1692 for (unsigned int k = 0; k < 18; ++k) {
1693 if (h_rbx[k] != 0) {
1694 TH1D *h1 = (TH1D *)h_rbx[k]->Clone();
1695 h1->Write();
1696 }
1697 }
1698 }
1699 std::cout << "All done" << std::endl;
1700 theFile->Close();
1701 }
1702
1703 void CalibMonitor::correctEnergy(double &eHcal, const Long64_t &entry) {
1704 bool debug(false);
1705 double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
1706 if ((ifDepth_ == 3) && (cFactor_ != nullptr)) {
1707 double cfac = cFactor_->getCorr(entry);
1708 eHcal *= cfac;
1709 if (debug)
1710 std::cout << "PU Factor for " << ifDepth_ << ":" << entry << " = " << cfac << ":" << eHcal << std::endl;
1711 } else if ((corrPU_ < 0) && (pmom > 0)) {
1712 double ediff = (t_eHcal30 - t_eHcal10);
1713 if (t_DetIds1 != 0 && t_DetIds3 != 0) {
1714 double Etot1(0), Etot3(0);
1715
1716 for (unsigned int idet = 0; idet < (*t_DetIds1).size(); idet++) {
1717 unsigned int id = truncateId((*t_DetIds1)[idet], truncateFlag_, false);
1718 double cfac = corrFactor_->getCorr(id);
1719 if (cFactor_ != 0)
1720 cfac *= cFactor_->getCorr(t_Run, (*t_DetIds1)[idet]);
1721 double hitEn = cfac * (*t_HitEnergies1)[idet];
1722 Etot1 += hitEn;
1723 }
1724 for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
1725 unsigned int id = truncateId((*t_DetIds3)[idet], truncateFlag_, false);
1726 double cfac = corrFactor_->getCorr(id);
1727 if (cFactor_ != 0)
1728 cfac *= cFactor_->getCorr(t_Run, (*t_DetIds3)[idet]);
1729 double hitEn = cfac * (*t_HitEnergies3)[idet];
1730 Etot3 += hitEn;
1731 }
1732 ediff = (Etot3 - Etot1);
1733 }
1734 double fac = puFactor(-corrPU_, t_ieta, pmom, eHcal, ediff, false);
1735 if (debug) {
1736 double fac1 = puFactor(-corrPU_, t_ieta, pmom, eHcal, ediff, true);
1737 double fac2 = puFactor(2, t_ieta, pmom, eHcal, ediff, true);
1738 std::cout << "PU Factor for " << -corrPU_ << " = " << fac1 << "; for 2 = " << fac2 << std::endl;
1739 }
1740 eHcal *= fac;
1741 } else if (corrPU_ > 0) {
1742 eHcal = puFactorRho(corrPU_, t_ieta, t_rhoh, eHcal);
1743 }
1744 }
1745
1746 class GetEntries {
1747 public:
1748 TTree *fChain;
1749 Int_t fCurrent;
1750
1751
1752 UInt_t t_RunNo;
1753 UInt_t t_EventNo;
1754 Int_t t_Tracks;
1755 Int_t t_TracksProp;
1756 Int_t t_TracksSaved;
1757 Int_t t_TracksLoose;
1758 Int_t t_TracksTight;
1759 Int_t t_allvertex;
1760 Bool_t t_TrigPass;
1761 Bool_t t_TrigPassSel;
1762 Bool_t t_L1Bit;
1763 std::vector<Bool_t> *t_hltbits;
1764 std::vector<int> *t_ietaAll;
1765 std::vector<int> *t_ietaGood;
1766 std::vector<int> *t_trackType;
1767
1768
1769 TBranch *b_t_RunNo;
1770 TBranch *b_t_EventNo;
1771 TBranch *b_t_Tracks;
1772 TBranch *b_t_TracksProp;
1773 TBranch *b_t_TracksSaved;
1774 TBranch *b_t_TracksLoose;
1775 TBranch *b_t_TracksTight;
1776 TBranch *b_t_allvertex;
1777 TBranch *b_t_TrigPass;
1778 TBranch *b_t_TrigPassSel;
1779 TBranch *b_t_L1Bit;
1780 TBranch *b_t_hltbits;
1781 TBranch *b_t_ietaAll;
1782 TBranch *b_t_ietaGood;
1783 TBranch *b_t_trackType;
1784
1785 GetEntries(const std::string &fname,
1786 const std::string &dirname,
1787 const char *dupFileName,
1788 const unsigned int bit1,
1789 const unsigned int bit2);
1790 virtual ~GetEntries();
1791 virtual Int_t Cut(Long64_t entry);
1792 virtual Int_t GetEntry(Long64_t entry);
1793 virtual Long64_t LoadTree(Long64_t entry);
1794 virtual void Init(TTree *tree, const char *dupFileName);
1795 virtual void Loop();
1796 virtual Bool_t Notify();
1797 virtual void Show(Long64_t entry = -1);
1798
1799 private:
1800 unsigned int bit_[2];
1801 std::vector<Long64_t> entries_;
1802 TH1I *h_tk[3], *h_eta[4], *h_pvx[3];
1803 TH1D *h_eff[3];
1804 };
1805
1806 GetEntries::GetEntries(const std::string &fname,
1807 const std::string &dirnm,
1808 const char *dupFileName,
1809 const unsigned int bit1,
1810 const unsigned int bit2) {
1811 TFile *file = new TFile(fname.c_str());
1812 TDirectory *dir = (TDirectory *)file->FindObjectAny(dirnm.c_str());
1813 std::cout << fname << " file " << file << " " << dirnm << " " << dir << std::endl;
1814 TTree *tree = (TTree *)dir->Get("EventInfo");
1815 std::cout << "CalibTree " << tree << std::endl;
1816 bit_[0] = bit1;
1817 bit_[1] = bit2;
1818 Init(tree, dupFileName);
1819 }
1820
1821 GetEntries::~GetEntries() {
1822 if (!fChain)
1823 return;
1824 delete fChain->GetCurrentFile();
1825 }
1826
1827 Int_t GetEntries::GetEntry(Long64_t entry) {
1828
1829 if (!fChain)
1830 return 0;
1831 return fChain->GetEntry(entry);
1832 }
1833
1834 Long64_t GetEntries::LoadTree(Long64_t entry) {
1835
1836 if (!fChain)
1837 return -5;
1838 Long64_t centry = fChain->LoadTree(entry);
1839 if (centry < 0)
1840 return centry;
1841 if (!fChain->InheritsFrom(TChain::Class()))
1842 return centry;
1843 TChain *chain = (TChain *)fChain;
1844 if (chain->GetTreeNumber() != fCurrent) {
1845 fCurrent = chain->GetTreeNumber();
1846 Notify();
1847 }
1848 return centry;
1849 }
1850
1851 void GetEntries::Init(TTree *tree, const char *dupFileName) {
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862 t_hltbits = 0;
1863 t_ietaAll = 0;
1864 t_ietaGood = 0;
1865 t_trackType = 0;
1866 t_L1Bit = false;
1867 fChain = tree;
1868 fCurrent = -1;
1869 if (!tree)
1870 return;
1871 fChain->SetMakeClass(1);
1872 fChain->SetBranchAddress("t_RunNo", &t_RunNo, &b_t_RunNo);
1873 fChain->SetBranchAddress("t_EventNo", &t_EventNo, &b_t_EventNo);
1874 fChain->SetBranchAddress("t_Tracks", &t_Tracks, &b_t_Tracks);
1875 fChain->SetBranchAddress("t_TracksProp", &t_TracksProp, &b_t_TracksProp);
1876 fChain->SetBranchAddress("t_TracksSaved", &t_TracksSaved, &b_t_TracksSaved);
1877 fChain->SetBranchAddress("t_TracksLoose", &t_TracksLoose, &b_t_TracksLoose);
1878 fChain->SetBranchAddress("t_TracksTight", &t_TracksTight, &b_t_TracksTight);
1879 fChain->SetBranchAddress("t_allvertex", &t_allvertex, &b_t_allvertex);
1880 fChain->SetBranchAddress("t_TrigPass", &t_TrigPass, &b_t_TrigPass);
1881 fChain->SetBranchAddress("t_TrigPassSel", &t_TrigPassSel, &b_t_TrigPassSel);
1882 fChain->SetBranchAddress("t_L1Bit", &t_L1Bit, &b_t_L1Bit);
1883 fChain->SetBranchAddress("t_hltbits", &t_hltbits, &b_t_hltbits);
1884 fChain->SetBranchAddress("t_ietaAll", &t_ietaAll, &b_t_ietaAll);
1885 fChain->SetBranchAddress("t_ietaGood", &t_ietaGood, &b_t_ietaGood);
1886 fChain->SetBranchAddress("t_trackType", &t_trackType, &b_t_trackType);
1887 Notify();
1888
1889 if (std::string(dupFileName) != "") {
1890 ifstream infile(dupFileName);
1891 if (!infile.is_open()) {
1892 std::cout << "Cannot open " << dupFileName << std::endl;
1893 } else {
1894 while (1) {
1895 Long64_t jentry;
1896 infile >> jentry;
1897 if (!infile.good())
1898 break;
1899 entries_.push_back(jentry);
1900 }
1901 infile.close();
1902 std::cout << "Reads a list of " << entries_.size() << " events from " << dupFileName << std::endl;
1903 }
1904 }
1905
1906 h_tk[0] = new TH1I("Track0", "# of tracks produced", 2000, 0, 2000);
1907 h_tk[1] = new TH1I("Track1", "# of tracks propagated", 2000, 0, 2000);
1908 h_tk[2] = new TH1I("Track2", "# of tracks saved in tree", 2000, 0, 2000);
1909 h_eta[0] = new TH1I("Eta0", "i#eta (All Tracks)", 60, -30, 30);
1910 h_eta[1] = new TH1I("Eta1", "i#eta (Good Tracks)", 60, -30, 30);
1911 h_eta[2] = new TH1I("Eta2", "i#eta (Loose Isolated Tracks)", 60, -30, 30);
1912 h_eta[3] = new TH1I("Eta3", "i#eta (Tight Isolated Tracks)", 60, -30, 30);
1913 h_eff[0] = new TH1D("Eff0", "i#eta (Selection Efficiency)", 60, -30, 30);
1914 h_eff[1] = new TH1D("Eff1", "i#eta (Loose Isolation Efficiency)", 60, -30, 30);
1915 h_eff[2] = new TH1D("Eff2", "i#eta (Tight Isolation Efficiency)", 60, -30, 30);
1916 h_pvx[0] = new TH1I("Pvx0", "Number of Good Vertex (All)", 100, 0, 100);
1917 h_pvx[1] = new TH1I("Pvx1", "Number of Good Vertex (Loose)", 100, 0, 100);
1918 h_pvx[2] = new TH1I("Pvx2", "Number of Good Vertex (Tight)", 100, 0, 100);
1919 }
1920
1921 Bool_t GetEntries::Notify() {
1922
1923
1924
1925
1926
1927
1928 return kTRUE;
1929 }
1930
1931 void GetEntries::Show(Long64_t entry) {
1932
1933
1934 if (!fChain)
1935 return;
1936 fChain->Show(entry);
1937 }
1938
1939 Int_t GetEntries::Cut(Long64_t) {
1940
1941
1942
1943 return 1;
1944 }
1945
1946 void GetEntries::Loop() {
1947
1948
1949
1950
1951
1952
1953
1954
1955
1956
1957
1958
1959
1960
1961
1962
1963
1964
1965
1966
1967
1968
1969
1970 if (fChain == 0)
1971 return;
1972
1973 Long64_t nentries = fChain->GetEntriesFast();
1974 Long64_t nbytes = 0, nb = 0;
1975 unsigned int kount(0), duplicate(0), selected(0);
1976 int l1(0), hlt(0), loose(0), tight(0);
1977 int allHLT[3] = {0, 0, 0};
1978 int looseHLT[3] = {0, 0, 0};
1979 int tightHLT[3] = {0, 0, 0};
1980 for (Long64_t jentry = 0; jentry < nentries; jentry++) {
1981 Long64_t ientry = LoadTree(jentry);
1982 if (ientry < 0)
1983 break;
1984 nb = fChain->GetEntry(jentry);
1985 nbytes += nb;
1986 bool select = (std::find(entries_.begin(), entries_.end(), jentry) == entries_.end());
1987 if (!select) {
1988 ++duplicate;
1989 continue;
1990 }
1991 h_tk[0]->Fill(t_Tracks);
1992 h_tk[1]->Fill(t_TracksProp);
1993 h_tk[2]->Fill(t_TracksSaved);
1994 h_pvx[0]->Fill(t_allvertex);
1995 if (t_TracksLoose > 0)
1996 h_pvx[1]->Fill(t_allvertex);
1997 if (t_TracksTight > 0)
1998 h_pvx[2]->Fill(t_allvertex);
1999 if (t_L1Bit) {
2000 ++l1;
2001 if (t_TracksLoose > 0)
2002 ++loose;
2003 if (t_TracksTight > 0)
2004 ++tight;
2005 if (t_TrigPass)
2006 ++hlt;
2007 }
2008 if (t_TrigPass) {
2009 ++kount;
2010 if (t_TrigPassSel)
2011 ++selected;
2012 }
2013 bool passt[2] = {false, false}, pass(false);
2014 for (unsigned k = 0; k < t_hltbits->size(); ++k) {
2015 if ((*t_hltbits)[k] > 0) {
2016 pass = true;
2017 for (int i = 0; i < 2; ++i)
2018 if (k == bit_[i])
2019 passt[i] = true;
2020 }
2021 }
2022 if (pass) {
2023 ++allHLT[0];
2024 for (int i = 0; i < 2; ++i)
2025 if (passt[i])
2026 ++allHLT[i + 1];
2027 if (t_TracksLoose > 0) {
2028 ++looseHLT[0];
2029 for (int i = 0; i < 2; ++i)
2030 if (passt[i])
2031 ++looseHLT[i + 1];
2032 }
2033 if (t_TracksTight > 0) {
2034 ++tightHLT[0];
2035 for (int i = 0; i < 2; ++i)
2036 if (passt[i])
2037 ++tightHLT[i + 1];
2038 }
2039 }
2040 for (unsigned int k = 0; k < t_ietaAll->size(); ++k)
2041 h_eta[0]->Fill((*t_ietaAll)[k]);
2042 for (unsigned int k = 0; k < t_ietaGood->size(); ++k) {
2043 h_eta[1]->Fill((*t_ietaGood)[k]);
2044 if (t_trackType->size() > 0) {
2045 if ((*t_trackType)[k] > 0)
2046 h_eta[2]->Fill((*t_ietaGood)[k]);
2047 if ((*t_trackType)[k] > 1)
2048 h_eta[3]->Fill((*t_ietaGood)[k]);
2049 }
2050 }
2051 }
2052 double ymaxk(0);
2053 for (int i = 1; i <= h_eff[0]->GetNbinsX(); ++i) {
2054 double rat(0), drat(0);
2055 if (h_eta[0]->GetBinContent(i) > ymaxk)
2056 ymaxk = h_eta[0]->GetBinContent(i);
2057 if ((h_eta[1]->GetBinContent(i) > 0) && (h_eta[0]->GetBinContent(i) > 0)) {
2058 rat = h_eta[1]->GetBinContent(i) / h_eta[0]->GetBinContent(i);
2059 drat = rat * std::sqrt(pow((h_eta[1]->GetBinError(i) / h_eta[1]->GetBinContent(i)), 2) +
2060 pow((h_eta[0]->GetBinError(i) / h_eta[0]->GetBinContent(i)), 2));
2061 }
2062 h_eff[0]->SetBinContent(i, rat);
2063 h_eff[0]->SetBinError(i, drat);
2064 if ((h_eta[1]->GetBinContent(i) > 0) && (h_eta[2]->GetBinContent(i) > 0)) {
2065 rat = h_eta[2]->GetBinContent(i) / h_eta[1]->GetBinContent(i);
2066 drat = rat * std::sqrt(pow((h_eta[2]->GetBinError(i) / h_eta[2]->GetBinContent(i)), 2) +
2067 pow((h_eta[1]->GetBinError(i) / h_eta[1]->GetBinContent(i)), 2));
2068 } else {
2069 rat = drat = 0;
2070 }
2071 h_eff[1]->SetBinContent(i, rat);
2072 h_eff[1]->SetBinError(i, drat);
2073 if ((h_eta[1]->GetBinContent(i) > 0) && (h_eta[3]->GetBinContent(i) > 0)) {
2074 rat = h_eta[3]->GetBinContent(i) / h_eta[1]->GetBinContent(i);
2075 drat = rat * std::sqrt(pow((h_eta[3]->GetBinError(i) / h_eta[3]->GetBinContent(i)), 2) +
2076 pow((h_eta[1]->GetBinError(i) / h_eta[1]->GetBinContent(i)), 2));
2077 } else {
2078 rat = drat = 0;
2079 }
2080 h_eff[1]->SetBinContent(i, rat);
2081 h_eff[1]->SetBinError(i, drat);
2082 }
2083 std::cout << "===== Remove " << duplicate << " events from " << nentries << "\n===== " << kount
2084 << " events passed trigger of which " << selected << " events get selected =====\n"
2085 << std::endl;
2086 std::cout << "===== " << l1 << " events passed L1 " << hlt << " events passed HLT and " << loose << ":" << tight
2087 << " events have at least 1 track candidate with loose:tight"
2088 << " isolation cut =====\n"
2089 << std::endl;
2090 for (int i = 0; i < 3; ++i) {
2091 char tbit[20];
2092 if (i == 0)
2093 sprintf(tbit, "Any");
2094 else
2095 sprintf(tbit, "%3d", bit_[i - 1]);
2096 std::cout << "Satisfying HLT trigger bit " << tbit << " Kount " << allHLT[i] << " Loose " << looseHLT[i]
2097 << " Tight " << tightHLT[i] << std::endl;
2098 }
2099
2100 gStyle->SetCanvasBorderMode(0);
2101 gStyle->SetCanvasColor(kWhite);
2102 gStyle->SetPadColor(kWhite);
2103 gStyle->SetFillColor(kWhite);
2104 gStyle->SetOptStat(1110);
2105 gStyle->SetOptTitle(0);
2106 int color[5] = {kBlack, kRed, kBlue, kMagenta, kCyan};
2107 int lines[5] = {1, 2, 3, 4, 5};
2108 TCanvas *pad1 = new TCanvas("c_track", "c_track", 500, 500);
2109 pad1->SetRightMargin(0.10);
2110 pad1->SetTopMargin(0.10);
2111 pad1->SetFillColor(kWhite);
2112 std::string titl1[3] = {"Reconstructed", "Propagated", "Saved"};
2113 TLegend *legend1 = new TLegend(0.11, 0.80, 0.50, 0.89);
2114 legend1->SetFillColor(kWhite);
2115 double ymax(0), xmax(0);
2116 for (int k = 0; k < 3; ++k) {
2117 int total(0), totaltk(0);
2118 for (int i = 1; i <= h_tk[k]->GetNbinsX(); ++i) {
2119 if (ymax < h_tk[k]->GetBinContent(i))
2120 ymax = h_tk[k]->GetBinContent(i);
2121 if (i > 1)
2122 total += (int)(h_tk[k]->GetBinContent(i));
2123 totaltk += (int)(h_tk[k]->GetBinContent(i)) * (i - 1);
2124 if (h_tk[k]->GetBinContent(i) > 0) {
2125 if (xmax < h_tk[k]->GetBinLowEdge(i) + h_tk[k]->GetBinWidth(i))
2126 xmax = h_tk[k]->GetBinLowEdge(i) + h_tk[k]->GetBinWidth(i);
2127 }
2128 }
2129 h_tk[k]->SetLineColor(color[k]);
2130 h_tk[k]->SetMarkerColor(color[k]);
2131 h_tk[k]->SetLineStyle(lines[k]);
2132 std::cout << h_tk[k]->GetTitle() << " Entries " << h_tk[k]->GetEntries() << " Events " << total << " Tracks "
2133 << totaltk << std::endl;
2134 legend1->AddEntry(h_tk[k], titl1[k].c_str(), "l");
2135 }
2136 int i1 = (int)(0.1 * xmax) + 1;
2137 xmax = 10.0 * i1;
2138 int i2 = (int)(0.01 * ymax) + 1;
2139
2140 ymax = 100.0 * i2;
2141 for (int k = 0; k < 3; ++k) {
2142 h_tk[k]->GetXaxis()->SetRangeUser(0, xmax);
2143 h_tk[k]->GetYaxis()->SetRangeUser(0.1, ymax);
2144 h_tk[k]->GetXaxis()->SetTitle("# Tracks");
2145 h_tk[k]->GetYaxis()->SetTitle("Events");
2146 h_tk[k]->GetYaxis()->SetLabelOffset(0.005);
2147 h_tk[k]->GetYaxis()->SetTitleOffset(1.20);
2148 if (k == 0)
2149 h_tk[k]->Draw("hist");
2150 else
2151 h_tk[k]->Draw("hist sames");
2152 }
2153 pad1->Update();
2154 pad1->SetLogy();
2155 ymax = 0.90;
2156 for (int k = 0; k < 3; ++k) {
2157 TPaveStats *st1 = (TPaveStats *)h_tk[k]->GetListOfFunctions()->FindObject("stats");
2158 if (st1 != NULL) {
2159 st1->SetLineColor(color[k]);
2160 st1->SetTextColor(color[k]);
2161 st1->SetY1NDC(ymax - 0.09);
2162 st1->SetY2NDC(ymax);
2163 st1->SetX1NDC(0.55);
2164 st1->SetX2NDC(0.90);
2165 ymax -= 0.09;
2166 }
2167 }
2168 pad1->Modified();
2169 pad1->Update();
2170 legend1->Draw("same");
2171 pad1->Update();
2172
2173 TCanvas *pad2 = new TCanvas("c_ieta", "c_ieta", 500, 500);
2174 pad2->SetRightMargin(0.10);
2175 pad2->SetTopMargin(0.10);
2176 pad2->SetFillColor(kWhite);
2177 pad2->SetLogy();
2178 std::string titl2[4] = {"All Tracks", "Selected Tracks", "Loose Isolation", "Tight Isolation"};
2179 TLegend *legend2 = new TLegend(0.11, 0.75, 0.50, 0.89);
2180 legend2->SetFillColor(kWhite);
2181 i2 = (int)(0.001 * ymaxk) + 1;
2182 ymax = 1000.0 * i2;
2183 for (int k = 0; k < 4; ++k) {
2184 h_eta[k]->GetYaxis()->SetRangeUser(1, ymax);
2185 h_eta[k]->SetLineColor(color[k]);
2186 h_eta[k]->SetMarkerColor(color[k]);
2187 h_eta[k]->SetLineStyle(lines[k]);
2188 h_eta[k]->GetXaxis()->SetTitle("i#eta");
2189 h_eta[k]->GetYaxis()->SetTitle("Tracks");
2190 h_eta[k]->GetYaxis()->SetLabelOffset(0.005);
2191 h_eta[k]->GetYaxis()->SetTitleOffset(1.20);
2192 legend2->AddEntry(h_eta[k], titl2[k].c_str(), "l");
2193 if (k == 0)
2194 h_eta[k]->Draw("hist");
2195 else
2196 h_eta[k]->Draw("hist sames");
2197 }
2198 pad2->Update();
2199 ymax = 0.90;
2200
2201 for (int k = 0; k < 4; ++k) {
2202 TPaveStats *st1 = (TPaveStats *)h_eta[k]->GetListOfFunctions()->FindObject("stats");
2203 if (st1 != NULL) {
2204 st1->SetLineColor(color[k]);
2205 st1->SetTextColor(color[k]);
2206 st1->SetY1NDC(ymax - 0.09);
2207 st1->SetY2NDC(ymax);
2208 st1->SetX1NDC(0.55);
2209 st1->SetX2NDC(0.90);
2210 ymax -= 0.09;
2211 }
2212 }
2213 pad2->Modified();
2214 pad2->Update();
2215 legend2->Draw("same");
2216 pad2->Update();
2217
2218 std::string titl3[3] = {"Selection", "Loose Isolation", "Tight Isolation"};
2219 TCanvas *pad3 = new TCanvas("c_effi", "c_effi", 500, 500);
2220 TLegend *legend3 = new TLegend(0.11, 0.785, 0.50, 0.89);
2221 pad3->SetRightMargin(0.10);
2222 pad3->SetTopMargin(0.10);
2223 pad3->SetFillColor(kWhite);
2224 pad3->SetLogy();
2225 for (int k = 0; k < 3; ++k) {
2226 h_eff[k]->SetStats(0);
2227 h_eff[k]->SetMarkerStyle(20);
2228 h_eff[k]->SetLineColor(color[k]);
2229 h_eff[k]->SetMarkerColor(color[k]);
2230 h_eff[k]->GetXaxis()->SetTitle("i#eta");
2231 h_eff[k]->GetYaxis()->SetTitle("Efficiency");
2232 h_eff[k]->GetYaxis()->SetLabelOffset(0.005);
2233 h_eff[k]->GetYaxis()->SetTitleOffset(1.20);
2234 if (k == 0)
2235 h_eff[k]->Draw("");
2236 else
2237 h_eff[k]->Draw("same");
2238 legend3->AddEntry(h_eff[k], titl3[k].c_str(), "l");
2239 }
2240 pad3->Modified();
2241 pad3->Update();
2242 legend3->Draw("same");
2243 pad3->Update();
2244
2245 std::string titl4[3] = {"All events", "Events with Loose Isolation", "Events with Tight Isolation"};
2246 TLegend *legend4 = new TLegend(0.11, 0.785, 0.50, 0.89);
2247 TCanvas *pad4 = new TCanvas("c_nvx", "c_nvx", 500, 500);
2248 pad4->SetRightMargin(0.10);
2249 pad4->SetTopMargin(0.10);
2250 pad4->SetFillColor(kWhite);
2251 pad4->SetLogy();
2252 for (int k = 0; k < 3; ++k) {
2253 h_pvx[k]->SetStats(1110);
2254 h_pvx[k]->SetMarkerStyle(20);
2255 h_pvx[k]->SetLineColor(color[k]);
2256 h_pvx[k]->SetMarkerColor(color[k]);
2257 h_pvx[k]->GetXaxis()->SetTitle("N_{PV}");
2258 h_pvx[k]->GetYaxis()->SetTitle("Events");
2259 h_pvx[k]->GetYaxis()->SetLabelOffset(0.005);
2260 h_pvx[k]->GetYaxis()->SetTitleOffset(1.20);
2261 if (k == 0)
2262 h_pvx[k]->Draw("");
2263 else
2264 h_pvx[k]->Draw("sames");
2265 legend4->AddEntry(h_pvx[k], titl4[k].c_str(), "l");
2266 }
2267 pad4->Update();
2268 ymax = 0.90;
2269 for (int k = 0; k < 3; ++k) {
2270 TPaveStats *st1 = (TPaveStats *)h_pvx[k]->GetListOfFunctions()->FindObject("stats");
2271 if (st1 != NULL) {
2272 st1->SetLineColor(color[k]);
2273 st1->SetTextColor(color[k]);
2274 st1->SetY1NDC(ymax - 0.09);
2275 st1->SetY2NDC(ymax);
2276 st1->SetX1NDC(0.55);
2277 st1->SetX2NDC(0.90);
2278 ymax -= 0.09;
2279 }
2280 }
2281 pad4->Modified();
2282 pad4->Update();
2283 legend4->Draw("same");
2284 pad4->Update();
2285 }
2286
2287 class CalibPlotProperties {
2288 public:
2289 TChain *fChain;
2290 Int_t fCurrent;
2291
2292
2293 Int_t t_Run;
2294 Int_t t_Event;
2295 Int_t t_DataType;
2296 Int_t t_ieta;
2297 Int_t t_iphi;
2298 Double_t t_EventWeight;
2299 Int_t t_nVtx;
2300 Int_t t_nTrk;
2301 Int_t t_goodPV;
2302 Double_t t_l1pt;
2303 Double_t t_l1eta;
2304 Double_t t_l1phi;
2305 Double_t t_l3pt;
2306 Double_t t_l3eta;
2307 Double_t t_l3phi;
2308 Double_t t_p;
2309 Double_t t_pt;
2310 Double_t t_phi;
2311 Double_t t_mindR1;
2312 Double_t t_mindR2;
2313 Double_t t_eMipDR;
2314 Double_t t_eHcal;
2315 Double_t t_eHcal10;
2316 Double_t t_eHcal30;
2317 Double_t t_hmaxNearP;
2318 Double_t t_rhoh;
2319 Bool_t t_selectTk;
2320 Bool_t t_qltyFlag;
2321 Bool_t t_qltyMissFlag;
2322 Bool_t t_qltyPVFlag;
2323 Double_t t_gentrackP;
2324 std::vector<unsigned int> *t_DetIds;
2325 std::vector<double> *t_HitEnergies;
2326 std::vector<bool> *t_trgbits;
2327 std::vector<unsigned int> *t_DetIds1;
2328 std::vector<unsigned int> *t_DetIds3;
2329 std::vector<double> *t_HitEnergies1;
2330 std::vector<double> *t_HitEnergies3;
2331
2332
2333 TBranch *b_t_Run;
2334 TBranch *b_t_Event;
2335 TBranch *b_t_DataType;
2336 TBranch *b_t_ieta;
2337 TBranch *b_t_iphi;
2338 TBranch *b_t_EventWeight;
2339 TBranch *b_t_nVtx;
2340 TBranch *b_t_nTrk;
2341 TBranch *b_t_goodPV;
2342 TBranch *b_t_l1pt;
2343 TBranch *b_t_l1eta;
2344 TBranch *b_t_l1phi;
2345 TBranch *b_t_l3pt;
2346 TBranch *b_t_l3eta;
2347 TBranch *b_t_l3phi;
2348 TBranch *b_t_p;
2349 TBranch *b_t_pt;
2350 TBranch *b_t_phi;
2351 TBranch *b_t_mindR1;
2352 TBranch *b_t_mindR2;
2353 TBranch *b_t_eMipDR;
2354 TBranch *b_t_eHcal;
2355 TBranch *b_t_eHcal10;
2356 TBranch *b_t_eHcal30;
2357 TBranch *b_t_hmaxNearP;
2358 TBranch *b_t_rhoh;
2359 TBranch *b_t_selectTk;
2360 TBranch *b_t_qltyFlag;
2361 TBranch *b_t_qltyMissFlag;
2362 TBranch *b_t_qltyPVFlag;
2363 TBranch *b_t_gentrackP;
2364 TBranch *b_t_DetIds;
2365 TBranch *b_t_HitEnergies;
2366 TBranch *b_t_trgbits;
2367 TBranch *b_t_DetIds1;
2368 TBranch *b_t_DetIds3;
2369 TBranch *b_t_HitEnergies1;
2370 TBranch *b_t_HitEnergies3;
2371
2372 CalibPlotProperties(const char *fname,
2373 const std::string &dirname,
2374 const char *dupFileName,
2375 const std::string &prefix = "",
2376 const char *corrFileName = "",
2377 const char *rcorFileName = "",
2378 int puCorr = -2,
2379 int flag = 1031,
2380 bool dataMC = true,
2381 int truncateFlag = 0,
2382 bool useGen = false,
2383 double scale = 1.0,
2384 int useScale = 0,
2385 int etalo = 0,
2386 int etahi = 30,
2387 int runlo = 0,
2388 int runhi = 99999999,
2389 int phimin = 1,
2390 int phimax = 72,
2391 int zside = 1,
2392 int nvxlo = 0,
2393 int nvxhi = 1000,
2394 int rbx = 0,
2395 bool exclude = false,
2396 bool etamax = false);
2397 virtual ~CalibPlotProperties();
2398 virtual Int_t Cut(Long64_t entry);
2399 virtual Int_t GetEntry(Long64_t entry);
2400 virtual Long64_t LoadTree(Long64_t entry);
2401 virtual void Init(TChain *, const char *);
2402 virtual void Loop();
2403 virtual Bool_t Notify();
2404 virtual void Show(Long64_t entry = -1);
2405 bool goodTrack(double &eHcal, double &cut, const Long64_t &entry, bool debug);
2406 bool selectPhi(bool debug);
2407 void savePlot(const std::string &theName, bool append = true, bool all = false);
2408 void correctEnergy(double &ener, const Long64_t &entry);
2409
2410 private:
2411 static const unsigned int npbin = 6, kp50 = 3, ndepth = 7;
2412 CalibCorrFactor *corrFactor_;
2413 CalibCorr *cFactor_;
2414 CalibSelectRBX *cSelect_;
2415 const std::string fname_, dirnm_, prefix_, outFileName_;
2416 const int corrPU_, flag_;
2417 const bool dataMC_, useGen_;
2418 const int truncateFlag_;
2419 const int etalo_, etahi_;
2420 int runlo_, runhi_;
2421 const int phimin_, phimax_, zside_, nvxlo_, nvxhi_, rbx_;
2422 int ifDepth_;
2423 bool exclude_, corrE_, cutL1T_;
2424 bool includeRun_, getHist_;
2425 int flexibleSelect_;
2426 bool plotBasic_, plotEnergy_, plotHists_;
2427 double log2by18_;
2428 std::ofstream fileout_;
2429 std::vector<Long64_t> entries_;
2430 std::vector<std::pair<int, int> > events_;
2431 std::vector<double> etas_, ps_, dl1_;
2432 std::vector<int> nvx_, ietas_;
2433 TH1D *h_p[5], *h_eta[5], *h_nvtx;
2434 std::vector<TH1D *> h_eta0, h_eta1, h_eta2, h_eta3, h_eta4;
2435 std::vector<TH1D *> h_dL1, h_vtx;
2436 std::vector<TH1D *> h_etaEH[npbin], h_etaEp[npbin], h_etaEE[npbin];
2437 std::vector<TH1D *> h_mom, h_eEcal, h_eHcal;
2438 std::vector<TH1F *> h_bvlist, h_bvlist2, h_evlist, h_evlist2;
2439 std::vector<TH1F *> h_bvlist3, h_evlist3;
2440 TH2F *h_etaE;
2441 };
2442
2443 CalibPlotProperties::CalibPlotProperties(const char *fname,
2444 const std::string &dirnm,
2445 const char *dupFileName,
2446 const std::string &prefix,
2447 const char *corrFileName,
2448 const char *rcorFileName,
2449 int puCorr,
2450 int flag,
2451 bool dataMC,
2452 int truncate,
2453 bool useGen,
2454 double scl,
2455 int useScale,
2456 int etalo,
2457 int etahi,
2458 int runlo,
2459 int runhi,
2460 int phimin,
2461 int phimax,
2462 int zside,
2463 int nvxlo,
2464 int nvxhi,
2465 int rbx,
2466 bool exc,
2467 bool etam)
2468 : corrFactor_(nullptr),
2469 cFactor_(nullptr),
2470 cSelect_(nullptr),
2471 fname_(fname),
2472 dirnm_(dirnm),
2473 prefix_(prefix),
2474 corrPU_(puCorr),
2475 flag_(flag),
2476 dataMC_(dataMC),
2477 useGen_(useGen),
2478 truncateFlag_(truncate),
2479 etalo_(etalo),
2480 etahi_(etahi),
2481 runlo_(runlo),
2482 runhi_(runhi),
2483 phimin_(phimin),
2484 phimax_(phimax),
2485 zside_(zside),
2486 nvxlo_(nvxlo),
2487 nvxhi_(nvxhi),
2488 rbx_(rbx),
2489 exclude_(exc),
2490 includeRun_(true) {
2491
2492
2493
2494 flexibleSelect_ = ((flag_ / 1) % 10);
2495 plotBasic_ = (((flag_ / 10) % 10) > 0);
2496 plotEnergy_ = (((flag_ / 10) % 10) > 0);
2497 int oneplace = ((flag_ / 1000) % 10);
2498 cutL1T_ = (oneplace % 2);
2499 bool marina = ((oneplace / 2) % 2);
2500 ifDepth_ = ((flag_ / 10000) % 10);
2501 plotHists_ = (((flag_ / 100000) % 10) > 0);
2502 log2by18_ = std::log(2.5) / 18.0;
2503 if (runlo_ < 0 || runhi_ < 0) {
2504 runlo_ = std::abs(runlo_);
2505 runhi_ = std::abs(runhi_);
2506 includeRun_ = false;
2507 }
2508 char treeName[400];
2509 sprintf(treeName, "%s/CalibTree", dirnm.c_str());
2510 TChain *chain = new TChain(treeName);
2511 std::cout << "Create a chain for " << treeName << " from " << fname << " flags " << flexibleSelect_ << "|"
2512 << plotBasic_ << "|"
2513 << "|" << plotEnergy_ << "|" << plotHists_ << "|" << corrPU_ << "\n cons " << log2by18_ << " eta range "
2514 << etalo_ << ":" << etahi_ << " run range " << runlo_ << ":" << runhi_ << " (inclusion flag " << includeRun_
2515 << ")\n Selection of RBX" << rbx << " Vertex Range " << nvxlo_ << ":" << nvxhi_
2516 << "\n corrFileName: " << corrFileName << " useScale " << useScale << ":" << scl << ":" << etam
2517 << "\n rcorFileName: " << rcorFileName << " flag " << ifDepth_ << std::endl;
2518 if (!fillChain(chain, fname)) {
2519 std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
2520 } else {
2521 std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
2522 Init(chain, dupFileName);
2523 corrFactor_ = new CalibCorrFactor(corrFileName, useScale, scl, etam, marina, false);
2524 if (std::string(rcorFileName) != "") {
2525 cFactor_ = new CalibCorr(rcorFileName, ifDepth_, false);
2526 } else {
2527 ifDepth_ = 0;
2528 }
2529 if (rbx != 0)
2530 cSelect_ = new CalibSelectRBX(rbx, false);
2531 }
2532 }
2533
2534 CalibPlotProperties::~CalibPlotProperties() {
2535 delete corrFactor_;
2536 delete cFactor_;
2537 delete cSelect_;
2538 if (!fChain)
2539 return;
2540 delete fChain->GetCurrentFile();
2541 }
2542
2543 Int_t CalibPlotProperties::GetEntry(Long64_t entry) {
2544
2545 if (!fChain)
2546 return 0;
2547 return fChain->GetEntry(entry);
2548 }
2549
2550 Long64_t CalibPlotProperties::LoadTree(Long64_t entry) {
2551
2552 if (!fChain)
2553 return -5;
2554 Long64_t centry = fChain->LoadTree(entry);
2555 if (centry < 0)
2556 return centry;
2557 if (!fChain->InheritsFrom(TChain::Class()))
2558 return centry;
2559 TChain *chain = (TChain *)fChain;
2560 if (chain->GetTreeNumber() != fCurrent) {
2561 fCurrent = chain->GetTreeNumber();
2562 Notify();
2563 }
2564 return centry;
2565 }
2566
2567 void CalibPlotProperties::Init(TChain *tree, const char *dupFileName) {
2568
2569
2570
2571
2572
2573
2574
2575
2576
2577 t_DetIds = 0;
2578 t_DetIds1 = 0;
2579 t_DetIds3 = 0;
2580 t_HitEnergies = 0;
2581 t_HitEnergies1 = 0;
2582 t_HitEnergies3 = 0;
2583 t_trgbits = 0;
2584
2585 fChain = tree;
2586 fCurrent = -1;
2587 if (!tree)
2588 return;
2589 fChain->SetMakeClass(1);
2590
2591 fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
2592 fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
2593 fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
2594 fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
2595 fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
2596 fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
2597 fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
2598 fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
2599 fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
2600 fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
2601 fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
2602 fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
2603 fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
2604 fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
2605 fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
2606 fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
2607 fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
2608 fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
2609 fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
2610 fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
2611 fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
2612 fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
2613 fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
2614 fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
2615 fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
2616 fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
2617 fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
2618 fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
2619 fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
2620 fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
2621 fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
2622 fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
2623 fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
2624 fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
2625 fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
2626 fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
2627 fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
2628 fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
2629 Notify();
2630
2631 if (std::string(dupFileName) != "") {
2632 ifstream infil1(dupFileName);
2633 if (!infil1.is_open()) {
2634 std::cout << "Cannot open duplicate file " << dupFileName << std::endl;
2635 } else {
2636 while (1) {
2637 Long64_t jentry;
2638 infil1 >> jentry;
2639 if (!infil1.good())
2640 break;
2641 entries_.push_back(jentry);
2642 }
2643 infil1.close();
2644 std::cout << "Reads a list of " << entries_.size() << " events from " << dupFileName << std::endl;
2645 }
2646 } else {
2647 std::cout << "No duplicate events in the input file" << std::endl;
2648 }
2649
2650 int ipbin[npbin] = {10, 20, 30, 40, 60, 100};
2651 for (unsigned int i = 0; i < npbin; ++i)
2652 ps_.push_back((double)(ipbin[i]));
2653 int ietas[4] = {0, 13, 18, 23};
2654 for (int i = 0; i < 4; ++i)
2655 ietas_.push_back(ietas[i]);
2656
2657 char name[20], title[200];
2658 unsigned int kp = ps_.size() - 1;
2659 unsigned int kk(0);
2660 std::string titl[5] = {
2661 "all tracks", "good quality tracks", "selected tracks", "isolated good tracks", "tracks having MIP in ECAL"};
2662
2663 if (plotBasic_) {
2664 std::cout << "Book Basic Histos" << std::endl;
2665 for (int k = 0; k < 5; ++k) {
2666 sprintf(name, "%sp%d", prefix_.c_str(), k);
2667 sprintf(title, "Momentum for %s", titl[k].c_str());
2668 h_p[k] = new TH1D(name, title, 100, 10.0, 110.0);
2669 sprintf(name, "%seta%d", prefix_.c_str(), k);
2670 sprintf(title, "#eta for %s", titl[k].c_str());
2671 h_eta[k] = new TH1D(name, title, 60, -30.0, 30.0);
2672 }
2673 for (unsigned int k = 0; k < kp; ++k) {
2674 sprintf(name, "%seta0%d", prefix_.c_str(), k);
2675 sprintf(title, "#eta for %s (p = %d:%d GeV)", titl[0].c_str(), ipbin[k], ipbin[k + 1]);
2676 h_eta0.push_back(new TH1D(name, title, 60, -30.0, 30.0));
2677 kk = h_eta0.size() - 1;
2678 h_eta0[kk]->Sumw2();
2679 sprintf(name, "%seta1%d", prefix_.c_str(), k);
2680 sprintf(title, "#eta for %s (p = %d:%d GeV)", titl[1].c_str(), ipbin[k], ipbin[k + 1]);
2681 h_eta1.push_back(new TH1D(name, title, 60, -30.0, 30.0));
2682 kk = h_eta1.size() - 1;
2683 h_eta1[kk]->Sumw2();
2684 sprintf(name, "%seta2%d", prefix_.c_str(), k);
2685 sprintf(title, "#eta for %s (p = %d:%d GeV)", titl[2].c_str(), ipbin[k], ipbin[k + 1]);
2686 h_eta2.push_back(new TH1D(name, title, 60, -30.0, 30.0));
2687 kk = h_eta2.size() - 1;
2688 h_eta2[kk]->Sumw2();
2689 sprintf(name, "%seta3%d", prefix_.c_str(), k);
2690 sprintf(title, "#eta for %s (p = %d:%d GeV)", titl[3].c_str(), ipbin[k], ipbin[k + 1]);
2691 h_eta3.push_back(new TH1D(name, title, 60, -30.0, 30.0));
2692 kk = h_eta3.size() - 1;
2693 h_eta3[kk]->Sumw2();
2694 sprintf(name, "%seta4%d", prefix_.c_str(), k);
2695 sprintf(title, "#eta for %s (p = %d:%d GeV)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
2696 h_eta4.push_back(new TH1D(name, title, 60, -30.0, 30.0));
2697 kk = h_eta4.size() - 1;
2698 h_eta4[kk]->Sumw2();
2699 sprintf(name, "%sdl1%d", prefix_.c_str(), k);
2700 sprintf(title, "Distance from L1 (p = %d:%d GeV)", ipbin[k], ipbin[k + 1]);
2701 h_dL1.push_back(new TH1D(name, title, 160, 0.0, 8.0));
2702 kk = h_dL1.size() - 1;
2703 h_dL1[kk]->Sumw2();
2704 sprintf(name, "%svtx%d", prefix_.c_str(), k);
2705 sprintf(title, "N_{Vertex} (p = %d:%d GeV)", ipbin[k], ipbin[k + 1]);
2706 h_vtx.push_back(new TH1D(name, title, 100, 0.0, 100.0));
2707 kk = h_vtx.size() - 1;
2708 h_vtx[kk]->Sumw2();
2709 }
2710 }
2711
2712 if (plotEnergy_) {
2713 std::cout << "Make plots for good tracks" << std::endl;
2714 for (unsigned int k = 0; k < kp; ++k) {
2715 for (int j = etalo_; j <= etahi_ + 1; ++j) {
2716 sprintf(name, "%senergyH%d%d", prefix_.c_str(), k, j);
2717 if (j > etahi_)
2718 sprintf(title,
2719 "HCAL energy for %s (p = %d:%d GeV |#eta| = %d:%d)",
2720 titl[3].c_str(),
2721 ipbin[k],
2722 ipbin[k + 1],
2723 etalo_,
2724 etahi_);
2725 else
2726 sprintf(title, "HCAL energy for %s (p = %d:%d GeV |#eta| = %d)", titl[3].c_str(), ipbin[k], ipbin[k + 1], j);
2727 h_etaEH[k].push_back(new TH1D(name, title, 200, 0, 200));
2728 kk = h_etaEH[k].size() - 1;
2729 h_etaEH[k][kk]->Sumw2();
2730 sprintf(name, "%senergyP%d%d", prefix_.c_str(), k, j);
2731 if (j > etahi_)
2732 sprintf(title,
2733 "momentum for %s (p = %d:%d GeV |#eta| = %d:%d)",
2734 titl[3].c_str(),
2735 ipbin[k],
2736 ipbin[k + 1],
2737 etalo_,
2738 etahi_);
2739 else
2740 sprintf(title, "momentum for %s (p = %d:%d GeV |#eta| = %d)", titl[3].c_str(), ipbin[k], ipbin[k + 1], j);
2741 h_etaEp[k].push_back(new TH1D(name, title, 100, 0, 100));
2742 kk = h_etaEp[k].size() - 1;
2743 h_etaEp[k][kk]->Sumw2();
2744 sprintf(name, "%senergyE%d%d", prefix_.c_str(), k, j);
2745 if (j > etahi_)
2746 sprintf(title,
2747 "ECAL energy for %s (p = %d:%d GeV |#eta| = %d:%d)",
2748 titl[3].c_str(),
2749 ipbin[k],
2750 ipbin[k + 1],
2751 etalo_,
2752 etahi_);
2753 else
2754 sprintf(title, "ECAL energy for %s (p = %d:%d GeV |#eta| = %d)", titl[3].c_str(), ipbin[k], ipbin[k + 1], j);
2755 h_etaEE[k].push_back(new TH1D(name, title, 100, 0, 10));
2756 kk = h_etaEE[k].size() - 1;
2757 h_etaEE[k][kk]->Sumw2();
2758 }
2759 }
2760
2761 for (unsigned int j = 0; j < ietas_.size(); ++j) {
2762 sprintf(name, "%senergyH%d", prefix_.c_str(), j);
2763 if (j == 0)
2764 sprintf(title, "HCAL energy for %s (All)", titl[3].c_str());
2765 else
2766 sprintf(title, "HCAL energy for %s (|#eta| = %d:%d)", titl[3].c_str(), ietas_[j - 1], ietas_[j]);
2767 h_eHcal.push_back(new TH1D(name, title, 200, 0, 200));
2768 kk = h_eHcal.size() - 1;
2769 h_eHcal[kk]->Sumw2();
2770 sprintf(name, "%senergyP%d", prefix_.c_str(), j);
2771 if (j == 0)
2772 sprintf(title, "Track momentum for %s (All)", titl[3].c_str());
2773 else
2774 sprintf(title, "Track momentum for %s (|#eta| = %d:%d)", titl[3].c_str(), ietas_[j - 1], ietas_[j]);
2775 h_mom.push_back(new TH1D(name, title, 100, 0, 100));
2776 kk = h_mom.size() - 1;
2777 h_mom[kk]->Sumw2();
2778 sprintf(name, "%senergyE%d", prefix_.c_str(), j);
2779 if (j == 0)
2780 sprintf(title, "ECAL energy for %s (All)", titl[3].c_str());
2781 else
2782 sprintf(title, "ECAL energy for %s (|#eta| = %d:%d)", titl[3].c_str(), ietas_[j - 1], ietas_[j]);
2783 h_eEcal.push_back(new TH1D(name, title, 100, 0, 10));
2784 kk = h_eEcal.size() - 1;
2785 h_eEcal[kk]->Sumw2();
2786 }
2787 }
2788
2789 if (plotHists_) {
2790 h_nvtx = new TH1D("hnvtx", "Number of vertices", 10, 0, 100);
2791 h_nvtx->Sumw2();
2792 for (unsigned int i = 0; i < ndepth; i++) {
2793 sprintf(name, "b_edepth%d", i);
2794 sprintf(title, "Total RecHit energy in depth %d (Barrel)", i + 1);
2795 h_bvlist.push_back(new TH1F(name, title, 1000, 0, 100));
2796 h_bvlist[i]->Sumw2();
2797 sprintf(name, "b_recedepth%d", i);
2798 sprintf(title, "RecHit energy in depth %d (Barrel)", i + 1);
2799 h_bvlist2.push_back(new TH1F(name, title, 1000, 0, 100));
2800 h_bvlist2[i]->Sumw2();
2801 sprintf(name, "b_nrecdepth%d", i);
2802 sprintf(title, "#RecHits in depth %d (Barrel)", i + 1);
2803 h_bvlist3.push_back(new TH1F(name, title, 1000, 0, 100));
2804 h_bvlist3[i]->Sumw2();
2805 sprintf(name, "e_edepth%d", i);
2806 sprintf(title, "Total RecHit energy in depth %d (Endcap)", i + 1);
2807 h_evlist.push_back(new TH1F(name, title, 1000, 0, 100));
2808 h_evlist[i]->Sumw2();
2809 sprintf(name, "e_recedepth%d", i);
2810 sprintf(title, "RecHit energy in depth %d (Endcap)", i + 1);
2811 h_evlist2.push_back(new TH1F(name, title, 1000, 0, 100));
2812 h_evlist2[i]->Sumw2();
2813 sprintf(name, "e_nrecdepth%d", i);
2814 sprintf(title, "#RecHits in depth %d (Endcap)", i + 1);
2815 h_evlist3.push_back(new TH1F(name, title, 1000, 0, 100));
2816 h_evlist3[i]->Sumw2();
2817 }
2818 h_etaE = new TH2F("heta", "", 50, -25, 25, 100, 0, 100);
2819 }
2820 }
2821
2822 Bool_t CalibPlotProperties::Notify() {
2823
2824
2825
2826
2827
2828
2829 return kTRUE;
2830 }
2831
2832 void CalibPlotProperties::Show(Long64_t entry) {
2833
2834
2835 if (!fChain)
2836 return;
2837 fChain->Show(entry);
2838 }
2839
2840 Int_t CalibPlotProperties::Cut(Long64_t) {
2841
2842
2843
2844 return 1;
2845 }
2846
2847 void CalibPlotProperties::Loop() {
2848
2849
2850
2851
2852
2853
2854
2855
2856
2857
2858
2859
2860
2861
2862
2863
2864
2865
2866
2867
2868
2869
2870
2871 const bool debug(false);
2872
2873 if (fChain == 0)
2874 return;
2875
2876
2877 Long64_t nentries = fChain->GetEntriesFast();
2878 std::cout << "Total entries " << nentries << std::endl;
2879 Long64_t nbytes(0), nb(0);
2880 unsigned int duplicate(0), good(0), kount(0);
2881 double sel(0), selHB(0), selHE(0);
2882 for (Long64_t jentry = 0; jentry < nentries; jentry++) {
2883 Long64_t ientry = LoadTree(jentry);
2884 if (ientry < 0)
2885 break;
2886 nb = fChain->GetEntry(jentry);
2887 nbytes += nb;
2888 if (jentry % 1000000 == 0)
2889 std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
2890 bool select = (std::find(entries_.begin(), entries_.end(), jentry) == entries_.end());
2891 if (!select) {
2892 ++duplicate;
2893 if (debug)
2894 std::cout << "Duplicate event " << t_Run << " " << t_Event << " " << t_p << std::endl;
2895 continue;
2896 }
2897 bool selRun = (includeRun_ ? ((t_Run >= runlo_) && (t_Run <= runhi_)) : ((t_Run < runlo_) || (t_Run > runhi_)));
2898 select =
2899 (selRun && (fabs(t_ieta) >= etalo_) && (fabs(t_ieta) <= etahi_) && (t_nVtx >= nvxlo_) && (t_nVtx <= nvxhi_));
2900 if (!select) {
2901 if (debug)
2902 std::cout << "Run # " << t_Run << " out of range of " << runlo_ << ":" << runhi_ << " or ieta " << t_ieta
2903 << " (" << etalo_ << ":" << etahi_ << ") or nvtx " << t_nVtx << " (" << nvxlo_ << ":" << nvxhi_
2904 << ") out of range" << std::endl;
2905 continue;
2906 }
2907 if (cSelect_ != nullptr) {
2908 if (exclude_) {
2909 if (cSelect_->isItRBX(t_DetIds))
2910 continue;
2911 } else {
2912 if (!(cSelect_->isItRBX(t_ieta, t_iphi)))
2913 continue;
2914 }
2915 }
2916 select = (!cutL1T_ || (t_mindR1 >= 0.5));
2917 if (!select) {
2918 if (debug)
2919 std::cout << "Reject Run # " << t_Run << " Event # " << t_Event << " too close to L1 trigger " << t_mindR1
2920 << std::endl;
2921 continue;
2922 }
2923 select = ((events_.size() == 0) ||
2924 (std::find(events_.begin(), events_.end(), std::pair<int, int>(t_Run, t_Event)) != events_.end()));
2925 if (!select) {
2926 if (debug)
2927 std::cout << "Reject Run # " << t_Run << " Event # " << t_Event << " not in the selection list" << std::endl;
2928 continue;
2929 }
2930
2931 if ((ifDepth_ == 3) && (cFactor_ != nullptr) && (cFactor_->absent(ientry)))
2932 continue;
2933
2934 unsigned int kp = ps_.size();
2935 double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
2936 for (unsigned int k = 1; k < ps_.size(); ++k) {
2937 if (pmom >= ps_[k - 1] && pmom < ps_[k]) {
2938 kp = k - 1;
2939 break;
2940 }
2941 }
2942 int jp1 = (((std::abs(t_ieta) >= etalo_) && (std::abs(t_ieta) <= etahi_)) ? (std::abs(t_ieta) - etalo_) : -1);
2943 int jp2 = (etahi_ - etalo_ + 1);
2944 unsigned int je1 = ietas_.size();
2945 for (unsigned int j = 1; j < ietas_.size(); ++j) {
2946 if (std::abs(t_ieta) > ietas_[j - 1] && std::abs(t_ieta) <= ietas_[j]) {
2947 je1 = j;
2948 break;
2949 }
2950 }
2951 int je2 = (je1 != ietas_.size()) ? 0 : -1;
2952 if (debug)
2953 std::cout << "Bin " << kp << ":" << je1 << ":" << je2 << ":" << jp1 << ":" << jp2 << std::endl;
2954 double cut = (pmom > 20) ? ((flexibleSelect_ == 0) ? 2.0 : 10.0) : 0.0;
2955 double rcut(-1000.0);
2956
2957
2958 double eHcal(t_eHcal);
2959 if (corrFactor_->doCorr() || (cFactor_ != 0)) {
2960 eHcal = 0;
2961 for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
2962
2963 double cfac(1.0);
2964 if (corrFactor_->doCorr()) {
2965 unsigned int id = truncateId((*t_DetIds)[k], truncateFlag_, false);
2966 cfac = corrFactor_->getCorr(id);
2967 }
2968 if (cFactor_ != 0)
2969 cfac *= cFactor_->getCorr(t_Run, (*t_DetIds)[k]);
2970 eHcal += (cfac * ((*t_HitEnergies)[k]));
2971 if (debug) {
2972 int subdet, zside, ieta, iphi, depth;
2973 unpackDetId((*t_DetIds)[k], subdet, zside, ieta, iphi, depth);
2974 std::cout << zside << ":" << ieta << ":" << depth << " Corr " << cfac << " " << (*t_HitEnergies)[k] << " Out "
2975 << eHcal << std::endl;
2976 }
2977 }
2978 }
2979 bool goodTk = goodTrack(eHcal, cut, jentry, debug);
2980 bool selPhi = selectPhi(debug);
2981 double rat = (pmom > 0) ? (eHcal / (pmom - t_eMipDR)) : 1.0;
2982 if (debug)
2983 std::cout << "Entry " << jentry << " p|eHcal|ratio " << pmom << "|" << t_eHcal << "|" << eHcal << "|" << rat
2984 << "|" << kp << " Cuts " << t_qltyFlag << "|" << t_selectTk << "|" << (t_hmaxNearP < cut) << "|"
2985 << (t_eMipDR < 1.0) << "|" << goodTk << "|" << (rat > rcut) << " Select Phi " << selPhi << std::endl;
2986 if (plotBasic_) {
2987 h_p[0]->Fill(pmom, t_EventWeight);
2988 h_eta[0]->Fill(t_ieta, t_EventWeight);
2989 if (kp < ps_.size())
2990 h_eta0[kp]->Fill(t_ieta, t_EventWeight);
2991 if (t_qltyFlag) {
2992 h_p[1]->Fill(pmom, t_EventWeight);
2993 h_eta[1]->Fill(t_ieta, t_EventWeight);
2994 if (kp < ps_.size())
2995 h_eta1[kp]->Fill(t_ieta, t_EventWeight);
2996 if (t_selectTk) {
2997 h_p[2]->Fill(pmom, t_EventWeight);
2998 h_eta[2]->Fill(t_ieta, t_EventWeight);
2999 if (kp < ps_.size())
3000 h_eta2[kp]->Fill(t_ieta, t_EventWeight);
3001 if (t_hmaxNearP < cut) {
3002 h_p[3]->Fill(pmom, t_EventWeight);
3003 h_eta[3]->Fill(t_ieta, t_EventWeight);
3004 if (kp < ps_.size())
3005 h_eta3[kp]->Fill(t_ieta, t_EventWeight);
3006 if (t_eMipDR < 1.0) {
3007 h_p[4]->Fill(pmom, t_EventWeight);
3008 h_eta[4]->Fill(t_ieta, t_EventWeight);
3009 if (kp < ps_.size()) {
3010 h_eta4[kp]->Fill(t_ieta, t_EventWeight);
3011 h_dL1[kp]->Fill(t_mindR1, t_EventWeight);
3012 h_vtx[kp]->Fill(t_goodPV, t_EventWeight);
3013 }
3014 }
3015 }
3016 }
3017 }
3018 }
3019
3020 if (goodTk && kp < ps_.size() && selPhi) {
3021 if (rat > rcut) {
3022 if (plotEnergy_) {
3023 if (jp1 >= 0) {
3024 h_etaEH[kp][jp1]->Fill(eHcal, t_EventWeight);
3025 h_etaEH[kp][jp2]->Fill(eHcal, t_EventWeight);
3026 h_etaEp[kp][jp1]->Fill(pmom, t_EventWeight);
3027 h_etaEp[kp][jp2]->Fill(pmom, t_EventWeight);
3028 h_etaEE[kp][jp1]->Fill(t_eMipDR, t_EventWeight);
3029 h_etaEE[kp][jp2]->Fill(t_eMipDR, t_EventWeight);
3030 }
3031 if (kp == kp50) {
3032 if (je1 != ietas_.size()) {
3033 h_eHcal[je1]->Fill(eHcal, t_EventWeight);
3034 h_eHcal[je2]->Fill(eHcal, t_EventWeight);
3035 h_mom[je1]->Fill(pmom, t_EventWeight);
3036 h_mom[je2]->Fill(pmom, t_EventWeight);
3037 h_eEcal[je1]->Fill(t_eMipDR, t_EventWeight);
3038 h_eEcal[je2]->Fill(t_eMipDR, t_EventWeight);
3039 }
3040 }
3041 }
3042
3043 if (plotHists_) {
3044 h_nvtx->Fill(t_nVtx);
3045 bool bad(false);
3046 for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
3047 unsigned int id = truncateId((*t_DetIds)[k], truncateFlag_, false);
3048 double cfac = corrFactor_->getCorr(id);
3049 if (cFactor_ != 0)
3050 cfac *= cFactor_->getCorr(t_Run, (*t_DetIds)[k]);
3051 double ener = cfac * (*t_HitEnergies)[k];
3052 if (corrPU_)
3053 correctEnergy(ener, jentry);
3054 if (ener < 0.001)
3055 bad = true;
3056 }
3057 if ((!bad) && (std::fabs(rat - 1) < 0.15) && (kp == kp50) &&
3058 ((std::abs(t_ieta) < 15) || (std::abs(t_ieta) > 17))) {
3059 float weight = (dataMC_ ? t_EventWeight : t_EventWeight * puweight(t_nVtx));
3060 h_etaE->Fill(t_ieta, eHcal, weight);
3061 sel += weight;
3062 std::vector<float> bv(7, 0.0f), ev(7, 0.0f);
3063 std::vector<int> bnrec(7, 0), enrec(7, 0);
3064 double eb(0), ee(0);
3065 for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
3066 unsigned int id = truncateId((*t_DetIds)[k], truncateFlag_, false);
3067 double cfac = corrFactor_->getCorr(id);
3068 if (cFactor_ != 0)
3069 cfac *= cFactor_->getCorr(t_Run, (*t_DetIds)[k]);
3070 double ener = cfac * (*t_HitEnergies)[k];
3071 if (corrPU_)
3072 correctEnergy(ener, jentry);
3073 unsigned int idx = (unsigned int)((*t_DetIds)[k]);
3074 int subdet, zside, ieta, iphi, depth;
3075 unpackDetId(idx, subdet, zside, ieta, iphi, depth);
3076 if (depth > 0 && depth <= (int)(ndepth)) {
3077 if (subdet == 1) {
3078 eb += ener;
3079 bv[depth - 1] += ener;
3080 h_bvlist2[depth - 1]->Fill(ener, weight);
3081 ++bnrec[depth - 1];
3082 } else if (subdet == 2) {
3083 ee += ener;
3084 ev[depth - 1] += ener;
3085 h_evlist2[depth - 1]->Fill(ener, weight);
3086 ++enrec[depth - 1];
3087 }
3088 }
3089 }
3090 bool barrel = (eb > ee);
3091 if (barrel)
3092 selHB += weight;
3093 else
3094 selHE += weight;
3095 for (unsigned int i = 0; i < ndepth; i++) {
3096 if (barrel) {
3097 h_bvlist[i]->Fill(bv[i], weight);
3098 h_bvlist3[i]->Fill((bnrec[i] + 0.001), weight);
3099 } else {
3100 h_evlist[i]->Fill(ev[i], weight);
3101 h_evlist3[i]->Fill((enrec[i] + 0.001), weight);
3102 }
3103 }
3104 }
3105 }
3106 }
3107 good++;
3108 }
3109 ++kount;
3110 }
3111 std::cout << "Finds " << duplicate << " Duplicate events out of " << kount << " events in this file and " << good
3112 << " selected events" << std::endl;
3113 if (plotHists_)
3114 std::cout << "Number of weighted selected events " << sel << " HB " << selHB << " HE " << selHE << std::endl;
3115 }
3116
3117 bool CalibPlotProperties::goodTrack(double &eHcal, double &cuti, const Long64_t &entry, bool debug) {
3118 bool select(true);
3119 double cut(cuti);
3120 if (debug) {
3121 std::cout << "goodTrack input " << eHcal << ":" << cut;
3122 }
3123 if (flexibleSelect_ > 1) {
3124 double eta = (t_ieta > 0) ? t_ieta : -t_ieta;
3125 cut = 8.0 * exp(eta * log2by18_);
3126 }
3127 correctEnergy(eHcal, entry);
3128 select = ((t_qltyFlag) && (t_selectTk) && (t_hmaxNearP < cut) && (t_eMipDR < 100.0) && (eHcal > 0.001));
3129 if (debug) {
3130 std::cout << " output " << eHcal << ":" << cut << ":" << select << std::endl;
3131 }
3132 return select;
3133 }
3134
3135 bool CalibPlotProperties::selectPhi(bool debug) {
3136 bool select(true);
3137 if (phimin_ > 1 || phimax_ < 72) {
3138 double eTotal(0), eSelec(0);
3139
3140 for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
3141 int iphi = ((*t_DetIds)[k]) & (0x3FF);
3142 int zside = ((*t_DetIds)[k] & 0x80000) ? (1) : (-1);
3143 eTotal += ((*t_HitEnergies)[k]);
3144 if (iphi >= phimin_ && iphi <= phimax_ && zside == zside_)
3145 eSelec += ((*t_HitEnergies)[k]);
3146 }
3147 if (eSelec < 0.9 * eTotal)
3148 select = false;
3149 if (debug) {
3150 std::cout << "Etotal " << eTotal << " and ESelec " << eSelec << " (phi " << phimin_ << ":" << phimax_ << " z "
3151 << zside_ << ") Selection " << select << std::endl;
3152 }
3153 }
3154 return select;
3155 }
3156
3157 void CalibPlotProperties::savePlot(const std::string &theName, bool append, bool all) {
3158 TFile *theFile(0);
3159 if (append) {
3160 theFile = new TFile(theName.c_str(), "UPDATE");
3161 } else {
3162 theFile = new TFile(theName.c_str(), "RECREATE");
3163 }
3164
3165 theFile->cd();
3166
3167 if (plotBasic_) {
3168 for (int k = 0; k < 5; ++k) {
3169 if (h_p[k] != 0) {
3170 TH1D *h1 = (TH1D *)h_p[k]->Clone();
3171 h1->Write();
3172 }
3173 if (h_eta[k] != 0) {
3174 TH1D *h2 = (TH1D *)h_eta[k]->Clone();
3175 h2->Write();
3176 }
3177 }
3178 for (unsigned int k = 0; k < h_eta0.size(); ++k) {
3179 if (h_eta0[k] != 0) {
3180 TH1D *h1 = (TH1D *)h_eta0[k]->Clone();
3181 h1->Write();
3182 }
3183 if (h_eta1[k] != 0) {
3184 TH1D *h2 = (TH1D *)h_eta1[k]->Clone();
3185 h2->Write();
3186 }
3187 if (h_eta2[k] != 0) {
3188 TH1D *h3 = (TH1D *)h_eta2[k]->Clone();
3189 h3->Write();
3190 }
3191 if (h_eta3[k] != 0) {
3192 TH1D *h4 = (TH1D *)h_eta3[k]->Clone();
3193 h4->Write();
3194 }
3195 if (h_eta4[k] != 0) {
3196 TH1D *h5 = (TH1D *)h_eta4[k]->Clone();
3197 h5->Write();
3198 }
3199 if (h_dL1[k] != 0) {
3200 TH1D *h6 = (TH1D *)h_dL1[k]->Clone();
3201 h6->Write();
3202 }
3203 if (h_vtx[k] != 0) {
3204 TH1D *h7 = (TH1D *)h_vtx[k]->Clone();
3205 h7->Write();
3206 }
3207 }
3208 }
3209
3210 if (plotEnergy_) {
3211 for (unsigned int k = 0; k < ps_.size() - 1; ++k) {
3212 for (unsigned int j = 0; j <= (unsigned int)(etahi_ - etalo_); ++j) {
3213 if (h_etaEH[k].size() > j && h_etaEH[k][j] != nullptr && (all || (k == kp50))) {
3214 TH1D *hist = (TH1D *)h_etaEH[k][j]->Clone();
3215 hist->Write();
3216 }
3217 if (h_etaEp[k].size() > j && h_etaEp[k][j] != nullptr && (all || (k == kp50))) {
3218 TH1D *hist = (TH1D *)h_etaEp[k][j]->Clone();
3219 hist->Write();
3220 }
3221 if (h_etaEE[k].size() > j && h_etaEE[k][j] != nullptr && (all || (k == kp50))) {
3222 TH1D *hist = (TH1D *)h_etaEE[k][j]->Clone();
3223 hist->Write();
3224 }
3225 }
3226 }
3227
3228 for (unsigned int j = 0; j < ietas_.size(); ++j) {
3229 if (h_eHcal.size() > j && (h_eHcal[j] != nullptr)) {
3230 TH1D *hist = (TH1D *)h_eHcal[j]->Clone();
3231 hist->Write();
3232 }
3233 if (h_mom.size() > j && (h_mom[j] != nullptr)) {
3234 TH1D *hist = (TH1D *)h_mom[j]->Clone();
3235 hist->Write();
3236 }
3237 if (h_eEcal.size() > j && (h_eEcal[j] != nullptr)) {
3238 TH1D *hist = (TH1D *)h_eEcal[j]->Clone();
3239 hist->Write();
3240 }
3241 }
3242 }
3243
3244 if (plotHists_) {
3245 if (h_nvtx != 0) {
3246 TH1D *h1 = (TH1D *)h_nvtx->Clone();
3247 h1->Write();
3248 }
3249 if (h_etaE != 0) {
3250 TH2D *h2 = (TH2D *)h_etaE->Clone();
3251 h2->Write();
3252 }
3253 for (unsigned int i = 0; i < ndepth; ++i) {
3254 if (h_bvlist[i] != 0) {
3255 TH1D *h1 = (TH1D *)h_bvlist[i]->Clone();
3256 h1->Write();
3257 }
3258 if (h_bvlist2[i] != 0) {
3259 TH1D *h2 = (TH1D *)h_bvlist2[i]->Clone();
3260 h2->Write();
3261 }
3262 if (h_bvlist3[i] != 0) {
3263 TH1D *h3 = (TH1D *)h_bvlist3[i]->Clone();
3264 h3->Write();
3265 }
3266 if (h_evlist[i] != 0) {
3267 TH1D *h4 = (TH1D *)h_evlist[i]->Clone();
3268 h4->Write();
3269 }
3270 if (h_evlist2[i] != 0) {
3271 TH1D *h5 = (TH1D *)h_evlist2[i]->Clone();
3272 h5->Write();
3273 }
3274 if (h_evlist3[i] != 0) {
3275 TH1D *h6 = (TH1D *)h_evlist3[i]->Clone();
3276 h6->Write();
3277 }
3278 }
3279 }
3280 std::cout << "All done" << std::endl;
3281 theFile->Close();
3282 }
3283
3284 void CalibPlotProperties::correctEnergy(double &eHcal, const Long64_t &entry) {
3285 double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
3286 if ((ifDepth_ == 3) && (cFactor_ != nullptr)) {
3287 double cfac = cFactor_->getCorr(entry);
3288 eHcal *= cfac;
3289 } else if ((corrPU_ < 0) && (pmom > 0)) {
3290 double ediff = (t_eHcal30 - t_eHcal10);
3291 if (t_DetIds1 != 0 && t_DetIds3 != 0) {
3292 double Etot1(0), Etot3(0);
3293
3294 for (unsigned int idet = 0; idet < (*t_DetIds1).size(); idet++) {
3295 unsigned int id = truncateId((*t_DetIds1)[idet], truncateFlag_, false);
3296 double cfac = corrFactor_->getCorr(id);
3297 if (cFactor_ != 0)
3298 cfac *= cFactor_->getCorr(t_Run, (*t_DetIds1)[idet]);
3299 double hitEn = cfac * (*t_HitEnergies1)[idet];
3300 Etot1 += hitEn;
3301 }
3302 for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
3303 unsigned int id = truncateId((*t_DetIds3)[idet], truncateFlag_, false);
3304 double cfac = corrFactor_->getCorr(id);
3305 if (cFactor_ != 0)
3306 cfac *= cFactor_->getCorr(t_Run, (*t_DetIds3)[idet]);
3307 double hitEn = cfac * (*t_HitEnergies3)[idet];
3308 Etot3 += hitEn;
3309 }
3310 ediff = (Etot3 - Etot1);
3311 }
3312 double fac = puFactor(-corrPU_, t_ieta, pmom, eHcal, ediff);
3313 eHcal *= fac;
3314 } else if (corrPU_ > 0) {
3315 eHcal = puFactorRho(corrPU_, t_ieta, t_rhoh, eHcal);
3316 }
3317 }