File indexing completed on 2025-01-04 00:29:12
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138
0139
0140
0141
0142
0143
0144
0145
0146
0147 #include <TROOT.h>
0148 #include <TChain.h>
0149 #include <TFile.h>
0150 #include <TH1D.h>
0151 #include <TH2F.h>
0152 #include <TProfile.h>
0153 #include <TFitResult.h>
0154 #include <TFitResultPtr.h>
0155 #include <TStyle.h>
0156 #include <TCanvas.h>
0157 #include <TLegend.h>
0158 #include <TPaveStats.h>
0159 #include <TPaveText.h>
0160
0161 #include <algorithm>
0162 #include <iomanip>
0163 #include <iostream>
0164 #include <fstream>
0165 #include <map>
0166 #include <vector>
0167 #include <string>
0168
0169 void unpackDetId(unsigned int, int &, int &, int &, int &, int &);
0170 #include "CalibCorr.C"
0171
0172 class CalibMonitor {
0173 public:
0174 TChain *fChain;
0175 Int_t fCurrent;
0176
0177
0178 Int_t t_Run;
0179 Int_t t_Event;
0180 Int_t t_DataType;
0181 Int_t t_ieta;
0182 Int_t t_iphi;
0183 Double_t t_EventWeight;
0184 Int_t t_nVtx;
0185 Int_t t_nTrk;
0186 Int_t t_goodPV;
0187 Double_t t_l1pt;
0188 Double_t t_l1eta;
0189 Double_t t_l1phi;
0190 Double_t t_l3pt;
0191 Double_t t_l3eta;
0192 Double_t t_l3phi;
0193 Double_t t_p;
0194 Double_t t_pt;
0195 Double_t t_phi;
0196 Double_t t_mindR1;
0197 Double_t t_mindR2;
0198 Double_t t_eMipDR;
0199 Double_t t_eHcal;
0200 Double_t t_eHcal10;
0201 Double_t t_eHcal30;
0202 Double_t t_hmaxNearP;
0203 Double_t t_rhoh;
0204 Bool_t t_selectTk;
0205 Bool_t t_qltyFlag;
0206 Bool_t t_qltyMissFlag;
0207 Bool_t t_qltyPVFlag;
0208 Double_t t_gentrackP;
0209 std::vector<unsigned int> *t_DetIds;
0210 std::vector<double> *t_HitEnergies;
0211 std::vector<bool> *t_trgbits;
0212 std::vector<unsigned int> *t_DetIds1;
0213 std::vector<unsigned int> *t_DetIds3;
0214 std::vector<double> *t_HitEnergies1;
0215 std::vector<double> *t_HitEnergies3;
0216
0217
0218 TBranch *b_t_Run;
0219 TBranch *b_t_Event;
0220 TBranch *b_t_DataType;
0221 TBranch *b_t_ieta;
0222 TBranch *b_t_iphi;
0223 TBranch *b_t_EventWeight;
0224 TBranch *b_t_nVtx;
0225 TBranch *b_t_nTrk;
0226 TBranch *b_t_goodPV;
0227 TBranch *b_t_l1pt;
0228 TBranch *b_t_l1eta;
0229 TBranch *b_t_l1phi;
0230 TBranch *b_t_l3pt;
0231 TBranch *b_t_l3eta;
0232 TBranch *b_t_l3phi;
0233 TBranch *b_t_p;
0234 TBranch *b_t_pt;
0235 TBranch *b_t_phi;
0236 TBranch *b_t_mindR1;
0237 TBranch *b_t_mindR2;
0238 TBranch *b_t_eMipDR;
0239 TBranch *b_t_eHcal;
0240 TBranch *b_t_eHcal10;
0241 TBranch *b_t_eHcal30;
0242 TBranch *b_t_hmaxNearP;
0243 TBranch *b_t_rhoh;
0244 TBranch *b_t_selectTk;
0245 TBranch *b_t_qltyFlag;
0246 TBranch *b_t_qltyMissFlag;
0247 TBranch *b_t_qltyPVFlag;
0248 TBranch *b_t_gentrackP;
0249 TBranch *b_t_DetIds;
0250 TBranch *b_t_HitEnergies;
0251 TBranch *b_t_trgbits;
0252 TBranch *b_t_DetIds1;
0253 TBranch *b_t_DetIds3;
0254 TBranch *b_t_HitEnergies1;
0255 TBranch *b_t_HitEnergies3;
0256
0257 struct counter {
0258 static const int npsize = 4;
0259 counter() {
0260 total = 0;
0261 for (int k = 0; k < npsize; ++k)
0262 count[k] = 0;
0263 };
0264 unsigned int total, count[npsize];
0265 };
0266
0267 CalibMonitor(const char *fname,
0268 const char *dirname,
0269 const char *dupFileName,
0270 const char *comFileName,
0271 const char *outFileName,
0272 const std::string &prefix = "",
0273 const char *corrFileName = "",
0274 const char *rcorFileName = "",
0275 int puCorr = -8,
0276 int flag = 1031,
0277 int numb = 50,
0278 bool datMC = true,
0279 int truncateFlag = 0,
0280 bool useGen = false,
0281 double scale = 1.0,
0282 int useScale = 0,
0283 int etalo = 0,
0284 int etahi = 30,
0285 int runlo = 0,
0286 int runhi = 99999999,
0287 int phimin = 1,
0288 int phimax = 72,
0289 int zside = 1,
0290 int nvxlo = 0,
0291 int nvxhi = 1000,
0292 const char *rbxFile = "",
0293 bool exclude = false,
0294 bool etamax = false);
0295 virtual ~CalibMonitor();
0296 virtual Int_t Cut(Long64_t entry);
0297 virtual Int_t GetEntry(Long64_t entry);
0298 virtual Long64_t LoadTree(Long64_t entry);
0299 virtual void Init(TChain *, const char *, const char *);
0300 virtual void Loop(Long64_t nmax = -1, bool debug = false);
0301 virtual Bool_t Notify();
0302 virtual void Show(Long64_t entry = -1);
0303 bool goodTrack(double &eHcal, double &cut, const Long64_t &entry, bool debug);
0304 bool selectPhi(bool debug);
0305 void plotHist(int type, int num, bool save = false);
0306 template <class Hist>
0307 void drawHist(Hist *, TCanvas *);
0308 void savePlot(const std::string &theName, bool append = true, bool all = false);
0309 void correctEnergy(double &ener, const Long64_t &entry);
0310
0311 private:
0312 static const unsigned int npbin = 7, kp50 = 3;
0313 CalibCorrFactor *corrFactor_;
0314 CalibCorr *cFactor_;
0315 CalibSelectRBX *cSelect_;
0316 CalibDuplicate *cDuplicate_;
0317 const std::string fname_, dirnm_, prefix_, outFileName_;
0318 const int corrPU_, flag_, numb_;
0319 const bool isRealData_, useGen_;
0320 const int truncateFlag_;
0321 const int etalo_, etahi_;
0322 int runlo_, runhi_;
0323 const int phimin_, phimax_, zside_, nvxlo_, nvxhi_;
0324 const char *rbxFile_;
0325 bool exclude_, cutL1T_, selRBX_;
0326 bool includeRun_;
0327 int coarseBin_, plotType_;
0328 int flexibleSelect_, ifDepth_, duplicate_, thrForm_;
0329 double log2by18_;
0330 std::ofstream fileout_;
0331 std::vector<std::pair<int, int> > events_;
0332 std::vector<double> etas_, ps_, dl1_;
0333 std::vector<int> nvx_, ietasL_, ietasH_;
0334 std::vector<TH1D *> h_rbx, h_etaF[npbin], h_etaB[npbin];
0335 std::vector<TProfile *> h_etaX[npbin];
0336 std::vector<TH1D *> h_etaR[npbin], h_nvxR[npbin], h_dL1R[npbin];
0337 std::vector<TH1D *> h_pp[npbin];
0338 std::vector<TH1D *> h_p;
0339 };
0340
0341 CalibMonitor::CalibMonitor(const char *fname,
0342 const char *dirnm,
0343 const char *dupFileName,
0344 const char *comFileName,
0345 const char *outFName,
0346 const std::string &prefix,
0347 const char *corrFileName,
0348 const char *rcorFileName,
0349 int puCorr,
0350 int flag,
0351 int numb,
0352 bool isRealData,
0353 int truncate,
0354 bool useGen,
0355 double scale,
0356 int useScale,
0357 int etalo,
0358 int etahi,
0359 int runlo,
0360 int runhi,
0361 int phimin,
0362 int phimax,
0363 int zside,
0364 int nvxlo,
0365 int nvxhi,
0366 const char *rbxFile,
0367 bool exc,
0368 bool etam)
0369 : corrFactor_(nullptr),
0370 cFactor_(nullptr),
0371 cSelect_(nullptr),
0372 cDuplicate_(nullptr),
0373 fname_(std::string(fname)),
0374 dirnm_(std::string(dirnm)),
0375 prefix_(prefix),
0376 outFileName_(std::string(outFName)),
0377 corrPU_(puCorr),
0378 flag_(flag),
0379 numb_(numb),
0380 isRealData_(isRealData),
0381 useGen_(useGen),
0382 truncateFlag_(truncate),
0383 etalo_(etalo),
0384 etahi_(etahi),
0385 runlo_(runlo),
0386 runhi_(runhi),
0387 phimin_(phimin),
0388 phimax_(phimax),
0389 zside_(zside),
0390 nvxlo_(nvxlo),
0391 nvxhi_(nvxhi),
0392 rbxFile_(rbxFile),
0393 exclude_(exc),
0394 includeRun_(true) {
0395
0396
0397
0398 plotType_ = ((flag_ / 10) % 10);
0399 if (plotType_ < 0 || plotType_ > 3)
0400 plotType_ = 3;
0401 flexibleSelect_ = (((flag_ / 1) % 10));
0402 int oneplace = ((flag_ / 1000) % 10);
0403 cutL1T_ = (oneplace % 2);
0404 bool marina = ((oneplace / 2) % 2);
0405 ifDepth_ = ((flag_ / 10000) % 10);
0406 selRBX_ = (((flag_ / 100000) % 10) > 0);
0407 duplicate_ = ((flag_ / 1000000) % 10);
0408 coarseBin_ = ((flag_ / 10000000) % 10);
0409 log2by18_ = std::log(2.5) / 18.0;
0410 if (runlo_ < 0 || runhi_ < 0) {
0411 runlo_ = std::abs(runlo_);
0412 runhi_ = std::abs(runhi_);
0413 includeRun_ = false;
0414 }
0415 int useScale0 = useScale % 10;
0416 thrForm_ = useScale / 10;
0417 char treeName[400];
0418 sprintf(treeName, "%s/CalibTree", dirnm_.c_str());
0419 TChain *chain = new TChain(treeName);
0420 std::cout << "Create a chain for " << treeName << " from " << fname << " flags " << flexibleSelect_ << "|"
0421 << plotType_ << "|" << corrPU_ << "\n cons " << log2by18_ << " eta range " << etalo_ << ":" << etahi_
0422 << " run range " << runlo_ << ":" << runhi_ << " (inclusion flag " << includeRun_ << ")\n Selection of RBX "
0423 << selRBX_ << " Vertex Range " << nvxlo_ << ":" << nvxhi_ << "\n corrFileName: " << corrFileName
0424 << " useScale " << useScale << ":" << scale << ":" << etam << "\n rcorFileName: " << rcorFileName
0425 << " flag " << ifDepth_ << ":" << cutL1T_ << ":" << marina << " Threshold Flag " << thrForm_ << std::endl;
0426 if (!fillChain(chain, fname)) {
0427 std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
0428 } else {
0429 std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
0430 corrFactor_ = new CalibCorrFactor(corrFileName, useScale0, scale, etam, marina, false);
0431 Init(chain, comFileName, outFName);
0432 if (std::string(dupFileName) != "")
0433 cDuplicate_ = new CalibDuplicate(dupFileName, duplicate_, false);
0434 if (std::string(rcorFileName) != "") {
0435 cFactor_ = new CalibCorr(rcorFileName, ifDepth_, false);
0436 if (cFactor_->absent())
0437 ifDepth_ = -1;
0438 } else {
0439 ifDepth_ = -1;
0440 }
0441 if (std::string(rbxFile) != "")
0442 cSelect_ = new CalibSelectRBX(rbxFile, false);
0443 }
0444 }
0445
0446 CalibMonitor::~CalibMonitor() {
0447 delete corrFactor_;
0448 delete cFactor_;
0449 delete cSelect_;
0450 delete cDuplicate_;
0451 if (!fChain)
0452 return;
0453 delete fChain->GetCurrentFile();
0454 }
0455
0456 Int_t CalibMonitor::GetEntry(Long64_t entry) {
0457
0458 if (!fChain)
0459 return 0;
0460 return fChain->GetEntry(entry);
0461 }
0462
0463 Long64_t CalibMonitor::LoadTree(Long64_t entry) {
0464
0465 if (!fChain)
0466 return -5;
0467 Long64_t centry = fChain->LoadTree(entry);
0468 if (centry < 0)
0469 return centry;
0470 if (!fChain->InheritsFrom(TChain::Class()))
0471 return centry;
0472 TChain *chain = (TChain *)fChain;
0473 if (chain->GetTreeNumber() != fCurrent) {
0474 fCurrent = chain->GetTreeNumber();
0475 Notify();
0476 }
0477 return centry;
0478 }
0479
0480 void CalibMonitor::Init(TChain *tree, const char *comFileName, const char *outFileName) {
0481
0482
0483
0484
0485
0486
0487
0488
0489
0490 t_DetIds = 0;
0491 t_DetIds1 = 0;
0492 t_DetIds3 = 0;
0493 t_HitEnergies = 0;
0494 t_HitEnergies1 = 0;
0495 t_HitEnergies3 = 0;
0496 t_trgbits = 0;
0497
0498 fChain = tree;
0499 fCurrent = -1;
0500 if (!tree)
0501 return;
0502 fChain->SetMakeClass(1);
0503
0504 fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
0505 fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
0506 fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
0507 fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
0508 fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
0509 fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
0510 fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
0511 fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
0512 fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
0513 fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
0514 fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
0515 fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
0516 fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
0517 fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
0518 fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
0519 fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
0520 fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
0521 fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
0522 fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
0523 fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
0524 fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
0525 fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
0526 fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
0527 fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
0528 fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
0529 fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
0530 fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
0531 fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
0532 fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
0533 fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
0534 fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
0535 fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
0536 fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
0537 fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
0538 fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
0539 fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
0540 fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
0541 fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
0542 Notify();
0543
0544 if (strcmp(comFileName, "") != 0) {
0545 std::ifstream infil2(comFileName);
0546 if (!infil2.is_open()) {
0547 std::cout << "Cannot open selection file " << comFileName << std::endl;
0548 } else {
0549 while (1) {
0550 int irun, ievt;
0551 infil2 >> irun >> ievt;
0552 if (!infil2.good())
0553 break;
0554 std::pair<int, int> key(irun, ievt);
0555 events_.push_back(key);
0556 }
0557 infil2.close();
0558 std::cout << "Reads a list of " << events_.size() << " events from " << comFileName << std::endl;
0559 }
0560 } else {
0561 std::cout << "No event list provided for selection" << std::endl;
0562 }
0563
0564 if (((flag_ / 100) % 10) > 0) {
0565 if (((flag_ / 100) % 10) == 2) {
0566 fileout_.open(outFileName, std::ofstream::out);
0567 std::cout << "Opens " << outFileName << " in output mode" << std::endl;
0568 } else {
0569 fileout_.open(outFileName, std::ofstream::app);
0570 std::cout << "Opens " << outFileName << " in append mode" << std::endl;
0571 }
0572 fileout_ << "Input file: " << fname_ << " Directory: " << dirnm_ << " Prefix: " << prefix_ << std::endl;
0573 }
0574
0575 double xbins[99];
0576 int nbins(-1);
0577 if (plotType_ == 0) {
0578 double xbin[9] = {-21.0, -16.0, -12.0, -6.0, 0.0, 6.0, 12.0, 16.0, 21.0};
0579 for (int i = 0; i < 9; ++i) {
0580 etas_.push_back(xbin[i]);
0581 xbins[i] = xbin[i];
0582 }
0583 nbins = 8;
0584 } else if (plotType_ == 1) {
0585 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};
0586 for (int i = 0; i < 11; ++i) {
0587 etas_.push_back(xbin[i]);
0588 xbins[i] = xbin[i];
0589 }
0590 nbins = 10;
0591 } else if (plotType_ == 2) {
0592 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,
0593 3.0, 5.0, 7.0, 9.0, 11.0, 13.0, 15.0, 17.0, 19.0, 21.0, 23.0};
0594 for (int i = 0; i < 23; ++i) {
0595 etas_.push_back(xbin[i]);
0596 xbins[i] = xbin[i];
0597 }
0598 nbins = 22;
0599 } else {
0600 double xbina[99];
0601 int neta = numb_ / 2;
0602 for (int k = 0; k < neta; ++k) {
0603 xbina[k] = (k - neta) - 0.5;
0604 xbina[numb_ - k] = (neta - k) + 0.5;
0605 }
0606 xbina[neta] = 0;
0607 for (int i = 0; i < numb_ + 1; ++i) {
0608 etas_.push_back(xbina[i]);
0609 xbins[i] = xbina[i];
0610 ++nbins;
0611 }
0612 }
0613 int ipbin[npbin] = {10, 20, 30, 40, 60, 100, 500};
0614 for (unsigned int i = 0; i < npbin; ++i)
0615 ps_.push_back((double)(ipbin[i]));
0616 int npvtx[6] = {0, 7, 10, 13, 16, 100};
0617 for (int i = 0; i < 6; ++i)
0618 nvx_.push_back(npvtx[i]);
0619 double dl1s[9] = {0, 0.10, 0.20, 0.50, 1.0, 2.0, 2.5, 3.0, 10.0};
0620 int ietasL[4] = {0, 13, 17, 0};
0621 int ietasH[4] = {14, 18, 25, 25};
0622 for (int i = 0; i < 4; ++i) {
0623 ietasL_.push_back(ietasL[i]);
0624 ietasH_.push_back(ietasH[i]);
0625 }
0626 int nxbin(100);
0627 double xlow(0.0), xhigh(5.0);
0628 if (coarseBin_ == 1) {
0629 xlow = 0.25;
0630 xhigh = 5.25;
0631 nxbin = 50;
0632 } else if (coarseBin_ > 1) {
0633 if (coarseBin_ == 2)
0634 nxbin = 500;
0635 else
0636 nxbin = 1000;
0637 }
0638
0639 char name[100], title[500];
0640 std::string titl[5] = {
0641 "All tracks", "Good quality tracks", "Selected tracks", "Tracks with charge isolation", "Tracks MIP in ECAL"};
0642 for (int i = 0; i < 9; ++i)
0643 dl1_.push_back(dl1s[i]);
0644 unsigned int kp = (ps_.size() - 1);
0645 for (unsigned int k = 0; k < kp; ++k) {
0646 for (unsigned int j = 0; j <= ietasL_.size(); ++j) {
0647 sprintf(name, "%spp%d%d", prefix_.c_str(), k, j);
0648 if (j == 0)
0649 sprintf(title, "E/p for %s (p = %d:%d GeV All)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0650 else
0651 sprintf(title,
0652 "E/p for %s (p = %d:%d GeV |#eta| %d:%d)",
0653 titl[4].c_str(),
0654 ipbin[k],
0655 ipbin[k + 1],
0656 (ietasL_[j - 1] + 1),
0657 ietasH_[j - 1]);
0658 h_pp[k].push_back(new TH1D(name, title, 100, 10.0, 110.0));
0659 int kk = h_pp[k].size() - 1;
0660 h_pp[k][kk]->Sumw2();
0661 }
0662 }
0663 if (plotType_ <= 1) {
0664 std::cout << "Book Histos for Standard" << std::endl;
0665 for (unsigned int k = 0; k < kp; ++k) {
0666 for (unsigned int i = 0; i < nvx_.size(); ++i) {
0667 sprintf(name, "%setaX%d%d", prefix_.c_str(), k, i);
0668 if (i == 0) {
0669 sprintf(title, "%s (p = %d:%d GeV all vertices)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0670 } else {
0671 sprintf(
0672 title, "%s (p = %d:%d GeV # Vtx %d:%d)", titl[4].c_str(), ipbin[k], ipbin[k + 1], nvx_[i - 1], nvx_[i]);
0673 }
0674 h_etaX[k].push_back(new TProfile(name, title, nbins, xbins));
0675 unsigned int kk = h_etaX[k].size() - 1;
0676 h_etaX[k][kk]->Sumw2();
0677 sprintf(name, "%snvxR%d%d", prefix_.c_str(), k, i);
0678 if (i == 0) {
0679 sprintf(title, "E/p for %s (p = %d:%d GeV all vertices)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0680 } else {
0681 sprintf(title,
0682 "E/p for %s (p = %d:%d GeV # Vtx %d:%d)",
0683 titl[4].c_str(),
0684 ipbin[k],
0685 ipbin[k + 1],
0686 nvx_[i - 1],
0687 nvx_[i]);
0688 }
0689 h_nvxR[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0690 kk = h_nvxR[k].size() - 1;
0691 h_nvxR[k][kk]->Sumw2();
0692 }
0693 for (unsigned int j = 0; j < etas_.size(); ++j) {
0694 sprintf(name, "%sratio%d%d", prefix_.c_str(), k, j);
0695 if (j == 0) {
0696 sprintf(title, "E/p for %s (p = %d:%d GeV)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0697 } else {
0698 sprintf(title,
0699 "E/p for %s (p = %d:%d GeV #eta %4.1f:%4.1f)",
0700 titl[4].c_str(),
0701 ipbin[k],
0702 ipbin[k + 1],
0703 etas_[j - 1],
0704 etas_[j]);
0705 }
0706 h_etaF[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0707 unsigned int kk = h_etaF[k].size() - 1;
0708 h_etaF[k][kk]->Sumw2();
0709 sprintf(name, "%setaR%d%d", prefix_.c_str(), k, j);
0710 h_etaR[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0711 kk = h_etaR[k].size() - 1;
0712 h_etaR[k][kk]->Sumw2();
0713 }
0714 for (unsigned int j = 1; j <= ietasL_.size(); ++j) {
0715 sprintf(name, "%setaB%d%d", prefix_.c_str(), k, j);
0716 sprintf(title,
0717 "E/p for %s (p = %d:%d GeV |#eta| %d:%d)",
0718 titl[4].c_str(),
0719 ipbin[k],
0720 ipbin[k + 1],
0721 (ietasL_[j - 1] + 1),
0722 ietasH_[j - 1]);
0723 h_etaB[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0724 unsigned int kk = h_etaB[k].size() - 1;
0725 h_etaB[k][kk]->Sumw2();
0726 }
0727 for (unsigned int j = 0; j < dl1_.size(); ++j) {
0728 sprintf(name, "%sdl1R%d%d", prefix_.c_str(), k, j);
0729 if (j == 0) {
0730 sprintf(title, "E/p for %s (p = %d:%d GeV All d_{L1})", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0731 } else {
0732 sprintf(title,
0733 "E/p for %s (p = %d:%d GeV d_{L1} %4.2f:%4.2f)",
0734 titl[4].c_str(),
0735 ipbin[k],
0736 ipbin[k + 1],
0737 dl1_[j - 1],
0738 dl1_[j]);
0739 }
0740 h_dL1R[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0741 unsigned int kk = h_dL1R[k].size() - 1;
0742 h_dL1R[k][kk]->Sumw2();
0743 }
0744 }
0745 for (unsigned int i = 0; i < nvx_.size(); ++i) {
0746 sprintf(name, "%setaX%d%d", prefix_.c_str(), kp, i);
0747 if (i == 0) {
0748 sprintf(title, "%s (All Momentum all vertices)", titl[4].c_str());
0749 } else {
0750 sprintf(title, "%s (All Momentum # Vtx %d:%d)", titl[4].c_str(), nvx_[i - 1], nvx_[i]);
0751 }
0752 h_etaX[npbin - 1].push_back(new TProfile(name, title, nbins, xbins));
0753 unsigned int kk = h_etaX[npbin - 1].size() - 1;
0754 h_etaX[npbin - 1][kk]->Sumw2();
0755 sprintf(name, "%snvxR%d%d", prefix_.c_str(), kp, i);
0756 if (i == 0) {
0757 sprintf(title, "E/p for %s (All Momentum all vertices)", titl[4].c_str());
0758 } else {
0759 sprintf(title, "E/p for %s (All Momentum # Vtx %d:%d)", titl[4].c_str(), nvx_[i - 1], nvx_[i]);
0760 }
0761 h_nvxR[npbin - 1].push_back(new TH1D(name, title, 200, 0., 10.));
0762 kk = h_nvxR[npbin - 1].size() - 1;
0763 h_nvxR[npbin - 1][kk]->Sumw2();
0764 }
0765 for (unsigned int j = 0; j < etas_.size(); ++j) {
0766 sprintf(name, "%sratio%d%d", prefix_.c_str(), kp, j);
0767 if (j == 0) {
0768 sprintf(title, "E/p for %s (All momentum)", titl[4].c_str());
0769 } else {
0770 sprintf(title, "E/p for %s (All momentum #eta %4.1f:%4.1f)", titl[4].c_str(), etas_[j - 1], etas_[j]);
0771 }
0772 h_etaF[npbin - 1].push_back(new TH1D(name, title, 200, 0., 10.));
0773 unsigned int kk = h_etaF[npbin - 1].size() - 1;
0774 h_etaF[npbin - 1][kk]->Sumw2();
0775 sprintf(name, "%setaR%d%d", prefix_.c_str(), kp, j);
0776 h_etaR[npbin - 1].push_back(new TH1D(name, title, 200, 0., 10.));
0777 kk = h_etaR[npbin - 1].size() - 1;
0778 h_etaR[npbin - 1][kk]->Sumw2();
0779 }
0780 for (unsigned int j = 1; j <= ietasL_.size(); ++j) {
0781 sprintf(name, "%setaB%d%d", prefix_.c_str(), kp, j);
0782 sprintf(title, "E/p for %s (All momentum |#eta| %d:%d)", titl[4].c_str(), (ietasL_[j - 1] + 1), ietasH_[j - 1]);
0783 h_etaB[npbin - 1].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0784 unsigned int kk = h_etaB[npbin - 1].size() - 1;
0785 h_etaB[npbin - 1][kk]->Sumw2();
0786 }
0787 for (unsigned int j = 0; j < dl1_.size(); ++j) {
0788 sprintf(name, "%sdl1R%d%d", prefix_.c_str(), kp, j);
0789 if (j == 0) {
0790 sprintf(title, "E/p for %s (All momentum)", titl[4].c_str());
0791 } else {
0792 sprintf(title, "E/p for %s (All momentum d_{L1} %4.2f:%4.2f)", titl[4].c_str(), dl1_[j - 1], dl1_[j]);
0793 }
0794 h_dL1R[npbin - 1].push_back(new TH1D(name, title, 200, 0., 10.));
0795 unsigned int kk = h_dL1R[npbin - 1].size() - 1;
0796 h_dL1R[npbin - 1][kk]->Sumw2();
0797 }
0798 } else {
0799 std::cout << "Book Histos for Non-Standard " << etas_.size() << ":" << kp50 << std::endl;
0800 unsigned int kp = (ps_.size() - 1);
0801 for (unsigned int k = 0; k < kp; ++k) {
0802 for (unsigned int j = 0; j < etas_.size(); ++j) {
0803 sprintf(name, "%sratio%d%d", prefix_.c_str(), k, j);
0804 if (j == 0) {
0805 sprintf(title, "E/p for %s (p = %d:%d GeV)", titl[4].c_str(), ipbin[k], ipbin[k + 1]);
0806 } else {
0807 sprintf(title,
0808 "E/p for %s (p = %d:%d GeV #eta %4.1f:%4.1f)",
0809 titl[4].c_str(),
0810 ipbin[k],
0811 ipbin[k + 1],
0812 etas_[j - 1],
0813 etas_[j]);
0814 }
0815 h_etaF[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0816 unsigned int kk = h_etaF[k].size() - 1;
0817 h_etaF[k][kk]->Sumw2();
0818 }
0819 for (unsigned int j = 1; j <= ietasL_.size(); ++j) {
0820 sprintf(name, "%setaB%d%d", prefix_.c_str(), k, j);
0821 sprintf(title,
0822 "E/p for %s (p = %d:%d GeV |#eta| %d:%d)",
0823 titl[4].c_str(),
0824 ipbin[k],
0825 ipbin[k + 1],
0826 (ietasL_[j - 1] + 1),
0827 ietasH_[j - 1]);
0828 h_etaB[k].push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0829 unsigned int kk = h_etaB[k].size() - 1;
0830 h_etaB[k][kk]->Sumw2();
0831 }
0832 }
0833 }
0834 if (selRBX_) {
0835 for (unsigned int j = 1; j <= 18; ++j) {
0836 sprintf(name, "%sRBX%d%d", prefix_.c_str(), kp50, j);
0837 sprintf(title, "E/p for RBX%d (p = %d:%d GeV |#eta| %d:%d)", j, ipbin[kp50], ipbin[kp50 + 1], etalo_, etahi_);
0838 h_rbx.push_back(new TH1D(name, title, nxbin, xlow, xhigh));
0839 h_rbx[j - 1]->Sumw2();
0840 }
0841 }
0842 for (unsigned int j = 1; j < npbin; ++j) {
0843 sprintf(name, "%sp%d", prefix_.c_str(), j);
0844 sprintf(title, "Momentum (GeV) of selected track (p = %d:%d GeV)", ipbin[j], ipbin[j + 1]);
0845 h_p.push_back(new TH1D(name, title, 100, ipbin[j], ipbin[j + 1]));
0846 h_p[j - 1]->Sumw2();
0847 }
0848 }
0849
0850 Bool_t CalibMonitor::Notify() {
0851
0852
0853
0854
0855
0856
0857 return kTRUE;
0858 }
0859
0860 void CalibMonitor::Show(Long64_t entry) {
0861
0862
0863 if (!fChain)
0864 return;
0865 fChain->Show(entry);
0866 }
0867
0868 Int_t CalibMonitor::Cut(Long64_t) {
0869
0870
0871
0872 return 1;
0873 }
0874
0875 void CalibMonitor::Loop(Long64_t nmax, bool debug) {
0876
0877
0878
0879
0880
0881
0882
0883
0884
0885
0886
0887
0888
0889
0890
0891
0892
0893
0894
0895
0896
0897
0898
0899
0900 if (fChain == 0)
0901 return;
0902 std::map<int, counter> runSum, runEn1, runEn2;
0903 if (debug) {
0904 std::cout << etas_.size() << " Eta Bins:";
0905 for (unsigned int j = 0; j < etas_.size(); ++j)
0906 std::cout << " " << etas_[j];
0907 std::cout << std::endl;
0908 }
0909
0910 Long64_t nentries = fChain->GetEntriesFast();
0911 Long64_t entries = (nmax > 0) ? nmax : nentries;
0912 std::cout << "Total entries " << nentries << ":" << entries << std::endl;
0913 Long64_t nbytes(0), nb(0);
0914 unsigned int duplicate(0), good(0), kount(0);
0915 unsigned int kp1 = ps_.size() - 1;
0916 unsigned int kv1 = 0;
0917 std::vector<int> kounts(kp1, 0);
0918 std::vector<int> kount50(20, 0);
0919 std::vector<int> kount0(20, 0);
0920 std::vector<int> kount1(20, 0);
0921 std::vector<int> kount2(20, 0);
0922 std::vector<int> kount3(20, 0);
0923 std::vector<int> kount4(20, 0);
0924 std::vector<int> kount5(20, 0);
0925 int32_t oddEven = (nmax == -2) ? 1 : ((nmax == -3) ? -1 : 0);
0926 for (Long64_t jentry = 0; jentry < entries; jentry++) {
0927
0928 Long64_t ientry = LoadTree(jentry);
0929 if (ientry < 0)
0930 break;
0931 nb = fChain->GetEntry(jentry);
0932 nbytes += nb;
0933 if (jentry % 1000000 == 0)
0934 std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
0935 if (oddEven != 0) {
0936 if ((oddEven < 0) && (jentry % 2 == 0))
0937 continue;
0938 else if ((oddEven > 0) && (jentry % 2 != 0))
0939 continue;
0940 }
0941 double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
0942 int kp(-1);
0943 for (unsigned int k = 1; k < ps_.size(); ++k) {
0944 if (pmom >= ps_[k - 1] && pmom < ps_[k]) {
0945 kp = k - 1;
0946 break;
0947 }
0948 }
0949 bool p4060 = ((pmom >= 40.0) && (pmom <= 60.0));
0950 if (p4060)
0951 ++kount50[0];
0952 if (kp == 0) {
0953 ++kount0[0];
0954 } else if (kp == 1) {
0955 ++kount1[0];
0956 } else if (kp == 2) {
0957 ++kount2[0];
0958 } else if (kp == 3) {
0959 ++kount3[0];
0960 } else if (kp == 4) {
0961 ++kount4[0];
0962 } else if (kp == 5) {
0963 ++kount5[0];
0964 }
0965 bool select = ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr(0))) ? (cDuplicate_->isDuplicate(jentry)) : true;
0966 if (!select) {
0967 ++duplicate;
0968 if (debug)
0969 std::cout << "Duplicate event " << t_Run << " " << t_Event << " " << t_p << std::endl;
0970 continue;
0971 }
0972 if (p4060)
0973 ++kount50[1];
0974 if (kp == 0) {
0975 ++kount0[1];
0976 } else if (kp == 1) {
0977 ++kount1[1];
0978 } else if (kp == 2) {
0979 ++kount2[1];
0980 } else if (kp == 3) {
0981 ++kount3[1];
0982 } else if (kp == 4) {
0983 ++kount4[1];
0984 } else if (kp == 5) {
0985 ++kount5[1];
0986 }
0987 bool selRun = (includeRun_ ? ((t_Run >= runlo_) && (t_Run <= runhi_)) : ((t_Run < runlo_) || (t_Run > runhi_)));
0988 if (select) {
0989 if (p4060)
0990 ++kount50[2];
0991 if (kp == 0) {
0992 ++kount0[2];
0993 } else if (kp == 1) {
0994 ++kount1[2];
0995 } else if (kp == 2) {
0996 ++kount2[2];
0997 } else if (kp == 3) {
0998 ++kount3[2];
0999 } else if (kp == 4) {
1000 ++kount4[2];
1001 } else if (kp == 5) {
1002 ++kount5[2];
1003 }
1004 }
1005 select =
1006 (selRun && (fabs(t_ieta) >= etalo_) && (fabs(t_ieta) <= etahi_) && (t_nVtx >= nvxlo_) && (t_nVtx <= nvxhi_));
1007 if (select) {
1008 if (p4060)
1009 ++kount50[3];
1010 if (kp == 0) {
1011 ++kount0[3];
1012 } else if (kp == 1) {
1013 ++kount1[3];
1014 } else if (kp == 2) {
1015 ++kount2[3];
1016 } else if (kp == 3) {
1017 ++kount3[3];
1018 } else if (kp == 4) {
1019 ++kount4[3];
1020 } else if (kp == 5) {
1021 ++kount5[3];
1022 }
1023 }
1024 if (!select) {
1025 if (debug)
1026 std::cout << "Run # " << t_Run << " out of range of " << runlo_ << ":" << runhi_ << " or ieta " << t_ieta
1027 << " (" << etalo_ << ":" << etahi_ << ") or nvtx " << t_nVtx << " (" << nvxlo_ << ":" << nvxhi_
1028 << ") out of range" << std::endl;
1029 continue;
1030 }
1031 if (cSelect_ != nullptr) {
1032 if (exclude_) {
1033 if (cSelect_->isItRBX(t_DetIds))
1034 continue;
1035 } else {
1036 if (!(cSelect_->isItRBX(t_ieta, t_iphi)))
1037 continue;
1038 }
1039 }
1040 if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr(2))) {
1041 if (cDuplicate_->select(t_ieta, t_iphi))
1042 continue;
1043 }
1044 if (p4060)
1045 ++kount50[4];
1046 if (kp == 0) {
1047 ++kount0[4];
1048 } else if (kp == 1) {
1049 ++kount1[4];
1050 } else if (kp == 2) {
1051 ++kount2[4];
1052 } else if (kp == 3) {
1053 ++kount3[4];
1054 } else if (kp == 4) {
1055 ++kount4[4];
1056 } else if (kp == 5) {
1057 ++kount5[4];
1058 }
1059 select = (!cutL1T_ || (t_mindR1 >= 0.5));
1060 if (!select) {
1061 if (debug)
1062 std::cout << "Reject Run # " << t_Run << " Event # " << t_Event << " too close to L1 trigger " << t_mindR1
1063 << std::endl;
1064 continue;
1065 }
1066 if (p4060)
1067 ++kount50[5];
1068 if (kp == 0) {
1069 ++kount0[5];
1070 } else if (kp == 1) {
1071 ++kount1[5];
1072 } else if (kp == 2) {
1073 ++kount2[5];
1074 } else if (kp == 3) {
1075 ++kount3[5];
1076 } else if (kp == 4) {
1077 ++kount4[5];
1078 } else if (kp == 5) {
1079 ++kount5[5];
1080 }
1081 select = ((events_.size() == 0) ||
1082 (std::find(events_.begin(), events_.end(), std::pair<int, int>(t_Run, t_Event)) != events_.end()));
1083 if (!select) {
1084 if (debug)
1085 std::cout << "Reject Run # " << t_Run << " Event # " << t_Event << " not in the selection list" << std::endl;
1086 continue;
1087 }
1088 if (p4060)
1089 ++kount50[6];
1090 if (kp == 0) {
1091 ++kount0[6];
1092 } else if (kp == 1) {
1093 ++kount1[6];
1094 } else if (kp == 2) {
1095 ++kount2[6];
1096 } else if (kp == 3) {
1097 ++kount3[6];
1098 } else if (kp == 4) {
1099 ++kount4[6];
1100 } else if (kp == 5) {
1101 ++kount5[6];
1102 }
1103 if ((ifDepth_ == 3) && (cFactor_ != nullptr) && (cFactor_->absent(ientry)))
1104 continue;
1105
1106
1107 int jp(-1), jp1(-1), jp2(-1);
1108 unsigned int kv = nvx_.size() - 1;
1109 for (unsigned int k = 1; k < nvx_.size(); ++k) {
1110 if (t_goodPV >= nvx_[k - 1] && t_goodPV < nvx_[k]) {
1111 kv = k;
1112 break;
1113 }
1114 }
1115 unsigned int kd1 = 0;
1116 unsigned int kd = dl1_.size() - 1;
1117 for (unsigned int k = 1; k < dl1_.size(); ++k) {
1118 if (t_mindR1 >= dl1_[k - 1] && t_mindR1 < dl1_[k]) {
1119 kd = k;
1120 break;
1121 }
1122 }
1123 double eta = (t_ieta > 0) ? ((double)(t_ieta)-0.001) : ((double)(t_ieta) + 0.001);
1124 for (unsigned int j = 1; j < etas_.size(); ++j) {
1125 if (eta > etas_[j - 1] && eta < etas_[j]) {
1126 jp = j;
1127 break;
1128 }
1129 }
1130 for (unsigned int j = 0; j < (ietasL_.size() - 1); ++j) {
1131 if (std::abs(t_ieta) > ietasL_[j] && std::abs(t_ieta) <= ietasH_[j]) {
1132 if (jp1 >= 0)
1133 jp2 = j;
1134 if (jp1 < 0)
1135 jp1 = j;
1136 }
1137 }
1138 int jp3 = (jp1 >= 0) ? (int)(ietasL_.size() - 1) : jp1;
1139 double cut = (pmom > 20) ? ((flexibleSelect_ == 0) ? 2.0 : 10.0) : 0.0;
1140
1141 double rcut(-1000.0);
1142 if (debug)
1143 std::cout << "Bin " << kp << ":" << kp1 << ":" << kv << ":" << kv1 << ":" << kd << ":" << kd1 << ":" << jp << ":"
1144 << jp1 << ":" << jp2 << ":" << jp3 << " pmom:ieta:pv:mindR " << pmom << ":" << std::abs(t_ieta) << ":"
1145 << t_ieta << ":" << t_goodPV << ":" << t_mindR1 << " Cuts " << cut << ":" << rcut << std::endl;
1146
1147
1148 double rat(1.0), eHcal(t_eHcal);
1149 if ((corrFactor_->doCorr()) || (cFactor_ != nullptr) || ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr()))) {
1150 eHcal = 0;
1151 for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
1152
1153 bool okcell = (thrForm_ == 0) || ((*t_HitEnergies)[k] > threshold((*t_DetIds)[k], thrForm_));
1154
1155 if (okcell) {
1156 double cfac(1.0);
1157 if (corrFactor_->doCorr()) {
1158 unsigned int id = truncateId((*t_DetIds)[k], truncateFlag_, false);
1159 cfac = corrFactor_->getCorr(id);
1160 }
1161 if ((cFactor_ != nullptr) && (ifDepth_ != 3) && (ifDepth_ > 0))
1162 cfac *= cFactor_->getCorr(t_Run, (*t_DetIds)[k]);
1163 if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr(1)))
1164 cfac *= cDuplicate_->getWeight((*t_DetIds)[k]);
1165 if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr(3))) {
1166 int subdet, zside, ieta, iphi, depth;
1167 unpackDetId((*t_DetIds)[k], subdet, zside, ieta, iphi, depth);
1168 cfac *= cDuplicate_->getCorr(t_Run, ieta, depth);
1169 }
1170 eHcal += (cfac * ((*t_HitEnergies)[k]));
1171 if (debug) {
1172 int subdet, zside, ieta, iphi, depth;
1173 unpackDetId((*t_DetIds)[k], subdet, zside, ieta, iphi, depth);
1174 std::cout << zside << ":" << ieta << ":" << depth << " Corr " << cfac << " " << (*t_HitEnergies)[k]
1175 << " Out " << eHcal << std::endl;
1176 }
1177 }
1178 }
1179 }
1180 bool goodTk = goodTrack(eHcal, cut, jentry, debug);
1181 bool selPhi = selectPhi(debug);
1182 if (t_qltyFlag) {
1183 if (p4060)
1184 ++kount50[7];
1185 if (kp == 0) {
1186 ++kount0[7];
1187 } else if (kp == 1) {
1188 ++kount1[7];
1189 } else if (kp == 2) {
1190 ++kount2[7];
1191 } else if (kp == 3) {
1192 ++kount3[7];
1193 } else if (kp == 4) {
1194 ++kount4[7];
1195 } else if (kp == 5) {
1196 ++kount5[7];
1197 }
1198 if (t_selectTk) {
1199 if (p4060)
1200 ++kount50[8];
1201 if (kp == 0) {
1202 ++kount0[8];
1203 } else if (kp == 1) {
1204 ++kount1[8];
1205 } else if (kp == 2) {
1206 ++kount2[8];
1207 } else if (kp == 3) {
1208 ++kount3[8];
1209 } else if (kp == 4) {
1210 ++kount4[8];
1211 } else if (kp == 5) {
1212 ++kount5[8];
1213 }
1214 if (t_hmaxNearP < cut) {
1215 if (p4060)
1216 ++kount50[9];
1217 if (kp == 0) {
1218 ++kount0[9];
1219 } else if (kp == 1) {
1220 ++kount1[9];
1221 } else if (kp == 2) {
1222 ++kount2[9];
1223 } else if (kp == 3) {
1224 ++kount3[9];
1225 } else if (kp == 4) {
1226 ++kount4[9];
1227 } else if (kp == 5) {
1228 ++kount5[9];
1229 }
1230 if (t_eMipDR < 1.0) {
1231 if (p4060)
1232 ++kount50[10];
1233 if (kp == 0) {
1234 ++kount0[10];
1235 } else if (kp == 1) {
1236 ++kount1[10];
1237 } else if (kp == 2) {
1238 ++kount2[10];
1239 } else if (kp == 3) {
1240 ++kount3[10];
1241 } else if (kp == 4) {
1242 ++kount4[10];
1243 } else if (kp == 5) {
1244 ++kount5[10];
1245 }
1246 if (eHcal > 0.001) {
1247 if (p4060)
1248 ++kount50[11];
1249 if (kp == 0) {
1250 ++kount0[11];
1251 } else if (kp == 1) {
1252 ++kount1[11];
1253 } else if (kp == 2) {
1254 ++kount2[11];
1255 } else if (kp == 3) {
1256 ++kount3[11];
1257 } else if (kp == 4) {
1258 ++kount4[11];
1259 } else if (kp == 5) {
1260 ++kount5[11];
1261 }
1262 if (selPhi) {
1263 if (p4060)
1264 ++kount50[12];
1265 if (kp == 0) {
1266 ++kount0[12];
1267 } else if (kp == 1) {
1268 ++kount1[12];
1269 } else if (kp == 2) {
1270 ++kount2[12];
1271 } else if (kp == 3) {
1272 ++kount3[12];
1273 } else if (kp == 4) {
1274 ++kount4[12];
1275 } else if (kp == 5) {
1276 ++kount5[12];
1277 }
1278 }
1279 }
1280 }
1281 }
1282 }
1283 }
1284 if (pmom > 0)
1285 rat = (eHcal / (pmom - t_eMipDR));
1286 if (debug) {
1287 std::cout << "Entry " << jentry << " p|eHcal|ratio " << pmom << "|" << t_eHcal << "|" << eHcal << "|" << rat
1288 << "|" << kp << "|" << kv << "|" << jp << " Cuts " << t_qltyFlag << "|" << t_selectTk << "|"
1289 << (t_hmaxNearP < cut) << "|" << (t_eMipDR < 1.0) << "|" << goodTk << "|" << (rat > rcut)
1290 << " Select Phi " << selPhi << " hmaxNearP " << t_hmaxNearP << " eMipDR " << t_eMipDR << std::endl;
1291 std::cout << "D1 : " << kp << ":" << kp1 << ":" << kv << ":" << kv1 << ":" << kd << ":" << kd1 << ":" << jp
1292 << std::endl;
1293 }
1294 if (goodTk && (kp >= 0) && selPhi) {
1295 if (p4060)
1296 ++kount50[13];
1297 if (kp == 0) {
1298 ++kount0[13];
1299 } else if (kp == 1) {
1300 ++kount1[13];
1301 } else if (kp == 2) {
1302 ++kount2[13];
1303 } else if (kp == 3) {
1304 ++kount3[13];
1305 } else if (kp == 4) {
1306 ++kount4[13];
1307 } else if (kp == 5) {
1308 ++kount5[13];
1309 }
1310 if (t_eHcal < 0.01) {
1311 std::map<int, counter>::const_iterator itr = runEn1.find(t_Run);
1312 if (itr == runEn1.end()) {
1313 counter knt;
1314 if ((kp >= 0) && (kp < counter::npsize))
1315 knt.count[kp] = 1;
1316 knt.total = 1;
1317 runEn1[t_Run] = knt;
1318 } else {
1319 counter knt = runEn1[t_Run];
1320 if ((kp >= 0) && (kp < counter::npsize))
1321 ++knt.count[kp];
1322 ++knt.total;
1323 runEn1[t_Run] = knt;
1324 }
1325 }
1326 if (t_eMipDR < 0.01 && t_eHcal < 0.01) {
1327 if (p4060)
1328 ++kount50[14];
1329 if (kp == 0) {
1330 ++kount0[14];
1331 } else if (kp == 1) {
1332 ++kount1[14];
1333 } else if (kp == 2) {
1334 ++kount2[14];
1335 } else if (kp == 3) {
1336 ++kount3[14];
1337 } else if (kp == 4) {
1338 ++kount4[14];
1339 } else if (kp == 5) {
1340 ++kount5[14];
1341 }
1342 std::map<int, counter>::const_iterator itr = runEn2.find(t_Run);
1343 if (itr == runEn2.end()) {
1344 counter knt;
1345 if ((kp >= 0) && (kp < counter::npsize))
1346 knt.count[kp] = 1;
1347 knt.total = 1;
1348 runEn2[t_Run] = knt;
1349 } else {
1350 counter knt = runEn2[t_Run];
1351 if ((kp >= 0) && (kp < counter::npsize))
1352 ++knt.count[kp];
1353 ++knt.total;
1354 runEn2[t_Run] = knt;
1355 }
1356 }
1357 if (rat > rcut) {
1358 if (debug)
1359 std::cout << "kp " << kp << " " << h_p[kp - 1]->GetName() << " p " << pmom << " wt " << t_EventWeight
1360 << std::endl;
1361 if (kp > 0)
1362 h_p[kp - 1]->Fill(pmom, t_EventWeight);
1363 if (p4060)
1364 ++kount50[15];
1365 if (kp == 0) {
1366 ++kount0[15];
1367 } else if (kp == 1) {
1368 ++kount1[15];
1369 } else if (kp == 2) {
1370 ++kount2[15];
1371 } else if (kp == 3) {
1372 ++kount3[15];
1373 } else if (kp == 4) {
1374 ++kount4[15];
1375 } else if (kp == 5) {
1376 ++kount5[15];
1377 }
1378 if (plotType_ <= 1) {
1379 h_etaX[kp][kv]->Fill(eta, rat, t_EventWeight);
1380 h_etaX[kp][kv1]->Fill(eta, rat, t_EventWeight);
1381 h_nvxR[kp][kv]->Fill(rat, t_EventWeight);
1382 h_nvxR[kp][kv1]->Fill(rat, t_EventWeight);
1383 h_dL1R[kp][kd]->Fill(rat, t_EventWeight);
1384 h_dL1R[kp][kd1]->Fill(rat, t_EventWeight);
1385 if (jp > 0)
1386 h_etaR[kp][jp]->Fill(rat, t_EventWeight);
1387 h_etaR[kp][0]->Fill(rat, t_EventWeight);
1388 }
1389 h_pp[kp][0]->Fill(pmom, t_EventWeight);
1390 if (jp1 >= 0) {
1391 h_pp[kp][jp1 + 1]->Fill(pmom, t_EventWeight);
1392 h_pp[kp][jp3 + 1]->Fill(pmom, t_EventWeight);
1393 }
1394 if (jp2 >= 0)
1395 h_pp[kp][jp2 + 1]->Fill(pmom, t_EventWeight);
1396 if (kp == (int)(kp50)) {
1397 std::map<int, counter>::const_iterator itr = runSum.find(t_Run);
1398 if (itr == runSum.end()) {
1399 counter knt;
1400 if ((kp >= 0) && (kp < counter::npsize))
1401 knt.count[kp] = 1;
1402 knt.total = 1;
1403 runSum[t_Run] = knt;
1404 } else {
1405 counter knt = runSum[t_Run];
1406 if ((kp >= 0) && (kp < counter::npsize))
1407 ++knt.count[kp];
1408 ++knt.total;
1409 runSum[t_Run] = knt;
1410 }
1411 }
1412 if ((!isRealData_) || (t_mindR1 > 0.5) || (t_DataType == 1)) {
1413 if (p4060)
1414 ++kount50[16];
1415 if (kp == 0) {
1416 ++kount0[16];
1417 } else if (kp == 1) {
1418 ++kount1[16];
1419 } else if (kp == 2) {
1420 ++kount2[16];
1421 } else if (kp == 3) {
1422 ++kount3[16];
1423 } else if (kp == 4) {
1424 ++kount4[16];
1425 } else if (kp == 5) {
1426 ++kount5[16];
1427 }
1428 ++kounts[kp];
1429 if (plotType_ <= 1) {
1430 if (jp > 0)
1431 h_etaF[kp][jp]->Fill(rat, t_EventWeight);
1432 h_etaF[kp][0]->Fill(rat, t_EventWeight);
1433 } else {
1434 if (debug) {
1435 std::cout << "kp " << kp << h_etaF[kp].size() << std::endl;
1436 }
1437 if (jp > 0)
1438 h_etaF[kp][jp]->Fill(rat, t_EventWeight);
1439 h_etaF[kp][0]->Fill(rat, t_EventWeight);
1440 if (jp1 >= 0) {
1441 h_etaB[kp][jp1]->Fill(rat, t_EventWeight);
1442 h_etaB[kp][jp3]->Fill(rat, t_EventWeight);
1443 }
1444 if (jp2 >= 0)
1445 h_etaB[kp][jp2]->Fill(rat, t_EventWeight);
1446 }
1447 if (selRBX_ && (kp == (int)(kp50)) && ((t_ieta * zside_) > 0)) {
1448 int rbx = (t_iphi > 70) ? 0 : (t_iphi + 1) / 4;
1449 h_rbx[rbx]->Fill(rat, t_EventWeight);
1450 }
1451 }
1452 if (pmom > 10.0) {
1453 if (plotType_ <= 1) {
1454 h_etaX[kp1][kv]->Fill(eta, rat, t_EventWeight);
1455 h_etaX[kp1][kv1]->Fill(eta, rat, t_EventWeight);
1456 h_nvxR[kp1][kv]->Fill(rat, t_EventWeight);
1457 h_nvxR[kp1][kv1]->Fill(rat, t_EventWeight);
1458 h_dL1R[kp1][kd]->Fill(rat, t_EventWeight);
1459 h_dL1R[kp1][kd1]->Fill(rat, t_EventWeight);
1460 if (jp > 0)
1461 h_etaR[kp1][jp]->Fill(rat, t_EventWeight);
1462 h_etaR[kp1][0]->Fill(rat, t_EventWeight);
1463 if (jp1 >= 0) {
1464 h_etaB[kp][jp1]->Fill(rat, t_EventWeight);
1465 h_etaB[kp][jp3]->Fill(rat, t_EventWeight);
1466 }
1467 if (jp2 >= 0)
1468 h_etaB[kp][jp2]->Fill(rat, t_EventWeight);
1469 }
1470 if (p4060)
1471 ++kount50[17];
1472 if (kp == 0) {
1473 ++kount0[17];
1474 } else if (kp == 1) {
1475 ++kount1[17];
1476 } else if (kp == 2) {
1477 ++kount2[17];
1478 } else if (kp == 3) {
1479 ++kount3[17];
1480 } else if (kp == 4) {
1481 ++kount4[17];
1482 } else if (kp == 5) {
1483 ++kount5[17];
1484 }
1485 }
1486 }
1487 }
1488 if (pmom > 10.0) {
1489 kount++;
1490 if (((flag_ / 100) % 10) != 0) {
1491 good++;
1492 fileout_ << good << " " << jentry << " " << t_Run << " " << t_Event << " " << t_ieta << " " << pmom
1493 << std::endl;
1494 }
1495 }
1496 }
1497 unsigned int k(0);
1498 std::cout << "\nSummary of entries with " << runSum.size() << " runs\n";
1499 for (std::map<int, counter>::iterator itr = runSum.begin(); itr != runSum.end(); ++itr, ++k)
1500 std::cout << "[" << k << "] Run " << itr->first << " Total " << (itr->second).total << " in p-bins "
1501 << (itr->second).count[0] << ":" << (itr->second).count[1] << ":" << (itr->second).count[2] << ":"
1502 << (itr->second).count[3] << std::endl;
1503 k = 0;
1504 std::cout << runEn1.size() << " runs with 0 energy in HCAL\n";
1505 for (std::map<int, counter>::iterator itr = runEn1.begin(); itr != runEn1.end(); ++itr, ++k)
1506 std::cout << "[" << k << "] Run " << itr->first << " Total " << (itr->second).total << " in p-bins "
1507 << (itr->second).count[0] << ":" << (itr->second).count[1] << ":" << (itr->second).count[2] << ":"
1508 << (itr->second).count[3] << std::endl;
1509 k = 0;
1510 std::cout << runEn2.size() << " runs with 0 energy in ECAL and HCAL\n";
1511 for (std::map<int, counter>::iterator itr = runEn2.begin(); itr != runEn2.end(); ++itr, ++k)
1512 std::cout << "[" << k << "] Run " << itr->first << " Total " << (itr->second).total << " in p-bins "
1513 << (itr->second).count[0] << ":" << (itr->second).count[1] << ":" << (itr->second).count[2] << ":"
1514 << (itr->second).count[3] << std::endl;
1515 if (((flag_ / 100) % 10) > 0) {
1516 fileout_.close();
1517 std::cout << "Writes " << good << " events in the file " << outFileName_ << std::endl;
1518 }
1519 std::cout << "Finds " << duplicate << " Duplicate events out of " << kount << " events in this file with p>10 Gev"
1520 << std::endl;
1521 std::cout << "Number of selected events:" << std::endl;
1522 for (unsigned int k = 1; k < ps_.size(); ++k)
1523 std::cout << ps_[k - 1] << ":" << ps_[k] << " " << kounts[k - 1] << std::endl;
1524 std::cout << "Number in each step for tracks of momentum 40-60 GeV: ";
1525 for (unsigned int k = 0; k < 18; ++k)
1526 std::cout << " [" << k << "] " << kount50[k];
1527 std::cout << std::endl;
1528 for (unsigned int k = 1; k < ps_.size(); ++k) {
1529 std::cout << "Number in each step for tracks of momentum " << ps_[k - 1] << "-" << ps_[k] << " Gev: ";
1530 for (unsigned int k1 = 0; k1 < 18; ++k1) {
1531 if (k == 1) {
1532 std::cout << " [" << k1 << "] " << kount0[k1];
1533 } else if (k == 2) {
1534 std::cout << " [" << k1 << "] " << kount1[k1];
1535 } else if (k == 3) {
1536 std::cout << " [" << k1 << "] " << kount2[k1];
1537 } else if (k == 4) {
1538 std::cout << " [" << k1 << "] " << kount3[k1];
1539 } else if (k == 5) {
1540 std::cout << " [" << k1 << "] " << kount4[k1];
1541 } else if (k == 6) {
1542 std::cout << " [" << k1 << "] " << kount5[k1];
1543 }
1544 }
1545 std::cout << std::endl;
1546 }
1547 }
1548
1549 bool CalibMonitor::goodTrack(double &eHcal, double &cuti, const Long64_t &entry, bool debug) {
1550 bool select(true);
1551 double cut(cuti);
1552 if (debug) {
1553 std::cout << "goodTrack input " << eHcal << ":" << cut;
1554 }
1555 if (flexibleSelect_ > 1) {
1556 double eta = (t_ieta > 0) ? t_ieta : -t_ieta;
1557 cut = 8.0 * exp(eta * log2by18_);
1558 }
1559 correctEnergy(eHcal, entry);
1560 select = ((t_qltyFlag) && (t_selectTk) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (eHcal > 0.001));
1561 if (debug) {
1562 std::cout << " output " << select << " Based on " << t_qltyFlag << ":" << t_selectTk << ":" << t_hmaxNearP << ":"
1563 << cut << ":" << t_eMipDR << ":" << eHcal << std::endl;
1564 }
1565 return select;
1566 }
1567
1568 bool CalibMonitor::selectPhi(bool debug) {
1569 bool select(true);
1570 if (phimin_ > 1 || phimax_ < 72) {
1571 double eTotal(0), eSelec(0);
1572
1573 for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
1574
1575 bool okcell = (thrForm_ == 0) || ((*t_HitEnergies)[k] > threshold((*t_DetIds)[k], thrForm_));
1576 if (okcell) {
1577 int iphi = ((*t_DetIds)[k]) & (0x3FF);
1578 int zside = ((*t_DetIds)[k] & 0x80000) ? (1) : (-1);
1579 eTotal += ((*t_HitEnergies)[k]);
1580 if (iphi >= phimin_ && iphi <= phimax_ && zside == zside_)
1581 eSelec += ((*t_HitEnergies)[k]);
1582 }
1583 }
1584 if (eSelec < 0.9 * eTotal)
1585 select = false;
1586 if (debug) {
1587 std::cout << "Etotal " << eTotal << " and ESelec " << eSelec << " (phi " << phimin_ << ":" << phimax_ << " z "
1588 << zside_ << ") Selection " << select << std::endl;
1589 }
1590 }
1591 return select;
1592 }
1593
1594 void CalibMonitor::plotHist(int itype, int inum, bool save) {
1595 gStyle->SetCanvasBorderMode(0);
1596 gStyle->SetCanvasColor(kWhite);
1597 gStyle->SetPadColor(kWhite);
1598 gStyle->SetFillColor(kWhite);
1599 gStyle->SetOptStat(111110);
1600 gStyle->SetOptFit(1);
1601 char name[100];
1602 int itmin = (itype >= 0 && itype < 4) ? itype : 0;
1603 int itmax = (itype >= 0 && itype < 4) ? itype : 3;
1604 std::string types[4] = {
1605 "E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})", "E_{HCAL}/(p-E_{ECAL})"};
1606 int nmax[4] = {npbin, npbin, npbin, npbin};
1607 for (int type = itmin; type <= itmax; ++type) {
1608 int inmin = (inum >= 0 && inum < nmax[type]) ? inum : 0;
1609 int inmax = (inum >= 0 && inum < nmax[type]) ? inum : nmax[type] - 1;
1610 int kmax = 1;
1611 if (type == 0)
1612 kmax = (int)(etas_.size());
1613 else if (type == 1)
1614 kmax = (int)(etas_.size());
1615 else if (type == 2)
1616 kmax = (int)(nvx_.size());
1617 else
1618 kmax = (int)(dl1_.size());
1619 for (int num = inmin; num <= inmax; ++num) {
1620 for (int k = 0; k < kmax; ++k) {
1621 sprintf(name, "c_%d%d%d", type, num, k);
1622 TCanvas *pad = new TCanvas(name, name, 700, 500);
1623 pad->SetRightMargin(0.10);
1624 pad->SetTopMargin(0.10);
1625 sprintf(name, "%s", types[type].c_str());
1626 if (type != 7) {
1627 TH1D *hist(0);
1628 if (type == 0)
1629 hist = (TH1D *)(h_etaR[num][k]->Clone());
1630 else if (type == 1)
1631 hist = (TH1D *)(h_etaF[num][k]->Clone());
1632 else if (type == 2)
1633 hist = (TH1D *)(h_nvxR[num][k]->Clone());
1634 else
1635 hist = (TH1D *)(h_dL1R[num][k]->Clone());
1636 hist->GetXaxis()->SetTitle(name);
1637 hist->GetYaxis()->SetTitle("Tracks");
1638 drawHist(hist, pad);
1639 if (save) {
1640 sprintf(name, "c_%s%d%d%d.gif", prefix_.c_str(), type, num, k);
1641 pad->Print(name);
1642 }
1643 } else {
1644 TProfile *hist = (TProfile *)(h_etaX[num][k]->Clone());
1645 hist->GetXaxis()->SetTitle(name);
1646 hist->GetYaxis()->SetTitle("<E_{HCAL}/(p-E_{ECAL})>");
1647 hist->GetYaxis()->SetRangeUser(0.4, 1.6);
1648 hist->Fit("pol0", "q");
1649 drawHist(hist, pad);
1650 if (save) {
1651 sprintf(name, "c_%s%d%d%d.gif", prefix_.c_str(), type, num, k);
1652 pad->Print(name);
1653 }
1654 }
1655 }
1656 }
1657 }
1658 }
1659
1660 template <class Hist>
1661 void CalibMonitor::drawHist(Hist *hist, TCanvas *pad) {
1662 hist->GetYaxis()->SetLabelOffset(0.005);
1663 hist->GetYaxis()->SetTitleOffset(1.20);
1664 hist->Draw();
1665 pad->Update();
1666 TPaveStats *st1 = (TPaveStats *)hist->GetListOfFunctions()->FindObject("stats");
1667 if (st1 != NULL) {
1668 st1->SetY1NDC(0.70);
1669 st1->SetY2NDC(0.90);
1670 st1->SetX1NDC(0.55);
1671 st1->SetX2NDC(0.90);
1672 }
1673 pad->Modified();
1674 pad->Update();
1675 }
1676
1677 void CalibMonitor::savePlot(const std::string &theName, bool append, bool all) {
1678 TFile *theFile(0);
1679 if (append) {
1680 theFile = new TFile(theName.c_str(), "UPDATE");
1681 } else {
1682 theFile = new TFile(theName.c_str(), "RECREATE");
1683 }
1684
1685 theFile->cd();
1686 for (unsigned int k = 0; k < ps_.size(); ++k) {
1687 for (unsigned int j = 0; j <= ietasL_.size(); ++j) {
1688 if ((all || k == kp50) && h_pp[k].size() > j && h_pp[k][j] != 0) {
1689 TH1D *hist = (TH1D *)h_pp[k][j]->Clone();
1690 hist->Write();
1691 }
1692 }
1693 if (plotType_ <= 1) {
1694 for (unsigned int i = 0; i < nvx_.size(); ++i) {
1695 if (h_etaX[k][i] != 0) {
1696 TProfile *hnew = (TProfile *)h_etaX[k][i]->Clone();
1697 hnew->Write();
1698 }
1699 if (h_nvxR[k].size() > i && h_nvxR[k][i] != 0) {
1700 TH1D *hist = (TH1D *)h_nvxR[k][i]->Clone();
1701 hist->Write();
1702 }
1703 }
1704 }
1705 for (unsigned int j = 0; j < etas_.size(); ++j) {
1706 if ((plotType_ <= 1) && (h_etaR[k][j] != 0)) {
1707 TH1D *hist = (TH1D *)h_etaR[k][j]->Clone();
1708 hist->Write();
1709 }
1710 if (h_etaF[k].size() > j && h_etaF[k][j] != 0 && (all || (k == kp50))) {
1711 TH1D *hist = (TH1D *)h_etaF[k][j]->Clone();
1712 hist->Write();
1713 }
1714 }
1715 if (plotType_ <= 1) {
1716 for (unsigned int j = 0; j < dl1_.size(); ++j) {
1717 if (h_dL1R[k][j] != 0) {
1718 TH1D *hist = (TH1D *)h_dL1R[k][j]->Clone();
1719 hist->Write();
1720 }
1721 }
1722 }
1723 for (unsigned int j = 0; j < ietasL_.size(); ++j) {
1724 if (h_etaB[k].size() > j && h_etaB[k][j] != 0 && (all || (k == kp50))) {
1725 TH1D *hist = (TH1D *)h_etaB[k][j]->Clone();
1726 hist->Write();
1727 }
1728 }
1729 }
1730 if (selRBX_) {
1731 for (unsigned int k = 0; k < 18; ++k) {
1732 if (h_rbx[k] != 0) {
1733 TH1D *h1 = (TH1D *)h_rbx[k]->Clone();
1734 h1->Write();
1735 }
1736 }
1737 }
1738 for (unsigned int k = 0; k < h_p.size(); ++k) {
1739 if (h_p[k] != 0) {
1740 TH1D *h1 = (TH1D *)h_p[k]->Clone();
1741 h1->Write();
1742 }
1743 }
1744 std::cout << "All done" << std::endl;
1745 theFile->Close();
1746 }
1747
1748 void CalibMonitor::correctEnergy(double &eHcal, const Long64_t &entry) {
1749 bool debug(false);
1750 double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
1751 if ((ifDepth_ == 3) && (cFactor_ != nullptr)) {
1752 double cfac = cFactor_->getCorr(entry);
1753 eHcal *= cfac;
1754 if (debug)
1755 std::cout << "PU Factor for " << ifDepth_ << ":" << entry << " = " << cfac << ":" << eHcal << std::endl;
1756 } else if ((corrPU_ < 0) && (pmom > 0)) {
1757 double ediff = (t_eHcal30 - t_eHcal10);
1758 if (t_DetIds1 != 0 && t_DetIds3 != 0) {
1759 double Etot1(0), Etot3(0);
1760
1761 for (unsigned int idet = 0; idet < (*t_DetIds1).size(); idet++) {
1762
1763 bool okcell = (thrForm_ == 0) || ((*t_HitEnergies1)[idet] > threshold((*t_DetIds1)[idet], thrForm_));
1764 if (okcell) {
1765 unsigned int id = truncateId((*t_DetIds1)[idet], truncateFlag_, false);
1766 double cfac = corrFactor_->getCorr(id);
1767 if ((cFactor_ != 0) && (ifDepth_ != 3) && (ifDepth_ > 0))
1768 cfac *= cFactor_->getCorr(t_Run, (*t_DetIds1)[idet]);
1769 if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr(1)))
1770 cfac *= cDuplicate_->getWeight((*t_DetIds1)[idet]);
1771 if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr(3))) {
1772 int subdet, zside, ieta, iphi, depth;
1773 unpackDetId((*t_DetIds1)[idet], subdet, zside, ieta, iphi, depth);
1774 cfac *= cDuplicate_->getCorr(t_Run, ieta, depth);
1775 }
1776 double hitEn = cfac * (*t_HitEnergies1)[idet];
1777 Etot1 += hitEn;
1778 }
1779 }
1780 for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
1781
1782 bool okcell = (thrForm_ == 0) || ((*t_HitEnergies3)[idet] > threshold((*t_DetIds3)[idet], thrForm_));
1783 if (okcell) {
1784 unsigned int id = truncateId((*t_DetIds3)[idet], truncateFlag_, false);
1785 double cfac = corrFactor_->getCorr(id);
1786 if ((cFactor_ != 0) && (ifDepth_ != 3) && (ifDepth_ > 0))
1787 cfac *= cFactor_->getCorr(t_Run, (*t_DetIds3)[idet]);
1788 if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr(1)))
1789 cfac *= cDuplicate_->getWeight((*t_DetIds3)[idet]);
1790 if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr(3))) {
1791 int subdet, zside, ieta, iphi, depth;
1792 unpackDetId((*t_DetIds3)[idet], subdet, zside, ieta, iphi, depth);
1793 cfac *= cDuplicate_->getCorr(t_Run, ieta, depth);
1794 }
1795 double hitEn = cfac * (*t_HitEnergies3)[idet];
1796 Etot3 += hitEn;
1797 }
1798 }
1799 ediff = (Etot3 - Etot1);
1800 }
1801 double fac = puFactor(-corrPU_, t_ieta, pmom, eHcal, ediff, false);
1802 if (debug) {
1803 double fac1 = puFactor(-corrPU_, t_ieta, pmom, eHcal, ediff, true);
1804 double fac2 = puFactor(2, t_ieta, pmom, eHcal, ediff, true);
1805 std::cout << "PU Factor for " << -corrPU_ << " = " << fac1 << "; for 2 = " << fac2 << std::endl;
1806 }
1807 eHcal *= fac;
1808 } else if (corrPU_ > 0) {
1809 eHcal = puFactorRho(corrPU_, t_ieta, t_rhoh, eHcal);
1810 }
1811 }
1812
1813 class GetEntries {
1814 public:
1815 TTree *fChain;
1816 Int_t fCurrent;
1817
1818
1819 UInt_t t_RunNo;
1820 UInt_t t_EventNo;
1821 Int_t t_Tracks;
1822 Int_t t_TracksProp;
1823 Int_t t_TracksSaved;
1824 Int_t t_TracksLoose;
1825 Int_t t_TracksTight;
1826 Int_t t_allvertex;
1827 Bool_t t_TrigPass;
1828 Bool_t t_TrigPassSel;
1829 Bool_t t_L1Bit;
1830 std::vector<Bool_t> *t_hltbits;
1831 std::vector<int> *t_ietaAll;
1832 std::vector<int> *t_ietaGood;
1833 std::vector<int> *t_trackType;
1834
1835
1836 TBranch *b_t_RunNo;
1837 TBranch *b_t_EventNo;
1838 TBranch *b_t_Tracks;
1839 TBranch *b_t_TracksProp;
1840 TBranch *b_t_TracksSaved;
1841 TBranch *b_t_TracksLoose;
1842 TBranch *b_t_TracksTight;
1843 TBranch *b_t_allvertex;
1844 TBranch *b_t_TrigPass;
1845 TBranch *b_t_TrigPassSel;
1846 TBranch *b_t_L1Bit;
1847 TBranch *b_t_hltbits;
1848 TBranch *b_t_ietaAll;
1849 TBranch *b_t_ietaGood;
1850 TBranch *b_t_trackType;
1851
1852 GetEntries(const char *fname,
1853 const char *dirname,
1854 const char *dupFileName,
1855 const unsigned int bit1,
1856 const unsigned int bit2);
1857 virtual ~GetEntries();
1858 virtual Int_t Cut(Long64_t entry);
1859 virtual Int_t GetEntry(Long64_t entry);
1860 virtual Long64_t LoadTree(Long64_t entry);
1861 virtual void Init(TTree *tree, const char *dupFileName);
1862 virtual void Loop(Long64_t nmax = -1);
1863 virtual Bool_t Notify();
1864 virtual void Show(Long64_t entry = -1);
1865
1866 private:
1867 unsigned int bit_[2];
1868 std::vector<Long64_t> entries_;
1869 TH1I *h_tk[3], *h_eta[4], *h_pvx[3];
1870 TH1D *h_eff[3];
1871 };
1872
1873 GetEntries::GetEntries(
1874 const char *fname, const char *dirnm, const char *dupFileName, const unsigned int bit1, const unsigned int bit2) {
1875 TFile *file = new TFile(fname);
1876 TDirectory *dir = (TDirectory *)file->FindObjectAny(dirnm);
1877 std::cout << fname << " file " << file << " " << dirnm << " " << dir << std::endl;
1878 TTree *tree = (TTree *)dir->Get("EventInfo");
1879 std::cout << "CalibTree " << tree << std::endl;
1880 bit_[0] = bit1;
1881 bit_[1] = bit2;
1882 Init(tree, dupFileName);
1883 }
1884
1885 GetEntries::~GetEntries() {
1886 if (!fChain)
1887 return;
1888 delete fChain->GetCurrentFile();
1889 }
1890
1891 Int_t GetEntries::GetEntry(Long64_t entry) {
1892
1893 if (!fChain)
1894 return 0;
1895 return fChain->GetEntry(entry);
1896 }
1897
1898 Long64_t GetEntries::LoadTree(Long64_t entry) {
1899
1900 if (!fChain)
1901 return -5;
1902 Long64_t centry = fChain->LoadTree(entry);
1903 if (centry < 0)
1904 return centry;
1905 if (!fChain->InheritsFrom(TChain::Class()))
1906 return centry;
1907 TChain *chain = (TChain *)fChain;
1908 if (chain->GetTreeNumber() != fCurrent) {
1909 fCurrent = chain->GetTreeNumber();
1910 Notify();
1911 }
1912 return centry;
1913 }
1914
1915 void GetEntries::Init(TTree *tree, const char *dupFileName) {
1916
1917
1918
1919
1920
1921
1922
1923
1924
1925
1926 t_hltbits = 0;
1927 t_ietaAll = 0;
1928 t_ietaGood = 0;
1929 t_trackType = 0;
1930 t_L1Bit = false;
1931 fChain = tree;
1932 fCurrent = -1;
1933 if (!tree)
1934 return;
1935 fChain->SetMakeClass(1);
1936 fChain->SetBranchAddress("t_RunNo", &t_RunNo, &b_t_RunNo);
1937 fChain->SetBranchAddress("t_EventNo", &t_EventNo, &b_t_EventNo);
1938 fChain->SetBranchAddress("t_Tracks", &t_Tracks, &b_t_Tracks);
1939 fChain->SetBranchAddress("t_TracksProp", &t_TracksProp, &b_t_TracksProp);
1940 fChain->SetBranchAddress("t_TracksSaved", &t_TracksSaved, &b_t_TracksSaved);
1941 fChain->SetBranchAddress("t_TracksLoose", &t_TracksLoose, &b_t_TracksLoose);
1942 fChain->SetBranchAddress("t_TracksTight", &t_TracksTight, &b_t_TracksTight);
1943 fChain->SetBranchAddress("t_allvertex", &t_allvertex, &b_t_allvertex);
1944 fChain->SetBranchAddress("t_TrigPass", &t_TrigPass, &b_t_TrigPass);
1945 fChain->SetBranchAddress("t_TrigPassSel", &t_TrigPassSel, &b_t_TrigPassSel);
1946 fChain->SetBranchAddress("t_L1Bit", &t_L1Bit, &b_t_L1Bit);
1947 fChain->SetBranchAddress("t_hltbits", &t_hltbits, &b_t_hltbits);
1948 fChain->SetBranchAddress("t_ietaAll", &t_ietaAll, &b_t_ietaAll);
1949 fChain->SetBranchAddress("t_ietaGood", &t_ietaGood, &b_t_ietaGood);
1950 fChain->SetBranchAddress("t_trackType", &t_trackType, &b_t_trackType);
1951 Notify();
1952
1953 if (std::string(dupFileName) != "") {
1954 std::ifstream infile(dupFileName);
1955 if (!infile.is_open()) {
1956 std::cout << "Cannot open " << dupFileName << std::endl;
1957 } else {
1958 while (1) {
1959 Long64_t jentry;
1960 infile >> jentry;
1961 if (!infile.good())
1962 break;
1963 entries_.push_back(jentry);
1964 }
1965 infile.close();
1966 std::cout << "Reads a list of " << entries_.size() << " events from " << dupFileName << std::endl;
1967 }
1968 }
1969
1970 h_tk[0] = new TH1I("Track0", "# of tracks produced", 2000, 0, 2000);
1971 h_tk[1] = new TH1I("Track1", "# of tracks propagated", 2000, 0, 2000);
1972 h_tk[2] = new TH1I("Track2", "# of tracks saved in tree", 2000, 0, 2000);
1973 h_eta[0] = new TH1I("Eta0", "i#eta (All Tracks)", 60, -30, 30);
1974 h_eta[1] = new TH1I("Eta1", "i#eta (Good Tracks)", 60, -30, 30);
1975 h_eta[2] = new TH1I("Eta2", "i#eta (Loose Isolated Tracks)", 60, -30, 30);
1976 h_eta[3] = new TH1I("Eta3", "i#eta (Tight Isolated Tracks)", 60, -30, 30);
1977 h_eff[0] = new TH1D("Eff0", "i#eta (Selection Efficiency)", 60, -30, 30);
1978 h_eff[1] = new TH1D("Eff1", "i#eta (Loose Isolation Efficiency)", 60, -30, 30);
1979 h_eff[2] = new TH1D("Eff2", "i#eta (Tight Isolation Efficiency)", 60, -30, 30);
1980 h_pvx[0] = new TH1I("Pvx0", "Number of Good Vertex (All)", 100, 0, 100);
1981 h_pvx[1] = new TH1I("Pvx1", "Number of Good Vertex (Loose)", 100, 0, 100);
1982 h_pvx[2] = new TH1I("Pvx2", "Number of Good Vertex (Tight)", 100, 0, 100);
1983 }
1984
1985 Bool_t GetEntries::Notify() {
1986
1987
1988
1989
1990
1991
1992 return kTRUE;
1993 }
1994
1995 void GetEntries::Show(Long64_t entry) {
1996
1997
1998 if (!fChain)
1999 return;
2000 fChain->Show(entry);
2001 }
2002
2003 Int_t GetEntries::Cut(Long64_t) {
2004
2005
2006
2007 return 1;
2008 }
2009
2010 void GetEntries::Loop(Long64_t nmax) {
2011
2012
2013
2014
2015
2016
2017
2018
2019
2020
2021
2022
2023
2024
2025
2026
2027
2028
2029
2030
2031
2032
2033
2034 if (fChain == 0)
2035 return;
2036
2037 Long64_t nentries = fChain->GetEntriesFast();
2038 Long64_t nbytes = 0, nb = 0;
2039 unsigned int kount(0), duplicate(0), selected(0);
2040 int l1(0), hlt(0), loose(0), tight(0);
2041 int allHLT[3] = {0, 0, 0};
2042 int looseHLT[3] = {0, 0, 0};
2043 int tightHLT[3] = {0, 0, 0};
2044 Long64_t entries = (nmax > 0) ? nmax : nentries;
2045 int32_t oddEven = (nmax == -2) ? 1 : ((nmax == -3) ? -1 : 0);
2046 for (Long64_t jentry = 0; jentry < entries; jentry++) {
2047 Long64_t ientry = LoadTree(jentry);
2048 if (ientry < 0)
2049 break;
2050 nb = fChain->GetEntry(jentry);
2051 nbytes += nb;
2052 if (oddEven != 0) {
2053 if ((oddEven < 0) && (jentry % 2 == 0))
2054 continue;
2055 else if ((oddEven > 0) && (jentry % 2 != 0))
2056 continue;
2057 }
2058 bool select = (std::find(entries_.begin(), entries_.end(), jentry) == entries_.end());
2059 if (!select) {
2060 ++duplicate;
2061 continue;
2062 }
2063 h_tk[0]->Fill(t_Tracks);
2064 h_tk[1]->Fill(t_TracksProp);
2065 h_tk[2]->Fill(t_TracksSaved);
2066 h_pvx[0]->Fill(t_allvertex);
2067 if (t_TracksLoose > 0)
2068 h_pvx[1]->Fill(t_allvertex);
2069 if (t_TracksTight > 0)
2070 h_pvx[2]->Fill(t_allvertex);
2071 if (t_L1Bit) {
2072 ++l1;
2073 if (t_TracksLoose > 0)
2074 ++loose;
2075 if (t_TracksTight > 0)
2076 ++tight;
2077 if (t_TrigPass)
2078 ++hlt;
2079 }
2080 if (t_TrigPass) {
2081 ++kount;
2082 if (t_TrigPassSel)
2083 ++selected;
2084 }
2085 bool passt[2] = {false, false}, pass(false);
2086 for (unsigned k = 0; k < t_hltbits->size(); ++k) {
2087 if ((*t_hltbits)[k] > 0) {
2088 pass = true;
2089 for (int i = 0; i < 2; ++i)
2090 if (k == bit_[i])
2091 passt[i] = true;
2092 }
2093 }
2094 if (pass) {
2095 ++allHLT[0];
2096 for (int i = 0; i < 2; ++i)
2097 if (passt[i])
2098 ++allHLT[i + 1];
2099 if (t_TracksLoose > 0) {
2100 ++looseHLT[0];
2101 for (int i = 0; i < 2; ++i)
2102 if (passt[i])
2103 ++looseHLT[i + 1];
2104 }
2105 if (t_TracksTight > 0) {
2106 ++tightHLT[0];
2107 for (int i = 0; i < 2; ++i)
2108 if (passt[i])
2109 ++tightHLT[i + 1];
2110 }
2111 }
2112 for (unsigned int k = 0; k < t_ietaAll->size(); ++k)
2113 h_eta[0]->Fill((*t_ietaAll)[k]);
2114 for (unsigned int k = 0; k < t_ietaGood->size(); ++k) {
2115 h_eta[1]->Fill((*t_ietaGood)[k]);
2116 if (t_trackType->size() > 0) {
2117 if ((*t_trackType)[k] > 0)
2118 h_eta[2]->Fill((*t_ietaGood)[k]);
2119 if ((*t_trackType)[k] > 1)
2120 h_eta[3]->Fill((*t_ietaGood)[k]);
2121 }
2122 }
2123 }
2124 double ymaxk(0);
2125 for (int i = 1; i <= h_eff[0]->GetNbinsX(); ++i) {
2126 double rat(0), drat(0);
2127 if (h_eta[0]->GetBinContent(i) > ymaxk)
2128 ymaxk = h_eta[0]->GetBinContent(i);
2129 if ((h_eta[1]->GetBinContent(i) > 0) && (h_eta[0]->GetBinContent(i) > 0)) {
2130 rat = h_eta[1]->GetBinContent(i) / h_eta[0]->GetBinContent(i);
2131 drat = rat * std::sqrt(pow((h_eta[1]->GetBinError(i) / h_eta[1]->GetBinContent(i)), 2) +
2132 pow((h_eta[0]->GetBinError(i) / h_eta[0]->GetBinContent(i)), 2));
2133 }
2134 h_eff[0]->SetBinContent(i, rat);
2135 h_eff[0]->SetBinError(i, drat);
2136 if ((h_eta[1]->GetBinContent(i) > 0) && (h_eta[2]->GetBinContent(i) > 0)) {
2137 rat = h_eta[2]->GetBinContent(i) / h_eta[1]->GetBinContent(i);
2138 drat = rat * std::sqrt(pow((h_eta[2]->GetBinError(i) / h_eta[2]->GetBinContent(i)), 2) +
2139 pow((h_eta[1]->GetBinError(i) / h_eta[1]->GetBinContent(i)), 2));
2140 } else {
2141 rat = drat = 0;
2142 }
2143 h_eff[1]->SetBinContent(i, rat);
2144 h_eff[1]->SetBinError(i, drat);
2145 if ((h_eta[1]->GetBinContent(i) > 0) && (h_eta[3]->GetBinContent(i) > 0)) {
2146 rat = h_eta[3]->GetBinContent(i) / h_eta[1]->GetBinContent(i);
2147 drat = rat * std::sqrt(pow((h_eta[3]->GetBinError(i) / h_eta[3]->GetBinContent(i)), 2) +
2148 pow((h_eta[1]->GetBinError(i) / h_eta[1]->GetBinContent(i)), 2));
2149 } else {
2150 rat = drat = 0;
2151 }
2152 h_eff[1]->SetBinContent(i, rat);
2153 h_eff[1]->SetBinError(i, drat);
2154 }
2155 std::cout << "===== Remove " << duplicate << " events from " << nentries << "\n===== " << kount
2156 << " events passed trigger of which " << selected << " events get selected =====\n"
2157 << std::endl;
2158 std::cout << "===== " << l1 << " events passed L1 " << hlt << " events passed HLT and " << loose << ":" << tight
2159 << " events have at least 1 track candidate with loose:tight"
2160 << " isolation cut =====\n"
2161 << std::endl;
2162 for (int i = 0; i < 3; ++i) {
2163 char tbit[20];
2164 if (i == 0)
2165 sprintf(tbit, "Any");
2166 else
2167 sprintf(tbit, "%3d", bit_[i - 1]);
2168 std::cout << "Satisfying HLT trigger bit " << tbit << " Kount " << allHLT[i] << " Loose " << looseHLT[i]
2169 << " Tight " << tightHLT[i] << std::endl;
2170 }
2171
2172 gStyle->SetCanvasBorderMode(0);
2173 gStyle->SetCanvasColor(kWhite);
2174 gStyle->SetPadColor(kWhite);
2175 gStyle->SetFillColor(kWhite);
2176 gStyle->SetOptStat(1110);
2177 gStyle->SetOptTitle(0);
2178 int color[5] = {kBlack, kRed, kBlue, kMagenta, kCyan};
2179 int lines[5] = {1, 2, 3, 4, 5};
2180 TCanvas *pad1 = new TCanvas("c_track", "c_track", 500, 500);
2181 pad1->SetRightMargin(0.10);
2182 pad1->SetTopMargin(0.10);
2183 pad1->SetFillColor(kWhite);
2184 std::string titl1[3] = {"Reconstructed", "Propagated", "Saved"};
2185 TLegend *legend1 = new TLegend(0.11, 0.80, 0.50, 0.89);
2186 legend1->SetFillColor(kWhite);
2187 double ymax(0), xmax(0);
2188 for (int k = 0; k < 3; ++k) {
2189 int total(0), totaltk(0);
2190 for (int i = 1; i <= h_tk[k]->GetNbinsX(); ++i) {
2191 if (ymax < h_tk[k]->GetBinContent(i))
2192 ymax = h_tk[k]->GetBinContent(i);
2193 if (i > 1)
2194 total += (int)(h_tk[k]->GetBinContent(i));
2195 totaltk += (int)(h_tk[k]->GetBinContent(i)) * (i - 1);
2196 if (h_tk[k]->GetBinContent(i) > 0) {
2197 if (xmax < h_tk[k]->GetBinLowEdge(i) + h_tk[k]->GetBinWidth(i))
2198 xmax = h_tk[k]->GetBinLowEdge(i) + h_tk[k]->GetBinWidth(i);
2199 }
2200 }
2201 h_tk[k]->SetLineColor(color[k]);
2202 h_tk[k]->SetMarkerColor(color[k]);
2203 h_tk[k]->SetLineStyle(lines[k]);
2204 std::cout << h_tk[k]->GetTitle() << " Entries " << h_tk[k]->GetEntries() << " Events " << total << " Tracks "
2205 << totaltk << std::endl;
2206 legend1->AddEntry(h_tk[k], titl1[k].c_str(), "l");
2207 }
2208 int i1 = (int)(0.1 * xmax) + 1;
2209 xmax = 10.0 * i1;
2210 int i2 = (int)(0.01 * ymax) + 1;
2211
2212 ymax = 100.0 * i2;
2213 for (int k = 0; k < 3; ++k) {
2214 h_tk[k]->GetXaxis()->SetRangeUser(0, xmax);
2215 h_tk[k]->GetYaxis()->SetRangeUser(0.1, ymax);
2216 h_tk[k]->GetXaxis()->SetTitle("# Tracks");
2217 h_tk[k]->GetYaxis()->SetTitle("Events");
2218 h_tk[k]->GetYaxis()->SetLabelOffset(0.005);
2219 h_tk[k]->GetYaxis()->SetTitleOffset(1.20);
2220 if (k == 0)
2221 h_tk[k]->Draw("hist");
2222 else
2223 h_tk[k]->Draw("hist sames");
2224 }
2225 pad1->Update();
2226 pad1->SetLogy();
2227 ymax = 0.90;
2228 for (int k = 0; k < 3; ++k) {
2229 TPaveStats *st1 = (TPaveStats *)h_tk[k]->GetListOfFunctions()->FindObject("stats");
2230 if (st1 != NULL) {
2231 st1->SetLineColor(color[k]);
2232 st1->SetTextColor(color[k]);
2233 st1->SetY1NDC(ymax - 0.09);
2234 st1->SetY2NDC(ymax);
2235 st1->SetX1NDC(0.55);
2236 st1->SetX2NDC(0.90);
2237 ymax -= 0.09;
2238 }
2239 }
2240 pad1->Modified();
2241 pad1->Update();
2242 legend1->Draw("same");
2243 pad1->Update();
2244
2245 TCanvas *pad2 = new TCanvas("c_ieta", "c_ieta", 500, 500);
2246 pad2->SetRightMargin(0.10);
2247 pad2->SetTopMargin(0.10);
2248 pad2->SetFillColor(kWhite);
2249 pad2->SetLogy();
2250 std::string titl2[4] = {"All Tracks", "Selected Tracks", "Loose Isolation", "Tight Isolation"};
2251 TLegend *legend2 = new TLegend(0.11, 0.75, 0.50, 0.89);
2252 legend2->SetFillColor(kWhite);
2253 i2 = (int)(0.001 * ymaxk) + 1;
2254 ymax = 1000.0 * i2;
2255 for (int k = 0; k < 4; ++k) {
2256 h_eta[k]->GetYaxis()->SetRangeUser(1, ymax);
2257 h_eta[k]->SetLineColor(color[k]);
2258 h_eta[k]->SetMarkerColor(color[k]);
2259 h_eta[k]->SetLineStyle(lines[k]);
2260 h_eta[k]->GetXaxis()->SetTitle("i#eta");
2261 h_eta[k]->GetYaxis()->SetTitle("Tracks");
2262 h_eta[k]->GetYaxis()->SetLabelOffset(0.005);
2263 h_eta[k]->GetYaxis()->SetTitleOffset(1.20);
2264 legend2->AddEntry(h_eta[k], titl2[k].c_str(), "l");
2265 if (k == 0)
2266 h_eta[k]->Draw("hist");
2267 else
2268 h_eta[k]->Draw("hist sames");
2269 }
2270 pad2->Update();
2271 ymax = 0.90;
2272
2273 for (int k = 0; k < 4; ++k) {
2274 TPaveStats *st1 = (TPaveStats *)h_eta[k]->GetListOfFunctions()->FindObject("stats");
2275 if (st1 != NULL) {
2276 st1->SetLineColor(color[k]);
2277 st1->SetTextColor(color[k]);
2278 st1->SetY1NDC(ymax - 0.09);
2279 st1->SetY2NDC(ymax);
2280 st1->SetX1NDC(0.55);
2281 st1->SetX2NDC(0.90);
2282 ymax -= 0.09;
2283 }
2284 }
2285 pad2->Modified();
2286 pad2->Update();
2287 legend2->Draw("same");
2288 pad2->Update();
2289
2290 std::string titl3[3] = {"Selection", "Loose Isolation", "Tight Isolation"};
2291 TCanvas *pad3 = new TCanvas("c_effi", "c_effi", 500, 500);
2292 TLegend *legend3 = new TLegend(0.11, 0.785, 0.50, 0.89);
2293 pad3->SetRightMargin(0.10);
2294 pad3->SetTopMargin(0.10);
2295 pad3->SetFillColor(kWhite);
2296 pad3->SetLogy();
2297 for (int k = 0; k < 3; ++k) {
2298 h_eff[k]->SetStats(0);
2299 h_eff[k]->SetMarkerStyle(20);
2300 h_eff[k]->SetLineColor(color[k]);
2301 h_eff[k]->SetMarkerColor(color[k]);
2302 h_eff[k]->GetXaxis()->SetTitle("i#eta");
2303 h_eff[k]->GetYaxis()->SetTitle("Efficiency");
2304 h_eff[k]->GetYaxis()->SetLabelOffset(0.005);
2305 h_eff[k]->GetYaxis()->SetTitleOffset(1.20);
2306 if (k == 0)
2307 h_eff[k]->Draw("");
2308 else
2309 h_eff[k]->Draw("same");
2310 legend3->AddEntry(h_eff[k], titl3[k].c_str(), "l");
2311 }
2312 pad3->Modified();
2313 pad3->Update();
2314 legend3->Draw("same");
2315 pad3->Update();
2316
2317 std::string titl4[3] = {"All events", "Events with Loose Isolation", "Events with Tight Isolation"};
2318 TLegend *legend4 = new TLegend(0.11, 0.785, 0.50, 0.89);
2319 TCanvas *pad4 = new TCanvas("c_nvx", "c_nvx", 500, 500);
2320 pad4->SetRightMargin(0.10);
2321 pad4->SetTopMargin(0.10);
2322 pad4->SetFillColor(kWhite);
2323 pad4->SetLogy();
2324 for (int k = 0; k < 3; ++k) {
2325 h_pvx[k]->SetStats(1110);
2326 h_pvx[k]->SetMarkerStyle(20);
2327 h_pvx[k]->SetLineColor(color[k]);
2328 h_pvx[k]->SetMarkerColor(color[k]);
2329 h_pvx[k]->GetXaxis()->SetTitle("N_{PV}");
2330 h_pvx[k]->GetYaxis()->SetTitle("Events");
2331 h_pvx[k]->GetYaxis()->SetLabelOffset(0.005);
2332 h_pvx[k]->GetYaxis()->SetTitleOffset(1.20);
2333 if (k == 0)
2334 h_pvx[k]->Draw("");
2335 else
2336 h_pvx[k]->Draw("sames");
2337 legend4->AddEntry(h_pvx[k], titl4[k].c_str(), "l");
2338 }
2339 pad4->Update();
2340 ymax = 0.90;
2341 for (int k = 0; k < 3; ++k) {
2342 TPaveStats *st1 = (TPaveStats *)h_pvx[k]->GetListOfFunctions()->FindObject("stats");
2343 if (st1 != NULL) {
2344 st1->SetLineColor(color[k]);
2345 st1->SetTextColor(color[k]);
2346 st1->SetY1NDC(ymax - 0.09);
2347 st1->SetY2NDC(ymax);
2348 st1->SetX1NDC(0.55);
2349 st1->SetX2NDC(0.90);
2350 ymax -= 0.09;
2351 }
2352 }
2353 pad4->Modified();
2354 pad4->Update();
2355 legend4->Draw("same");
2356 pad4->Update();
2357 }