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