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