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