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