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