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