File indexing completed on 2022-09-08 03:21:21
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 #include <TCanvas.h>
0115 #include <TChain.h>
0116 #include <TF1.h>
0117 #include <TFile.h>
0118 #include <TFitResult.h>
0119 #include <TFitResultPtr.h>
0120 #include <TGraph.h>
0121 #include <TH1D.h>
0122 #include <TH2D.h>
0123 #include <TLegend.h>
0124 #include <TPaveStats.h>
0125 #include <TPaveText.h>
0126 #include <TProfile.h>
0127 #include <TROOT.h>
0128 #include <TStyle.h>
0129
0130 #include <algorithm>
0131 #include <fstream>
0132 #include <iomanip>
0133 #include <iostream>
0134 #include <map>
0135 #include <sstream>
0136 #include <string>
0137 #include <vector>
0138
0139 #include "CalibCorr.C"
0140
0141 struct record {
0142 record(int ser = 0, int ent = 0, int r = 0, int ev = 0, int ie = 0, double p = 0)
0143 : serial_(ser), entry_(ent), run_(r), event_(ev), ieta_(ie), p_(p) {}
0144
0145 int serial_, entry_, run_, event_, ieta_;
0146 double p_;
0147 };
0148
0149 struct recordLess {
0150 bool operator()(const record &a, const record &b) {
0151 return ((a.run_ < b.run_) || ((a.run_ == b.run_) && (a.event_ < b.event_)) ||
0152 ((a.run_ == b.run_) && (a.event_ == b.event_) && (a.ieta_ < b.ieta_)));
0153 }
0154 };
0155
0156 struct recordEvent {
0157 recordEvent(unsigned int ser = 0, unsigned int ent = 0, unsigned int r = 0, unsigned int ev = 0)
0158 : serial_(ser), entry_(ent), run_(r), event_(ev) {}
0159
0160 unsigned int serial_, entry_, run_, event_;
0161 };
0162
0163 struct recordEventLess {
0164 bool operator()(const recordEvent &a, const recordEvent &b) {
0165 return ((a.run_ < b.run_) || ((a.run_ == b.run_) && (a.event_ < b.event_)));
0166 }
0167 };
0168
0169 class CalibSort {
0170 public:
0171 CalibSort(const char *fname,
0172 std::string dirname = "HcalIsoTrkAnalyzer",
0173 std::string prefix = "",
0174 bool allEvent = false,
0175 int flag = 0,
0176 double mipCut = 2.0);
0177 virtual ~CalibSort();
0178 virtual Int_t Cut(Long64_t entry);
0179 virtual Int_t GetEntry(Long64_t entry);
0180 virtual Long64_t LoadTree(Long64_t entry);
0181 virtual void Init(TChain *);
0182 virtual void Loop(const char *);
0183 virtual Bool_t Notify();
0184 virtual void Show(Long64_t entry = -1);
0185
0186 private:
0187 TChain *fChain;
0188 Int_t fCurrent;
0189
0190
0191 Int_t t_Run;
0192 Int_t t_Event;
0193 Int_t t_DataType;
0194 Int_t t_ieta;
0195 Int_t t_iphi;
0196 Double_t t_EventWeight;
0197 Int_t t_nVtx;
0198 Int_t t_nTrk;
0199 Int_t t_goodPV;
0200 Double_t t_l1pt;
0201 Double_t t_l1eta;
0202 Double_t t_l1phi;
0203 Double_t t_l3pt;
0204 Double_t t_l3eta;
0205 Double_t t_l3phi;
0206 Double_t t_p;
0207 Double_t t_pt;
0208 Double_t t_phi;
0209 Double_t t_mindR1;
0210 Double_t t_mindR2;
0211 Double_t t_eMipDR;
0212 Double_t t_eHcal;
0213 Double_t t_eHcal10;
0214 Double_t t_eHcal30;
0215 Double_t t_hmaxNearP;
0216 Double_t t_rhoh;
0217 Bool_t t_selectTk;
0218 Bool_t t_qltyFlag;
0219 Bool_t t_qltyMissFlag;
0220 Bool_t t_qltyPVFlag;
0221 Double_t t_gentrackP;
0222 std::vector<unsigned int> *t_DetIds;
0223 std::vector<double> *t_HitEnergies;
0224 std::vector<bool> *t_trgbits;
0225 std::vector<unsigned int> *t_DetIds1;
0226 std::vector<unsigned int> *t_DetIds3;
0227 std::vector<double> *t_HitEnergies1;
0228 std::vector<double> *t_HitEnergies3;
0229
0230
0231 TBranch *b_t_Run;
0232 TBranch *b_t_Event;
0233 TBranch *b_t_DataType;
0234 TBranch *b_t_ieta;
0235 TBranch *b_t_iphi;
0236 TBranch *b_t_EventWeight;
0237 TBranch *b_t_nVtx;
0238 TBranch *b_t_nTrk;
0239 TBranch *b_t_goodPV;
0240 TBranch *b_t_l1pt;
0241 TBranch *b_t_l1eta;
0242 TBranch *b_t_l1phi;
0243 TBranch *b_t_l3pt;
0244 TBranch *b_t_l3eta;
0245 TBranch *b_t_l3phi;
0246 TBranch *b_t_p;
0247 TBranch *b_t_pt;
0248 TBranch *b_t_phi;
0249 TBranch *b_t_mindR1;
0250 TBranch *b_t_mindR2;
0251 TBranch *b_t_eMipDR;
0252 TBranch *b_t_eHcal;
0253 TBranch *b_t_eHcal10;
0254 TBranch *b_t_eHcal30;
0255 TBranch *b_t_hmaxNearP;
0256 TBranch *b_t_rhoh;
0257 TBranch *b_t_selectTk;
0258 TBranch *b_t_qltyFlag;
0259 TBranch *b_t_qltyMissFlag;
0260 TBranch *b_t_qltyPVFlag;
0261 TBranch *b_t_gentrackP;
0262 TBranch *b_t_DetIds;
0263 TBranch *b_t_HitEnergies;
0264 TBranch *b_t_trgbits;
0265 TBranch *b_t_DetIds1;
0266 TBranch *b_t_DetIds3;
0267 TBranch *b_t_HitEnergies1;
0268 TBranch *b_t_HitEnergies3;
0269
0270 std::string fname_, dirnm_, prefix_;
0271 bool allEvent_;
0272 int flag_;
0273 double mipCut_;
0274 };
0275
0276 CalibSort::CalibSort(const char *fname, std::string dirnm, std::string prefix, bool allEvent, int flag, double mipCut)
0277 : fname_(fname), dirnm_(dirnm), prefix_(prefix), allEvent_(allEvent), flag_(flag), mipCut_(mipCut) {
0278 char treeName[400];
0279 sprintf(treeName, "%s/CalibTree", dirnm.c_str());
0280 TChain *chain = new TChain(treeName);
0281 std::cout << "Create a chain for " << treeName << " from " << fname << std::endl;
0282 if (!fillChain(chain, fname)) {
0283 std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
0284 } else {
0285 std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
0286 Init(chain);
0287 }
0288 }
0289
0290 CalibSort::~CalibSort() {
0291 if (!fChain)
0292 return;
0293 delete fChain->GetCurrentFile();
0294 }
0295
0296 Int_t CalibSort::GetEntry(Long64_t entry) {
0297
0298 if (!fChain)
0299 return 0;
0300 return fChain->GetEntry(entry);
0301 }
0302
0303 Long64_t CalibSort::LoadTree(Long64_t entry) {
0304
0305 if (!fChain)
0306 return -5;
0307 Long64_t centry = fChain->LoadTree(entry);
0308 if (centry < 0)
0309 return centry;
0310 if (!fChain->InheritsFrom(TChain::Class()))
0311 return centry;
0312 TChain *chain = (TChain *)fChain;
0313 if (chain->GetTreeNumber() != fCurrent) {
0314 fCurrent = chain->GetTreeNumber();
0315 Notify();
0316 }
0317 return centry;
0318 }
0319
0320 void CalibSort::Init(TChain *tree) {
0321
0322
0323
0324
0325
0326
0327
0328
0329
0330 t_DetIds = 0;
0331 t_DetIds1 = 0;
0332 t_DetIds3 = 0;
0333 t_HitEnergies = 0;
0334 t_HitEnergies1 = 0;
0335 t_HitEnergies3 = 0;
0336 t_trgbits = 0;
0337
0338 fChain = tree;
0339 fCurrent = -1;
0340 if (!tree)
0341 return;
0342 fChain->SetMakeClass(1);
0343
0344 fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
0345 fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
0346 fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
0347 fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
0348 fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
0349 fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
0350 fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
0351 fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
0352 fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
0353 fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
0354 fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
0355 fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
0356 fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
0357 fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
0358 fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
0359 fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
0360 fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
0361 fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
0362 fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
0363 fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
0364 fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
0365 fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
0366 fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
0367 fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
0368 fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
0369 fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
0370 fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
0371 fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
0372 fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
0373 fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
0374 fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
0375 fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
0376 fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
0377 fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
0378 fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
0379 fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
0380 fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
0381 fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
0382
0383 Notify();
0384 }
0385
0386 Bool_t CalibSort::Notify() {
0387
0388
0389
0390
0391
0392
0393 return kTRUE;
0394 }
0395
0396 void CalibSort::Show(Long64_t entry) {
0397
0398
0399 if (!fChain)
0400 return;
0401 fChain->Show(entry);
0402 }
0403
0404 Int_t CalibSort::Cut(Long64_t) {
0405
0406
0407
0408 return 1;
0409 }
0410
0411 void CalibSort::Loop(const char *outFile) {
0412
0413
0414
0415
0416
0417
0418
0419
0420
0421
0422
0423
0424
0425
0426
0427
0428
0429
0430
0431
0432
0433
0434
0435 if (fChain == 0)
0436 return;
0437
0438 std::ofstream fileout;
0439 if ((flag_ % 10) == 1) {
0440 fileout.open(outFile, std::ofstream::out);
0441 std::cout << "Opens " << outFile << " in output mode" << std::endl;
0442 } else {
0443 fileout.open(outFile, std::ofstream::app);
0444 std::cout << "Opens " << outFile << " in append mode" << std::endl;
0445 }
0446 if (!allEvent_)
0447 fileout << "Input file: " << fname_ << " Directory: " << dirnm_ << " Prefix: " << prefix_ << std::endl;
0448 Int_t runLow(99999999), runHigh(0);
0449 Long64_t nbytes(0), nb(0), good(0);
0450 Long64_t nentries = fChain->GetEntriesFast();
0451 for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0452 Long64_t ientry = LoadTree(jentry);
0453 if (ientry < 0)
0454 break;
0455 nb = fChain->GetEntry(jentry);
0456 nbytes += nb;
0457 if (t_Run > 200000 && t_Run < 800000) {
0458 if (t_Run < runLow)
0459 runLow = t_Run;
0460 if (t_Run > runHigh)
0461 runHigh = t_Run;
0462 }
0463 double cut = (t_p > 20) ? 10.0 : 0.0;
0464 if ((flag_ / 10) % 10 > 0)
0465 std::cout << "Entry " << jentry << " p " << t_p << " Cuts " << t_qltyFlag << "|" << t_selectTk << "|"
0466 << (t_hmaxNearP < cut) << "|" << (t_eMipDR < mipCut_) << std::endl;
0467 if ((t_qltyFlag && t_selectTk && (t_hmaxNearP < cut) && (t_eMipDR < mipCut_)) || allEvent_) {
0468 good++;
0469 fileout << good << " " << jentry << " " << t_Run << " " << t_Event << " " << t_ieta << " " << t_p << std::endl;
0470 }
0471 }
0472 fileout.close();
0473 std::cout << "Writes " << good << " events in the file " << outFile << " from " << nentries
0474 << " entries in run range " << runLow << ":" << runHigh << std::endl;
0475 }
0476
0477 class CalibSortEvent {
0478 public:
0479 TChain *fChain;
0480 Int_t fCurrent;
0481
0482
0483 UInt_t t_RunNo;
0484 UInt_t t_EventNo;
0485 Int_t t_Tracks;
0486 Int_t t_TracksProp;
0487 Int_t t_TracksSaved;
0488 Int_t t_TracksLoose;
0489 Int_t t_TracksTight;
0490 Int_t t_allvertex;
0491 Bool_t t_TrigPass;
0492 Bool_t t_TrigPassSel;
0493 Bool_t t_L1Bit;
0494 std::vector<Bool_t> *t_hltbits;
0495 std::vector<int> *t_ietaAll;
0496 std::vector<int> *t_ietaGood;
0497 std::vector<int> *t_trackType;
0498
0499
0500 TBranch *b_t_RunNo;
0501 TBranch *b_t_EventNo;
0502 TBranch *b_t_Tracks;
0503 TBranch *b_t_TracksProp;
0504 TBranch *b_t_TracksSaved;
0505 TBranch *b_t_TracksLoose;
0506 TBranch *b_t_TracksTight;
0507 TBranch *b_t_allvertex;
0508 TBranch *b_t_TrigPass;
0509 TBranch *b_t_TrigPassSel;
0510 TBranch *b_t_L1Bit;
0511 TBranch *b_t_hltbits;
0512 TBranch *b_t_ietaAll;
0513 TBranch *b_t_ietaGood;
0514 TBranch *b_t_trackType;
0515
0516 CalibSortEvent(const char *fname, std::string dirname, std::string prefix = "", bool append = false);
0517 virtual ~CalibSortEvent();
0518 virtual Int_t Cut(Long64_t entry);
0519 virtual Int_t GetEntry(Long64_t entry);
0520 virtual Long64_t LoadTree(Long64_t entry);
0521 virtual void Init(TChain *tree);
0522 virtual void Loop();
0523 virtual Bool_t Notify();
0524 virtual void Show(Long64_t entry = -1);
0525
0526 private:
0527 std::string fname_, dirnm_, prefix_;
0528 bool append_;
0529 };
0530
0531 CalibSortEvent::CalibSortEvent(const char *fname, std::string dirnm, std::string prefix, bool append)
0532 : fname_(fname), dirnm_(dirnm), prefix_(prefix), append_(append) {
0533 char treeName[400];
0534 sprintf(treeName, "%s/EventInfo", dirnm.c_str());
0535 TChain *chain = new TChain(treeName);
0536 std::cout << "Create a chain for " << treeName << " from " << fname << std::endl;
0537 if (!fillChain(chain, fname)) {
0538 std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
0539 } else {
0540 std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
0541 Init(chain);
0542 }
0543 }
0544
0545 CalibSortEvent::~CalibSortEvent() {
0546 if (!fChain)
0547 return;
0548 delete fChain->GetCurrentFile();
0549 }
0550
0551 Int_t CalibSortEvent::GetEntry(Long64_t entry) {
0552
0553 if (!fChain)
0554 return 0;
0555 return fChain->GetEntry(entry);
0556 }
0557
0558 Long64_t CalibSortEvent::LoadTree(Long64_t entry) {
0559
0560 if (!fChain)
0561 return -5;
0562 Long64_t centry = fChain->LoadTree(entry);
0563 if (centry < 0)
0564 return centry;
0565 if (!fChain->InheritsFrom(TChain::Class()))
0566 return centry;
0567 TChain *chain = (TChain *)fChain;
0568 if (chain->GetTreeNumber() != fCurrent) {
0569 fCurrent = chain->GetTreeNumber();
0570 Notify();
0571 }
0572 return centry;
0573 }
0574
0575 void CalibSortEvent::Init(TChain *tree) {
0576
0577
0578
0579
0580
0581
0582
0583
0584
0585
0586 t_hltbits = 0;
0587 t_ietaAll = 0;
0588 t_ietaGood = 0;
0589 t_trackType = 0;
0590 t_L1Bit = false;
0591 fChain = tree;
0592 fCurrent = -1;
0593 if (!tree)
0594 return;
0595 fChain->SetMakeClass(1);
0596 fChain->SetBranchAddress("t_RunNo", &t_RunNo, &b_t_RunNo);
0597 fChain->SetBranchAddress("t_EventNo", &t_EventNo, &b_t_EventNo);
0598 fChain->SetBranchAddress("t_Tracks", &t_Tracks, &b_t_Tracks);
0599 fChain->SetBranchAddress("t_TracksProp", &t_TracksProp, &b_t_TracksProp);
0600 fChain->SetBranchAddress("t_TracksSaved", &t_TracksSaved, &b_t_TracksSaved);
0601 fChain->SetBranchAddress("t_TracksLoose", &t_TracksLoose, &b_t_TracksLoose);
0602 fChain->SetBranchAddress("t_TracksTight", &t_TracksTight, &b_t_TracksTight);
0603 fChain->SetBranchAddress("t_allvertex", &t_allvertex, &b_t_allvertex);
0604 fChain->SetBranchAddress("t_TrigPass", &t_TrigPass, &b_t_TrigPass);
0605 fChain->SetBranchAddress("t_TrigPassSel", &t_TrigPassSel, &b_t_TrigPassSel);
0606 fChain->SetBranchAddress("t_L1Bit", &t_L1Bit, &b_t_L1Bit);
0607 fChain->SetBranchAddress("t_hltbits", &t_hltbits, &b_t_hltbits);
0608 fChain->SetBranchAddress("t_ietaAll", &t_ietaAll, &b_t_ietaAll);
0609 fChain->SetBranchAddress("t_ietaGood", &t_ietaGood, &b_t_ietaGood);
0610 fChain->SetBranchAddress("t_trackType", &t_trackType, &b_t_trackType);
0611 Notify();
0612 }
0613
0614 Bool_t CalibSortEvent::Notify() {
0615
0616
0617
0618
0619
0620
0621 return kTRUE;
0622 }
0623
0624 void CalibSortEvent::Show(Long64_t entry) {
0625
0626
0627 if (!fChain)
0628 return;
0629 fChain->Show(entry);
0630 }
0631
0632 Int_t CalibSortEvent::Cut(Long64_t) {
0633
0634
0635
0636 return 1;
0637 }
0638
0639 void CalibSortEvent::Loop() {
0640
0641
0642
0643
0644
0645
0646
0647
0648
0649
0650
0651
0652
0653
0654
0655
0656
0657
0658
0659
0660
0661
0662
0663 if (fChain == 0)
0664 return;
0665
0666 std::ofstream fileout;
0667 if (!append_) {
0668 fileout.open("runevents.txt", std::ofstream::out);
0669 std::cout << "Opens runevents.txt in output mode" << std::endl;
0670 } else {
0671 fileout.open("runevents.txt", std::ofstream::app);
0672 std::cout << "Opens runevents.txt in append mode" << std::endl;
0673 }
0674 fileout << "Input file: " << fname_ << " Directory: " << dirnm_ << " Prefix: " << prefix_ << std::endl;
0675 UInt_t runLow(99999999), runHigh(0);
0676 Long64_t nbytes(0), nb(0), good(0);
0677 Long64_t nentries = fChain->GetEntriesFast();
0678 for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0679 Long64_t ientry = LoadTree(jentry);
0680 if (ientry < 0)
0681 break;
0682 nb = fChain->GetEntry(jentry);
0683 nbytes += nb;
0684 if (t_RunNo > 200000 && t_RunNo < 800000) {
0685 if (t_RunNo < runLow)
0686 runLow = t_RunNo;
0687 if (t_RunNo > runHigh)
0688 runHigh = t_RunNo;
0689 }
0690 good++;
0691 fileout << good << " " << jentry << " " << t_RunNo << " " << t_EventNo << std::endl;
0692 }
0693 fileout.close();
0694 std::cout << "Writes " << good << " events in the file events.txt from " << nentries << " entries in run range "
0695 << runLow << ":" << runHigh << std::endl;
0696 }
0697
0698 void readRecords(std::string fname, std::vector<record> &records, bool debug) {
0699 records.clear();
0700 ifstream infile(fname.c_str());
0701 if (!infile.is_open()) {
0702 std::cout << "Cannot open " << fname << std::endl;
0703 } else {
0704 while (1) {
0705 int ser, ent, r, ev, ie;
0706 double p;
0707 infile >> ser >> ent >> r >> ev >> ie >> p;
0708 if (!infile.good())
0709 break;
0710 record rec(ser, ent, r, ev, ie, p);
0711 records.push_back(rec);
0712 }
0713 infile.close();
0714 }
0715 std::cout << "Reads " << records.size() << " records from " << fname << std::endl;
0716 if (debug) {
0717 for (unsigned int k = 0; k < records.size(); ++k) {
0718 if (k % 100 == 0)
0719 std::cout << "[" << records[k].serial_ << ":" << records[k].entry_ << "] " << records[k].run_ << ":"
0720 << records[k].event_ << " " << records[k].ieta_ << " " << records[k].p_ << std::endl;
0721 }
0722 }
0723 }
0724
0725 void readMap(std::string fname, std::map<std::pair<int, int>, int> &records, bool debug) {
0726 records.clear();
0727 ifstream infile(fname.c_str());
0728 if (!infile.is_open()) {
0729 std::cout << "Cannot open " << fname << std::endl;
0730 } else {
0731 while (1) {
0732 int ser, ent, r, ev, ie;
0733 double p;
0734 infile >> ser >> ent >> r >> ev >> ie >> p;
0735 if (!infile.good())
0736 break;
0737 std::pair<int, int> key(r, ev);
0738 if (records.find(key) == records.end())
0739 records[key] = ent;
0740 }
0741 infile.close();
0742 }
0743 std::cout << "Reads " << records.size() << " records from " << fname << std::endl;
0744 if (debug) {
0745 unsigned k(0);
0746 for (std::map<std::pair<int, int>, int>::iterator itr = records.begin(); itr != records.end(); ++itr, ++k) {
0747 if (k % 100 == 0)
0748 std::cout << "[" << k << "] " << itr->second << ":" << (itr->first).first << ":" << (itr->first).second << "\n";
0749 }
0750 }
0751 }
0752
0753 void sort(std::vector<record> &records, bool debug) {
0754
0755 std::sort(records.begin(), records.end(), recordLess());
0756 if (debug) {
0757 for (unsigned int k = 0; k < records.size(); ++k) {
0758 std::cout << "[" << k << ":" << records[k].serial_ << ":" << records[k].entry_ << "] " << records[k].run_ << ":"
0759 << records[k].event_ << " " << records[k].ieta_ << " " << records[k].p_ << std::endl;
0760 }
0761 }
0762 }
0763
0764 void duplicate(std::string fname, std::vector<record> &records, bool debug) {
0765 std::ofstream file;
0766 file.open(fname.c_str(), std::ofstream::out);
0767 std::cout << "List of entry names of duplicate events" << std::endl;
0768 int duplicate(0), dupl40(0);
0769 for (unsigned int k = 1; k < records.size(); ++k) {
0770 if ((records[k].run_ == records[k - 1].run_) && (records[k].event_ == records[k - 1].event_) &&
0771 (records[k].ieta_ == records[k - 1].ieta_) && (fabs(records[k].p_ - records[k - 1].p_) < 0.0001)) {
0772
0773 if (records[k].entry_ < records[k - 1].entry_) {
0774 record swap = records[k - 1];
0775 records[k - 1] = records[k];
0776 records[k] = swap;
0777 }
0778 if (debug) {
0779 std::cout << "Serial " << records[k - 1].serial_ << ":" << records[k].serial_ << " Entry "
0780 << records[k - 1].entry_ << ":" << records[k].entry_ << " Run " << records[k - 1].run_ << ":"
0781 << records[k].run_ << " Event " << records[k - 1].event_ << " " << records[k].event_ << " Eta "
0782 << records[k - 1].ieta_ << " " << records[k].ieta_ << " p " << records[k - 1].p_ << ":"
0783 << records[k].p_ << std::endl;
0784 }
0785 file << records[k].entry_ << std::endl;
0786 duplicate++;
0787 if (records[k].p_ >= 40.0 && records[k].p_ <= 60.0)
0788 dupl40++;
0789 }
0790 }
0791 file.close();
0792 std::cout << "Total # of duplcate events " << duplicate << " (" << dupl40 << " with p 40:60)" << std::endl;
0793 }
0794
0795 void findDuplicate(std::string infile, std::string outfile, bool debug = false) {
0796 std::vector<record> records;
0797 readRecords(infile, records, debug);
0798 sort(records, debug);
0799 duplicate(outfile, records, debug);
0800 }
0801
0802 void findCommon(std::string infile1, std::string infile2, std::string infile3, std::string outfile, bool debug = false) {
0803 std::map<std::pair<int, int>, int> map1, map2, map3;
0804 readMap(infile1, map1, debug);
0805 readMap(infile2, map2, debug);
0806 readMap(infile3, map3, debug);
0807 bool check3 = (map3.size() > 0);
0808 std::ofstream file;
0809 file.open(outfile.c_str(), std::ofstream::out);
0810 unsigned int k(0), good(0);
0811 for (std::map<std::pair<int, int>, int>::iterator itr = map1.begin(); itr != map1.end(); ++itr, ++k) {
0812 std::pair<int, int> key = itr->first;
0813 bool ok = (map2.find(key) != map2.end());
0814 if (ok && check3)
0815 ok = (map3.find(key) != map3.end());
0816 if (debug && k % 100 == 0)
0817 std::cout << "[" << k << "] Run " << key.first << " Event " << key.second << " Flag " << ok << std::endl;
0818 if (ok) {
0819 ++good;
0820 file << key.first << " " << key.second << std::endl;
0821 }
0822 }
0823 file.close();
0824 std::cout << "Total # of common events " << good << " written to o/p file " << outfile << std::endl;
0825 }
0826
0827 void readRecordEvents(std::string fname, std::vector<recordEvent> &records, bool debug) {
0828 records.clear();
0829 ifstream infile(fname.c_str());
0830 if (!infile.is_open()) {
0831 std::cout << "Cannot open " << fname << std::endl;
0832 } else {
0833 while (1) {
0834 unsigned int ser, ent, r, ev;
0835 infile >> ser >> ent >> r >> ev;
0836 if (!infile.good())
0837 break;
0838 recordEvent rec(ser, ent, r, ev);
0839 records.push_back(rec);
0840 }
0841 infile.close();
0842 }
0843 std::cout << "Reads " << records.size() << " records from " << fname << std::endl;
0844 if (debug) {
0845 for (unsigned int k = 0; k < records.size(); ++k) {
0846 if (k % 100 == 0)
0847 std::cout << "[" << records[k].serial_ << ":" << records[k].entry_ << "] " << records[k].run_ << ":"
0848 << records[k].event_ << std::endl;
0849 }
0850 }
0851 }
0852
0853 void sortEvent(std::vector<recordEvent> &records, bool debug) {
0854
0855 std::sort(records.begin(), records.end(), recordEventLess());
0856 if (debug) {
0857 for (unsigned int k = 0; k < records.size(); ++k) {
0858 std::cout << "[" << k << ":" << records[k].serial_ << ":" << records[k].entry_ << "] " << records[k].run_ << ":"
0859 << records[k].event_ << std::endl;
0860 }
0861 }
0862 }
0863
0864 void duplicateEvent(std::string fname, std::vector<recordEvent> &records, bool debug) {
0865 std::ofstream file;
0866 file.open(fname.c_str(), std::ofstream::out);
0867 std::cout << "List of entry names of duplicate events" << std::endl;
0868 int duplicate(0);
0869 for (unsigned int k = 1; k < records.size(); ++k) {
0870 if ((records[k].run_ == records[k - 1].run_) && (records[k].event_ == records[k - 1].event_)) {
0871
0872 if (records[k].entry_ < records[k - 1].entry_) {
0873 recordEvent swap = records[k - 1];
0874 records[k - 1] = records[k];
0875 records[k] = swap;
0876 }
0877 if (debug) {
0878 std::cout << "Serial " << records[k - 1].serial_ << ":" << records[k].serial_ << " Entry "
0879 << records[k - 1].entry_ << ":" << records[k].entry_ << " Run " << records[k - 1].run_ << ":"
0880 << records[k].run_ << " Event " << records[k - 1].event_ << " " << records[k].event_ << std::endl;
0881 }
0882 file << records[k].entry_ << std::endl;
0883 duplicate++;
0884 }
0885 }
0886 file.close();
0887 std::cout << "Total # of duplcate events " << duplicate << std::endl;
0888 }
0889
0890 void findDuplicateEvent(std::string infile, std::string outfile, bool debug = false) {
0891 std::vector<recordEvent> records;
0892 readRecordEvents(infile, records, debug);
0893 sortEvent(records, debug);
0894 duplicateEvent(outfile, records, debug);
0895 }
0896
0897 void combine(const std::string fname1,
0898 const std::string fname2,
0899 const std::string fname,
0900 const std::string outfile = "",
0901 int debug = 0) {
0902 ifstream infile1(fname1.c_str());
0903 ifstream infile2(fname2.c_str());
0904 if ((!infile1.is_open()) || (!infile2.is_open())) {
0905 std::cout << "Cannot open " << fname1 << " or " << fname2 << std::endl;
0906 } else {
0907 std::ofstream file;
0908 file.open(fname.c_str(), std::ofstream::out);
0909 TH1D *hist0 = new TH1D("W0", "Correction Factor (All)", 100, 0.0, 10.0);
0910 TH1D *hist1 = new TH1D("W1", "Correction Factor (Barrel)", 100, 0.0, 10.0);
0911 TH1D *hist2 = new TH1D("W2", "Correction Factor (Endcap)", 100, 0.0, 10.0);
0912 TProfile *prof = new TProfile("P", "Correction vs i#eta", 60, -30, 30, 0, 100, "");
0913 unsigned int kount(0), kout(0);
0914 double wmaxb(0), wmaxe(0);
0915 while (1) {
0916 Long64_t entry1, entry2;
0917 int ieta1, ieta2;
0918 double wt1, wt2;
0919 infile1 >> entry1 >> ieta1 >> wt1;
0920 infile2 >> entry2 >> ieta2 >> wt2;
0921 if ((!infile1.good()) || (!infile2.good()))
0922 break;
0923 ++kount;
0924 if (debug > 1) {
0925 std::cout << kount << " Enrty No. " << entry1 << ":" << entry2 << " eta " << ieta1 << ":" << ieta2 << " wt "
0926 << wt1 << ":" << wt2 << std::endl;
0927 }
0928 if (entry1 == entry2) {
0929 int ieta = fabs(ieta1);
0930 double wt = (ieta >= 16) ? wt1 : wt2;
0931 file << std::setw(8) << entry1 << " " << std::setw(12) << std::setprecision(8) << wt << std::endl;
0932 hist0->Fill(wt);
0933 if (ieta >= 16) {
0934 hist2->Fill(wt);
0935 if (wt > wmaxe)
0936 wmaxe = wt;
0937 } else {
0938 hist1->Fill(wt);
0939 if (wt > wmaxb)
0940 wmaxb = wt;
0941 }
0942 prof->Fill(ieta1, wt);
0943 ++kout;
0944 } else if (debug > 0) {
0945 std::cout << kount << " Entry " << entry1 << ":" << entry2 << " eta " << ieta1 << ":" << ieta2 << " wt " << wt1
0946 << ":" << wt2 << " mismatch" << std::endl;
0947 }
0948 }
0949 infile1.close();
0950 infile2.close();
0951 file.close();
0952 std::cout << "Writes " << kout << " entries to " << fname << " from " << kount << " events from " << fname1
0953 << " and " << fname2 << " maximum correction factor " << wmaxb << " (Barrel) " << wmaxe << " (Endcap)"
0954 << std::endl;
0955 if (outfile != "") {
0956 TFile *theFile = new TFile(outfile.c_str(), "RECREATE");
0957 theFile->cd();
0958 hist0->Write();
0959 hist1->Write();
0960 hist2->Write();
0961 prof->Write();
0962 theFile->Close();
0963 }
0964 }
0965 }
0966
0967 void plotCombine(const char *infile, int flag = -1, bool drawStatBox = true, bool save = true) {
0968 int flag1 = (flag >= 0 ? (flag / 10) % 10 : 0);
0969 int flag2 = (flag >= 0 ? (flag % 10) : 0);
0970 std::string name[4] = {"W0", "W1", "W2", "P"};
0971 std::string xtitl[4] = {
0972 "Correction Factor (All)", "Correction Factor (Barrel)", "Correction Factor (Endcap)", "i#eta"};
0973 std::string ytitl[4] = {"Track", "Track", "Track", "Correction Factor"};
0974 char title[100];
0975 std::string mom[2] = {"all momentum", "p = 40:60 GeV"};
0976 std::string iso[2] = {"loose", "tight"};
0977 if (flag < 0) {
0978 sprintf(title, "All tracks");
0979 } else {
0980 sprintf(title, "Good tracks of %s with %s isolation", mom[flag1].c_str(), iso[flag2].c_str());
0981 }
0982
0983 gStyle->SetCanvasBorderMode(0);
0984 gStyle->SetCanvasColor(kWhite);
0985 gStyle->SetPadColor(kWhite);
0986 gStyle->SetFillColor(kWhite);
0987 gStyle->SetOptTitle(0);
0988 gStyle->SetOptStat(111110);
0989 gStyle->SetOptFit(0);
0990 TFile *file = new TFile(infile);
0991 for (int k = 0; k < 4; ++k) {
0992 char nameh[40], namep[50];
0993 if (flag >= 0)
0994 sprintf(nameh, "%s%d%d", name[k].c_str(), flag1, flag2);
0995 else
0996 sprintf(nameh, "%s", name[k].c_str());
0997 sprintf(namep, "c_%s", nameh);
0998 TCanvas *pad(nullptr);
0999 if (k == 3) {
1000 TProfile *hist = (TProfile *)file->FindObjectAny(nameh);
1001 if (hist != nullptr) {
1002 pad = new TCanvas(namep, namep, 700, 500);
1003 pad->SetRightMargin(0.10);
1004 pad->SetTopMargin(0.10);
1005 hist->GetXaxis()->SetTitleSize(0.04);
1006 hist->GetXaxis()->SetTitle(xtitl[k].c_str());
1007 hist->GetYaxis()->SetTitle(ytitl[k].c_str());
1008 hist->GetYaxis()->SetLabelOffset(0.005);
1009 hist->GetYaxis()->SetTitleSize(0.04);
1010 hist->GetYaxis()->SetLabelSize(0.035);
1011 hist->GetYaxis()->SetTitleOffset(1.10);
1012 hist->SetMarkerStyle(20);
1013 hist->Draw();
1014 pad->Update();
1015 TPaveStats *st1 = (TPaveStats *)hist->GetListOfFunctions()->FindObject("stats");
1016 if (st1 != nullptr && drawStatBox) {
1017 st1->SetY1NDC(0.70);
1018 st1->SetY2NDC(0.90);
1019 st1->SetX1NDC(0.65);
1020 st1->SetX2NDC(0.90);
1021 }
1022 }
1023 } else {
1024 TH1D *hist = (TH1D *)file->FindObjectAny(nameh);
1025 if (hist != nullptr) {
1026 pad = new TCanvas(namep, namep, 700, 500);
1027 pad->SetRightMargin(0.10);
1028 pad->SetTopMargin(0.10);
1029 pad->SetLogy();
1030 hist->GetXaxis()->SetTitleSize(0.04);
1031 hist->GetXaxis()->SetTitle(xtitl[k].c_str());
1032 hist->GetYaxis()->SetTitle(ytitl[k].c_str());
1033 hist->GetYaxis()->SetLabelOffset(0.005);
1034 hist->GetYaxis()->SetTitleSize(0.04);
1035 hist->GetYaxis()->SetLabelSize(0.035);
1036 hist->GetYaxis()->SetTitleOffset(1.10);
1037 hist->SetMarkerStyle(20);
1038 hist->Draw();
1039 pad->Update();
1040 TPaveStats *st1 = (TPaveStats *)hist->GetListOfFunctions()->FindObject("stats");
1041 if (st1 != nullptr && drawStatBox) {
1042 st1->SetY1NDC(0.70);
1043 st1->SetY2NDC(0.90);
1044 st1->SetX1NDC(0.65);
1045 st1->SetX2NDC(0.90);
1046 }
1047 }
1048 TPaveText *txt0 = new TPaveText(0.10, 0.91, 0.90, 0.96, "blNDC");
1049 txt0->SetFillColor(0);
1050 txt0->AddText(title);
1051 txt0->Draw("same");
1052 pad->Update();
1053 }
1054 if (pad != nullptr) {
1055 pad->Modified();
1056 pad->Update();
1057 if (save) {
1058 sprintf(namep, "%s.pdf", pad->GetName());
1059 pad->Print(namep);
1060 }
1061 }
1062 }
1063 }
1064
1065 class CalibPlotCombine {
1066 public:
1067 TChain *fChain;
1068 Int_t fCurrent;
1069
1070
1071 Int_t t_Run;
1072 Int_t t_Event;
1073 Int_t t_DataType;
1074 Int_t t_ieta;
1075 Int_t t_iphi;
1076 Double_t t_EventWeight;
1077 Int_t t_nVtx;
1078 Int_t t_nTrk;
1079 Int_t t_goodPV;
1080 Double_t t_l1pt;
1081 Double_t t_l1eta;
1082 Double_t t_l1phi;
1083 Double_t t_l3pt;
1084 Double_t t_l3eta;
1085 Double_t t_l3phi;
1086 Double_t t_p;
1087 Double_t t_pt;
1088 Double_t t_phi;
1089 Double_t t_mindR1;
1090 Double_t t_mindR2;
1091 Double_t t_eMipDR;
1092 Double_t t_eHcal;
1093 Double_t t_eHcal10;
1094 Double_t t_eHcal30;
1095 Double_t t_hmaxNearP;
1096 Double_t t_rhoh;
1097 Bool_t t_selectTk;
1098 Bool_t t_qltyFlag;
1099 Bool_t t_qltyMissFlag;
1100 Bool_t t_qltyPVFlag;
1101 Double_t t_gentrackP;
1102 std::vector<unsigned int> *t_DetIds;
1103 std::vector<double> *t_HitEnergies;
1104 std::vector<bool> *t_trgbits;
1105 std::vector<unsigned int> *t_DetIds1;
1106 std::vector<unsigned int> *t_DetIds3;
1107 std::vector<double> *t_HitEnergies1;
1108 std::vector<double> *t_HitEnergies3;
1109
1110
1111 TBranch *b_t_Run;
1112 TBranch *b_t_Event;
1113 TBranch *b_t_DataType;
1114 TBranch *b_t_ieta;
1115 TBranch *b_t_iphi;
1116 TBranch *b_t_EventWeight;
1117 TBranch *b_t_nVtx;
1118 TBranch *b_t_nTrk;
1119 TBranch *b_t_goodPV;
1120 TBranch *b_t_l1pt;
1121 TBranch *b_t_l1eta;
1122 TBranch *b_t_l1phi;
1123 TBranch *b_t_l3pt;
1124 TBranch *b_t_l3eta;
1125 TBranch *b_t_l3phi;
1126 TBranch *b_t_p;
1127 TBranch *b_t_pt;
1128 TBranch *b_t_phi;
1129 TBranch *b_t_mindR1;
1130 TBranch *b_t_mindR2;
1131 TBranch *b_t_eMipDR;
1132 TBranch *b_t_eHcal;
1133 TBranch *b_t_eHcal10;
1134 TBranch *b_t_eHcal30;
1135 TBranch *b_t_hmaxNearP;
1136 TBranch *b_t_rhoh;
1137 TBranch *b_t_selectTk;
1138 TBranch *b_t_qltyFlag;
1139 TBranch *b_t_qltyMissFlag;
1140 TBranch *b_t_qltyPVFlag;
1141 TBranch *b_t_gentrackP;
1142 TBranch *b_t_DetIds;
1143 TBranch *b_t_HitEnergies;
1144 TBranch *b_t_trgbits;
1145 TBranch *b_t_DetIds1;
1146 TBranch *b_t_DetIds3;
1147 TBranch *b_t_HitEnergies1;
1148 TBranch *b_t_HitEnergies3;
1149
1150 CalibPlotCombine(const char *fname,
1151 const char *corrFileName,
1152 const std::string dirname = "HcalIsoTrkAnalyzer",
1153 int flag = 10);
1154 virtual ~CalibPlotCombine();
1155 virtual Int_t Cut(Long64_t entry);
1156 virtual Int_t GetEntry(Long64_t entry);
1157 virtual Long64_t LoadTree(Long64_t entry);
1158 virtual void Init(TChain *);
1159 virtual void Loop();
1160 virtual Bool_t Notify();
1161 virtual void Show(Long64_t entry = -1);
1162 bool goodTrack(double eHcal, double cut);
1163 void savePlot(const std::string &theName, bool append = false);
1164
1165 private:
1166 CalibCorr *cFactor_;
1167 int flag_, ifDepth_;
1168 bool selectP_, tightSelect_;
1169 TH1D *hist0_, *hist1_, *hist2_;
1170 TProfile *prof_;
1171 };
1172
1173 CalibPlotCombine::CalibPlotCombine(const char *fname, const char *corrFileName, const std::string dirnm, int flag)
1174 : cFactor_(nullptr), flag_(flag), ifDepth_(3) {
1175 selectP_ = (((flag_ / 10) % 10) > 0);
1176 tightSelect_ = ((flag_ % 10) > 0);
1177 flag_ = 0;
1178 if (selectP_)
1179 flag_ += 10;
1180 if (tightSelect_)
1181 ++flag_;
1182 std::cout << "Selection on momentum range: " << selectP_ << std::endl
1183 << "Tight isolation flag: " << tightSelect_ << std::endl
1184 << "Flag: " << flag_ << std::endl;
1185 char treeName[400];
1186 sprintf(treeName, "%s/CalibTree", dirnm.c_str());
1187 TChain *chain = new TChain(treeName);
1188 std::cout << "Create a chain for " << treeName << " from " << fname << std::endl;
1189 if (!fillChain(chain, fname)) {
1190 std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
1191 } else {
1192 std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
1193 Init(chain);
1194 }
1195 cFactor_ = new CalibCorr(corrFileName, ifDepth_, false);
1196 }
1197
1198 CalibPlotCombine::~CalibPlotCombine() {
1199 delete cFactor_;
1200 if (!fChain)
1201 return;
1202 delete fChain->GetCurrentFile();
1203 }
1204
1205 Int_t CalibPlotCombine::GetEntry(Long64_t entry) {
1206
1207 if (!fChain)
1208 return 0;
1209 return fChain->GetEntry(entry);
1210 }
1211
1212 Long64_t CalibPlotCombine::LoadTree(Long64_t entry) {
1213
1214 if (!fChain)
1215 return -5;
1216 Long64_t centry = fChain->LoadTree(entry);
1217 if (centry < 0)
1218 return centry;
1219 if (!fChain->InheritsFrom(TChain::Class()))
1220 return centry;
1221 TChain *chain = (TChain *)fChain;
1222 if (chain->GetTreeNumber() != fCurrent) {
1223 fCurrent = chain->GetTreeNumber();
1224 Notify();
1225 }
1226 return centry;
1227 }
1228
1229 void CalibPlotCombine::Init(TChain *tree) {
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239 t_DetIds = 0;
1240 t_DetIds1 = 0;
1241 t_DetIds3 = 0;
1242 t_HitEnergies = 0;
1243 t_HitEnergies1 = 0;
1244 t_HitEnergies3 = 0;
1245 t_trgbits = 0;
1246
1247 fChain = tree;
1248 fCurrent = -1;
1249 if (!tree)
1250 return;
1251 fChain->SetMakeClass(1);
1252
1253 fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
1254 fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
1255 fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
1256 fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
1257 fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
1258 fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
1259 fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
1260 fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
1261 fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
1262 fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
1263 fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
1264 fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
1265 fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
1266 fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
1267 fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
1268 fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
1269 fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
1270 fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
1271 fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
1272 fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
1273 fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
1274 fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
1275 fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
1276 fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
1277 fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
1278 fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
1279 fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
1280 fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
1281 fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
1282 fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
1283 fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
1284 fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
1285 fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
1286 fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
1287 fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
1288 fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
1289 fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
1290 fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
1291 Notify();
1292
1293
1294 char name[20], title[100];
1295 std::string mom[2] = {"all momentum", "p = 40:60 GeV"};
1296 std::string iso[2] = {"loose", "tight"};
1297 std::string rng[3] = {"All", "Barrel", "Endcap"};
1298 int flag1 = (flag_ / 10) % 10;
1299 int flag2 = (flag_ % 10);
1300 sprintf(name, "W0%d%d", flag1, flag2);
1301 sprintf(title,
1302 "Correction Factor (%s) for tracks with %s and %s isolation",
1303 rng[0].c_str(),
1304 mom[flag1].c_str(),
1305 iso[flag2].c_str());
1306 hist0_ = new TH1D(name, title, 100, 0.0, 10.0);
1307 sprintf(name, "W1%d%d", flag1, flag2);
1308 sprintf(title,
1309 "Correction Factor (%s) for tracks with %s and %s isolation",
1310 rng[1].c_str(),
1311 mom[flag1].c_str(),
1312 iso[flag2].c_str());
1313 hist1_ = new TH1D(name, title, 100, 0.0, 10.0);
1314 sprintf(name, "W2%d%d", flag1, flag2);
1315 sprintf(title,
1316 "Correction Factor (%s) for tracks with %s and %s isolation",
1317 rng[2].c_str(),
1318 mom[flag1].c_str(),
1319 iso[flag2].c_str());
1320 hist2_ = new TH1D(name, title, 100, 0.0, 10.0);
1321 sprintf(name, "P%d%d", flag1, flag2);
1322 sprintf(
1323 title, "Correction Factor vs i#eta for tracks with %s and %s isolation", mom[flag1].c_str(), iso[flag2].c_str());
1324 prof_ = new TProfile(name, title, 60, -30, 30, 0, 100, "");
1325 }
1326
1327 Bool_t CalibPlotCombine::Notify() {
1328
1329
1330
1331
1332
1333
1334 return kTRUE;
1335 }
1336
1337 void CalibPlotCombine::Show(Long64_t entry) {
1338
1339
1340 if (!fChain)
1341 return;
1342 fChain->Show(entry);
1343 }
1344
1345 Int_t CalibPlotCombine::Cut(Long64_t) {
1346
1347
1348
1349 return 1;
1350 }
1351
1352 void CalibPlotCombine::Loop() {
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377 if (fChain == 0)
1378 return;
1379 Long64_t nentries = fChain->GetEntriesFast();
1380 std::cout << "Total entries " << nentries << std::endl;
1381 Long64_t nbytes(0), nb(0);
1382 double wminb(99999999), wmaxb(0), wmine(99999999), wmaxe(0);
1383 for (Long64_t jentry = 0; jentry < nentries; jentry++) {
1384 Long64_t ientry = LoadTree(jentry);
1385 if (ientry < 0)
1386 break;
1387 nb = fChain->GetEntry(jentry);
1388 nbytes += nb;
1389 if (jentry % 1000000 == 0)
1390 std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
1391 double cut = (t_p > 20) ? (tightSelect_ ? 2.0 : 10.0) : 0.0;
1392 double eHcal(t_eHcal);
1393 if (goodTrack(eHcal, cut)) {
1394 if ((!selectP_) || ((t_p > 40.0) && (t_p < 60.0))) {
1395 int ieta = fabs(t_ieta);
1396 double wt = cFactor_->getTrueCorr(ientry);
1397 hist0_->Fill(wt);
1398 if (ieta >= 16) {
1399 hist2_->Fill(wt);
1400 wmaxe = std::max(wmaxe, wt);
1401 wmine = std::min(wmine, wt);
1402 } else {
1403 hist1_->Fill(wt);
1404 wmaxb = std::max(wmaxb, wt);
1405 wminb = std::min(wminb, wt);
1406 }
1407 prof_->Fill(t_ieta, wt);
1408 }
1409 }
1410 }
1411 std::cout << "Minimum and maximum correction factors " << wminb << ":" << wmaxb << " (Barrel) " << wmine << ":"
1412 << wmaxe << " (Endcap)" << std::endl;
1413 }
1414
1415 bool CalibPlotCombine::goodTrack(double eHcal, double cut) {
1416 return ((t_qltyFlag) && t_selectTk && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (eHcal > 0.001));
1417 }
1418
1419 void CalibPlotCombine::savePlot(const std::string &theName, bool append) {
1420 TFile *theFile = (append ? (new TFile(theName.c_str(), "UPDATE")) : (new TFile(theName.c_str(), "RECREATE")));
1421 theFile->cd();
1422 TH1D *hist;
1423 hist = (TH1D *)(hist0_->Clone());
1424 hist->Write();
1425 hist = (TH1D *)(hist1_->Clone());
1426 hist->Write();
1427 hist = (TH1D *)(hist2_->Clone());
1428 hist->Write();
1429 TProfile *prof = (TProfile *)(prof_->Clone());
1430 prof->Write();
1431 std::cout << "All done" << std::endl;
1432 theFile->Close();
1433 }
1434
1435 class CalibFitPU {
1436 private:
1437 TTree *fChain;
1438 Int_t fCurrent;
1439
1440
1441
1442
1443 Int_t t_Event;
1444 Double_t t_p_PU;
1445 Double_t t_eHcal_PU;
1446 Double_t t_delta_PU;
1447 Double_t t_p_NoPU;
1448 Double_t t_eHcal_noPU;
1449 Double_t t_delta_NoPU;
1450 Int_t t_ieta;
1451
1452
1453 TBranch *b_t_Event;
1454 TBranch *b_t_p_PU;
1455 TBranch *b_t_eHcal_PU;
1456 TBranch *b_t_delta_PU;
1457 TBranch *b_t_p_NoPU;
1458 TBranch *b_t_eHcal_noPU;
1459 TBranch *b_t_delta_NoPU;
1460 TBranch *b_t_ieta;
1461
1462 public:
1463 CalibFitPU(const char *fname = "isotrackRelval.root");
1464 virtual ~CalibFitPU();
1465 virtual Int_t Cut(Long64_t entry);
1466 virtual Int_t GetEntry(Long64_t entry);
1467 virtual Long64_t LoadTree(Long64_t entry);
1468 virtual void Init(TTree *tree);
1469 virtual void Loop(bool extract_PU_parameters, std::string fileName);
1470 virtual Bool_t Notify();
1471 virtual void Show(Long64_t entry = -1);
1472 };
1473
1474 CalibFitPU::CalibFitPU(const char *fname) : fChain(0) {
1475
1476
1477 TFile *f = new TFile(fname);
1478 TTree *tree = new TTree();
1479 f->GetObject("tree", tree);
1480 std::cout << "Find tree Tree in " << tree << " from " << fname << std::endl;
1481 Init(tree);
1482 }
1483
1484 CalibFitPU::~CalibFitPU() {
1485 if (!fChain)
1486 return;
1487 delete fChain->GetCurrentFile();
1488 }
1489
1490 Int_t CalibFitPU::GetEntry(Long64_t entry) {
1491
1492 if (!fChain)
1493 return 0;
1494 return fChain->GetEntry(entry);
1495 }
1496
1497 Long64_t CalibFitPU::LoadTree(Long64_t entry) {
1498
1499 if (!fChain)
1500 return -5;
1501 Long64_t centry = fChain->LoadTree(entry);
1502 if (centry < 0)
1503 return centry;
1504 if (fChain->GetTreeNumber() != fCurrent) {
1505 fCurrent = fChain->GetTreeNumber();
1506 Notify();
1507 }
1508 return centry;
1509 }
1510
1511 void CalibFitPU::Init(TTree *tree) {
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521 if (!tree)
1522 return;
1523 fChain = tree;
1524 fCurrent = -1;
1525 fChain->SetMakeClass(1);
1526
1527 fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
1528 fChain->SetBranchAddress("t_p_PU", &t_p_PU, &b_t_p_PU);
1529 fChain->SetBranchAddress("t_eHcal_PU", &t_eHcal_PU, &b_t_eHcal_PU);
1530 fChain->SetBranchAddress("t_delta_PU", &t_delta_PU, &b_t_delta_PU);
1531 fChain->SetBranchAddress("t_p_NoPU", &t_p_NoPU, &b_t_p_NoPU);
1532 fChain->SetBranchAddress("t_eHcal_noPU", &t_eHcal_noPU, &b_t_eHcal_noPU);
1533 fChain->SetBranchAddress("t_delta_NoPU", &t_delta_NoPU, &b_t_delta_NoPU);
1534 fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
1535 Notify();
1536 }
1537
1538 Bool_t CalibFitPU::Notify() {
1539
1540
1541
1542
1543
1544
1545 return kTRUE;
1546 }
1547
1548 void CalibFitPU::Show(Long64_t entry) {
1549
1550
1551 if (!fChain)
1552 return;
1553 fChain->Show(entry);
1554 }
1555
1556 Int_t CalibFitPU::Cut(Long64_t) {
1557
1558
1559
1560 return 1;
1561 }
1562
1563 void CalibFitPU::Loop(bool extract_PU_parameters, std::string fileName) {
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587 if (fChain == 0)
1588 return;
1589
1590 Long64_t nentries = fChain->GetEntriesFast();
1591 Long64_t nbytes = 0, nb = 0;
1592
1593 char filename[100];
1594
1595 gStyle->SetCanvasBorderMode(0);
1596 gStyle->SetCanvasColor(kWhite);
1597 gStyle->SetPadColor(kWhite);
1598 gStyle->SetFillColor(kWhite);
1599
1600 const int n = 7;
1601 int ieta_grid[n] = {7, 16, 25, 26, 27, 28, 29};
1602 double a00[n], a10[n], a20[n], a01[n], a11[n], a21[n], a02[n], a12[n], a22[n];
1603 const int nbin1 = 15;
1604 double bins1[nbin1 + 1] = {0.0, 0.1, 0.2, 0.3, 0.4, 0.6, 0.8, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 6.0, 9.0};
1605 const int nbin2 = 7;
1606 double bins2[nbin1 + 1] = {1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0};
1607 int nbins[n] = {4, 5, 12, nbin2, nbin2, nbin2, nbin2};
1608 char name[20], title[100];
1609
1610 std::vector<TH2D *> vec_h2;
1611 std::vector<TGraph *> vec_gr;
1612 std::vector<TProfile *> vec_hp;
1613
1614 if (extract_PU_parameters) {
1615 for (int k = 0; k < n; k++) {
1616 sprintf(name, "h2_ieta%d", k);
1617 int ieta1 = (k == 0) ? 1 : ieta_grid[k - 1];
1618 sprintf(title, "PU Energy vs #Delta/p (i#eta = %d:%d)", ieta1, (ieta_grid[k] - 1));
1619 vec_h2.push_back(new TH2D(name, title, 100, 0, 10, 100, 0, 2));
1620 vec_gr.push_back(new TGraph());
1621 sprintf(name, "hp_ieta%d", k);
1622 if (k < 3)
1623 vec_hp.push_back(new TProfile(name, title, nbins[k], bins1));
1624 else
1625 vec_hp.push_back(new TProfile(name, title, nbins[k], bins2));
1626 }
1627
1628 int points[7] = {0, 0, 0, 0, 0, 0, 0};
1629
1630
1631 for (Long64_t jentry = 0; jentry < nentries; jentry++) {
1632 Long64_t ientry = LoadTree(jentry);
1633 if (ientry < 0)
1634 break;
1635 nb = fChain->GetEntry(jentry);
1636 nbytes += nb;
1637
1638 double deltaOvP = t_delta_PU / t_p_PU;
1639 double diffEpuEnopuOvP = t_eHcal_noPU / t_eHcal_PU;
1640
1641 for (int k = 0; k < n; k++) {
1642 if (std::abs(t_ieta) < ieta_grid[k]) {
1643 points[k]++;
1644 vec_h2[k]->Fill(deltaOvP, diffEpuEnopuOvP);
1645 vec_gr[k]->SetPoint(points[k], deltaOvP, diffEpuEnopuOvP);
1646 vec_hp[k]->Fill(deltaOvP, diffEpuEnopuOvP);
1647 break;
1648 }
1649 }
1650
1651 }
1652
1653 for (int k = 0; k < n; k++)
1654 std::cout << "Bin " << k << " for 2D Hist " << vec_h2[k] << ":" << vec_h2[k]->GetEntries() << " Graph "
1655 << vec_gr[k] << ":" << points[k] << " Profile " << vec_hp[k] << ":" << vec_hp[k]->GetEntries()
1656 << std::endl;
1657
1658 std::ofstream myfile0, myfile1, myfile2;
1659 sprintf(filename, "%s_par2d.txt", fileName.c_str());
1660 myfile0.open(filename);
1661 sprintf(filename, "%s_parProf.txt", fileName.c_str());
1662 myfile1.open(filename);
1663 sprintf(filename, "%s_parGraph.txt", fileName.c_str());
1664 myfile2.open(filename);
1665
1666 char namepng[20];
1667 for (int k = 0; k < n; k++) {
1668 gStyle->SetOptStat(1100);
1669 gStyle->SetOptFit(1);
1670
1671 TF1 *f1 = ((k < 2) ? (new TF1("f1", "[0]+[1]*x", 0, 5)) : (new TF1("f1", "[0]+[1]*x+[2]*x*x", 0, 5)));
1672 sprintf(name, "c_ieta%d2D", k);
1673 TCanvas *pad1 = new TCanvas(name, name, 500, 500);
1674 pad1->SetLeftMargin(0.10);
1675 pad1->SetRightMargin(0.10);
1676 pad1->SetTopMargin(0.10);
1677 vec_h2[k]->GetXaxis()->SetLabelSize(0.035);
1678 vec_h2[k]->GetYaxis()->SetLabelSize(0.035);
1679 vec_h2[k]->GetXaxis()->SetTitle("#Delta/p");
1680 vec_h2[k]->GetYaxis()->SetTitle("E_{NoPU} / E_{PU}");
1681 vec_h2[k]->Draw();
1682 vec_h2[k]->Fit(f1);
1683
1684 a00[k] = f1->GetParameter(0);
1685 a10[k] = f1->GetParameter(1);
1686 a20[k] = (k < 2) ? 0 : f1->GetParameter(2);
1687 myfile0 << k << "\t" << a00[k] << "\t" << a10[k] << "\t" << a20[k] << "\n";
1688 pad1->Update();
1689 TPaveStats *st1 = (TPaveStats *)vec_h2[k]->GetListOfFunctions()->FindObject("stats");
1690 if (st1 != nullptr) {
1691 st1->SetY1NDC(0.70);
1692 st1->SetY2NDC(0.90);
1693 st1->SetX1NDC(0.65);
1694 st1->SetX2NDC(0.90);
1695 }
1696 pad1->Update();
1697 sprintf(namepng, "%s.png", pad1->GetName());
1698 pad1->Print(namepng);
1699
1700 TF1 *f2 = ((k < 2) ? (new TF1("f2", "[0]+[1]*x", 0, 5)) : (new TF1("f2", "[0]+[1]*x+[2]*x*x", 0, 5)));
1701 sprintf(name, "c_ieta%dPr", k);
1702 TCanvas *pad2 = new TCanvas(name, name, 500, 500);
1703 pad2->SetLeftMargin(0.10);
1704 pad2->SetRightMargin(0.10);
1705 pad2->SetTopMargin(0.10);
1706 vec_hp[k]->GetXaxis()->SetLabelSize(0.035);
1707 vec_hp[k]->GetYaxis()->SetLabelSize(0.035);
1708 vec_hp[k]->GetXaxis()->SetTitle("#Delta/p");
1709 vec_hp[k]->GetYaxis()->SetTitle("E_{NoPU} / E_{PU}");
1710 vec_hp[k]->Draw();
1711 vec_hp[k]->Fit(f2);
1712
1713 a01[k] = f2->GetParameter(0);
1714 a11[k] = f2->GetParameter(1);
1715 a21[k] = (k < 2) ? 0 : f2->GetParameter(2);
1716 myfile1 << k << "\t" << a01[k] << "\t" << a11[k] << "\t" << a21[k] << "\n";
1717 pad2->Update();
1718 TPaveStats *st2 = (TPaveStats *)vec_hp[k]->GetListOfFunctions()->FindObject("stats");
1719 if (st2 != nullptr) {
1720 st2->SetY1NDC(0.70);
1721 st2->SetY2NDC(0.90);
1722 st2->SetX1NDC(0.65);
1723 st2->SetX2NDC(0.90);
1724 }
1725 pad2->Update();
1726 sprintf(namepng, "%s.png", pad2->GetName());
1727 pad2->Print(namepng);
1728
1729 TF1 *f3 = ((k < 2) ? (new TF1("f3", "[0]+[1]*x", 0, 5)) : (new TF1("f3", "[0]+[1]*x+[2]*x*x", 0, 5)));
1730 sprintf(name, "c_ieta%dGr", k);
1731 TCanvas *pad3 = new TCanvas(name, name, 500, 500);
1732 pad3->SetLeftMargin(0.10);
1733 pad3->SetRightMargin(0.10);
1734 pad3->SetTopMargin(0.10);
1735 gStyle->SetOptFit(1111);
1736 vec_gr[k]->GetXaxis()->SetLabelSize(0.035);
1737 vec_gr[k]->GetYaxis()->SetLabelSize(0.035);
1738 vec_gr[k]->GetXaxis()->SetTitle("#Delta/p");
1739 vec_gr[k]->GetYaxis()->SetTitle("E_{PU} - E_{NoPU}/p");
1740 vec_gr[k]->Fit(f3, "R");
1741 vec_gr[k]->Draw("Ap");
1742 f3->Draw("same");
1743
1744 a02[k] = f3->GetParameter(0);
1745 a12[k] = f3->GetParameter(1);
1746 a22[k] = (k < 2) ? 0 : f3->GetParameter(2);
1747 myfile2 << k << "\t" << a02[k] << "\t" << a12[k] << "\t" << a22[k] << "\n";
1748 pad3->Update();
1749 TPaveStats *st3 = (TPaveStats *)vec_gr[k]->GetListOfFunctions()->FindObject("stats");
1750 if (st3 != nullptr) {
1751 st3->SetY1NDC(0.70);
1752 st3->SetY2NDC(0.90);
1753 st3->SetX1NDC(0.65);
1754 st3->SetX2NDC(0.90);
1755 }
1756 pad3->Update();
1757 pad3->Modified();
1758 sprintf(namepng, "%s.png", pad3->GetName());
1759 pad3->Print(namepng);
1760 }
1761 } else {
1762 std::string line;
1763 double number;
1764 sprintf(filename, "%s_par2d.txt", fileName.c_str());
1765 std::ifstream myfile0(filename);
1766 if (myfile0.is_open()) {
1767 int iii = 0;
1768 while (getline(myfile0, line)) {
1769 std::istringstream iss2(line);
1770 int ii = 0;
1771 while (iss2 >> number) {
1772 if (ii == 0)
1773 a00[iii] = number;
1774 else if (ii == 1)
1775 a10[iii] = number;
1776 else if (ii == 2)
1777 a20[iii] = number;
1778 ++ii;
1779 }
1780 ++iii;
1781 }
1782 }
1783 sprintf(filename, "%s_parProf.txt", fileName.c_str());
1784 std::ifstream myfile1(filename);
1785 if (myfile1.is_open()) {
1786 int iii = 0;
1787 while (getline(myfile1, line)) {
1788 std::istringstream iss2(line);
1789 int ii = 0;
1790 while (iss2 >> number) {
1791 if (ii == 0)
1792 a01[iii] = number;
1793 else if (ii == 1)
1794 a11[iii] = number;
1795 else if (ii == 2)
1796 a21[iii] = number;
1797 ++ii;
1798 }
1799 ++iii;
1800 }
1801 }
1802 sprintf(filename, "%s_parGraph.txt", fileName.c_str());
1803 std::ifstream myfile2(filename);
1804 if (myfile2.is_open()) {
1805 int iii = 0;
1806 while (getline(myfile2, line)) {
1807 std::istringstream iss2(line);
1808 int ii = 0;
1809 while (iss2 >> number) {
1810 if (ii == 0)
1811 a02[iii] = number;
1812 else if (ii == 1)
1813 a12[iii] = number;
1814 else if (ii == 2)
1815 a22[iii] = number;
1816 ++ii;
1817 }
1818 ++iii;
1819 }
1820 }
1821 }
1822
1823 std::cout << "\nParameter Values:\n";
1824 for (int i = 0; i < n; i++) {
1825 std::cout << "[" << i << "] (" << a00[i] << ", " << a10[i] << ", " << a20[i] << ") (" << a01[i] << ", " << a11[i]
1826 << ", " << a21[i] << ") (" << a02[i] << ", " << a12[i] << ", " << a22[i] << ")\n";
1827 }
1828 std::cout << "\n\n";
1829 std::vector<TH1F *> vec_EovP_UnCorr_PU, vec_EovP_UnCorr_NoPU;
1830 std::vector<TH1F *> vec_EovP_Corr0_PU, vec_EovP_Corr0_NoPU;
1831 std::vector<TH1F *> vec_EovP_Corr1_PU, vec_EovP_Corr1_NoPU;
1832 std::vector<TH1F *> vec_EovP_Corr2_PU, vec_EovP_Corr2_NoPU;
1833 TH1F *h_current;
1834
1835 for (int k = 0; k < (2 * 28 + 1); k++) {
1836 if (k != 28) {
1837 sprintf(name, "EovP_ieta%dUnPU", k - 28);
1838 sprintf(title, "E/p (Uncorrected PU) for i#eta = %d", k - 28);
1839 h_current = new TH1F(name, title, 100, 0, 5);
1840 h_current->GetXaxis()->SetTitle(title);
1841 h_current->GetYaxis()->SetTitle("Tracks");
1842 vec_EovP_UnCorr_PU.push_back(h_current);
1843 sprintf(name, "EovP_ieta%dUnNoPU", k - 28);
1844 sprintf(title, "E/p (Uncorrected No PU) for i#eta = %d", k - 28);
1845 h_current = new TH1F(name, title, 100, 0, 5);
1846 h_current->GetXaxis()->SetTitle(title);
1847 h_current->GetYaxis()->SetTitle("Tracks");
1848 vec_EovP_UnCorr_NoPU.push_back(h_current);
1849 sprintf(name, "EovP_ieta%dCor0NoPU", k - 28);
1850 sprintf(title, "E/p (Corrected using 2D No PU) for i#eta = %d", k - 28);
1851 h_current = new TH1F(name, title, 100, 0, 5);
1852 h_current->GetXaxis()->SetTitle(title);
1853 h_current->GetYaxis()->SetTitle("Tracks");
1854 vec_EovP_Corr0_NoPU.push_back(h_current);
1855 sprintf(name, "EovP_ieta%dCor0PU", k - 28);
1856 sprintf(title, "E/p (Corrected using 2D PU) for i#eta = %d", k - 28);
1857 h_current = new TH1F(name, title, 100, 0, 5);
1858 h_current->GetXaxis()->SetTitle(title);
1859 h_current->GetYaxis()->SetTitle("Tracks");
1860 vec_EovP_Corr0_PU.push_back(h_current);
1861 sprintf(name, "EovP_ieta%dCor1NoPU", k - 28);
1862 sprintf(title, "E/p (Corrected using profile No PU) for i#eta = %d", k - 28);
1863 h_current = new TH1F(name, title, 100, 0, 5);
1864 h_current->GetXaxis()->SetTitle(title);
1865 h_current->GetYaxis()->SetTitle("Tracks");
1866 vec_EovP_Corr1_NoPU.push_back(h_current);
1867 sprintf(name, "EovP_ieta%dCor1PU", k - 28);
1868 sprintf(title, "E/p (Corrected using profile PU) for i#eta = %d", k - 28);
1869 h_current = new TH1F(name, title, 100, 0, 5);
1870 h_current->GetXaxis()->SetTitle(title);
1871 h_current->GetYaxis()->SetTitle("Tracks");
1872 vec_EovP_Corr1_PU.push_back(h_current);
1873 sprintf(name, "EovP_ieta%dCor2NoPU", k - 28);
1874 sprintf(title, "E/p (Corrected using graph No PU) for i#eta = %d", k - 28);
1875 h_current = new TH1F(name, title, 100, 0, 5);
1876 h_current->GetXaxis()->SetTitle(title);
1877 h_current->GetYaxis()->SetTitle("Tracks");
1878 vec_EovP_Corr2_NoPU.push_back(h_current);
1879 sprintf(name, "EovP_ieta%dCor2PU", k - 28);
1880 sprintf(title, "E/p (Corrected using graph PU) for i#eta = %d", k - 28);
1881 h_current = new TH1F(name, title, 100, 0, 5);
1882 h_current->GetXaxis()->SetTitle(title);
1883 h_current->GetYaxis()->SetTitle("Tracks");
1884 vec_EovP_Corr2_PU.push_back(h_current);
1885 }
1886 }
1887 std::cout << "Book " << (8 * vec_EovP_UnCorr_PU.size()) << " histograms\n";
1888
1889
1890
1891 nbytes = nb = 0;
1892 std::cout << nentries << " entries in the root file" << std::endl;
1893 for (Long64_t jentry = 0; jentry < nentries; jentry++) {
1894 Long64_t ientry = LoadTree(jentry);
1895 if (ientry < 0)
1896 break;
1897 nb = fChain->GetEntry(jentry);
1898 nbytes += nb;
1899
1900 int i1 = n;
1901 for (int k = 0; k < n; k++) {
1902 if (std::abs(t_ieta) < ieta_grid[k]) {
1903 i1 = k;
1904 break;
1905 }
1906 }
1907 if ((t_ieta == 0) || (i1 >= n))
1908 continue;
1909 int i2 = (t_ieta < 0) ? (t_ieta + 28) : (t_ieta + 27);
1910
1911 double EpuOvP = t_eHcal_PU / t_p_PU;
1912 double EnopuOvP = t_eHcal_noPU / t_p_PU;
1913 double deltaOvP = t_delta_PU / t_p_PU;
1914 double deltaNopuOvP = t_delta_NoPU / t_p_PU;
1915
1916 vec_EovP_UnCorr_PU[i2]->Fill(EpuOvP);
1917 vec_EovP_UnCorr_NoPU[i2]->Fill(EnopuOvP);
1918
1919 double c0p =
1920 (((i1 < 3) || (deltaOvP > 1.0)) ? (a00[i1] + a10[i1] * deltaOvP + a20[i1] * deltaOvP * deltaOvP) : 1.0);
1921 double c1p =
1922 (((i1 < 3) || (deltaOvP > 1.0)) ? (a01[i1] + a11[i1] * deltaOvP + a21[i1] * deltaOvP * deltaOvP) : 1.0);
1923 double c2p =
1924 (((i1 < 3) || (deltaOvP > 1.0)) ? (a02[i1] + a12[i1] * deltaOvP + a22[i1] * deltaOvP * deltaOvP) : 1.0);
1925 double c0np =
1926 (((i1 < 3) || (deltaNopuOvP > 1.0)) ? (a00[i1] + a10[i1] * deltaNopuOvP + a12[i1] * deltaNopuOvP * deltaNopuOvP)
1927 : 1.0);
1928 double c1np =
1929 (((i1 < 3) || (deltaNopuOvP > 1.0)) ? (a01[i1] + a11[i1] * deltaNopuOvP + a21[i1] * deltaNopuOvP * deltaNopuOvP)
1930 : 1.0);
1931 double c2np =
1932 (((i1 < 3) || (deltaNopuOvP > 1.0)) ? (a02[i1] + a12[i1] * deltaNopuOvP + a22[i1] * deltaNopuOvP * deltaNopuOvP)
1933 : 1.0);
1934
1935 vec_EovP_Corr0_PU[i2]->Fill(EpuOvP * c0p);
1936 vec_EovP_Corr0_NoPU[i2]->Fill(EnopuOvP * c0np);
1937 vec_EovP_Corr1_PU[i2]->Fill(EpuOvP * c1p);
1938 vec_EovP_Corr1_NoPU[i2]->Fill(EnopuOvP * c1np);
1939 vec_EovP_Corr2_PU[i2]->Fill(EpuOvP * c2p);
1940 vec_EovP_Corr2_NoPU[i2]->Fill(EnopuOvP * c2np);
1941 }
1942
1943 sprintf(filename, "%s.root", fileName.c_str());
1944 TFile *f1 = new TFile(filename, "RECREATE");
1945 f1->cd();
1946 for (unsigned int k = 0; k < vec_EovP_UnCorr_PU.size(); k++) {
1947 vec_EovP_UnCorr_PU[k]->Write();
1948 vec_EovP_UnCorr_NoPU[k]->Write();
1949 vec_EovP_Corr0_PU[k]->Write();
1950 vec_EovP_Corr0_NoPU[k]->Write();
1951 vec_EovP_Corr1_PU[k]->Write();
1952 vec_EovP_Corr1_NoPU[k]->Write();
1953 vec_EovP_Corr2_PU[k]->Write();
1954 vec_EovP_Corr2_NoPU[k]->Write();
1955 }
1956 for (unsigned int k = 0; k < vec_hp.size(); ++k) {
1957 vec_h2[k]->Write();
1958 vec_hp[k]->Write();
1959 vec_gr[k]->Write();
1960 }
1961 }
1962
1963 class CalibMerge {
1964 public:
1965 struct isotk {
1966 isotk(int pv = 0, double rho = 0, double pp = 0, double eH = 0, double eH10 = 0, double eH30 = 0, double del = 0)
1967 : goodPV(pv), rhoh(rho), p(pp), eHcal(eH), eHcal10(eH10), eHcal30(eH30), delta(del) {}
1968
1969 int goodPV;
1970 double rhoh, p, eHcal, eHcal10, eHcal30, delta;
1971 };
1972
1973 CalibMerge();
1974 ~CalibMerge();
1975 Int_t Cut(Long64_t entry);
1976 Int_t GetEntry(Long64_t entry);
1977 Long64_t LoadTree(Long64_t entry);
1978 Bool_t Notify();
1979 void Show(Long64_t entry = -1);
1980 void LoopNoPU(const char *fnameNoPU, std::string dirnm = "HcalIsoTrkAnalyzer");
1981 void LoopPU(const char *fnamePU, const char *fnameOut, std::string dirnm = "HcalIsoTrkAnalyzer");
1982 bool Init(const char *fname, const std::string &dirnm);
1983 void bookTree(const char *fname);
1984 void close();
1985
1986 private:
1987 TChain *fChain;
1988 Int_t fCurrent;
1989
1990 std::map<std::pair<int, int>, std::vector<isotk> > isotks_;
1991
1992
1993 Int_t t_Run;
1994 Int_t t_Event;
1995 Int_t t_DataType;
1996 Int_t t_ieta;
1997 Int_t t_iphi;
1998 Double_t t_EventWeight;
1999 Int_t t_nVtx;
2000 Int_t t_nTrk;
2001 Int_t t_goodPV;
2002 Double_t t_l1pt;
2003 Double_t t_l1eta;
2004 Double_t t_l1phi;
2005 Double_t t_l3pt;
2006 Double_t t_l3eta;
2007 Double_t t_l3phi;
2008 Double_t t_p;
2009 Double_t t_pt;
2010 Double_t t_phi;
2011 Double_t t_mindR1;
2012 Double_t t_mindR2;
2013 Double_t t_eMipDR;
2014 Double_t t_eHcal;
2015 Double_t t_eHcal10;
2016 Double_t t_eHcal30;
2017 Double_t t_hmaxNearP;
2018 Double_t t_rhoh;
2019 Bool_t t_selectTk;
2020 Bool_t t_qltyFlag;
2021 Bool_t t_qltyMissFlag;
2022 Bool_t t_qltyPVFlag;
2023 Double_t t_gentrackP;
2024 std::vector<unsigned int> *t_DetIds;
2025 std::vector<double> *t_HitEnergies;
2026 std::vector<bool> *t_trgbits;
2027 std::vector<unsigned int> *t_DetIds1;
2028 std::vector<unsigned int> *t_DetIds3;
2029 std::vector<double> *t_HitEnergies1;
2030 std::vector<double> *t_HitEnergies3;
2031
2032
2033 TBranch *b_t_Run;
2034 TBranch *b_t_Event;
2035 TBranch *b_t_DataType;
2036 TBranch *b_t_ieta;
2037 TBranch *b_t_iphi;
2038 TBranch *b_t_EventWeight;
2039 TBranch *b_t_nVtx;
2040 TBranch *b_t_nTrk;
2041 TBranch *b_t_goodPV;
2042 TBranch *b_t_l1pt;
2043 TBranch *b_t_l1eta;
2044 TBranch *b_t_l1phi;
2045 TBranch *b_t_l3pt;
2046 TBranch *b_t_l3eta;
2047 TBranch *b_t_l3phi;
2048 TBranch *b_t_p;
2049 TBranch *b_t_pt;
2050 TBranch *b_t_phi;
2051 TBranch *b_t_mindR1;
2052 TBranch *b_t_mindR2;
2053 TBranch *b_t_eMipDR;
2054 TBranch *b_t_eHcal;
2055 TBranch *b_t_eHcal10;
2056 TBranch *b_t_eHcal30;
2057 TBranch *b_t_hmaxNearP;
2058 TBranch *b_t_rhoh;
2059 TBranch *b_t_selectTk;
2060 TBranch *b_t_qltyFlag;
2061 TBranch *b_t_qltyMissFlag;
2062 TBranch *b_t_qltyPVFlag;
2063 TBranch *b_t_gentrackP;
2064 TBranch *b_t_DetIds;
2065 TBranch *b_t_HitEnergies;
2066 TBranch *b_t_trgbits;
2067 TBranch *b_t_DetIds1;
2068 TBranch *b_t_DetIds3;
2069 TBranch *b_t_HitEnergies1;
2070 TBranch *b_t_HitEnergies3;
2071
2072 TFile *outputFile_;
2073 TTree *outtree_;
2074 Int_t event_;
2075 Int_t ieta_;
2076 Int_t nvtx_;
2077 Double_t rhoh_;
2078 Double_t p_NoPU_;
2079 Double_t p_PU_;
2080 Double_t eHcal_NoPU_;
2081 Double_t eHcal_PU_;
2082 Double_t eHcal10_NoPU_;
2083 Double_t eHcal10_PU_;
2084 Double_t eHcal30_NoPU_;
2085 Double_t eHcal30_PU_;
2086 Double_t delta_NoPU_;
2087 Double_t delta_PU_;
2088 };
2089
2090 CalibMerge::CalibMerge() : outputFile_(nullptr), outtree_(nullptr) {}
2091
2092 CalibMerge::~CalibMerge() { close(); }
2093
2094 Int_t CalibMerge::GetEntry(Long64_t entry) {
2095
2096 if (!fChain)
2097 return 0;
2098 return fChain->GetEntry(entry);
2099 }
2100
2101 Long64_t CalibMerge::LoadTree(Long64_t entry) {
2102
2103 if (!fChain)
2104 return -5;
2105 Long64_t centry = fChain->LoadTree(entry);
2106 if (centry < 0)
2107 return centry;
2108 if (fChain->GetTreeNumber() != fCurrent) {
2109 fCurrent = fChain->GetTreeNumber();
2110 Notify();
2111 }
2112 return centry;
2113 }
2114
2115 Bool_t CalibMerge::Notify() {
2116
2117
2118
2119
2120
2121
2122 return kTRUE;
2123 }
2124
2125 void CalibMerge::Show(Long64_t entry) {
2126
2127
2128 if (!fChain)
2129 return;
2130 fChain->Show(entry);
2131 }
2132
2133 Int_t CalibMerge::Cut(Long64_t) {
2134
2135
2136
2137 return 1;
2138 }
2139
2140 void CalibMerge::LoopNoPU(const char *fname, std::string dirnm) {
2141 bool ok = Init(fname, dirnm);
2142 std::cout << "Opens no PU file " << fname << " with flag " << ok << std::endl;
2143 isotks_.clear();
2144 if (ok) {
2145 Long64_t nentries = fChain->GetEntriesFast();
2146 Long64_t nbytes = 0, nb = 0;
2147 for (Long64_t jentry = 0; jentry < nentries; jentry++) {
2148 Long64_t ientry = LoadTree(jentry);
2149 if (ientry < 0)
2150 break;
2151 nb = fChain->GetEntry(jentry);
2152 nbytes += nb;
2153 std::pair<int, int> key(t_Event, t_ieta);
2154 isotk tk(t_goodPV, t_rhoh, t_p, t_eHcal, t_eHcal10, t_eHcal30, (t_eHcal30 - t_eHcal10));
2155 std::map<std::pair<int, int>, std::vector<isotk> >::iterator itr = isotks_.find(key);
2156 if (itr == isotks_.end()) {
2157 std::vector<isotk> v;
2158 v.push_back(tk);
2159 isotks_[key] = v;
2160 } else {
2161 (itr->second).push_back(tk);
2162 }
2163 }
2164 std::cout << "Finds " << isotks_.size() << " tracks from " << nentries << " entries with " << nbytes << " bytes"
2165 << std::endl;
2166 }
2167 }
2168
2169 void CalibMerge::LoopPU(const char *fname, const char *fout, std::string dirnm) {
2170 bool ok = Init(fname, dirnm);
2171 std::cout << "Opens PU file " << fname << " with flag " << ok << std::endl;
2172 if (ok && (isotks_.size() > 0)) {
2173 bookTree(fout);
2174 Long64_t nentries = fChain->GetEntriesFast();
2175 Long64_t nbytes(0), nb(0), merge(0);
2176 for (Long64_t jentry = 0; jentry < nentries; jentry++) {
2177 Long64_t ientry = LoadTree(jentry);
2178 if (ientry < 0)
2179 break;
2180 nb = fChain->GetEntry(jentry);
2181 nbytes += nb;
2182 std::pair<int, int> key(t_Event, t_ieta);
2183 std::map<std::pair<int, int>, std::vector<isotk> >::iterator itr = isotks_.find(key);
2184 if (itr != isotks_.end()) {
2185 for (unsigned int k = 0; (itr->second).size(); ++k) {
2186 if (std::abs(t_p - (itr->second)[k].p) < 0.01) {
2187 event_ = (itr->first).first;
2188 ieta_ = (itr->first).second;
2189 nvtx_ = (itr->second)[k].goodPV;
2190 rhoh_ = (itr->second)[k].rhoh;
2191 p_NoPU_ = (itr->second)[k].p;
2192 p_PU_ = t_p;
2193 eHcal_NoPU_ = (itr->second)[k].eHcal;
2194 eHcal_PU_ = t_eHcal;
2195 eHcal10_NoPU_ = (itr->second)[k].eHcal10;
2196 eHcal10_PU_ = t_eHcal10;
2197 eHcal30_NoPU_ = (itr->second)[k].eHcal30;
2198 eHcal30_PU_ = t_eHcal30;
2199 delta_NoPU_ = (itr->second)[k].delta;
2200 delta_PU_ = (t_eHcal30 - t_eHcal10);
2201 outtree_->Fill();
2202 ++merge;
2203 break;
2204 }
2205 }
2206 }
2207 }
2208 std::cout << "Use " << isotks_.size() << ":" << nentries << " from noPU|PU files to produce " << merge
2209 << " merged tracks";
2210 }
2211 }
2212
2213 bool CalibMerge::Init(const char *fname, const std::string &dirnm) {
2214 char treeName[400];
2215 sprintf(treeName, "%s/CalibTree", dirnm.c_str());
2216 TChain *fChain = new TChain(treeName);
2217 if (!fillChain(fChain, fname)) {
2218 std::cout << "Creation of the chain for " << treeName << " from " << fname << " failed " << std::endl;
2219 return false;
2220 } else {
2221 std::cout << "Proceed with a tree chain with " << fChain->GetEntries() << " entries" << std::endl;
2222 fCurrent = -1;
2223 fChain->SetMakeClass(1);
2224
2225 t_DetIds = 0;
2226 t_DetIds1 = 0;
2227 t_DetIds3 = 0;
2228 t_HitEnergies = 0;
2229 t_HitEnergies1 = 0;
2230 t_HitEnergies3 = 0;
2231 t_trgbits = 0;
2232 fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
2233 fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
2234 fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
2235 fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
2236 fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
2237 fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
2238 fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
2239 fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
2240 fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
2241 fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
2242 fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
2243 fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
2244 fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
2245 fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
2246 fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
2247 fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
2248 fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
2249 fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
2250 fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
2251 fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
2252 fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
2253 fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
2254 fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
2255 fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
2256 fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
2257 fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
2258 fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
2259 fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
2260 fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
2261 fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
2262 fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
2263 fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
2264 fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
2265 fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
2266 fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
2267 fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
2268 fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
2269 fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
2270 Notify();
2271 return true;
2272 }
2273 }
2274
2275 void CalibMerge::bookTree(const char *fname) {
2276 std::cout << "BookTree in " << fname << std::endl;
2277 outputFile_ = TFile::Open(fname, "RECREATE");
2278 outputFile_->cd();
2279 outtree_ = new TTree("Tree", "Tree");
2280 outtree_->Branch("t_Event", &event_);
2281 outtree_->Branch("t_ieta", &ieta_);
2282 outtree_->Branch("t_nvtx", &nvtx_);
2283 outtree_->Branch("t_rhoh", &rhoh_);
2284 outtree_->Branch("t_p_noPU", &p_NoPU_);
2285 outtree_->Branch("t_p_PU", &p_PU_);
2286 outtree_->Branch("t_eHcal_noPU", &eHcal_NoPU_);
2287 outtree_->Branch("t_eHcal_PU", &eHcal_PU_);
2288 outtree_->Branch("t_eHcal10_noPU", &eHcal10_NoPU_);
2289 outtree_->Branch("t_eHcal10_PU", &eHcal10_PU_);
2290 outtree_->Branch("t_eHcal30_noPU", &eHcal30_NoPU_);
2291 outtree_->Branch("t_eHcal30_PU", &eHcal30_PU_);
2292 outtree_->Branch("t_delta_noPU", &delta_NoPU_);
2293 outtree_->Branch("t_delta_PU", &delta_PU_);
2294 }
2295
2296 void CalibMerge::close() {
2297 if (outputFile_ != nullptr) {
2298 outputFile_->cd();
2299 std::cout << "file yet to be Written" << std::endl;
2300 outtree_->Write();
2301 std::cout << "file Written" << std::endl;
2302 outputFile_->Close();
2303 }
2304 outputFile_ = nullptr;
2305 std::cout << "now doing return" << std::endl;
2306 }