File indexing completed on 2023-03-17 10:42:55
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118 #include <TROOT.h>
0119 #include <TChain.h>
0120 #include <TFile.h>
0121 #include <TH1D.h>
0122 #include <TH2F.h>
0123 #include <TProfile.h>
0124 #include <TFitResult.h>
0125 #include <TFitResultPtr.h>
0126 #include <TStyle.h>
0127 #include <TCanvas.h>
0128 #include <TLegend.h>
0129 #include <TPaveStats.h>
0130 #include <TPaveText.h>
0131
0132 #include <algorithm>
0133 #include <iomanip>
0134 #include <iostream>
0135 #include <fstream>
0136 #include <sstream>
0137 #include <map>
0138 #include <vector>
0139 #include <string>
0140
0141 void unpackDetId(unsigned int, int &, int &, int &, int &, int &);
0142 #include "CalibCorr.C"
0143
0144 namespace CalibPlots {
0145 static const int npbin = 4;
0146 static const int netabin = 4;
0147 static const int ndepth = 7;
0148 static const int ntitles = 5;
0149 static const int npbin0 = (npbin + 1);
0150 int getP(int k) {
0151 const int ipbin[npbin0] = {20, 30, 40, 60, 100};
0152 return ((k >= 0 && k < npbin0) ? ipbin[k] : 0);
0153 }
0154 double getMomentum(int k) { return (double)(getP(k)); }
0155 int getEta(int k) {
0156 int ietas[netabin] = {0, 13, 18, 23};
0157 return ((k >= 0 && k < netabin) ? ietas[k] : -1);
0158 }
0159 std::string getTitle(int k) {
0160 std::string titl[ntitles] = {
0161 "all tracks", "good quality tracks", "selected tracks", "isolated good tracks", "tracks having MIP in ECAL"};
0162 return ((k >= 0 && k < ntitles) ? titl[k] : "");
0163 }
0164 }
0165
0166 class CalibPlotProperties {
0167 public:
0168 TChain *fChain;
0169 Int_t fCurrent;
0170
0171
0172 Int_t t_Run;
0173 Int_t t_Event;
0174 Int_t t_DataType;
0175 Int_t t_ieta;
0176 Int_t t_iphi;
0177 Double_t t_EventWeight;
0178 Int_t t_nVtx;
0179 Int_t t_nTrk;
0180 Int_t t_goodPV;
0181 Double_t t_l1pt;
0182 Double_t t_l1eta;
0183 Double_t t_l1phi;
0184 Double_t t_l3pt;
0185 Double_t t_l3eta;
0186 Double_t t_l3phi;
0187 Double_t t_p;
0188 Double_t t_pt;
0189 Double_t t_phi;
0190 Double_t t_mindR1;
0191 Double_t t_mindR2;
0192 Double_t t_eMipDR;
0193 Double_t t_eHcal;
0194 Double_t t_eHcal10;
0195 Double_t t_eHcal30;
0196 Double_t t_hmaxNearP;
0197 Double_t t_rhoh;
0198 Bool_t t_selectTk;
0199 Bool_t t_qltyFlag;
0200 Bool_t t_qltyMissFlag;
0201 Bool_t t_qltyPVFlag;
0202 Double_t t_gentrackP;
0203 std::vector<unsigned int> *t_DetIds;
0204 std::vector<double> *t_HitEnergies;
0205 std::vector<bool> *t_trgbits;
0206 std::vector<unsigned int> *t_DetIds1;
0207 std::vector<unsigned int> *t_DetIds3;
0208 std::vector<double> *t_HitEnergies1;
0209 std::vector<double> *t_HitEnergies3;
0210
0211
0212 TBranch *b_t_Run;
0213 TBranch *b_t_Event;
0214 TBranch *b_t_DataType;
0215 TBranch *b_t_ieta;
0216 TBranch *b_t_iphi;
0217 TBranch *b_t_EventWeight;
0218 TBranch *b_t_nVtx;
0219 TBranch *b_t_nTrk;
0220 TBranch *b_t_goodPV;
0221 TBranch *b_t_l1pt;
0222 TBranch *b_t_l1eta;
0223 TBranch *b_t_l1phi;
0224 TBranch *b_t_l3pt;
0225 TBranch *b_t_l3eta;
0226 TBranch *b_t_l3phi;
0227 TBranch *b_t_p;
0228 TBranch *b_t_pt;
0229 TBranch *b_t_phi;
0230 TBranch *b_t_mindR1;
0231 TBranch *b_t_mindR2;
0232 TBranch *b_t_eMipDR;
0233 TBranch *b_t_eHcal;
0234 TBranch *b_t_eHcal10;
0235 TBranch *b_t_eHcal30;
0236 TBranch *b_t_hmaxNearP;
0237 TBranch *b_t_rhoh;
0238 TBranch *b_t_selectTk;
0239 TBranch *b_t_qltyFlag;
0240 TBranch *b_t_qltyMissFlag;
0241 TBranch *b_t_qltyPVFlag;
0242 TBranch *b_t_gentrackP;
0243 TBranch *b_t_DetIds;
0244 TBranch *b_t_HitEnergies;
0245 TBranch *b_t_trgbits;
0246 TBranch *b_t_DetIds1;
0247 TBranch *b_t_DetIds3;
0248 TBranch *b_t_HitEnergies1;
0249 TBranch *b_t_HitEnergies3;
0250
0251 CalibPlotProperties(const char *fname,
0252 const std::string &dirname,
0253 const char *dupFileName,
0254 const std::string &prefix = "",
0255 const char *corrFileName = "",
0256 const char *rcorFileName = "",
0257 int puCorr = -8,
0258 int flag = 101111,
0259 bool dataMC = true,
0260 int truncateFlag = 0,
0261 bool useGen = false,
0262 double scale = 1.0,
0263 int useScale = 0,
0264 int etalo = 0,
0265 int etahi = 30,
0266 int runlo = 0,
0267 int runhi = 99999999,
0268 int phimin = 1,
0269 int phimax = 72,
0270 int zside = 1,
0271 int nvxlo = 0,
0272 int nvxhi = 1000,
0273 int rbx = 0,
0274 bool exclude = false,
0275 bool etamax = false);
0276 virtual ~CalibPlotProperties();
0277 virtual Int_t Cut(Long64_t entry);
0278 virtual Int_t GetEntry(Long64_t entry);
0279 virtual Long64_t LoadTree(Long64_t entry);
0280 virtual void Init(TChain *);
0281 virtual void Loop(Long64_t nentries = -1);
0282 virtual Bool_t Notify();
0283 virtual void Show(Long64_t entry = -1);
0284 bool goodTrack(double &eHcal, double &cut, bool debug);
0285 bool selectPhi(bool debug);
0286 void savePlot(const std::string &theName, bool append, bool all = false, bool debug = false);
0287 void correctEnergy(double &ener);
0288
0289 private:
0290 static const int kp50 = 2;
0291 CalibCorrFactor *corrFactor_;
0292 CalibCorr *cFactor_;
0293 CalibSelectRBX *cSelect_;
0294 CalibDuplicate *cDuplicate_;
0295 const std::string fname_, dirnm_, prefix_, outFileName_;
0296 const int corrPU_, flag_;
0297 const bool dataMC_, useGen_;
0298 const int truncateFlag_;
0299 const int etalo_, etahi_;
0300 int runlo_, runhi_;
0301 const int phimin_, phimax_, zside_, nvxlo_, nvxhi_, rbx_;
0302 bool exclude_, corrE_, cutL1T_;
0303 bool includeRun_, getHist_;
0304 int flexibleSelect_, ifDepth_, duplicate_;
0305 bool plotBasic_, plotEnergy_, plotHists_;
0306 double log2by18_;
0307 std::ofstream fileout_;
0308 std::vector<std::pair<int, int> > events_;
0309 TH1D *h_p[CalibPlots::ntitles];
0310 TH1D *h_eta[CalibPlots::ntitles], *h_nvtx, *h_nvtxEv, *h_nvtxTk;
0311 std::vector<TH1D *> h_eta0, h_eta1, h_eta2, h_eta3, h_eta4;
0312 std::vector<TH1D *> h_dL1, h_vtx;
0313 std::vector<TH1D *> h_etaEH[CalibPlots::npbin0];
0314 std::vector<TH1D *> h_etaEp[CalibPlots::npbin0];
0315 std::vector<TH1D *> h_etaEE[CalibPlots::npbin0];
0316 std::vector<TH1D *> h_etaEE0[CalibPlots::npbin0];
0317 std::vector<TH1D *> h_mom, h_eEcal, h_eHcal;
0318 std::vector<TH1F *> h_bvlist, h_bvlist2, h_evlist, h_evlist2;
0319 std::vector<TH1F *> h_bvlist3, h_evlist3;
0320 TH2F *h_etaE;
0321 };
0322
0323 CalibPlotProperties::CalibPlotProperties(const char *fname,
0324 const std::string &dirnm,
0325 const char *dupFileName,
0326 const std::string &prefix,
0327 const char *corrFileName,
0328 const char *rcorFileName,
0329 int puCorr,
0330 int flag,
0331 bool dataMC,
0332 int truncate,
0333 bool useGen,
0334 double scl,
0335 int useScale,
0336 int etalo,
0337 int etahi,
0338 int runlo,
0339 int runhi,
0340 int phimin,
0341 int phimax,
0342 int zside,
0343 int nvxlo,
0344 int nvxhi,
0345 int rbx,
0346 bool exc,
0347 bool etam)
0348 : corrFactor_(nullptr),
0349 cFactor_(nullptr),
0350 cSelect_(nullptr),
0351 cDuplicate_(nullptr),
0352 fname_(fname),
0353 dirnm_(dirnm),
0354 prefix_(prefix),
0355 corrPU_(puCorr),
0356 flag_(flag),
0357 dataMC_(dataMC),
0358 useGen_(useGen),
0359 truncateFlag_(truncate),
0360 etalo_(etalo),
0361 etahi_(etahi),
0362 runlo_(runlo),
0363 runhi_(runhi),
0364 phimin_(phimin),
0365 phimax_(phimax),
0366 zside_(zside),
0367 nvxlo_(nvxlo),
0368 nvxhi_(nvxhi),
0369 rbx_(rbx),
0370 exclude_(exc),
0371 includeRun_(true) {
0372
0373
0374
0375 flexibleSelect_ = ((flag_ / 1) % 10);
0376 plotBasic_ = (((flag_ / 10) % 10) > 0);
0377 plotEnergy_ = (((flag_ / 10) % 10) > 0);
0378 int oneplace = ((flag_ / 1000) % 10);
0379 cutL1T_ = (oneplace % 2);
0380 bool marina = ((oneplace / 2) % 2);
0381 ifDepth_ = ((flag_ / 10000) % 10);
0382 plotHists_ = (((flag_ / 100000) % 10) > 0);
0383 duplicate_ = ((flag_ / 1000000) % 10);
0384 log2by18_ = std::log(2.5) / 18.0;
0385 if (runlo_ < 0 || runhi_ < 0) {
0386 runlo_ = std::abs(runlo_);
0387 runhi_ = std::abs(runhi_);
0388 includeRun_ = false;
0389 }
0390 char treeName[400];
0391 sprintf(treeName, "%s/CalibTree", dirnm.c_str());
0392 TChain *chain = new TChain(treeName);
0393 std::cout << "Create a chain for " << treeName << " from " << fname << " flags " << flexibleSelect_ << "|"
0394 << plotBasic_ << "|"
0395 << "|" << plotEnergy_ << "|" << plotHists_ << "|" << corrPU_ << " cons " << log2by18_ << " eta range "
0396 << etalo_ << ":" << etahi_ << " run range " << runlo_ << ":" << runhi_ << " (inclusion flag " << includeRun_
0397 << ") Vertex Range " << nvxlo_ << ":" << nvxhi_ << std::endl;
0398 corrFactor_ = new CalibCorrFactor(corrFileName, useScale, scl, etam, marina, false);
0399 if (!fillChain(chain, fname)) {
0400 std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
0401 } else {
0402 std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
0403 Init(chain);
0404 if (std::string(rcorFileName) != "") {
0405 cFactor_ = new CalibCorr(rcorFileName, ifDepth_, false);
0406 if (cFactor_->absent())
0407 ifDepth_ = -1;
0408 } else {
0409 ifDepth_ = -1;
0410 }
0411 if (std::string(dupFileName) != "")
0412 cDuplicate_ = new CalibDuplicate(dupFileName, duplicate_, false);
0413 if (rbx != 0)
0414 cSelect_ = new CalibSelectRBX(rbx, false);
0415 }
0416 }
0417
0418 CalibPlotProperties::~CalibPlotProperties() {
0419 delete corrFactor_;
0420 delete cFactor_;
0421 delete cSelect_;
0422 if (!fChain)
0423 return;
0424 delete fChain->GetCurrentFile();
0425 }
0426
0427 Int_t CalibPlotProperties::GetEntry(Long64_t entry) {
0428
0429 if (!fChain)
0430 return 0;
0431 return fChain->GetEntry(entry);
0432 }
0433
0434 Long64_t CalibPlotProperties::LoadTree(Long64_t entry) {
0435
0436 if (!fChain)
0437 return -5;
0438 Long64_t centry = fChain->LoadTree(entry);
0439 if (centry < 0)
0440 return centry;
0441 if (!fChain->InheritsFrom(TChain::Class()))
0442 return centry;
0443 TChain *chain = (TChain *)fChain;
0444 if (chain->GetTreeNumber() != fCurrent) {
0445 fCurrent = chain->GetTreeNumber();
0446 Notify();
0447 }
0448 return centry;
0449 }
0450
0451 void CalibPlotProperties::Init(TChain *tree) {
0452
0453
0454
0455
0456
0457
0458
0459
0460
0461 t_DetIds = 0;
0462 t_DetIds1 = 0;
0463 t_DetIds3 = 0;
0464 t_HitEnergies = 0;
0465 t_HitEnergies1 = 0;
0466 t_HitEnergies3 = 0;
0467 t_trgbits = 0;
0468
0469 fChain = tree;
0470 fCurrent = -1;
0471 if (!tree)
0472 return;
0473 fChain->SetMakeClass(1);
0474
0475 fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
0476 fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
0477 fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
0478 fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
0479 fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
0480 fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
0481 fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
0482 fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
0483 fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
0484 fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
0485 fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
0486 fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
0487 fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
0488 fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
0489 fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
0490 fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
0491 fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
0492 fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
0493 fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
0494 fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
0495 fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
0496 fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
0497 fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
0498 fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
0499 fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
0500 fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
0501 fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
0502 fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
0503 fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
0504 fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
0505 fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
0506 fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
0507 fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
0508 fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
0509 fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
0510 fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
0511 fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
0512 fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
0513 Notify();
0514
0515 char name[20], title[200];
0516 unsigned int kk(0);
0517
0518 if (plotBasic_) {
0519 std::cout << "Book Basic Histos" << std::endl;
0520 h_nvtx = new TH1D("hnvtx", "Number of vertices (selected entries)", 10, 0, 100);
0521 h_nvtx->Sumw2();
0522 h_nvtxEv = new TH1D("hnvtxEv", "Number of vertices (selected events)", 10, 0, 100);
0523 h_nvtxEv->Sumw2();
0524 h_nvtxTk = new TH1D("hnvtxTk", "Number of vertices (selected tracks)", 10, 0, 100);
0525 h_nvtxTk->Sumw2();
0526 for (int k = 0; k < CalibPlots::ntitles; ++k) {
0527 sprintf(name, "%sp%d", prefix_.c_str(), k);
0528 sprintf(title, "Momentum for %s", CalibPlots::getTitle(k).c_str());
0529 h_p[k] = new TH1D(name, title, 150, 0.0, 150.0);
0530 sprintf(name, "%seta%d", prefix_.c_str(), k);
0531 sprintf(title, "#eta for %s", CalibPlots::getTitle(k).c_str());
0532 h_eta[k] = new TH1D(name, title, 60, -30.0, 30.0);
0533 }
0534 for (int k = 0; k < CalibPlots::npbin; ++k) {
0535 sprintf(name, "%seta0%d", prefix_.c_str(), k);
0536 sprintf(title,
0537 "#eta for %s (p = %d:%d GeV)",
0538 CalibPlots::getTitle(0).c_str(),
0539 CalibPlots::getP(k),
0540 CalibPlots::getP(k + 1));
0541 h_eta0.push_back(new TH1D(name, title, 60, -30.0, 30.0));
0542 kk = h_eta0.size() - 1;
0543 h_eta0[kk]->Sumw2();
0544 sprintf(name, "%seta1%d", prefix_.c_str(), k);
0545 sprintf(title,
0546 "#eta for %s (p = %d:%d GeV)",
0547 CalibPlots::getTitle(1).c_str(),
0548 CalibPlots::getP(k),
0549 CalibPlots::getP(k + 1));
0550 h_eta1.push_back(new TH1D(name, title, 60, -30.0, 30.0));
0551 kk = h_eta1.size() - 1;
0552 h_eta1[kk]->Sumw2();
0553 sprintf(name, "%seta2%d", prefix_.c_str(), k);
0554 sprintf(title,
0555 "#eta for %s (p = %d:%d GeV)",
0556 CalibPlots::getTitle(2).c_str(),
0557 CalibPlots::getP(k),
0558 CalibPlots::getP(k + 1));
0559 h_eta2.push_back(new TH1D(name, title, 60, -30.0, 30.0));
0560 kk = h_eta2.size() - 1;
0561 h_eta2[kk]->Sumw2();
0562 sprintf(name, "%seta3%d", prefix_.c_str(), k);
0563 sprintf(title,
0564 "#eta for %s (p = %d:%d GeV)",
0565 CalibPlots::getTitle(3).c_str(),
0566 CalibPlots::getP(k),
0567 CalibPlots::getP(k + 1));
0568 h_eta3.push_back(new TH1D(name, title, 60, -30.0, 30.0));
0569 kk = h_eta3.size() - 1;
0570 h_eta3[kk]->Sumw2();
0571 sprintf(name, "%seta4%d", prefix_.c_str(), k);
0572 sprintf(title,
0573 "#eta for %s (p = %d:%d GeV)",
0574 CalibPlots::getTitle(4).c_str(),
0575 CalibPlots::getP(k),
0576 CalibPlots::getP(k + 1));
0577 h_eta4.push_back(new TH1D(name, title, 60, -30.0, 30.0));
0578 kk = h_eta4.size() - 1;
0579 h_eta4[kk]->Sumw2();
0580 sprintf(name, "%sdl1%d", prefix_.c_str(), k);
0581 sprintf(title, "Distance from L1 (p = %d:%d GeV)", CalibPlots::getP(k), CalibPlots::getP(k + 1));
0582 h_dL1.push_back(new TH1D(name, title, 160, 0.0, 8.0));
0583 kk = h_dL1.size() - 1;
0584 h_dL1[kk]->Sumw2();
0585 sprintf(name, "%svtx%d", prefix_.c_str(), k);
0586 sprintf(title, "N_{Vertex} (p = %d:%d GeV)", CalibPlots::getP(k), CalibPlots::getP(k + 1));
0587 h_vtx.push_back(new TH1D(name, title, 100, 0.0, 100.0));
0588 kk = h_vtx.size() - 1;
0589 h_vtx[kk]->Sumw2();
0590 }
0591 }
0592
0593 if (plotEnergy_) {
0594 std::cout << "Make plots for good tracks" << std::endl;
0595 for (int k = 0; k < CalibPlots::npbin; ++k) {
0596 for (int j = etalo_; j <= etahi_ + 1; ++j) {
0597 sprintf(name, "%senergyH%d%d", prefix_.c_str(), k, j);
0598 if (j > etahi_)
0599 sprintf(title,
0600 "HCAL energy for %s (p = %d:%d GeV |#eta| = %d:%d)",
0601 CalibPlots::getTitle(3).c_str(),
0602 CalibPlots::getP(k),
0603 CalibPlots::getP(k + 1),
0604 etalo_,
0605 etahi_);
0606 else
0607 sprintf(title,
0608 "HCAL energy for %s (p = %d:%d GeV |#eta| = %d)",
0609 CalibPlots::getTitle(3).c_str(),
0610 CalibPlots::getP(k),
0611 CalibPlots::getP(k + 1),
0612 j);
0613 h_etaEH[k].push_back(new TH1D(name, title, 200, 0, 200));
0614 kk = h_etaEH[k].size() - 1;
0615 h_etaEH[k][kk]->Sumw2();
0616 sprintf(name, "%senergyP%d%d", prefix_.c_str(), k, j);
0617 if (j > etahi_)
0618 sprintf(title,
0619 "momentum for %s (p = %d:%d GeV |#eta| = %d:%d)",
0620 CalibPlots::getTitle(3).c_str(),
0621 CalibPlots::getP(k),
0622 CalibPlots::getP(k + 1),
0623 etalo_,
0624 etahi_);
0625 else
0626 sprintf(title,
0627 "momentum for %s (p = %d:%d GeV |#eta| = %d)",
0628 CalibPlots::getTitle(3).c_str(),
0629 CalibPlots::getP(k),
0630 CalibPlots::getP(k + 1),
0631 j);
0632 h_etaEp[k].push_back(new TH1D(name, title, 100, 0, 100));
0633 kk = h_etaEp[k].size() - 1;
0634 h_etaEp[k][kk]->Sumw2();
0635 sprintf(name, "%senergyE%d%d", prefix_.c_str(), k, j);
0636 if (j > etahi_)
0637 sprintf(title,
0638 "ECAL energy for %s (p = %d:%d GeV |#eta| = %d:%d)",
0639 CalibPlots::getTitle(3).c_str(),
0640 CalibPlots::getP(k),
0641 CalibPlots::getP(k + 1),
0642 etalo_,
0643 etahi_);
0644 else
0645 sprintf(title,
0646 "ECAL energy for %s (p = %d:%d GeV |#eta| = %d)",
0647 CalibPlots::getTitle(3).c_str(),
0648 CalibPlots::getP(k),
0649 CalibPlots::getP(k + 1),
0650 j);
0651 h_etaEE[k].push_back(new TH1D(name, title, 100, 0, 10));
0652 kk = h_etaEE[k].size() - 1;
0653 h_etaEE[k][kk]->Sumw2();
0654 sprintf(name, "%senergyER%d%d", prefix_.c_str(), k, j);
0655 h_etaEE0[k].push_back(new TH1D(name, title, 100, 0, 1));
0656 kk = h_etaEE0[k].size() - 1;
0657 h_etaEE0[k][kk]->Sumw2();
0658 }
0659 }
0660
0661 for (int j = 0; j < CalibPlots::netabin; ++j) {
0662 sprintf(name, "%senergyH%d", prefix_.c_str(), j);
0663 if (j == 0)
0664 sprintf(title, "HCAL energy for %s (All)", CalibPlots::getTitle(3).c_str());
0665 else
0666 sprintf(title,
0667 "HCAL energy for %s (|#eta| = %d:%d)",
0668 CalibPlots::getTitle(3).c_str(),
0669 CalibPlots::getEta(j - 1),
0670 CalibPlots::getEta(j));
0671 h_eHcal.push_back(new TH1D(name, title, 200, 0, 200));
0672 kk = h_eHcal.size() - 1;
0673 h_eHcal[kk]->Sumw2();
0674 sprintf(name, "%senergyP%d", prefix_.c_str(), j);
0675 if (j == 0)
0676 sprintf(title, "Track momentum for %s (All)", CalibPlots::getTitle(3).c_str());
0677 else
0678 sprintf(title,
0679 "Track momentum for %s (|#eta| = %d:%d)",
0680 CalibPlots::getTitle(3).c_str(),
0681 CalibPlots::getEta(j - 1),
0682 CalibPlots::getEta(j));
0683 h_mom.push_back(new TH1D(name, title, 100, 0, 100));
0684 kk = h_mom.size() - 1;
0685 h_mom[kk]->Sumw2();
0686 sprintf(name, "%senergyE%d", prefix_.c_str(), j);
0687 if (j == 0)
0688 sprintf(title, "ECAL energy for %s (All)", CalibPlots::getTitle(3).c_str());
0689 else
0690 sprintf(title,
0691 "ECAL energy for %s (|#eta| = %d:%d)",
0692 CalibPlots::getTitle(3).c_str(),
0693 CalibPlots::getEta(j - 1),
0694 CalibPlots::getEta(j));
0695 h_eEcal.push_back(new TH1D(name, title, 100, 0, 10));
0696 kk = h_eEcal.size() - 1;
0697 h_eEcal[kk]->Sumw2();
0698 }
0699 }
0700
0701 if (plotHists_) {
0702 for (int i = 0; i < CalibPlots::ndepth; i++) {
0703 sprintf(name, "b_edepth%d", i);
0704 sprintf(title, "Total RecHit energy in depth %d (Barrel)", i + 1);
0705 h_bvlist.push_back(new TH1F(name, title, 1000, 0, 100));
0706 h_bvlist[i]->Sumw2();
0707 sprintf(name, "b_recedepth%d", i);
0708 sprintf(title, "RecHit energy in depth %d (Barrel)", i + 1);
0709 h_bvlist2.push_back(new TH1F(name, title, 1000, 0, 100));
0710 h_bvlist2[i]->Sumw2();
0711 sprintf(name, "b_nrecdepth%d", i);
0712 sprintf(title, "#RecHits in depth %d (Barrel)", i + 1);
0713 h_bvlist3.push_back(new TH1F(name, title, 1000, 0, 100));
0714 h_bvlist3[i]->Sumw2();
0715 sprintf(name, "e_edepth%d", i);
0716 sprintf(title, "Total RecHit energy in depth %d (Endcap)", i + 1);
0717 h_evlist.push_back(new TH1F(name, title, 1000, 0, 100));
0718 h_evlist[i]->Sumw2();
0719 sprintf(name, "e_recedepth%d", i);
0720 sprintf(title, "RecHit energy in depth %d (Endcap)", i + 1);
0721 h_evlist2.push_back(new TH1F(name, title, 1000, 0, 100));
0722 h_evlist2[i]->Sumw2();
0723 sprintf(name, "e_nrecdepth%d", i);
0724 sprintf(title, "#RecHits in depth %d (Endcap)", i + 1);
0725 h_evlist3.push_back(new TH1F(name, title, 1000, 0, 100));
0726 h_evlist3[i]->Sumw2();
0727 }
0728 h_etaE = new TH2F("heta", "", 60, -30, 30, 100, 0, 100);
0729 }
0730 }
0731
0732 Bool_t CalibPlotProperties::Notify() {
0733
0734
0735
0736
0737
0738
0739 return kTRUE;
0740 }
0741
0742 void CalibPlotProperties::Show(Long64_t entry) {
0743
0744
0745 if (!fChain)
0746 return;
0747 fChain->Show(entry);
0748 }
0749
0750 Int_t CalibPlotProperties::Cut(Long64_t) {
0751
0752
0753
0754 return 1;
0755 }
0756
0757 void CalibPlotProperties::Loop(Long64_t nentries) {
0758
0759
0760
0761
0762
0763
0764
0765
0766
0767
0768
0769
0770
0771
0772
0773
0774
0775
0776
0777
0778
0779
0780
0781 const bool debug(false);
0782
0783 if (fChain == 0)
0784 return;
0785
0786
0787 if (nentries < 0)
0788 nentries = fChain->GetEntriesFast();
0789 std::cout << "Total entries " << nentries << ":" << fChain->GetEntriesFast() << std::endl;
0790 std::vector<std::pair<int, int> > runEvent;
0791 Long64_t nbytes(0), nb(0);
0792 unsigned int duplicate(0), good(0), kount(0);
0793 double sel(0), selHB(0), selHE(0);
0794 for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0795 Long64_t ientry = LoadTree(jentry);
0796 if (ientry < 0)
0797 break;
0798 nb = fChain->GetEntry(jentry);
0799 nbytes += nb;
0800 if (jentry % 1000000 == 0)
0801 std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
0802 bool select = ((cDuplicate_ != nullptr) && (duplicate_ == 1)) ? (cDuplicate_->isDuplicate(jentry)) : true;
0803 if (!select) {
0804 ++duplicate;
0805 if (debug)
0806 std::cout << "Duplicate event " << t_Run << " " << t_Event << " " << t_p << std::endl;
0807 continue;
0808 }
0809 bool selRun = (includeRun_ ? ((t_Run >= runlo_) && (t_Run <= runhi_)) : ((t_Run < runlo_) || (t_Run > runhi_)));
0810 select =
0811 (selRun && (fabs(t_ieta) >= etalo_) && (fabs(t_ieta) <= etahi_) && (t_nVtx >= nvxlo_) && (t_nVtx <= nvxhi_));
0812 if (!select) {
0813 if (debug)
0814 std::cout << "Run # " << t_Run << " out of range of " << runlo_ << ":" << runhi_ << " or ieta " << t_ieta
0815 << " (" << etalo_ << ":" << etahi_ << ") or nvtx " << t_nVtx << " (" << nvxlo_ << ":" << nvxhi_
0816 << ") out of range" << std::endl;
0817 continue;
0818 }
0819 if (cSelect_ != nullptr) {
0820 if (exclude_) {
0821 if (cSelect_->isItRBX(t_DetIds))
0822 continue;
0823 } else {
0824 if (!(cSelect_->isItRBX(t_ieta, t_iphi)))
0825 continue;
0826 }
0827 }
0828 select = (!cutL1T_ || (t_mindR1 >= 0.5));
0829 if (!select) {
0830 if (debug)
0831 std::cout << "Reject Run # " << t_Run << " Event # " << t_Event << " too close to L1 trigger " << t_mindR1
0832 << std::endl;
0833 continue;
0834 }
0835 select = ((events_.size() == 0) ||
0836 (std::find(events_.begin(), events_.end(), std::pair<int, int>(t_Run, t_Event)) != events_.end()));
0837 if (!select) {
0838 if (debug)
0839 std::cout << "Reject Run # " << t_Run << " Event # " << t_Event << " not in the selection list" << std::endl;
0840 continue;
0841 }
0842
0843 if (plotBasic_) {
0844 h_nvtx->Fill(t_nVtx);
0845 std::pair<int, int> runEv(t_Run, t_Event);
0846 if (std::find(runEvent.begin(), runEvent.end(), runEv) == runEvent.end()) {
0847 h_nvtxEv->Fill(t_nVtx);
0848 runEvent.push_back(runEv);
0849 }
0850 }
0851
0852
0853 int kp = CalibPlots::npbin0;
0854 double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
0855 for (int k = 1; k < CalibPlots::npbin0; ++k) {
0856 if (pmom >= CalibPlots::getMomentum(k - 1) && pmom < CalibPlots::getMomentum(k)) {
0857 kp = k - 1;
0858 break;
0859 }
0860 }
0861 int jp1 = (((std::abs(t_ieta) >= etalo_) && (std::abs(t_ieta) <= etahi_)) ? (std::abs(t_ieta) - etalo_) : -1);
0862 int jp2 = (etahi_ - etalo_ + 1);
0863 int je1 = CalibPlots::netabin;
0864 for (int j = 1; j < CalibPlots::netabin; ++j) {
0865 if (std::abs(t_ieta) > CalibPlots::getEta(j - 1) && std::abs(t_ieta) <= CalibPlots::getEta(j)) {
0866 je1 = j;
0867 break;
0868 }
0869 }
0870 int je2 = (je1 != CalibPlots::netabin) ? 0 : -1;
0871 if (debug)
0872 std::cout << "Bin " << kp << ":" << je1 << ":" << je2 << ":" << jp1 << ":" << jp2 << std::endl;
0873 double cut = (pmom > 20) ? ((flexibleSelect_ == 0) ? 2.0 : 10.0) : 0.0;
0874 double rcut(-1000.0);
0875
0876
0877 double eHcal(t_eHcal);
0878 if (corrFactor_->doCorr()) {
0879 eHcal = 0;
0880 for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
0881
0882 unsigned int id = truncateId((*t_DetIds)[k], truncateFlag_, false);
0883 double cfac = corrFactor_->getCorr(id);
0884 if ((cFactor_ != 0) && (ifDepth_ != 3) && (ifDepth_ > 0))
0885 cfac *= cFactor_->getCorr(t_Run, (*t_DetIds)[k]);
0886 if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr()))
0887 cfac *= cDuplicate_->getWeight((*t_DetIds)[k]);
0888 eHcal += (cfac * ((*t_HitEnergies)[k]));
0889 if (debug) {
0890 int subdet, zside, ieta, iphi, depth;
0891 unpackDetId(id, subdet, zside, ieta, iphi, depth);
0892 std::cout << zside << ":" << ieta << ":" << depth << " Corr " << cfac << " " << (*t_HitEnergies)[k] << " Out "
0893 << eHcal << std::endl;
0894 }
0895 }
0896 }
0897 bool goodTk = goodTrack(eHcal, cut, debug);
0898 bool selPhi = selectPhi(debug);
0899 double rat = (pmom > 0) ? (eHcal / (pmom - t_eMipDR)) : 1.0;
0900 if (debug)
0901 std::cout << "Entry " << jentry << " p|eHcal|ratio " << pmom << "|" << t_eHcal << "|" << eHcal << "|" << rat
0902 << "|" << kp << " Cuts " << t_qltyFlag << "|" << t_selectTk << "|" << (t_hmaxNearP < cut) << "|"
0903 << (t_eMipDR < 1.0) << "|" << goodTk << "|" << (rat > rcut) << " Select Phi " << selPhi << std::endl;
0904 if (plotBasic_) {
0905 h_nvtxTk->Fill(t_nVtx);
0906 h_p[0]->Fill(pmom, t_EventWeight);
0907 h_eta[0]->Fill(t_ieta, t_EventWeight);
0908 if (kp < CalibPlots::npbin0)
0909 h_eta0[kp]->Fill(t_ieta, t_EventWeight);
0910 if (t_qltyFlag) {
0911 h_p[1]->Fill(pmom, t_EventWeight);
0912 h_eta[1]->Fill(t_ieta, t_EventWeight);
0913 if (kp < CalibPlots::npbin0)
0914 h_eta1[kp]->Fill(t_ieta, t_EventWeight);
0915 if (t_selectTk) {
0916 h_p[2]->Fill(pmom, t_EventWeight);
0917 h_eta[2]->Fill(t_ieta, t_EventWeight);
0918 if (kp < CalibPlots::npbin0)
0919 h_eta2[kp]->Fill(t_ieta, t_EventWeight);
0920 if (t_hmaxNearP < cut) {
0921 h_p[3]->Fill(pmom, t_EventWeight);
0922 h_eta[3]->Fill(t_ieta, t_EventWeight);
0923 if (kp < CalibPlots::npbin0)
0924 h_eta3[kp]->Fill(t_ieta, t_EventWeight);
0925 if (t_eMipDR < 1.0) {
0926 h_p[4]->Fill(pmom, t_EventWeight);
0927 h_eta[4]->Fill(t_ieta, t_EventWeight);
0928 if (kp < CalibPlots::npbin0) {
0929 h_eta4[kp]->Fill(t_ieta, t_EventWeight);
0930 h_dL1[kp]->Fill(t_mindR1, t_EventWeight);
0931 h_vtx[kp]->Fill(t_goodPV, t_EventWeight);
0932 }
0933 }
0934 }
0935 }
0936 }
0937 }
0938
0939 if (goodTk && kp != CalibPlots::npbin0 && selPhi) {
0940 if (rat > rcut) {
0941 if (plotEnergy_) {
0942 if (jp1 >= 0) {
0943 h_etaEH[kp][jp1]->Fill(eHcal, t_EventWeight);
0944 h_etaEH[kp][jp2]->Fill(eHcal, t_EventWeight);
0945 h_etaEp[kp][jp1]->Fill(pmom, t_EventWeight);
0946 h_etaEp[kp][jp2]->Fill(pmom, t_EventWeight);
0947 h_etaEE[kp][jp1]->Fill(t_eMipDR, t_EventWeight);
0948 h_etaEE[kp][jp2]->Fill(t_eMipDR, t_EventWeight);
0949 h_etaEE0[kp][jp1]->Fill(t_eMipDR, t_EventWeight);
0950 h_etaEE0[kp][jp2]->Fill(t_eMipDR, t_EventWeight);
0951 }
0952 if (kp == kp50) {
0953 if (je1 != CalibPlots::netabin) {
0954 h_eHcal[je1]->Fill(eHcal, t_EventWeight);
0955 h_eHcal[je2]->Fill(eHcal, t_EventWeight);
0956 h_mom[je1]->Fill(pmom, t_EventWeight);
0957 h_mom[je2]->Fill(pmom, t_EventWeight);
0958 h_eEcal[je1]->Fill(t_eMipDR, t_EventWeight);
0959 h_eEcal[je2]->Fill(t_eMipDR, t_EventWeight);
0960 }
0961 }
0962 }
0963
0964 if (plotHists_) {
0965 if ((std::fabs(rat - 1) < 0.15) && (kp == kp50) && ((std::abs(t_ieta) < 15) || (std::abs(t_ieta) > 17))) {
0966 float weight = (dataMC_ ? t_EventWeight : t_EventWeight * puweight(t_nVtx));
0967 h_etaE->Fill(t_ieta, eHcal, weight);
0968 sel += weight;
0969 std::vector<float> bv(7, 0.0f), ev(7, 0.0f);
0970 std::vector<int> bnrec(7, 0), enrec(7, 0);
0971 double eb(0), ee(0);
0972 for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
0973 unsigned int id = truncateId((*t_DetIds)[k], truncateFlag_, false);
0974 double cfac = corrFactor_->getCorr(id);
0975 if ((cFactor_ != 0) && (ifDepth_ != 3) && (ifDepth_ > 0))
0976 cfac *= cFactor_->getCorr(t_Run, (*t_DetIds)[k]);
0977 if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr()))
0978 cfac *= cDuplicate_->getWeight((*t_DetIds)[k]);
0979 double ener = cfac * (*t_HitEnergies)[k];
0980 if (corrPU_)
0981 correctEnergy(ener);
0982 unsigned int idx = (unsigned int)((*t_DetIds)[k]);
0983 int subdet, zside, ieta, iphi, depth;
0984 unpackDetId(idx, subdet, zside, ieta, iphi, depth);
0985 if (depth > 0 && depth <= CalibPlots::ndepth) {
0986 if (subdet == 1) {
0987 eb += ener;
0988 bv[depth - 1] += ener;
0989 h_bvlist2[depth - 1]->Fill(ener, weight);
0990 ++bnrec[depth - 1];
0991 } else if (subdet == 2) {
0992 ee += ener;
0993 ev[depth - 1] += ener;
0994 h_evlist2[depth - 1]->Fill(ener, weight);
0995 ++enrec[depth - 1];
0996 }
0997 }
0998 }
0999 bool barrel = (eb > ee);
1000 if (barrel)
1001 selHB += weight;
1002 else
1003 selHE += weight;
1004 for (int i = 0; i < CalibPlots::ndepth; i++) {
1005 if (barrel) {
1006 h_bvlist[i]->Fill(bv[i], weight);
1007 h_bvlist3[i]->Fill((bnrec[i] + 0.001), weight);
1008 } else {
1009 h_evlist[i]->Fill(ev[i], weight);
1010 h_evlist3[i]->Fill((enrec[i] + 0.001), weight);
1011 }
1012 }
1013 }
1014 }
1015 }
1016 good++;
1017 }
1018 ++kount;
1019 }
1020 std::cout << "Finds " << duplicate << " Duplicate events out of " << kount << " events in this file and " << good
1021 << " selected events" << std::endl;
1022 if (plotHists_)
1023 std::cout << "Number of weighted selected events " << sel << " HB " << selHB << " HE " << selHE << std::endl;
1024 }
1025
1026 bool CalibPlotProperties::goodTrack(double &eHcal, double &cuti, bool debug) {
1027 bool select(true);
1028 double cut(cuti);
1029 if (debug) {
1030 std::cout << "goodTrack input " << eHcal << ":" << cut;
1031 }
1032 if (flexibleSelect_ > 1) {
1033 double eta = (t_ieta > 0) ? t_ieta : -t_ieta;
1034 cut = 8.0 * exp(eta * log2by18_);
1035 }
1036 correctEnergy(eHcal);
1037 select = ((t_qltyFlag) && (t_selectTk) && (t_hmaxNearP < cut) && (t_eMipDR < 100.0));
1038 if (debug) {
1039 std::cout << " output " << eHcal << ":" << cut << ":" << select << std::endl;
1040 }
1041 return select;
1042 }
1043
1044 bool CalibPlotProperties::selectPhi(bool debug) {
1045 bool select(true);
1046 if (phimin_ > 1 || phimax_ < 72) {
1047 double eTotal(0), eSelec(0);
1048
1049 for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
1050 int iphi = ((*t_DetIds)[k]) & (0x3FF);
1051 int zside = ((*t_DetIds)[k] & 0x80000) ? (1) : (-1);
1052 eTotal += ((*t_HitEnergies)[k]);
1053 if (iphi >= phimin_ && iphi <= phimax_ && zside == zside_)
1054 eSelec += ((*t_HitEnergies)[k]);
1055 }
1056 if (eSelec < 0.9 * eTotal)
1057 select = false;
1058 if (debug) {
1059 std::cout << "Etotal " << eTotal << " and ESelec " << eSelec << " (phi " << phimin_ << ":" << phimax_ << " z "
1060 << zside_ << ") Selection " << select << std::endl;
1061 }
1062 }
1063 return select;
1064 }
1065
1066 void CalibPlotProperties::savePlot(const std::string &theName, bool append, bool all, bool debug) {
1067 TFile *theFile(0);
1068 if (append) {
1069 theFile = new TFile(theName.c_str(), "UPDATE");
1070 } else {
1071 theFile = new TFile(theName.c_str(), "RECREATE");
1072 }
1073
1074 theFile->cd();
1075
1076 if (plotBasic_) {
1077 if (debug)
1078 std::cout << "nvtx " << h_nvtx << ":" << h_nvtxEv << ":" << h_nvtxTk << std::endl;
1079 h_nvtx->Write();
1080 h_nvtxEv->Write();
1081 h_nvtxTk->Write();
1082 for (int k = 0; k < CalibPlots::ntitles; ++k) {
1083 if (debug)
1084 std::cout << "[" << k << "] p " << h_p[k] << " eta " << h_eta[k] << std::endl;
1085 h_p[k]->Write();
1086 h_eta[k]->Write();
1087 }
1088 for (int k = 0; k < CalibPlots::npbin; ++k) {
1089 if (debug)
1090 std::cout << "[" << k << "] eta " << h_eta0[k] << ":" << h_eta1[k] << ":" << h_eta2[k] << ":" << h_eta3[k]
1091 << ":" << h_eta4[k] << " dl " << h_dL1[k] << " vtx " << h_vtx[k] << std::endl;
1092 if (h_eta0[k] != 0) {
1093 TH1D *h1 = (TH1D *)h_eta0[k]->Clone();
1094 h1->Write();
1095 }
1096 if (h_eta1[k] != 0) {
1097 TH1D *h2 = (TH1D *)h_eta1[k]->Clone();
1098 h2->Write();
1099 }
1100 if (h_eta2[k] != 0) {
1101 TH1D *h3 = (TH1D *)h_eta2[k]->Clone();
1102 h3->Write();
1103 }
1104 if (h_eta3[k] != 0) {
1105 TH1D *h4 = (TH1D *)h_eta3[k]->Clone();
1106 h4->Write();
1107 }
1108 if (h_eta4[k] != 0) {
1109 TH1D *h5 = (TH1D *)h_eta4[k]->Clone();
1110 h5->Write();
1111 }
1112 if (h_dL1[k] != 0) {
1113 TH1D *h6 = (TH1D *)h_dL1[k]->Clone();
1114 h6->Write();
1115 }
1116 if (h_vtx[k] != 0) {
1117 TH1D *h7 = (TH1D *)h_vtx[k]->Clone();
1118 h7->Write();
1119 }
1120 }
1121 }
1122
1123 if (plotEnergy_) {
1124 for (int k = 0; k < CalibPlots::npbin; ++k) {
1125 if (debug)
1126 std::cout << "Energy[" << k << "] "
1127 << " eta " << etalo_ << ":" << etahi_ << ":" << CalibPlots::netabin << " etaEH " << h_etaEH[k].size()
1128 << " etaEp " << h_etaEp[k].size() << " etaEE " << h_etaEE[k].size() << " eHcal " << h_eHcal.size()
1129 << " mom " << h_mom.size() << " eEcal " << h_eEcal.size() << std::endl;
1130 for (unsigned int j = 0; j <= (unsigned int)(etahi_ - etalo_); ++j) {
1131 if (h_etaEH[k].size() > j && h_etaEH[k][j] != nullptr && (all || (k == kp50))) {
1132 TH1D *hist = (TH1D *)h_etaEH[k][j]->Clone();
1133 hist->Write();
1134 }
1135 if (h_etaEp[k].size() > j && h_etaEp[k][j] != nullptr && (all || (k == kp50))) {
1136 TH1D *hist = (TH1D *)h_etaEp[k][j]->Clone();
1137 hist->Write();
1138 }
1139 if (h_etaEE[k].size() > j && h_etaEE[k][j] != nullptr && (all || (k == kp50))) {
1140 TH1D *hist = (TH1D *)h_etaEE[k][j]->Clone();
1141 hist->Write();
1142 }
1143 if (h_etaEE0[k].size() > j && h_etaEE0[k][j] != nullptr && (all || (k == kp50))) {
1144 TH1D *hist = (TH1D *)h_etaEE0[k][j]->Clone();
1145 hist->Write();
1146 }
1147 }
1148 }
1149
1150 for (int j = 0; j < CalibPlots::netabin; ++j) {
1151 if (h_eHcal.size() > (unsigned int)(j) && (h_eHcal[j] != nullptr)) {
1152 TH1D *hist = (TH1D *)h_eHcal[j]->Clone();
1153 hist->Write();
1154 }
1155 if (h_mom.size() > (unsigned int)(j) && (h_mom[j] != nullptr)) {
1156 TH1D *hist = (TH1D *)h_mom[j]->Clone();
1157 hist->Write();
1158 }
1159 if (h_eEcal.size() > (unsigned int)(j) && (h_eEcal[j] != nullptr)) {
1160 TH1D *hist = (TH1D *)h_eEcal[j]->Clone();
1161 hist->Write();
1162 }
1163 }
1164 }
1165
1166 if (plotHists_) {
1167 if (debug)
1168 std::cout << "etaE " << h_etaE << std::endl;
1169 h_etaE->Write();
1170 for (int i = 0; i < CalibPlots::ndepth; ++i) {
1171 if (debug)
1172 std::cout << "Depth[" << i << "] bv " << h_bvlist[i] << ":" << h_bvlist2[i] << ":" << h_bvlist3[i] << " ev "
1173 << h_evlist[i] << ":" << h_evlist2[i] << ":" << h_evlist3[i] << std::endl;
1174 h_bvlist[i]->Write();
1175 h_bvlist2[i]->Write();
1176 h_bvlist3[i]->Write();
1177 h_evlist[i]->Write();
1178 h_evlist2[i]->Write();
1179 h_evlist3[i]->Write();
1180 }
1181 }
1182 std::cout << "All done" << std::endl;
1183 theFile->Close();
1184 }
1185
1186 void CalibPlotProperties::correctEnergy(double &eHcal) {
1187 double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
1188 if ((corrPU_ < 0) && (pmom > 0)) {
1189 double ediff = (t_eHcal30 - t_eHcal10);
1190 if (t_DetIds1 != 0 && t_DetIds3 != 0) {
1191 double Etot1(0), Etot3(0);
1192
1193 for (unsigned int idet = 0; idet < (*t_DetIds1).size(); idet++) {
1194 unsigned int id = truncateId((*t_DetIds1)[idet], truncateFlag_, false);
1195 double cfac = corrFactor_->getCorr(id);
1196 if ((cFactor_ != 0) && (ifDepth_ != 3) && (ifDepth_ > 0))
1197 cfac *= cFactor_->getCorr(t_Run, (*t_DetIds1)[idet]);
1198 if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr()))
1199 cfac *= cDuplicate_->getWeight((*t_DetIds1)[idet]);
1200 double hitEn = cfac * (*t_HitEnergies1)[idet];
1201 Etot1 += hitEn;
1202 }
1203 for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
1204 unsigned int id = truncateId((*t_DetIds3)[idet], truncateFlag_, false);
1205 double cfac = corrFactor_->getCorr(id);
1206 if ((cFactor_ != 0) && (ifDepth_ != 3) && (ifDepth_ > 0))
1207 cfac *= cFactor_->getCorr(t_Run, (*t_DetIds3)[idet]);
1208 if ((cDuplicate_ != nullptr) && (cDuplicate_->doCorr()))
1209 cfac *= cDuplicate_->getWeight((*t_DetIds)[idet]);
1210 double hitEn = cfac * (*t_HitEnergies3)[idet];
1211 Etot3 += hitEn;
1212 }
1213 ediff = (Etot3 - Etot1);
1214 }
1215 double fac = puFactor(-corrPU_, t_ieta, pmom, eHcal, ediff);
1216 eHcal *= fac;
1217 } else if (corrPU_ > 1) {
1218 eHcal = puFactorRho(corrPU_, t_ieta, t_rhoh, eHcal);
1219 }
1220 }
1221
1222 void PlotThisHist(TH1D *hist, const std::string &text, int save) {
1223 char namep[120];
1224 sprintf(namep, "c_%s", hist->GetName());
1225 TCanvas *pad = new TCanvas(namep, namep, 700, 500);
1226 pad->SetRightMargin(0.10);
1227 pad->SetTopMargin(0.10);
1228 hist->GetXaxis()->SetTitleSize(0.04);
1229 hist->GetYaxis()->SetTitle("Tracks");
1230 hist->GetYaxis()->SetLabelOffset(0.005);
1231 hist->GetYaxis()->SetTitleSize(0.04);
1232 hist->GetYaxis()->SetLabelSize(0.035);
1233 hist->GetYaxis()->SetTitleOffset(1.10);
1234 hist->SetMarkerStyle(20);
1235 hist->SetMarkerColor(2);
1236 hist->SetLineColor(2);
1237 hist->Draw("Hist");
1238 pad->Modified();
1239 pad->Update();
1240 TPaveStats *st1 = (TPaveStats *)hist->GetListOfFunctions()->FindObject("stats");
1241 TPaveText *txt0 = new TPaveText(0.12, 0.91, 0.49, 0.96, "blNDC");
1242 txt0->SetFillColor(0);
1243 char txt[100];
1244 sprintf(txt, "CMS Simulation Preliminary");
1245 txt0->AddText(txt);
1246 txt0->Draw("same");
1247 TPaveText *txt1 = new TPaveText(0.51, 0.91, 0.90, 0.96, "blNDC");
1248 txt1->SetFillColor(0);
1249 sprintf(txt, "%s", text.c_str());
1250 txt1->AddText(txt);
1251 txt1->Draw("same");
1252 if (st1 != nullptr) {
1253 st1->SetY1NDC(0.70);
1254 st1->SetY2NDC(0.90);
1255 st1->SetX1NDC(0.65);
1256 st1->SetX2NDC(0.90);
1257 }
1258 pad->Update();
1259 if (save != 0) {
1260 if (save > 0)
1261 sprintf(namep, "%s.pdf", pad->GetName());
1262 else
1263 sprintf(namep, "%s.jpg", pad->GetName());
1264 pad->Print(namep);
1265 }
1266 }
1267
1268 void PlotHist(const char *hisFileName,
1269 const std::string &prefix = "",
1270 const std::string &text = "",
1271 int flagC = 111,
1272 int etalo = 0,
1273 int etahi = 30,
1274 int save = 0) {
1275 gStyle->SetCanvasBorderMode(0);
1276 gStyle->SetCanvasColor(kWhite);
1277 gStyle->SetPadColor(kWhite);
1278 gStyle->SetFillColor(kWhite);
1279 gStyle->SetOptTitle(0);
1280 gStyle->SetOptStat(1110);
1281
1282 bool plotBasic = (((flagC / 1) % 10) > 0);
1283 bool plotEnergy = (((flagC / 10) % 10) > 0);
1284 bool plotHists = (((flagC / 100) % 10) > 0);
1285
1286 TFile *file = new TFile(hisFileName);
1287 char name[100], title[200];
1288 TH1D *hist;
1289 if ((file != nullptr) && plotBasic) {
1290 std::cout << "Plot Basic Histos" << std::endl;
1291 hist = (TH1D *)(file->FindObjectAny("hnvtx"));
1292 if (hist != nullptr) {
1293 hist->GetXaxis()->SetTitle("Number of vertices (selected entries)");
1294 PlotThisHist(hist, text, save);
1295 }
1296 hist = (TH1D *)(file->FindObjectAny("hnvtxEv"));
1297 if (hist != nullptr) {
1298 hist->GetXaxis()->SetTitle("Number of vertices (selected events)");
1299 PlotThisHist(hist, text, save);
1300 }
1301 hist = (TH1D *)(file->FindObjectAny("hnvtxTk"));
1302 if (hist != nullptr) {
1303 hist->GetXaxis()->SetTitle("Number of vertices (selected tracks)");
1304 PlotThisHist(hist, text, save);
1305 }
1306 for (int k = 0; k < CalibPlots::ntitles; ++k) {
1307 sprintf(name, "%sp%d", prefix.c_str(), k);
1308 hist = (TH1D *)(file->FindObjectAny(name));
1309 if (hist != nullptr) {
1310 sprintf(title, "Momentum for %s (GeV)", CalibPlots::getTitle(k).c_str());
1311 hist->GetXaxis()->SetTitle(title);
1312 PlotThisHist(hist, text, save);
1313 }
1314 sprintf(name, "%seta%d", prefix.c_str(), k);
1315 hist = (TH1D *)(file->FindObjectAny(name));
1316 if (hist != nullptr) {
1317 sprintf(title, "#eta for %s", CalibPlots::getTitle(k).c_str());
1318 hist->GetXaxis()->SetTitle(title);
1319 PlotThisHist(hist, text, save);
1320 }
1321 }
1322 for (int k = 0; k < CalibPlots::npbin; ++k) {
1323 sprintf(name, "%seta0%d", prefix.c_str(), k);
1324 hist = (TH1D *)(file->FindObjectAny(name));
1325 if (hist != nullptr) {
1326 sprintf(title,
1327 "#eta for %s (p = %d:%d GeV)",
1328 CalibPlots::getTitle(0).c_str(),
1329 CalibPlots::getP(k),
1330 CalibPlots::getP(k + 1));
1331 hist->GetXaxis()->SetTitle(title);
1332 PlotThisHist(hist, text, save);
1333 }
1334 sprintf(name, "%seta1%d", prefix.c_str(), k);
1335 hist = (TH1D *)(file->FindObjectAny(name));
1336 if (hist != nullptr) {
1337 sprintf(title,
1338 "#eta for %s (p = %d:%d GeV)",
1339 CalibPlots::getTitle(1).c_str(),
1340 CalibPlots::getP(k),
1341 CalibPlots::getP(k + 1));
1342 hist->GetXaxis()->SetTitle(title);
1343 PlotThisHist(hist, text, save);
1344 }
1345 sprintf(name, "%seta2%d", prefix.c_str(), k);
1346 hist = (TH1D *)(file->FindObjectAny(name));
1347 if (hist != nullptr) {
1348 sprintf(title,
1349 "#eta for %s (p = %d:%d GeV)",
1350 CalibPlots::getTitle(2).c_str(),
1351 CalibPlots::getP(k),
1352 CalibPlots::getP(k + 1));
1353 hist->GetXaxis()->SetTitle(title);
1354 PlotThisHist(hist, text, save);
1355 }
1356 sprintf(name, "%seta3%d", prefix.c_str(), k);
1357 hist = (TH1D *)(file->FindObjectAny(name));
1358 if (hist != nullptr) {
1359 sprintf(title,
1360 "#eta for %s (p = %d:%d GeV)",
1361 CalibPlots::getTitle(3).c_str(),
1362 CalibPlots::getP(k),
1363 CalibPlots::getP(k + 1));
1364 hist->GetXaxis()->SetTitle(title);
1365 PlotThisHist(hist, text, save);
1366 }
1367 sprintf(name, "%seta4%d", prefix.c_str(), k);
1368 hist = (TH1D *)(file->FindObjectAny(name));
1369 if (hist != nullptr) {
1370 sprintf(title,
1371 "#eta for %s (p = %d:%d GeV)",
1372 CalibPlots::getTitle(4).c_str(),
1373 CalibPlots::getP(k),
1374 CalibPlots::getP(k + 1));
1375 hist->GetXaxis()->SetTitle(title);
1376 PlotThisHist(hist, text, save);
1377 }
1378 sprintf(name, "%sdl1%d", prefix.c_str(), k);
1379 hist = (TH1D *)(file->FindObjectAny(name));
1380 if (hist != nullptr) {
1381 sprintf(title, "Distance from L1 (p = %d:%d GeV)", CalibPlots::getP(k), CalibPlots::getP(k + 1));
1382 hist->GetXaxis()->SetTitle(title);
1383 PlotThisHist(hist, text, save);
1384 }
1385 sprintf(name, "%svtx%d", prefix.c_str(), k);
1386 hist = (TH1D *)(file->FindObjectAny(name));
1387 if (hist != nullptr) {
1388 sprintf(title, "N_{Vertex} (p = %d:%d GeV)", CalibPlots::getP(k), CalibPlots::getP(k + 1));
1389 hist->GetXaxis()->SetTitle(title);
1390 PlotThisHist(hist, text, save);
1391 }
1392 }
1393 }
1394
1395 if ((file != nullptr) && plotEnergy) {
1396 std::cout << "Make plots for good tracks" << std::endl;
1397 for (int k = 0; k < CalibPlots::npbin; ++k) {
1398 for (int j = etalo; j <= etahi + 1; ++j) {
1399 sprintf(name, "%senergyH%d%d", prefix.c_str(), k, j);
1400 hist = (TH1D *)(file->FindObjectAny(name));
1401 if (hist != nullptr) {
1402 if (j > etahi)
1403 sprintf(title,
1404 "HCAL energy for %s (p = %d:%d GeV |#eta| = %d:%d)",
1405 CalibPlots::getTitle(3).c_str(),
1406 CalibPlots::getP(k),
1407 CalibPlots::getP(k + 1),
1408 etalo,
1409 etahi);
1410 else
1411 sprintf(title,
1412 "HCAL energy for %s (p = %d:%d GeV |#eta| = %d)",
1413 CalibPlots::getTitle(3).c_str(),
1414 CalibPlots::getP(k),
1415 CalibPlots::getP(k + 1),
1416 j);
1417 hist->GetXaxis()->SetTitle(title);
1418 PlotThisHist(hist, text, save);
1419 }
1420 sprintf(name, "%senergyP%d%d", prefix.c_str(), k, j);
1421 hist = (TH1D *)(file->FindObjectAny(name));
1422 if (hist != nullptr) {
1423 if (j > etahi)
1424 sprintf(title,
1425 "momentum for %s (p = %d:%d GeV |#eta| = %d:%d)",
1426 CalibPlots::getTitle(3).c_str(),
1427 CalibPlots::getP(k),
1428 CalibPlots::getP(k + 1),
1429 etalo,
1430 etahi);
1431 else
1432 sprintf(title,
1433 "momentum for %s (p = %d:%d GeV |#eta| = %d)",
1434 CalibPlots::getTitle(3).c_str(),
1435 CalibPlots::getP(k),
1436 CalibPlots::getP(k + 1),
1437 j);
1438 hist->GetXaxis()->SetTitle(title);
1439 PlotThisHist(hist, text, save);
1440 }
1441 sprintf(name, "%senergyE%d%d", prefix.c_str(), k, j);
1442 hist = (TH1D *)(file->FindObjectAny(name));
1443 if (hist != nullptr) {
1444 if (j > etahi)
1445 sprintf(title,
1446 "ECAL energy for %s (p = %d:%d GeV |#eta| = %d:%d)",
1447 CalibPlots::getTitle(3).c_str(),
1448 CalibPlots::getP(k),
1449 CalibPlots::getP(k + 1),
1450 etalo,
1451 etahi);
1452 else
1453 sprintf(title,
1454 "ECAL energy for %s (p = %d:%d GeV |#eta| = %d)",
1455 CalibPlots::getTitle(3).c_str(),
1456 CalibPlots::getP(k),
1457 CalibPlots::getP(k + 1),
1458 j);
1459 hist->GetXaxis()->SetTitle(title);
1460 PlotThisHist(hist, text, save);
1461 }
1462 sprintf(name, "%senergyER%d%d", prefix.c_str(), k, j);
1463 hist = (TH1D *)(file->FindObjectAny(name));
1464 if (hist != nullptr) {
1465 std::cout << name << " Mean " << hist->GetMean() << " +- " << hist->GetMeanError() << " Entries "
1466 << hist->GetEntries() << " RMS " << hist->GetRMS() << std::endl;
1467 }
1468 }
1469 }
1470
1471 for (int j = 0; j < CalibPlots::netabin; ++j) {
1472 sprintf(name, "%senergyH%d", prefix.c_str(), j);
1473 hist = (TH1D *)(file->FindObjectAny(name));
1474 if (hist != nullptr) {
1475 if (j == 0)
1476 sprintf(title, "HCAL energy for %s (All)", CalibPlots::getTitle(3).c_str());
1477 else
1478 sprintf(title,
1479 "HCAL energy for %s (|#eta| = %d:%d)",
1480 CalibPlots::getTitle(3).c_str(),
1481 CalibPlots::getEta(j - 1),
1482 CalibPlots::getEta(j));
1483 hist->GetXaxis()->SetTitle(title);
1484 PlotThisHist(hist, text, save);
1485 }
1486 sprintf(name, "%senergyP%d", prefix.c_str(), j);
1487 hist = (TH1D *)(file->FindObjectAny(name));
1488 if (hist != nullptr) {
1489 if (j == 0)
1490 sprintf(title, "Track momentum for %s (All)", CalibPlots::getTitle(3).c_str());
1491 else
1492 sprintf(title,
1493 "Track momentum for %s (|#eta| = %d:%d)",
1494 CalibPlots::getTitle(3).c_str(),
1495 CalibPlots::getEta(j - 1),
1496 CalibPlots::getEta(j));
1497 hist->GetXaxis()->SetTitle(title);
1498 PlotThisHist(hist, text, save);
1499 }
1500 sprintf(name, "%senergyE%d", prefix.c_str(), j);
1501 hist = (TH1D *)(file->FindObjectAny(name));
1502 if (hist != nullptr) {
1503 if (j == 0)
1504 sprintf(title, "ECAL energy for %s (All)", CalibPlots::getTitle(3).c_str());
1505 else
1506 sprintf(title,
1507 "ECAL energy for %s (|#eta| = %d:%d)",
1508 CalibPlots::getTitle(3).c_str(),
1509 CalibPlots::getEta(j - 1),
1510 CalibPlots::getEta(j));
1511 hist->GetXaxis()->SetTitle(title);
1512 PlotThisHist(hist, text, save);
1513 }
1514 }
1515 }
1516
1517 if (plotHists) {
1518 for (int i = 0; i < CalibPlots::ndepth; i++) {
1519 sprintf(name, "b_edepth%d", i);
1520 hist = (TH1D *)(file->FindObjectAny(name));
1521 if (hist != nullptr) {
1522 sprintf(title, "Total RecHit energy in depth %d (Barrel)", i + 1);
1523 hist->GetXaxis()->SetTitle(title);
1524 PlotThisHist(hist, text, save);
1525 }
1526 sprintf(name, "b_recedepth%d", i);
1527 hist = (TH1D *)(file->FindObjectAny(name));
1528 if (hist != nullptr) {
1529 sprintf(title, "RecHit energy in depth %d (Barrel)", i + 1);
1530 hist->GetXaxis()->SetTitle(title);
1531 PlotThisHist(hist, text, save);
1532 }
1533 sprintf(name, "b_nrecdepth%d", i);
1534 hist = (TH1D *)(file->FindObjectAny(name));
1535 if (hist != nullptr) {
1536 sprintf(title, "#RecHits in depth %d (Barrel)", i + 1);
1537 hist->GetXaxis()->SetTitle(title);
1538 PlotThisHist(hist, text, save);
1539 }
1540 sprintf(name, "e_edepth%d", i);
1541 hist = (TH1D *)(file->FindObjectAny(name));
1542 if (hist != nullptr) {
1543 sprintf(title, "Total RecHit energy in depth %d (Endcap)", i + 1);
1544 hist->GetXaxis()->SetTitle(title);
1545 PlotThisHist(hist, text, save);
1546 }
1547 sprintf(name, "e_recedepth%d", i);
1548 hist = (TH1D *)(file->FindObjectAny(name));
1549 if (hist != nullptr) {
1550 sprintf(title, "RecHit energy in depth %d (Endcap)", i + 1);
1551 hist->GetXaxis()->SetTitle(title);
1552 PlotThisHist(hist, text, save);
1553 }
1554 sprintf(name, "e_nrecdepth%d", i);
1555 hist = (TH1D *)(file->FindObjectAny(name));
1556 if (hist != nullptr) {
1557 sprintf(title, "#RecHits in depth %d (Endcap)", i + 1);
1558 hist->GetXaxis()->SetTitle(title);
1559 PlotThisHist(hist, text, save);
1560 }
1561 }
1562 TH2F *h_etaE = (TH2F *)(file->FindObjectAny("heta"));
1563 if (h_etaE != nullptr) {
1564 h_etaE->GetXaxis()->SetTitle("i#eta");
1565 h_etaE->GetYaxis()->SetTitle("Energy (GeV)");
1566 char namep[120];
1567 sprintf(namep, "c_%s", h_etaE->GetName());
1568 TCanvas *pad = new TCanvas(namep, namep, 700, 700);
1569 pad->SetRightMargin(0.10);
1570 pad->SetTopMargin(0.10);
1571 h_etaE->GetXaxis()->SetTitleSize(0.04);
1572 h_etaE->GetYaxis()->SetLabelOffset(0.005);
1573 h_etaE->GetYaxis()->SetTitleSize(0.04);
1574 h_etaE->GetYaxis()->SetLabelSize(0.035);
1575 h_etaE->GetYaxis()->SetTitleOffset(1.10);
1576 h_etaE->SetMarkerStyle(20);
1577 h_etaE->SetMarkerColor(2);
1578 h_etaE->SetLineColor(2);
1579 h_etaE->Draw();
1580 pad->Update();
1581 TPaveStats *st1 = (TPaveStats *)h_etaE->GetListOfFunctions()->FindObject("stats");
1582 if (st1 != nullptr) {
1583 st1->SetY1NDC(0.70);
1584 st1->SetY2NDC(0.90);
1585 st1->SetX1NDC(0.65);
1586 st1->SetX2NDC(0.90);
1587 }
1588 if (save != 0) {
1589 if (save > 0)
1590 sprintf(namep, "%s.pdf", pad->GetName());
1591 else
1592 sprintf(namep, "%s.jpg", pad->GetName());
1593 pad->Print(namep);
1594 }
1595 }
1596 }
1597 }
1598
1599 class CalibSplit {
1600 public:
1601 TChain *fChain;
1602 Int_t fCurrent;
1603
1604
1605 Int_t t_Run;
1606 Int_t t_Event;
1607 Int_t t_DataType;
1608 Int_t t_ieta;
1609 Int_t t_iphi;
1610 Double_t t_EventWeight;
1611 Int_t t_nVtx;
1612 Int_t t_nTrk;
1613 Int_t t_goodPV;
1614 Double_t t_l1pt;
1615 Double_t t_l1eta;
1616 Double_t t_l1phi;
1617 Double_t t_l3pt;
1618 Double_t t_l3eta;
1619 Double_t t_l3phi;
1620 Double_t t_p;
1621 Double_t t_pt;
1622 Double_t t_phi;
1623 Double_t t_mindR1;
1624 Double_t t_mindR2;
1625 Double_t t_eMipDR;
1626 Double_t t_eHcal;
1627 Double_t t_eHcal10;
1628 Double_t t_eHcal30;
1629 Double_t t_hmaxNearP;
1630 Double_t t_rhoh;
1631 Bool_t t_selectTk;
1632 Bool_t t_qltyFlag;
1633 Bool_t t_qltyMissFlag;
1634 Bool_t t_qltyPVFlag;
1635 Double_t t_gentrackP;
1636 std::vector<unsigned int> *t_DetIds;
1637 std::vector<double> *t_HitEnergies;
1638 std::vector<bool> *t_trgbits;
1639 std::vector<unsigned int> *t_DetIds1;
1640 std::vector<unsigned int> *t_DetIds3;
1641 std::vector<double> *t_HitEnergies1;
1642 std::vector<double> *t_HitEnergies3;
1643
1644
1645 TBranch *b_t_Run;
1646 TBranch *b_t_Event;
1647 TBranch *b_t_DataType;
1648 TBranch *b_t_ieta;
1649 TBranch *b_t_iphi;
1650 TBranch *b_t_EventWeight;
1651 TBranch *b_t_nVtx;
1652 TBranch *b_t_nTrk;
1653 TBranch *b_t_goodPV;
1654 TBranch *b_t_l1pt;
1655 TBranch *b_t_l1eta;
1656 TBranch *b_t_l1phi;
1657 TBranch *b_t_l3pt;
1658 TBranch *b_t_l3eta;
1659 TBranch *b_t_l3phi;
1660 TBranch *b_t_p;
1661 TBranch *b_t_pt;
1662 TBranch *b_t_phi;
1663 TBranch *b_t_mindR1;
1664 TBranch *b_t_mindR2;
1665 TBranch *b_t_eMipDR;
1666 TBranch *b_t_eHcal;
1667 TBranch *b_t_eHcal10;
1668 TBranch *b_t_eHcal30;
1669 TBranch *b_t_hmaxNearP;
1670 TBranch *b_t_rhoh;
1671 TBranch *b_t_selectTk;
1672 TBranch *b_t_qltyFlag;
1673 TBranch *b_t_qltyMissFlag;
1674 TBranch *b_t_qltyPVFlag;
1675 TBranch *b_t_gentrackP;
1676 TBranch *b_t_DetIds;
1677 TBranch *b_t_HitEnergies;
1678 TBranch *b_t_trgbits;
1679 TBranch *b_t_DetIds1;
1680 TBranch *b_t_DetIds3;
1681 TBranch *b_t_HitEnergies1;
1682 TBranch *b_t_HitEnergies3;
1683
1684
1685 Int_t tout_Run;
1686 Int_t tout_Event;
1687 Int_t tout_DataType;
1688 Int_t tout_ieta;
1689 Int_t tout_iphi;
1690 Double_t tout_EventWeight;
1691 Int_t tout_nVtx;
1692 Int_t tout_nTrk;
1693 Int_t tout_goodPV;
1694 Double_t tout_l1pt;
1695 Double_t tout_l1eta;
1696 Double_t tout_l1phi;
1697 Double_t tout_l3pt;
1698 Double_t tout_l3eta;
1699 Double_t tout_l3phi;
1700 Double_t tout_p;
1701 Double_t tout_pt;
1702 Double_t tout_phi;
1703 Double_t tout_mindR1;
1704 Double_t tout_mindR2;
1705 Double_t tout_eMipDR;
1706 Double_t tout_eHcal;
1707 Double_t tout_eHcal10;
1708 Double_t tout_eHcal30;
1709 Double_t tout_hmaxNearP;
1710 Double_t tout_rhoh;
1711 Bool_t tout_selectTk;
1712 Bool_t tout_qltyFlag;
1713 Bool_t tout_qltyMissFlag;
1714 Bool_t tout_qltyPVFlag;
1715 Double_t tout_gentrackP;
1716 std::vector<unsigned int> *tout_DetIds;
1717 std::vector<double> *tout_HitEnergies;
1718 std::vector<bool> *tout_trgbits;
1719 std::vector<unsigned int> *tout_DetIds1;
1720 std::vector<unsigned int> *tout_DetIds3;
1721 std::vector<double> *tout_HitEnergies1;
1722 std::vector<double> *tout_HitEnergies3;
1723
1724 CalibSplit(const char *fname,
1725 const std::string &dirname,
1726 const std::string &outFileName,
1727 double pmin = 40.0,
1728 double pmax = 60.0,
1729 bool debug = false);
1730 virtual ~CalibSplit();
1731 virtual Int_t Cut(Long64_t entry);
1732 virtual Int_t GetEntry(Long64_t entry);
1733 virtual Long64_t LoadTree(Long64_t entry);
1734 virtual void Init(TChain *);
1735 virtual void Loop(Long64_t nentries = -1);
1736 virtual Bool_t Notify();
1737 virtual void Show(Long64_t entry = -1);
1738 void copyTree();
1739 void close();
1740
1741 private:
1742 const std::string fname_, dirnm_, outFileName_;
1743 const double pmin_, pmax_;
1744 const bool debug_;
1745 TFile *outputFile_;
1746 TDirectoryFile *outputDir_;
1747 TTree *outputTree_;
1748 };
1749
1750 CalibSplit::CalibSplit(
1751 const char *fname, const std::string &dirnm, const std::string &outFileName, double pmin, double pmax, bool debug)
1752 : fname_(fname),
1753 dirnm_(dirnm),
1754 outFileName_(outFileName),
1755 pmin_(pmin),
1756 pmax_(pmax),
1757 debug_(debug),
1758 outputFile_(nullptr),
1759 outputDir_(nullptr),
1760 outputTree_(nullptr) {
1761 char treeName[400];
1762 sprintf(treeName, "%s/CalibTree", dirnm.c_str());
1763 TChain *chain = new TChain(treeName);
1764 std::cout << "Create a chain for " << treeName << " from " << fname << " to write tracs of momentum in the range "
1765 << pmin_ << ":" << pmax_ << " to file " << outFileName_ << std::endl;
1766 if (!fillChain(chain, fname)) {
1767 std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
1768 } else {
1769 std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
1770 Init(chain);
1771 }
1772 }
1773
1774 CalibSplit::~CalibSplit() {
1775 if (!fChain)
1776 return;
1777 delete fChain->GetCurrentFile();
1778 }
1779
1780 Int_t CalibSplit::GetEntry(Long64_t entry) {
1781
1782 if (!fChain)
1783 return 0;
1784 return fChain->GetEntry(entry);
1785 }
1786
1787 Long64_t CalibSplit::LoadTree(Long64_t entry) {
1788
1789 if (!fChain)
1790 return -5;
1791 Long64_t centry = fChain->LoadTree(entry);
1792 if (centry < 0)
1793 return centry;
1794 if (!fChain->InheritsFrom(TChain::Class()))
1795 return centry;
1796 TChain *chain = (TChain *)fChain;
1797 if (chain->GetTreeNumber() != fCurrent) {
1798 fCurrent = chain->GetTreeNumber();
1799 Notify();
1800 }
1801 return centry;
1802 }
1803
1804 void CalibSplit::Init(TChain *tree) {
1805
1806
1807
1808
1809
1810
1811
1812
1813
1814 t_DetIds = nullptr;
1815 t_DetIds1 = nullptr;
1816 t_DetIds3 = nullptr;
1817 t_HitEnergies = nullptr;
1818 t_HitEnergies1 = nullptr;
1819 t_HitEnergies3 = nullptr;
1820 t_trgbits = nullptr;
1821
1822 fChain = tree;
1823 fCurrent = -1;
1824 if (!tree)
1825 return;
1826 fChain->SetMakeClass(1);
1827
1828 fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
1829 fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
1830 fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
1831 fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
1832 fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
1833 fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
1834 fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
1835 fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
1836 fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
1837 fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
1838 fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
1839 fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
1840 fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
1841 fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
1842 fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
1843 fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
1844 fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
1845 fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
1846 fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
1847 fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
1848 fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
1849 fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
1850 fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
1851 fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
1852 fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
1853 fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
1854 fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
1855 fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
1856 fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
1857 fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
1858 fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
1859 fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
1860 fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
1861 fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
1862 fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
1863 fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
1864 fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
1865 fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
1866 Notify();
1867
1868 tout_DetIds = new std::vector<unsigned int>();
1869 tout_HitEnergies = new std::vector<double>();
1870 tout_trgbits = new std::vector<bool>();
1871 tout_DetIds1 = new std::vector<unsigned int>();
1872 tout_DetIds3 = new std::vector<unsigned int>();
1873 tout_HitEnergies1 = new std::vector<double>();
1874 tout_HitEnergies3 = new std::vector<double>();
1875
1876 outputFile_ = TFile::Open(outFileName_.c_str(), "RECREATE");
1877 outputFile_->cd();
1878 outputDir_ = new TDirectoryFile(dirnm_.c_str(), dirnm_.c_str());
1879 outputDir_->cd();
1880 outputTree_ = new TTree("CalibTree", "CalibTree");
1881 outputTree_->Branch("t_Run", &tout_Run);
1882 outputTree_->Branch("t_Event", &tout_Event);
1883 outputTree_->Branch("t_DataType", &tout_DataType);
1884 outputTree_->Branch("t_ieta", &tout_ieta);
1885 outputTree_->Branch("t_iphi", &tout_iphi);
1886 outputTree_->Branch("t_EventWeight", &tout_EventWeight);
1887 outputTree_->Branch("t_nVtx", &tout_nVtx);
1888 outputTree_->Branch("t_nTrk", &tout_nTrk);
1889 outputTree_->Branch("t_goodPV", &tout_goodPV);
1890 outputTree_->Branch("t_l1pt", &tout_l1pt);
1891 outputTree_->Branch("t_l1eta", &tout_l1eta);
1892 outputTree_->Branch("t_l1phi", &tout_l1phi);
1893 outputTree_->Branch("t_l3pt", &tout_l3pt);
1894 outputTree_->Branch("t_l3eta", &tout_l3eta);
1895 outputTree_->Branch("t_l3phi", &tout_l3phi);
1896 outputTree_->Branch("t_p", &tout_p);
1897 outputTree_->Branch("t_pt", &tout_pt);
1898 outputTree_->Branch("t_phi", &tout_phi);
1899 outputTree_->Branch("t_mindR1", &tout_mindR1);
1900 outputTree_->Branch("t_mindR2", &tout_mindR2);
1901 outputTree_->Branch("t_eMipDR", &tout_eMipDR);
1902 outputTree_->Branch("t_eHcal", &tout_eHcal);
1903 outputTree_->Branch("t_eHcal10", &tout_eHcal10);
1904 outputTree_->Branch("t_eHcal30", &tout_eHcal30);
1905 outputTree_->Branch("t_hmaxNearP", &tout_hmaxNearP);
1906 outputTree_->Branch("t_rhoh", &tout_rhoh);
1907 outputTree_->Branch("t_selectTk", &tout_selectTk);
1908 outputTree_->Branch("t_qltyFlag", &tout_qltyFlag);
1909 outputTree_->Branch("t_qltyMissFlag", &tout_qltyMissFlag);
1910 outputTree_->Branch("t_qltyPVFlag", &tout_qltyPVFlag);
1911 outputTree_->Branch("t_gentrackP", &tout_gentrackP);
1912 outputTree_->Branch("t_DetIds", &tout_DetIds);
1913 outputTree_->Branch("t_HitEnergies", &tout_HitEnergies);
1914 outputTree_->Branch("t_trgbits", &tout_trgbits);
1915 outputTree_->Branch("t_DetIds1", &tout_DetIds1);
1916 outputTree_->Branch("t_DetIds3", &tout_DetIds3);
1917 outputTree_->Branch("t_HitEnergies1", &tout_HitEnergies1);
1918 outputTree_->Branch("t_HitEnergies3", &tout_HitEnergies3);
1919
1920 std::cout << "Output CalibTree is created in directory " << dirnm_ << " of " << outFileName_ << std::endl;
1921 }
1922
1923 Bool_t CalibSplit::Notify() {
1924
1925
1926
1927
1928
1929
1930 return kTRUE;
1931 }
1932
1933 void CalibSplit::Show(Long64_t entry) {
1934
1935
1936 if (!fChain)
1937 return;
1938 fChain->Show(entry);
1939 }
1940
1941 Int_t CalibSplit::Cut(Long64_t) {
1942
1943
1944
1945 return 1;
1946 }
1947
1948 void CalibSplit::Loop(Long64_t nentries) {
1949
1950
1951
1952
1953
1954
1955
1956
1957
1958
1959
1960
1961
1962
1963
1964
1965
1966
1967
1968
1969
1970
1971
1972
1973 if (fChain == 0)
1974 return;
1975
1976
1977 if (nentries < 0)
1978 nentries = fChain->GetEntriesFast();
1979 std::cout << "Total entries " << nentries << ":" << fChain->GetEntriesFast() << std::endl;
1980 Long64_t nbytes(0), nb(0);
1981 unsigned int good(0), reject(0), kount(0);
1982 for (Long64_t jentry = 0; jentry < nentries; jentry++) {
1983 Long64_t ientry = LoadTree(jentry);
1984 if (ientry < 0)
1985 break;
1986 nb = fChain->GetEntry(jentry);
1987 nbytes += nb;
1988 if (jentry % 1000000 == 0)
1989 std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
1990 ++kount;
1991 bool select = ((t_p >= pmin_) && (t_p < pmax_));
1992 if (!select) {
1993 ++reject;
1994 if (debug_)
1995 std::cout << "Rejected event " << t_Run << " " << t_Event << " " << t_p << std::endl;
1996 continue;
1997 }
1998 ++good;
1999 copyTree();
2000 outputTree_->Fill();
2001 }
2002 std::cout << "\nSelect " << good << " Reject " << reject << " entries out of a total of " << kount << " counts"
2003 << std::endl;
2004 close();
2005 }
2006
2007 void CalibSplit::copyTree() {
2008 tout_Run = t_Run;
2009 tout_Event = t_Event;
2010 tout_DataType = t_DataType;
2011 tout_ieta = t_ieta;
2012 tout_iphi = t_iphi;
2013 tout_EventWeight = t_EventWeight;
2014 tout_nVtx = t_nVtx;
2015 tout_nTrk = t_nTrk;
2016 tout_goodPV = t_goodPV;
2017 tout_l1pt = t_l1pt;
2018 tout_l1eta = t_l1eta;
2019 tout_l1phi = t_l1phi;
2020 tout_l3pt = t_l3pt;
2021 tout_l3eta = t_l3eta;
2022 tout_l3phi = t_l3phi;
2023 tout_p = t_p;
2024 tout_pt = t_pt;
2025 tout_phi = t_phi;
2026 tout_mindR1 = t_mindR1;
2027 tout_mindR2 = t_mindR2;
2028 tout_eMipDR = t_eMipDR;
2029 tout_eHcal = t_eHcal;
2030 tout_eHcal10 = t_eHcal10;
2031 tout_eHcal30 = t_eHcal30;
2032 tout_hmaxNearP = t_hmaxNearP;
2033 tout_rhoh = t_rhoh;
2034 tout_selectTk = t_selectTk;
2035 tout_qltyFlag = t_qltyFlag;
2036 tout_qltyMissFlag = t_qltyMissFlag;
2037 tout_qltyPVFlag = t_qltyPVFlag;
2038 tout_gentrackP = t_gentrackP;
2039 tout_DetIds->clear();
2040 if (t_DetIds != nullptr) {
2041 tout_DetIds->reserve(t_DetIds->size());
2042 for (unsigned int i = 0; i < t_DetIds->size(); ++i)
2043 tout_DetIds->push_back((*t_DetIds)[i]);
2044 }
2045 tout_HitEnergies->clear();
2046 if (t_HitEnergies != nullptr) {
2047 tout_HitEnergies->reserve(t_HitEnergies->size());
2048 for (unsigned int i = 0; i < t_HitEnergies->size(); ++i)
2049 tout_HitEnergies->push_back((*t_HitEnergies)[i]);
2050 }
2051 tout_trgbits->clear();
2052 if (t_trgbits != nullptr) {
2053 tout_trgbits->reserve(t_trgbits->size());
2054 for (unsigned int i = 0; i < t_trgbits->size(); ++i)
2055 tout_trgbits->push_back((*t_trgbits)[i]);
2056 }
2057 tout_DetIds1->clear();
2058 if (t_DetIds1 != nullptr) {
2059 tout_DetIds1->reserve(t_DetIds1->size());
2060 for (unsigned int i = 0; i < t_DetIds1->size(); ++i)
2061 tout_DetIds1->push_back((*t_DetIds1)[i]);
2062 }
2063 tout_DetIds3->clear();
2064 if (t_DetIds3 != nullptr) {
2065 tout_DetIds3->reserve(t_DetIds3->size());
2066 for (unsigned int i = 0; i < t_DetIds3->size(); ++i)
2067 tout_DetIds3->push_back((*t_DetIds3)[i]);
2068 }
2069 tout_HitEnergies1->clear();
2070 if (t_HitEnergies1 != nullptr) {
2071 tout_HitEnergies1->reserve(t_HitEnergies1->size());
2072 for (unsigned int i = 0; i < t_HitEnergies1->size(); ++i)
2073 tout_HitEnergies1->push_back((*t_HitEnergies1)[i]);
2074 }
2075 tout_HitEnergies3->clear();
2076 if (t_HitEnergies1 != nullptr) {
2077 tout_HitEnergies3->reserve(t_HitEnergies3->size());
2078 for (unsigned int i = 0; i < t_HitEnergies3->size(); ++i)
2079 tout_HitEnergies3->push_back((*t_HitEnergies3)[i]);
2080 }
2081 }
2082
2083 void CalibSplit::close() {
2084 if (outputFile_) {
2085 outputDir_->cd();
2086 std::cout << "file yet to be Written" << std::endl;
2087 outputTree_->Write();
2088 std::cout << "file Written" << std::endl;
2089 outputFile_->Close();
2090 std::cout << "now doing return" << std::endl;
2091 }
2092 outputFile_ = nullptr;
2093 }