File indexing completed on 2021-04-09 02:50:18
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 #include <TStyle.h>
0091 #include <TCanvas.h>
0092 #include <TROOT.h>
0093 #include <TChain.h>
0094 #include <TFile.h>
0095 #include <TFitResult.h>
0096 #include <TFitResultPtr.h>
0097 #include <TTree.h>
0098 #include <TH1.h>
0099 #include <TGraph.h>
0100 #include <TProfile.h>
0101 #include <algorithm>
0102 #include <vector>
0103 #include <string>
0104 #include <iomanip>
0105 #include <iostream>
0106 #include <fstream>
0107 #include <sstream>
0108
0109 #include "CalibCorr.C"
0110
0111 void Run(const char *inFileName = "Silver.root",
0112 const char *dirName = "HcalIsoTrkAnalyzer",
0113 const char *treeName = "CalibTree",
0114 const char *outFileName = "Silver_out.root",
0115 const char *corrFileName = "Silver_corr.txt",
0116 const char *dupFileName = "events_DXS2.txt",
0117 const char *rcorFileName = "",
0118 bool useweight = true,
0119 bool useMean = false,
0120 int nMin = 0,
0121 bool inverse = true,
0122 double ratMin = 0.25,
0123 double ratMax = 3.,
0124 int ietaMax = 25,
0125 int ietaTrack = -1,
0126 int sysmode = -1,
0127 int puCorr = -5,
0128 int applyL1Cut = 1,
0129 double l1Cut = 0.5,
0130 int truncateFlag = 0,
0131 int maxIter = 30,
0132 int rcorForm = 0,
0133 bool useGen = false,
0134 int runlo = 0,
0135 int runhi = 99999999,
0136 int phimin = 1,
0137 int phimax = 72,
0138 int zside = 0,
0139 int nvxlo = 0,
0140 int nvxhi = 1000,
0141 int rbx = 0,
0142 bool exclude = true,
0143 int higheta = 1,
0144 double fraction = 1.0,
0145 bool writeDebugHisto = false,
0146 bool debug = false);
0147
0148
0149
0150 class CalibTree {
0151 public:
0152 struct myEntry {
0153 myEntry(int k = 0, double f0 = 0, double f1 = 0, double f2 = 0) : kount(k), fact0(f0), fact1(f1), fact2(f2) {}
0154 int kount;
0155 double fact0, fact1, fact2;
0156 };
0157
0158 struct energyCalor {
0159 energyCalor(double e1 = 0, double e2 = 0, double e3 = 0) : Etot(e1), Etot2(e2), ehcal(e3) {}
0160 double Etot, Etot2, ehcal;
0161 };
0162
0163 CalibTree(const char *dupFileName,
0164 const char *rcorFileName,
0165 int truncateFlag,
0166 bool useMean,
0167 int runlo,
0168 int runhi,
0169 int phimin,
0170 int phimax,
0171 int zside,
0172 int nvxlo,
0173 int nvxhi,
0174 int sysmode,
0175 int rbx,
0176 int puCorr,
0177 int rcorForm,
0178 bool useGen,
0179 bool exclude,
0180 int higheta,
0181 TChain *tree);
0182 virtual ~CalibTree();
0183 virtual Int_t Cut(Long64_t entry);
0184 virtual Int_t GetEntry(Long64_t entry);
0185 virtual Long64_t LoadTree(Long64_t entry);
0186 virtual void Init(TChain *tree, const char *dupFileName);
0187 virtual Double_t Loop(int k,
0188 TFile *fout,
0189 bool useweight,
0190 int nMin,
0191 bool inverse,
0192 double rMin,
0193 double rMax,
0194 int ietaMax,
0195 int ietaTrack,
0196 int applyL1Cut,
0197 double l1Cut,
0198 bool last,
0199 double fraction,
0200 bool writeHisto,
0201 bool debug);
0202 virtual Bool_t Notify();
0203 virtual void Show(Long64_t entry = -1);
0204 bool goodTrack();
0205 void writeCorrFactor(const char *corrFileName, int ietaMax);
0206 bool selectPhi(unsigned int detId);
0207 std::pair<double, double> fitMean(TH1D *, int);
0208 void makeplots(double rmin, double rmax, int ietaMax, bool useWeight, double fraction, bool debug);
0209 void fitPol0(TH1D *hist, bool debug);
0210 void highEtaFactors(int ietaMax, bool debug);
0211 energyCalor energyHcal(double pmom, const Long64_t &entry, bool final);
0212
0213 TChain *fChain;
0214 Int_t fCurrent;
0215 TH1D *h_pbyE, *h_cvg;
0216 TProfile *h_Ebyp_bfr, *h_Ebyp_aftr;
0217
0218 private:
0219
0220 Int_t t_Run;
0221 Int_t t_Event;
0222 Int_t t_DataType;
0223 Int_t t_ieta;
0224 Int_t t_iphi;
0225 Double_t t_EventWeight;
0226 Int_t t_nVtx;
0227 Int_t t_nTrk;
0228 Int_t t_goodPV;
0229 Double_t t_l1pt;
0230 Double_t t_l1eta;
0231 Double_t t_l1phi;
0232 Double_t t_l3pt;
0233 Double_t t_l3eta;
0234 Double_t t_l3phi;
0235 Double_t t_p;
0236 Double_t t_pt;
0237 Double_t t_phi;
0238 Double_t t_mindR1;
0239 Double_t t_mindR2;
0240 Double_t t_eMipDR;
0241 Double_t t_eHcal;
0242 Double_t t_eHcal10;
0243 Double_t t_eHcal30;
0244 Double_t t_hmaxNearP;
0245 Double_t t_rhoh;
0246 Bool_t t_selectTk;
0247 Bool_t t_qltyFlag;
0248 Bool_t t_qltyMissFlag;
0249 Bool_t t_qltyPVFlag;
0250 Double_t t_gentrackP;
0251 std::vector<unsigned int> *t_DetIds;
0252 std::vector<double> *t_HitEnergies;
0253 std::vector<bool> *t_trgbits;
0254 std::vector<unsigned int> *t_DetIds1;
0255 std::vector<unsigned int> *t_DetIds3;
0256 std::vector<double> *t_HitEnergies1;
0257 std::vector<double> *t_HitEnergies3;
0258
0259
0260 TBranch *b_t_Run;
0261 TBranch *b_t_Event;
0262 TBranch *b_t_DataType;
0263 TBranch *b_t_ieta;
0264 TBranch *b_t_iphi;
0265 TBranch *b_t_EventWeight;
0266 TBranch *b_t_nVtx;
0267 TBranch *b_t_nTrk;
0268 TBranch *b_t_goodPV;
0269 TBranch *b_t_l1pt;
0270 TBranch *b_t_l1eta;
0271 TBranch *b_t_l1phi;
0272 TBranch *b_t_l3pt;
0273 TBranch *b_t_l3eta;
0274 TBranch *b_t_l3phi;
0275 TBranch *b_t_p;
0276 TBranch *b_t_pt;
0277 TBranch *b_t_phi;
0278 TBranch *b_t_mindR1;
0279 TBranch *b_t_mindR2;
0280 TBranch *b_t_eMipDR;
0281 TBranch *b_t_eHcal;
0282 TBranch *b_t_eHcal10;
0283 TBranch *b_t_eHcal30;
0284 TBranch *b_t_hmaxNearP;
0285 TBranch *b_t_rhoh;
0286 TBranch *b_t_selectTk;
0287 TBranch *b_t_qltyFlag;
0288 TBranch *b_t_qltyMissFlag;
0289 TBranch *b_t_qltyPVFlag;
0290 TBranch *b_t_gentrackP;
0291 TBranch *b_t_DetIds;
0292 TBranch *b_t_HitEnergies;
0293 TBranch *b_t_trgbits;
0294 TBranch *b_t_DetIds1;
0295 TBranch *b_t_DetIds3;
0296 TBranch *b_t_HitEnergies1;
0297 TBranch *b_t_HitEnergies3;
0298
0299 CalibCorr *cFactor_;
0300 CalibSelectRBX *cSelect_;
0301 const int truncateFlag_;
0302 const bool useMean_;
0303 int runlo_, runhi_;
0304 const int phimin_, phimax_;
0305 const int zside_, nvxlo_, nvxhi_;
0306 const int sysmode_, rbx_, puCorr_;
0307 int rcorForm_;
0308 const bool useGen_, exclude_;
0309 const int higheta_;
0310 bool includeRun_;
0311 double log2by18_, eHcalDelta_;
0312 std::vector<Long64_t> entries;
0313 std::vector<unsigned int> detIds;
0314 std::map<unsigned int, TH1D *> histos;
0315 std::map<unsigned int, std::pair<double, double> > Cprev;
0316 };
0317
0318 void doIt(const char *infile, const char *dup) {
0319 char outf1[100], outf2[100];
0320 double lumt(1.0), fac(0.5);
0321 for (int k = 0; k < 5; ++k) {
0322 sprintf(outf1, "%s_%d.root", infile, k);
0323 sprintf(outf2, "%s_%d.txt", infile, k);
0324 double lumi = (k == 0) ? -1 : lumt;
0325 lumt *= fac;
0326 Run(infile,
0327 "HcalIsoTrkAnalyzer",
0328 "CalibTree",
0329 outf1,
0330 outf2,
0331 dup,
0332 "",
0333 true,
0334 false,
0335 0,
0336 true,
0337 0.25,
0338 3.0,
0339 25,
0340 -1,
0341 -1,
0342 -5,
0343 1,
0344 0.5,
0345 0,
0346 30,
0347 0,
0348 false,
0349 0,
0350 99999999,
0351 1,
0352 72,
0353 0,
0354 1000,
0355 0,
0356 true,
0357 1,
0358 lumi,
0359 false,
0360 false);
0361 }
0362 }
0363
0364 void Run(const char *inFileName,
0365 const char *dirName,
0366 const char *treeName,
0367 const char *outFileName,
0368 const char *corrFileName,
0369 const char *dupFileName,
0370 const char *rcorFileName,
0371 bool useweight,
0372 bool useMean,
0373 int nMin,
0374 bool inverse,
0375 double ratMin,
0376 double ratMax,
0377 int ietaMax,
0378 int ietaTrack,
0379 int sysmode,
0380 int puCorr,
0381 int applyL1Cut,
0382 double l1Cut,
0383 int truncateFlag,
0384 int maxIter,
0385 int rcorForm,
0386 bool useGen,
0387 int runlo,
0388 int runhi,
0389 int phimin,
0390 int phimax,
0391 int zside,
0392 int nvxlo,
0393 int nvxhi,
0394 int rbx,
0395 bool exclude,
0396 int higheta,
0397 double fraction,
0398 bool writeHisto,
0399 bool debug) {
0400 char name[500];
0401 sprintf(name, "%s/%s", dirName, treeName);
0402 TChain *chain = new TChain(name);
0403 std::cout << "Create a chain for " << name << " from " << inFileName << std::endl;
0404 if (!fillChain(chain, inFileName)) {
0405 std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
0406 } else {
0407 std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
0408 Long64_t nentryTot = chain->GetEntries();
0409 Long64_t nentries = (fraction > 0.01 && fraction < 0.99) ? (Long64_t)(fraction * nentryTot) : nentryTot;
0410 static const int maxIterMax = 100;
0411 if (maxIter > maxIterMax)
0412 maxIter = maxIterMax;
0413 std::cout << "Tree " << name << " " << chain << " in directory " << dirName << " from file " << inFileName
0414 << " with nentries (tracks): " << nentries << std::endl;
0415 unsigned int k(0), kmax(maxIter);
0416 CalibTree t(dupFileName,
0417 rcorFileName,
0418 truncateFlag,
0419 useMean,
0420 runlo,
0421 runhi,
0422 phimin,
0423 phimax,
0424 zside,
0425 nvxlo,
0426 nvxhi,
0427 sysmode,
0428 rbx,
0429 puCorr,
0430 rcorForm,
0431 useGen,
0432 exclude,
0433 higheta,
0434 chain);
0435 t.h_pbyE = new TH1D("pbyE", "pbyE", 100, -1.0, 9.0);
0436 t.h_Ebyp_bfr = new TProfile("Ebyp_bfr", "Ebyp_bfr", 60, -30, 30, 0, 10);
0437 t.h_Ebyp_aftr = new TProfile("Ebyp_aftr", "Ebyp_aftr", 60, -30, 30, 0, 10);
0438 t.h_cvg = new TH1D("Cvg0", "Convergence", kmax, 0, kmax);
0439 t.h_cvg->SetMarkerStyle(7);
0440 t.h_cvg->SetMarkerSize(5.0);
0441
0442 TFile *fout = new TFile(outFileName, "RECREATE");
0443 std::cout << "Output file: " << outFileName << " opened in recreate mode" << std::endl;
0444 fout->cd();
0445
0446 double cvgs[maxIterMax], itrs[maxIterMax];
0447 for (; k <= kmax; ++k) {
0448 std::cout << "Calling Loop() " << k << "th time" << std::endl;
0449 double cvg = t.Loop(k,
0450 fout,
0451 useweight,
0452 nMin,
0453 inverse,
0454 ratMin,
0455 ratMax,
0456 ietaMax,
0457 ietaTrack,
0458 applyL1Cut,
0459 l1Cut,
0460 k == kmax,
0461 fraction,
0462 writeHisto,
0463 debug);
0464 itrs[k] = k;
0465 cvgs[k] = cvg;
0466 if (cvg < 0.00001)
0467 break;
0468 }
0469
0470 t.writeCorrFactor(corrFileName, ietaMax);
0471
0472 fout->cd();
0473 TGraph *g_cvg;
0474 g_cvg = new TGraph(k, itrs, cvgs);
0475 g_cvg->SetMarkerStyle(7);
0476 g_cvg->SetMarkerSize(5.0);
0477 g_cvg->Draw("AP");
0478 g_cvg->Write("Cvg");
0479 std::cout << "Finish looping after " << k << " iterations" << std::endl;
0480 t.makeplots(ratMin, ratMax, ietaMax, useweight, fraction, debug);
0481 fout->Close();
0482 }
0483 }
0484
0485 CalibTree::CalibTree(const char *dupFileName,
0486 const char *rcorFileName,
0487 int flag,
0488 bool useMean,
0489 int runlo,
0490 int runhi,
0491 int phimin,
0492 int phimax,
0493 int zside,
0494 int nvxlo,
0495 int nvxhi,
0496 int mode,
0497 int rbx,
0498 int pu,
0499 int rForm,
0500 bool gen,
0501 bool excl,
0502 int heta,
0503 TChain *tree)
0504 : fChain(nullptr),
0505 cFactor_(nullptr),
0506 cSelect_(nullptr),
0507 truncateFlag_(flag),
0508 useMean_(useMean),
0509 runlo_(runlo),
0510 runhi_(runhi),
0511 phimin_(phimin),
0512 phimax_(phimax),
0513 zside_(zside),
0514 nvxlo_(nvxlo),
0515 nvxhi_(nvxhi),
0516 sysmode_(mode),
0517 rbx_(rbx),
0518 puCorr_(pu),
0519 rcorForm_(rForm),
0520 useGen_(gen),
0521 exclude_(excl),
0522 higheta_(heta),
0523 includeRun_(true) {
0524 if (runlo_ < 0 || runhi_ < 0) {
0525 runlo_ = std::abs(runlo_);
0526 runhi_ = std::abs(runhi_);
0527 includeRun_ = false;
0528 }
0529 log2by18_ = std::log(2.5) / 18.0;
0530 eHcalDelta_ = 0;
0531 std::cout << "Initialize CalibTree with TruncateFlag " << truncateFlag_ << " UseMean " << useMean_ << " Run Range "
0532 << runlo_ << ":" << runhi_ << " Phi Range " << phimin_ << ":" << phimax_ << ":" << zside_
0533 << " Vertex Range " << nvxlo_ << ":" << nvxhi_ << " Mode " << sysmode_ << " PU " << puCorr_
0534 << " rcorFormat " << rcorForm_ << " Gen " << useGen_ << " High Eta " << higheta_ << std::endl;
0535 std::cout << "Duplicate events read from " << dupFileName << " RadDam Corrections read from " << rcorFileName
0536 << " Treat RBX " << rbx_ << " with exclusion mode " << exclude_ << std::endl;
0537 Init(tree, dupFileName);
0538 if (std::string(rcorFileName) != "") {
0539 cFactor_ = new CalibCorr(rcorFileName, rcorForm_, false);
0540 } else {
0541 rcorForm_ = 0;
0542 }
0543 if (rbx != 0)
0544 cSelect_ = new CalibSelectRBX(rbx);
0545 }
0546
0547 CalibTree::~CalibTree() {
0548 delete cFactor_;
0549 delete cSelect_;
0550 if (!fChain)
0551 return;
0552 delete fChain->GetCurrentFile();
0553 }
0554
0555 Int_t CalibTree::GetEntry(Long64_t entry) {
0556
0557 if (!fChain)
0558 return 0;
0559 return fChain->GetEntry(entry);
0560 }
0561
0562 Long64_t CalibTree::LoadTree(Long64_t entry) {
0563
0564 if (!fChain)
0565 return -5;
0566 Long64_t centry = fChain->LoadTree(entry);
0567 if (centry < 0)
0568 return centry;
0569 if (fChain->GetTreeNumber() != fCurrent) {
0570 fCurrent = fChain->GetTreeNumber();
0571 Notify();
0572 }
0573 return centry;
0574 }
0575
0576 void CalibTree::Init(TChain *tree, const char *dupFileName) {
0577
0578
0579
0580
0581
0582
0583
0584
0585
0586 t_DetIds = 0;
0587 t_HitEnergies = 0;
0588 t_trgbits = 0;
0589 t_DetIds1 = 0;
0590 t_DetIds3 = 0;
0591 t_HitEnergies1 = 0;
0592 t_HitEnergies3 = 0;
0593
0594 fChain = tree;
0595 if (!tree)
0596 return;
0597 fCurrent = -1;
0598 fChain->SetMakeClass(1);
0599
0600 fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
0601 fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
0602 fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
0603 fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
0604 fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
0605 fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
0606 fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
0607 fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
0608 fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
0609 fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
0610 fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
0611 fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
0612 fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
0613 fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
0614 fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
0615 fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
0616 fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
0617 fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
0618 fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
0619 fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
0620 fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
0621 fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
0622 fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
0623 fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
0624 fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
0625 fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
0626 fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
0627 fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
0628 fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
0629 fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
0630 fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
0631 fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
0632 fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
0633 fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
0634 fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
0635 fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
0636 fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
0637 fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
0638 Notify();
0639
0640 if (std::string(dupFileName) != "") {
0641 ifstream infil1(dupFileName);
0642 if (!infil1.is_open()) {
0643 std::cout << "Cannot open " << dupFileName << std::endl;
0644 } else {
0645 while (1) {
0646 Long64_t jentry;
0647 infil1 >> jentry;
0648 if (!infil1.good())
0649 break;
0650 entries.push_back(jentry);
0651 }
0652 infil1.close();
0653 std::cout << "Reads a list of " << entries.size() << " events from " << dupFileName << std::endl;
0654 }
0655 } else {
0656 std::cout << "No duplicate events in the input file" << std::endl;
0657 }
0658 }
0659
0660 Bool_t CalibTree::Notify() {
0661
0662
0663
0664
0665
0666
0667 return kTRUE;
0668 }
0669
0670 void CalibTree::Show(Long64_t entry) {
0671
0672
0673 if (!fChain)
0674 return;
0675 fChain->Show(entry);
0676 }
0677
0678 Int_t CalibTree::Cut(Long64_t) {
0679
0680
0681
0682 return 1;
0683 }
0684
0685 Double_t CalibTree::Loop(int loop,
0686 TFile *fout,
0687 bool useweight,
0688 int nMin,
0689 bool inverse,
0690 double rmin,
0691 double rmax,
0692 int ietaMax,
0693 int ietaTrack,
0694 int applyL1Cut,
0695 double l1Cut,
0696 bool last,
0697 double fraction,
0698 bool writeHisto,
0699 bool debug) {
0700 if (fChain == 0)
0701 return 0;
0702 Long64_t nbytes(0), nb(0), kprint(0);
0703 Long64_t nentryTot = fChain->GetEntriesFast();
0704 Long64_t nentries = (fraction > 0.01 && fraction < 0.99) ? (Long64_t)(fraction * nentryTot) : nentryTot;
0705 if (detIds.size() == 0) {
0706 for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0707 Long64_t ientry = LoadTree(jentry);
0708 if (ientry < 0)
0709 break;
0710 nb = fChain->GetEntry(jentry);
0711 nbytes += nb;
0712 if (jentry % 1000000 == 0)
0713 std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
0714
0715 bool selRun = (includeRun_ ? ((t_Run >= runlo_) && (t_Run <= runhi_)) : ((t_Run < runlo_) || (t_Run > runhi_)));
0716 bool selTrack = ((ietaTrack <= 0) || (abs(t_ieta) <= ietaTrack));
0717 if (selRun && (t_nVtx >= nvxlo_) && (t_nVtx <= nvxhi_) && selTrack) {
0718 bool isItRBX(false);
0719 if (cSelect_ != nullptr) {
0720 bool temp = cSelect_->isItRBX(t_DetIds);
0721 if (exclude_)
0722 isItRBX = temp;
0723 else
0724 isItRBX = !(temp);
0725 }
0726 ++kprint;
0727 if (!(isItRBX)) {
0728 for (unsigned int idet = 0; idet < (*t_DetIds).size(); idet++) {
0729 if (selectPhi((*t_DetIds)[idet])) {
0730 unsigned int detid = truncateId((*t_DetIds)[idet], truncateFlag_, debug);
0731 if (debug && (kprint <= 10)) {
0732 std::cout << "DetId[" << idet << "] Original " << std::hex << (*t_DetIds)[idet] << " truncated "
0733 << detid << std::dec;
0734 }
0735 if (std::find(detIds.begin(), detIds.end(), detid) == detIds.end()) {
0736 detIds.push_back(detid);
0737 if (debug && (kprint <= 10))
0738 std::cout << " new";
0739 }
0740 if (debug && (kprint <= 10))
0741 std::cout << std::endl;
0742 }
0743 }
0744
0745 if (t_DetIds3 != 0) {
0746 for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
0747 if (selectPhi((*t_DetIds3)[idet])) {
0748 unsigned int detid = truncateId((*t_DetIds3)[idet], truncateFlag_, debug);
0749 if (std::find(detIds.begin(), detIds.end(), detid) == detIds.end()) {
0750 detIds.push_back(detid);
0751 }
0752 }
0753 }
0754 }
0755 }
0756 }
0757 }
0758 }
0759 if (debug) {
0760 std::cout << "Total of " << detIds.size() << " detIds and " << histos.size() << " histos found" << std::endl;
0761
0762 for (unsigned int k = 0; k < detIds.size(); ++k) {
0763 int subdet, depth, zside, ieta, iphi;
0764 unpackDetId(detIds[k], subdet, zside, ieta, iphi, depth);
0765 std::cout << "DetId[" << k << "] " << subdet << ":" << zside * ieta << ":" << depth << ":" << iphi << " "
0766 << std::hex << detIds[k] << std::dec << std::endl;
0767 }
0768 }
0769 unsigned int k(0);
0770 for (std::map<unsigned int, TH1D *>::const_iterator itr = histos.begin(); itr != histos.end(); ++itr, ++k) {
0771 if (debug) {
0772 std::cout << "histos[" << k << "] " << std::hex << itr->first << std::dec << " " << itr->second;
0773 if (itr->second != 0)
0774 std::cout << " " << itr->second->GetTitle();
0775 std::cout << std::endl;
0776 }
0777 if (itr->second != 0)
0778 itr->second->Delete();
0779 }
0780
0781 for (unsigned int k = 0; k < detIds.size(); ++k) {
0782 char name[20], title[100];
0783 sprintf(name, "Hist%d_%d", detIds[k], loop);
0784 int subdet, depth, zside, ieta, iphi;
0785 unpackDetId(detIds[k], subdet, zside, ieta, iphi, depth);
0786 sprintf(title, "Correction for Subdet %d #eta %d depth %d (Loop %d)", subdet, zside * ieta, depth, loop);
0787 TH1D *hist = new TH1D(name, title, 100, 0.0, 5.0);
0788 hist->Sumw2();
0789 if (debug)
0790 std::cout << "Book Histo " << k << " " << title << std::endl;
0791 histos[detIds[k]] = hist;
0792 }
0793 std::cout << "Total of " << detIds.size() << " detIds and " << histos.size() << " found in " << nentries << std::endl;
0794
0795 nbytes = nb = 0;
0796 std::map<unsigned int, myEntry> SumW;
0797 std::map<unsigned int, double> nTrks;
0798
0799 int ntkgood(0);
0800 for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0801 Long64_t ientry = LoadTree(jentry);
0802 if (ientry < 0)
0803 break;
0804 nb = fChain->GetEntry(jentry);
0805 nbytes += nb;
0806 if (jentry % 1000000 == 0)
0807 std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
0808 if (std::find(entries.begin(), entries.end(), jentry) != entries.end())
0809 continue;
0810 bool selRun = (includeRun_ ? ((t_Run >= runlo_) && (t_Run <= runhi_)) : ((t_Run < runlo_) || (t_Run > runhi_)));
0811 if (!selRun)
0812 continue;
0813 if ((t_nVtx < nvxlo_) || (t_nVtx > nvxhi_))
0814 continue;
0815 if (cSelect_ != nullptr) {
0816 if (exclude_) {
0817 if (cSelect_->isItRBX(t_DetIds))
0818 continue;
0819 } else {
0820 if (!(cSelect_->isItRBX(t_ieta, t_iphi)))
0821 continue;
0822 }
0823 }
0824 bool selTrack = ((ietaTrack <= 0) || (abs(t_ieta) <= ietaTrack));
0825 if (!selTrack)
0826 continue;
0827 if ((rcorForm_ == 3) && (cFactor_ != nullptr) && (cFactor_->absent(ientry)))
0828 continue;
0829
0830 if (debug) {
0831 std::cout << "***Entry (Track) Number : " << ientry << std::endl;
0832 std::cout << "p/eHCal/eMipDR/nDets : " << t_p << "/" << t_eHcal << "/" << t_eMipDR << "/" << (*t_DetIds).size()
0833 << std::endl;
0834 }
0835 double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
0836 if (goodTrack()) {
0837 ++ntkgood;
0838 CalibTree::energyCalor en = energyHcal(pmom, jentry, true);
0839 double evWt = (useweight) ? t_EventWeight : 1.0;
0840 if (en.ehcal > 0.001) {
0841 double pufac = (en.Etot > 0) ? (en.ehcal / en.Etot) : 1.0;
0842 double ratio = en.ehcal / (pmom - t_eMipDR);
0843 if (debug)
0844 std::cout << " Weights " << evWt << ":" << pufac << " Energy " << en.Etot2 << ":" << en.Etot << ":" << pmom
0845 << ":" << t_eMipDR << ":" << t_eHcal << ":" << en.ehcal << " ratio " << ratio << std::endl;
0846 if (loop == 0) {
0847 h_pbyE->Fill(ratio, evWt);
0848 h_Ebyp_bfr->Fill(t_ieta, ratio, evWt);
0849 }
0850 if (last) {
0851 h_Ebyp_aftr->Fill(t_ieta, ratio, evWt);
0852 }
0853 bool l1c(true);
0854 if (applyL1Cut != 0)
0855 l1c = ((t_mindR1 >= l1Cut) || ((applyL1Cut == 1) && (t_DataType == 1)));
0856 if ((rmin >= 0 && ratio > rmin) && (rmax >= 0 && ratio < rmax) && l1c) {
0857 for (unsigned int idet = 0; idet < (*t_DetIds).size(); idet++) {
0858 if (selectPhi((*t_DetIds)[idet])) {
0859 unsigned int id = (*t_DetIds)[idet];
0860 unsigned int detid = truncateId(id, truncateFlag_, false);
0861 double hitEn = 0.0;
0862 if (debug) {
0863 std::cout << "idet " << idet << " detid/hitenergy : " << std::hex << (*t_DetIds)[idet] << ":" << detid
0864 << "/" << (*t_HitEnergies)[idet] << std::endl;
0865 }
0866 if (Cprev.find(detid) != Cprev.end())
0867 hitEn = Cprev[detid].first * (*t_HitEnergies)[idet];
0868 else
0869 hitEn = (*t_HitEnergies)[idet];
0870 if (cFactor_)
0871 hitEn *= cFactor_->getCorr(t_Run, id);
0872 double Wi = evWt * hitEn / en.Etot;
0873 double Fac = (inverse) ? (en.ehcal / (pmom - t_eMipDR)) : ((pmom - t_eMipDR) / en.ehcal);
0874 double Fac2 = Wi * Fac * Fac;
0875 TH1D *hist(0);
0876 std::map<unsigned int, TH1D *>::const_iterator itr = histos.find(detid);
0877 if (itr != histos.end())
0878 hist = itr->second;
0879 if (debug) {
0880 std::cout << "Det Id " << std::hex << detid << std::dec << " " << hist << std::endl;
0881 }
0882 if (hist != 0)
0883 hist->Fill(Fac, Wi);
0884 Fac *= Wi;
0885 if (SumW.find(detid) != SumW.end()) {
0886 Wi += SumW[detid].fact0;
0887 Fac += SumW[detid].fact1;
0888 Fac2 += SumW[detid].fact2;
0889 int kount = SumW[detid].kount + 1;
0890 SumW[detid] = myEntry(kount, Wi, Fac, Fac2);
0891 nTrks[detid] += evWt;
0892 } else {
0893 SumW.insert(std::pair<unsigned int, myEntry>(detid, myEntry(1, Wi, Fac, Fac2)));
0894 nTrks.insert(std::pair<unsigned int, unsigned int>(detid, evWt));
0895 }
0896 }
0897 }
0898 }
0899 }
0900 }
0901 }
0902 if (debug)
0903 std::cout << "# of Good Tracks " << ntkgood << " out of " << nentries << std::endl;
0904 if (loop == 0) {
0905 h_pbyE->Write("h_pbyE");
0906 h_Ebyp_bfr->Write("h_Ebyp_bfr");
0907 }
0908 if (last) {
0909 h_Ebyp_aftr->Write("h_Ebyp_aftr");
0910 }
0911
0912 std::map<unsigned int, std::pair<double, double> > cfactors;
0913 unsigned int kount(0), kountus(0);
0914 double sumfactor(0);
0915 for (std::map<unsigned int, TH1D *>::const_iterator itr = histos.begin(); itr != histos.end(); ++itr) {
0916 if (writeHisto) {
0917 std::pair<double, double> result_write = fitMean(itr->second, 0);
0918 (itr->second)->Write();
0919 }
0920
0921 int subdet, depth, zside, ieta, iphi;
0922 unpackDetId(itr->first, subdet, zside, ieta, iphi, depth);
0923 if (debug) {
0924 std::cout << "DETID :" << subdet << " IETA :" << ieta << " HIST ENTRIES :" << (itr->second)->GetEntries()
0925 << std::endl;
0926 }
0927 }
0928
0929 if (debug)
0930 std::cout << "Histos with " << histos.size() << " entries\n";
0931 for (std::map<unsigned int, TH1D *>::const_iterator itr = histos.begin(); itr != histos.end(); ++itr, ++kount) {
0932 std::pair<double, double> result = fitMean(itr->second, 0);
0933 double factor = (inverse) ? (2. - result.first) : result.first;
0934 if (debug) {
0935 int subdet, depth, zside, ieta, iphi;
0936 unpackDetId(itr->first, subdet, zside, ieta, iphi, depth);
0937 std::cout << "DetId[" << kount << "] " << subdet << ":" << zside * ieta << ":" << depth << " Factor " << factor
0938 << " +- " << result.second << std::endl;
0939 }
0940 if (!useMean_) {
0941 cfactors[itr->first] = std::pair<double, double>(factor, result.second);
0942 if (itr->second->GetEntries() > nMin) {
0943 kountus++;
0944 if (factor > 1)
0945 sumfactor += (1 - 1 / factor);
0946 else
0947 sumfactor += (1 - factor);
0948 }
0949 }
0950 }
0951
0952 if (debug)
0953 std::cout << "SumW with " << SumW.size() << " entries\n";
0954 std::map<unsigned int, myEntry>::const_iterator SumWItr = SumW.begin();
0955 for (; SumWItr != SumW.end(); SumWItr++) {
0956 unsigned int detid = SumWItr->first;
0957 int subdet, depth, zside, ieta, iphi;
0958 unpackDetId(detid, subdet, zside, ieta, iphi, depth);
0959 if (debug) {
0960 std::cout << "Detid|kount|SumWi|SumFac|myId : " << subdet << ":" << zside * ieta << ":" << depth << " | "
0961 << (SumWItr->second).kount << " | " << (SumWItr->second).fact0 << "|" << (SumWItr->second).fact1 << "|"
0962 << (SumWItr->second).fact2 << std::endl;
0963 }
0964 double factor = (SumWItr->second).fact1 / (SumWItr->second).fact0;
0965 double dfac1 = ((SumWItr->second).fact2 / (SumWItr->second).fact0 - factor * factor);
0966 if (dfac1 < 0)
0967 dfac1 = 0;
0968 double dfac = sqrt(dfac1 / (SumWItr->second).kount);
0969 if (debug) {
0970 std::cout << "Factor " << factor << " " << dfac1 << " " << dfac << std::endl;
0971 }
0972 if (inverse)
0973 factor = 2. - factor;
0974 if (useMean_) {
0975 cfactors[detid] = std::pair<double, double>(factor, dfac);
0976 if ((SumWItr->second).kount > nMin) {
0977 kountus++;
0978 if (factor > 1)
0979 sumfactor += (1 - 1 / factor);
0980 else
0981 sumfactor += (1 - factor);
0982 }
0983 }
0984 }
0985
0986 static const int maxch = 500;
0987 double dets[maxch], cfacs[maxch], wfacs[maxch], myId[maxch], nTrk[maxch];
0988 std::cout << "cafctors: " << cfactors.size() << ":" << maxch << std::endl;
0989 kount = 0;
0990 std::map<unsigned int, std::pair<double, double> >::const_iterator itr = cfactors.begin();
0991 for (; itr != cfactors.end(); ++itr, ++kount) {
0992 unsigned int detid = itr->first;
0993 int subdet, depth, zside, ieta, iphi;
0994 unpackDetId(detid, subdet, zside, ieta, iphi, depth);
0995 double id = ieta * zside + 0.25 * (depth - 1);
0996 double factor = (itr->second).first;
0997 double dfac = (itr->second).second;
0998 if (ieta > ietaMax) {
0999 factor = 1;
1000 dfac = 0;
1001 }
1002 std::pair<double, double> cfac(factor, dfac);
1003 if (Cprev.find(detid) != Cprev.end()) {
1004 dfac /= factor;
1005 factor *= Cprev[detid].first;
1006 dfac *= factor;
1007 Cprev[detid] = std::pair<double, double>(factor, dfac);
1008 cfacs[kount] = factor;
1009 } else {
1010 Cprev[detid] = std::pair<double, double>(factor, dfac);
1011 cfacs[kount] = factor;
1012 }
1013 wfacs[kount] = factor;
1014 dets[kount] = detid;
1015 myId[kount] = id;
1016 nTrk[kount] = nTrks[detid];
1017 }
1018 if (higheta_ > 0)
1019 highEtaFactors(ietaMax, debug);
1020
1021 std::cout << kountus << " detids out of " << kount << " have tracks > " << nMin << std::endl;
1022
1023 char fname[50];
1024 fout->cd();
1025 TGraph *g_fac1 = new TGraph(kount, dets, cfacs);
1026 sprintf(fname, "Cfacs%d", loop);
1027 g_fac1->SetMarkerStyle(7);
1028 g_fac1->SetMarkerSize(5.0);
1029 g_fac1->Draw("AP");
1030 g_fac1->Write(fname);
1031 TGraph *g_fac2 = new TGraph(kount, dets, wfacs);
1032 sprintf(fname, "Wfacs%d", loop);
1033 g_fac2->SetMarkerStyle(7);
1034 g_fac2->SetMarkerSize(5.0);
1035 g_fac2->Draw("AP");
1036 g_fac2->Write(fname);
1037 TGraph *g_fac3 = new TGraph(kount, myId, cfacs);
1038 sprintf(fname, "CfacsVsMyId%d", loop);
1039 g_fac3->SetMarkerStyle(7);
1040 g_fac3->SetMarkerSize(5.0);
1041 g_fac3->Draw("AP");
1042 g_fac3->Write(fname);
1043 TGraph *g_fac4 = new TGraph(kount, myId, wfacs);
1044 sprintf(fname, "WfacsVsMyId%d", loop);
1045 g_fac4->SetMarkerStyle(7);
1046 g_fac4->SetMarkerSize(5.0);
1047 g_fac4->Draw("AP");
1048 g_fac4->Write(fname);
1049 TGraph *g_nTrk = new TGraph(kount, myId, nTrk);
1050 sprintf(fname, "nTrk");
1051 if (loop == 0) {
1052 g_nTrk->SetMarkerStyle(7);
1053 g_nTrk->SetMarkerSize(5.0);
1054 g_nTrk->Draw("AP");
1055 g_nTrk->Write(fname);
1056 }
1057 std::cout << "The new factors are :" << std::endl;
1058 std::map<unsigned int, std::pair<double, double> >::const_iterator CprevItr = Cprev.begin();
1059 unsigned int indx(0);
1060 for (; CprevItr != Cprev.end(); CprevItr++, indx++) {
1061 unsigned int detid = CprevItr->first;
1062 int subdet, depth, zside, ieta, iphi;
1063 unpackDetId(detid, subdet, zside, ieta, iphi, depth);
1064 std::cout << "DetId[" << indx << "] " << std::hex << detid << std::dec << "(" << ieta * zside << "," << depth
1065 << ") (nTrks:" << nTrks[detid] << ") : " << CprevItr->second.first << " +- " << CprevItr->second.second
1066 << std::endl;
1067 }
1068 double mean = (kountus > 0) ? (sumfactor / kountus) : 0;
1069 std::cout << "Mean deviation " << mean << " from 1 for " << kountus << " DetIds" << std::endl;
1070 h_cvg->SetBinContent(loop + 1, mean);
1071 if (last)
1072 h_cvg->Write("Cvg0");
1073 return mean;
1074 }
1075
1076 bool CalibTree::goodTrack() {
1077 bool ok(true);
1078 double cut(2.0);
1079 double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
1080 if (sysmode_ == 1) {
1081 ok =
1082 ((t_qltyFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) && (pmom > 40.0) && (pmom < 60.0));
1083 } else if (sysmode_ == 2) {
1084 ok = ((t_qltyFlag) && (t_qltyPVFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) &&
1085 (pmom > 40.0) && (pmom < 60.0));
1086 } else if (sysmode_ == 3) {
1087 ok =
1088 ((t_selectTk) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) && (pmom > 40.0) && (pmom < 60.0));
1089 } else if (sysmode_ == 4) {
1090 ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < 0.0) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) &&
1091 (pmom > 40.0) && (pmom < 60.0));
1092 } else if (sysmode_ == 5) {
1093 ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 0.5) && (t_mindR1 > 1.0) &&
1094 (pmom > 40.0) && (pmom < 60.0));
1095 } else if (sysmode_ == 6) {
1096 ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 2.0) && (t_mindR1 > 1.0) &&
1097 (pmom > 40.0) && (pmom < 60.0));
1098 } else if (sysmode_ == 7) {
1099 ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 0.5) &&
1100 (pmom > 40.0) && (pmom < 60.0));
1101 } else {
1102 if (sysmode_ < 0) {
1103 double eta = (t_ieta > 0) ? t_ieta : -t_ieta;
1104 if (sysmode_ == -2)
1105 cut = 8.0 * exp(eta * log2by18_);
1106 else
1107 cut = 10.0;
1108 }
1109 ok = ((t_selectTk) && (t_qltyMissFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0) && (t_mindR1 > 1.0) &&
1110 (pmom > 40.0) && (pmom < 60.0));
1111 }
1112 return ok;
1113 }
1114
1115 void CalibTree::writeCorrFactor(const char *corrFileName, int ietaMax) {
1116 ofstream myfile;
1117 myfile.open(corrFileName);
1118 if (!myfile.is_open()) {
1119 std::cout << "** ERROR: Can't open '" << corrFileName << std::endl;
1120 } else {
1121 myfile << "#" << std::setprecision(4) << std::setw(10) << "detId" << std::setw(10) << "ieta" << std::setw(10)
1122 << "depth" << std::setw(15) << "corrFactor" << std::endl;
1123 std::map<unsigned int, std::pair<double, double> >::const_iterator itr;
1124 for (itr = Cprev.begin(); itr != Cprev.end(); ++itr) {
1125 unsigned int detId = itr->first;
1126 int subdet, depth, zside, ieta, iphi;
1127 unpackDetId(detId, subdet, zside, ieta, iphi, depth);
1128 if (ieta <= ietaMax) {
1129 myfile << std::setw(10) << std::hex << detId << std::setw(10) << std::dec << zside * ieta << std::setw(10)
1130 << depth << std::setw(10) << itr->second.first << " " << std::setw(10) << itr->second.second
1131 << std::endl;
1132 std::cout << itr->second.first << ",";
1133 }
1134 }
1135 myfile.close();
1136 std::cout << std::endl;
1137 }
1138 }
1139
1140 bool CalibTree::selectPhi(unsigned int detId) {
1141 bool flag(true);
1142
1143 if ((phimin_ > 1) || (phimax_ < 72) || (zside_ != 0)) {
1144 int subdet, depth, zside, ieta, iphi;
1145 unpackDetId(detId, subdet, zside, ieta, iphi, depth);
1146 if (phimin_ > 1 || phimax_ < 72) {
1147 if ((iphi < phimin_) || (iphi > phimax_))
1148 flag = false;
1149 }
1150 if (zside_ != 0) {
1151 if (zside != zside_)
1152 flag = false;
1153 }
1154 }
1155 return flag;
1156 }
1157
1158 std::pair<double, double> CalibTree::fitMean(TH1D *hist, int mode) {
1159 std::pair<double, double> results = std::pair<double, double>(1, 0);
1160 if (hist != 0) {
1161 double mean = hist->GetMean(), rms = hist->GetRMS();
1162 double LowEdge(0.7), HighEdge(1.3);
1163 char option[20];
1164 if (mode == 1) {
1165 LowEdge = mean - 1.5 * rms;
1166 HighEdge = mean + 1.5 * rms;
1167 int nbin = hist->GetNbinsX();
1168 if (LowEdge < hist->GetBinLowEdge(1))
1169 LowEdge = hist->GetBinLowEdge(1);
1170 if (HighEdge > hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin))
1171 HighEdge = hist->GetBinLowEdge(nbin) + hist->GetBinWidth(nbin);
1172 }
1173 if (hist->GetEntries() > 100)
1174 sprintf(option, "+QRLS");
1175 else
1176 sprintf(option, "+QRWLS");
1177 double value(mean);
1178 double error = rms / sqrt(hist->GetEntries());
1179 if (hist->GetEntries() > 20) {
1180 TFitResultPtr Fit = hist->Fit("gaus", option, "", LowEdge, HighEdge);
1181 value = Fit->Value(1);
1182 error = Fit->FitResult::Error(1);
1183
1184
1185
1186
1187
1188
1189
1190 }
1191 results = std::pair<double, double>(value, error);
1192 }
1193 return results;
1194 }
1195
1196 void CalibTree::makeplots(double rmin, double rmax, int ietaMax, bool useweight, double fraction, bool debug) {
1197 if (fChain == 0)
1198 return;
1199 Long64_t nentryTot = fChain->GetEntriesFast();
1200 Long64_t nentries = (fraction > 0.01 && fraction < 0.99) ? (Long64_t)(fraction * nentryTot) : nentryTot;
1201
1202
1203 std::map<int, std::pair<TH1D *, TH1D *> > histos;
1204 for (int ieta = -ietaMax; ieta <= ietaMax; ++ieta) {
1205 char name[20], title[100];
1206 sprintf(name, "begin%d", ieta);
1207 if (ieta == 0)
1208 sprintf(title, "Ratio at start");
1209 else
1210 sprintf(title, "Ratio at start for i#eta=%d", ieta);
1211 TH1D *h1 = new TH1D(name, title, 50, rmin, rmax);
1212 h1->Sumw2();
1213 sprintf(name, "end%d", ieta);
1214 if (ieta == 0)
1215 sprintf(title, "Ratio at the end");
1216 else
1217 sprintf(title, "Ratio at the end for i#eta=%d", ieta);
1218 TH1D *h2 = new TH1D(name, title, 50, rmin, rmax);
1219 h2->Sumw2();
1220 histos[ieta] = std::pair<TH1D *, TH1D *>(h1, h2);
1221 }
1222
1223 Long64_t nbytes(0), nb(0);
1224 for (Long64_t jentry = 0; jentry < nentries; jentry++) {
1225 Long64_t ientry = LoadTree(jentry);
1226 if (ientry < 0)
1227 break;
1228 if (std::find(entries.begin(), entries.end(), jentry) != entries.end())
1229 break;
1230 nb = fChain->GetEntry(jentry);
1231 nbytes += nb;
1232 if (goodTrack()) {
1233 double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
1234 CalibTree::energyCalor en1 = energyHcal(pmom, jentry, false);
1235 CalibTree::energyCalor en2 = energyHcal(pmom, jentry, true);
1236 if ((en1.ehcal > 0.001) && (en2.ehcal > 0.001)) {
1237 double evWt = (useweight) ? t_EventWeight : 1.0;
1238 double ratioi = en1.ehcal / (pmom - t_eMipDR);
1239 double ratiof = en2.ehcal / (pmom - t_eMipDR);
1240 if (t_ieta >= -ietaMax && t_ieta <= ietaMax && t_ieta != 0) {
1241 if (ratioi >= rmin && ratioi <= rmax) {
1242 histos[0].first->Fill(ratioi, evWt);
1243 histos[t_ieta].first->Fill(ratioi, evWt);
1244 }
1245 if (ratiof >= rmin && ratiof <= rmax) {
1246 histos[0].second->Fill(ratiof, evWt);
1247 histos[t_ieta].second->Fill(ratiof, evWt);
1248 }
1249 }
1250 }
1251 }
1252 }
1253
1254
1255 TH1D *hbef1 = new TH1D("Eta1Bf", "Mean vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1256 TH1D *hbef2 = new TH1D("Eta2Bf", "Median vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1257 TH1D *haft1 = new TH1D("Eta1Af", "Mean vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1258 TH1D *haft2 = new TH1D("Eta2Af", "Median vs i#eta", 2 * ietaMax, -ietaMax, ietaMax);
1259 for (int ieta = -ietaMax; ieta <= ietaMax; ++ieta) {
1260 int bin = (ieta < 0) ? (ieta + ietaMax + 1) : (ieta + ietaMax);
1261 TH1D *h1 = histos[ieta].first;
1262 double mean1 = h1->GetMean();
1263 double err1 = h1->GetMeanError();
1264 std::pair<double, double> fit1 = fitMean(h1, 1);
1265 if (debug) {
1266 std::cout << ieta << " " << h1->GetName() << " " << mean1 << " +- " << err1 << " and " << fit1.first << " +- "
1267 << fit1.second << std::endl;
1268 }
1269 if (ieta != 0) {
1270 hbef1->SetBinContent(bin, mean1);
1271 hbef1->SetBinError(bin, err1);
1272 hbef2->SetBinContent(bin, fit1.first);
1273 hbef2->SetBinError(bin, fit1.second);
1274 }
1275 h1->Write();
1276 TH1D *h2 = histos[ieta].second;
1277 double mean2 = h2->GetMean();
1278 double err2 = h2->GetMeanError();
1279 std::pair<double, double> fit2 = fitMean(h2, 1);
1280 if (debug) {
1281 std::cout << ieta << " " << h2->GetName() << " " << mean2 << " +- " << err2 << " and " << fit2.first << " +- "
1282 << fit2.second << std::endl;
1283 }
1284 if (ieta != 0) {
1285 haft1->SetBinContent(bin, mean2);
1286 haft1->SetBinError(bin, err2);
1287 haft2->SetBinContent(bin, fit2.first);
1288 haft2->SetBinError(bin, fit2.second);
1289 }
1290 h2->Write();
1291 }
1292 fitPol0(hbef1, debug);
1293 fitPol0(hbef2, debug);
1294 fitPol0(haft1, debug);
1295 fitPol0(haft2, debug);
1296 }
1297
1298 void CalibTree::fitPol0(TH1D *hist, bool debug) {
1299 hist->GetXaxis()->SetTitle("i#eta");
1300 hist->GetYaxis()->SetTitle("<E_{HCAL}/(p-E_{ECAL})>");
1301 hist->GetYaxis()->SetRangeUser(0.4, 1.6);
1302 TFitResultPtr Fit = hist->Fit("pol0", "+QRWLS");
1303 if (debug) {
1304 std::cout << "Fit to Pol0 to " << hist->GetTitle() << ": " << Fit->Value(0) << " +- " << Fit->FitResult::Error(0)
1305 << std::endl;
1306 }
1307 hist->Write();
1308 }
1309
1310 void CalibTree::highEtaFactors(int ietaMax, bool debug) {
1311 std::map<unsigned int, std::pair<double, double> >::const_iterator itr;
1312 std::pair<double, double> cfacp, cfacn;
1313 cfacp = cfacn = std::pair<double, double>(1.0, 0.0);
1314 for (itr = Cprev.begin(); itr != Cprev.end(); ++itr) {
1315 unsigned int detid = itr->first;
1316 int subdet, depth, zside, ieta, iphi;
1317 unpackDetId(detid, subdet, zside, ieta, iphi, depth);
1318 if ((ieta == ietaMax) && (depth == 1)) {
1319 if (zside > 0)
1320 cfacp = itr->second;
1321 else
1322 cfacn = itr->second;
1323 }
1324 }
1325 if (debug) {
1326 std::cout << "Correction factor for (" << ietaMax << ",1) = " << cfacp.first << ":" << cfacp.second << " and ("
1327 << -ietaMax << ",1) = " << cfacn.first << ":" << cfacn.second << std::endl;
1328 }
1329 for (itr = Cprev.begin(); itr != Cprev.end(); ++itr) {
1330 unsigned int detid = itr->first;
1331 int subdet, depth, zside, ieta, iphi;
1332 unpackDetId(detid, subdet, zside, ieta, iphi, depth);
1333 if (ieta > ietaMax) {
1334 Cprev[detid] = (zside > 0) ? cfacp : cfacn;
1335 if (debug) {
1336 std::cout << "Set correction factor for (" << zside * ieta << "," << depth << ") = " << Cprev[detid].first
1337 << ":" << Cprev[detid].second << std::endl;
1338 }
1339 }
1340 }
1341 }
1342
1343 CalibTree::energyCalor CalibTree::energyHcal(double pmom, const Long64_t &entry, bool final) {
1344 double etot = t_eHcal;
1345 double etot2 = t_eHcal;
1346 double ediff = (t_eHcal30 - t_eHcal10);
1347 if (final) {
1348 etot = etot2 = 0;
1349 for (unsigned int idet = 0; idet < (*t_DetIds).size(); idet++) {
1350 if (selectPhi((*t_DetIds)[idet])) {
1351 unsigned int id = (*t_DetIds)[idet];
1352 double hitEn(0);
1353 unsigned int detid = truncateId(id, truncateFlag_, false);
1354 if (Cprev.find(detid) != Cprev.end())
1355 hitEn = Cprev[detid].first * (*t_HitEnergies)[idet];
1356 else
1357 hitEn = (*t_HitEnergies)[idet];
1358 if (cFactor_)
1359 hitEn *= cFactor_->getCorr(t_Run, id);
1360 etot += hitEn;
1361 etot2 += ((*t_HitEnergies)[idet]);
1362 }
1363 }
1364
1365 double etot1(0), etot3(0);
1366 if (t_DetIds1 != 0 && t_DetIds3 != 0) {
1367 for (unsigned int idet = 0; idet < (*t_DetIds1).size(); idet++) {
1368 if (selectPhi((*t_DetIds1)[idet])) {
1369 unsigned int id = (*t_DetIds1)[idet];
1370 unsigned int detid = truncateId(id, truncateFlag_, false);
1371 double hitEn(0);
1372 if (Cprev.find(detid) != Cprev.end())
1373 hitEn = Cprev[detid].first * (*t_HitEnergies1)[idet];
1374 else
1375 hitEn = (*t_HitEnergies1)[idet];
1376 if (cFactor_)
1377 hitEn *= cFactor_->getCorr(t_Run, id);
1378 etot1 += hitEn;
1379 }
1380 }
1381 for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
1382 if (selectPhi((*t_DetIds3)[idet])) {
1383 unsigned int id = (*t_DetIds3)[idet];
1384 unsigned int detid = truncateId(id, truncateFlag_, false);
1385 double hitEn(0);
1386 if (Cprev.find(detid) != Cprev.end())
1387 hitEn = Cprev[detid].first * (*t_HitEnergies3)[idet];
1388 else
1389 hitEn = (*t_HitEnergies3)[idet];
1390 if (cFactor_)
1391 hitEn *= cFactor_->getCorr(t_Run, id);
1392 etot3 += hitEn;
1393 }
1394 }
1395 }
1396 ediff = etot3 - etot1;
1397 }
1398
1399 double ehcal = (((rcorForm_ == 3) && (cFactor_ != nullptr))
1400 ? (etot * cFactor_->getCorr(entry))
1401 : ((puCorr_ == 0) ? etot
1402 : ((puCorr_ < 0) ? (etot * puFactor(-puCorr_, t_ieta, pmom, etot, ediff))
1403 : puFactorRho(puCorr_, t_ieta, t_rhoh, etot))));
1404 return CalibTree::energyCalor(etot, etot2, ehcal);
1405 }