From 8aa9ff4fda4178f8f937d9b0b04f9fcf8bb84cc2 Mon Sep 17 00:00:00 2001 From: Sunanda Date: Mon, 8 Apr 2019 19:05:07 +0200 Subject: [PATCH] Correct the calibration algos --- Calibration/HcalCalibAlgos/macros/CalibCorr.C | 104 ++++--- .../HcalCalibAlgos/macros/CalibFitPlots.C | 19 +- .../HcalCalibAlgos/macros/CalibMonitor.C | 33 ++- Calibration/HcalCalibAlgos/macros/CalibTree.C | 262 +++++++++--------- 4 files changed, 240 insertions(+), 178 deletions(-) diff --git a/Calibration/HcalCalibAlgos/macros/CalibCorr.C b/Calibration/HcalCalibAlgos/macros/CalibCorr.C index 25641e2ab879a..6c8a2c3e12f2d 100644 --- a/Calibration/HcalCalibAlgos/macros/CalibCorr.C +++ b/Calibration/HcalCalibAlgos/macros/CalibCorr.C @@ -64,9 +64,12 @@ unsigned int truncateId(unsigned int detId, int truncateFlag, bool debug=false){ return id; } -double puFactor(int type, int ieta, double pmom, double eHcal, double ediff) { +double puFactor(int type, int ieta, double pmom, double eHcal, double ediff, + bool debug=false) { double fac(1.0); + if (debug) std::cout << "Input Type " << type << " ieta " << ieta + << " pmon " << pmom << " E " << eHcal << ":" << ediff; if (type <=2) { double frac = (type == 1) ? 0.02 : 0.03; if (pmom > 0 && ediff > frac*pmom) { @@ -87,50 +90,65 @@ double puFactor(int type, int ieta, double pmom, double eHcal, double ediff) { } } fac = (1.0+a1*(eHcal/pmom)*(ediff/pmom)*(1+a2*(ediff/pmom))); - } else { - int jeta = std::abs(ieta); - double d2p = (ediff/pmom); - const double DELTA_CUT = 0.03; - const int PU_IETA_3 = 25; - if (type == 3) { // 16pu - const double CONST_COR_COEF[4] = { 0.971, 1.008, 0.985, 1.086 }; - const double LINEAR_COR_COEF[4] = { 0, -0.359, -0.251, -0.535 }; - const double SQUARE_COR_COEF[4] = { 0, 0, 0.048, 0.143 }; - const int PU_IETA_1 = 9; - const int PU_IETA_2 = 16; - unsigned icor = (unsigned(jeta >= PU_IETA_1) + - unsigned(jeta >= PU_IETA_2) + - unsigned(jeta >= PU_IETA_3)); - if (d2p > DELTA_CUT) fac = (CONST_COR_COEF[icor] + - LINEAR_COR_COEF[icor]*d2p + - SQUARE_COR_COEF[icor]*d2p*d2p); - } else if (type == 4) { // 17pu - const double CONST_COR_COEF[4] = { 0.974, 1.023, 0.989, 1.077 }; - const double LINEAR_COR_COEF[4] = { 0, -0.524, -0.268, -0.584 }; - const double SQUARE_COR_COEF[4] = { 0, 0, 0.053, 0.170 }; - const int PU_IETA_1 = 9; - const int PU_IETA_2 = 18; - unsigned icor = (unsigned(jeta >= PU_IETA_1) + - unsigned(jeta >= PU_IETA_2) + - unsigned(jeta >= PU_IETA_3)); - if (d2p > DELTA_CUT) fac = (CONST_COR_COEF[icor] + - LINEAR_COR_COEF[icor]*d2p + - SQUARE_COR_COEF[icor]*d2p*d2p); - } else { // 18pu - const double CONST_COR_COEF[4] = { 0.973, 0.998, 0.992, 0.965 }; - const double LINEAR_COR_COEF[4] = { 0, -0.318, -0.261, -0.406 }; - const double SQUARE_COR_COEF[4] = { 0, 0, 0.047, 0.089 }; - const int PU_IETA_1 = 7; - const int PU_IETA_2 = 16; - unsigned icor = (unsigned(jeta >= PU_IETA_1) + - unsigned(jeta >= PU_IETA_2) + - unsigned(jeta >= PU_IETA_3)); - if (d2p > DELTA_CUT) fac = (CONST_COR_COEF[icor] + - LINEAR_COR_COEF[icor]*d2p + - SQUARE_COR_COEF[icor]*d2p*d2p); - } + if (debug) std::cout << " coeff " << a1 << ":" << a2 << " Fac " << fac; + } + } else { + int jeta = std::abs(ieta); + double d2p = (ediff/pmom); + const double DELTA_CUT = 0.03; + const int PU_IETA_3 = 25; + if (type == 3) { // 16pu + const double CONST_COR_COEF[4] = { 0.971, 1.008, 0.985, 1.086 }; + const double LINEAR_COR_COEF[4] = { 0, -0.359, -0.251, -0.535 }; + const double SQUARE_COR_COEF[4] = { 0, 0, 0.048, 0.143 }; + const int PU_IETA_1 = 9; + const int PU_IETA_2 = 16; + unsigned icor = (unsigned(jeta >= PU_IETA_1) + + unsigned(jeta >= PU_IETA_2) + + unsigned(jeta >= PU_IETA_3)); + if (d2p > DELTA_CUT) fac = (CONST_COR_COEF[icor] + + LINEAR_COR_COEF[icor]*d2p + + SQUARE_COR_COEF[icor]*d2p*d2p); + if (debug) std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " + << icor << ":" << CONST_COR_COEF[icor] << ":" + << LINEAR_COR_COEF[icor] << ":" + << SQUARE_COR_COEF[icor] << " Fac " << fac; + } else if (type == 4) { // 17pu + const double CONST_COR_COEF[4] = { 0.974, 1.023, 0.989, 1.077 }; + const double LINEAR_COR_COEF[4] = { 0, -0.524, -0.268, -0.584 }; + const double SQUARE_COR_COEF[4] = { 0, 0, 0.053, 0.170 }; + const int PU_IETA_1 = 9; + const int PU_IETA_2 = 18; + unsigned icor = (unsigned(jeta >= PU_IETA_1) + + unsigned(jeta >= PU_IETA_2) + + unsigned(jeta >= PU_IETA_3)); + if (d2p > DELTA_CUT) fac = (CONST_COR_COEF[icor] + + LINEAR_COR_COEF[icor]*d2p + + SQUARE_COR_COEF[icor]*d2p*d2p); + if (debug) std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " + << icor << ":" << CONST_COR_COEF[icor] << ":" + << LINEAR_COR_COEF[icor] << ":" + << SQUARE_COR_COEF[icor] << " Fac " << fac; + } else { // 18pu + const double CONST_COR_COEF[4] = { 0.973, 0.998, 0.992, 0.965 }; + const double LINEAR_COR_COEF[4] = { 0, -0.318, -0.261, -0.406 }; + const double SQUARE_COR_COEF[4] = { 0, 0, 0.047, 0.089 }; + const int PU_IETA_1 = 7; + const int PU_IETA_2 = 16; + unsigned icor = (unsigned(jeta >= PU_IETA_1) + + unsigned(jeta >= PU_IETA_2) + + unsigned(jeta >= PU_IETA_3)); + if (d2p > DELTA_CUT) fac = (CONST_COR_COEF[icor] + + LINEAR_COR_COEF[icor]*d2p + + SQUARE_COR_COEF[icor]*d2p*d2p); + if (debug) std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " + << icor << ":" << CONST_COR_COEF[icor] << ":" + << LINEAR_COR_COEF[icor] << ":" + << SQUARE_COR_COEF[icor] << " Fac " << fac; } } + if (fac < 0 || fac > 1) fac = 0; + if (debug) std::cout << " Final factor " << fac << std::endl; return fac; } diff --git a/Calibration/HcalCalibAlgos/macros/CalibFitPlots.C b/Calibration/HcalCalibAlgos/macros/CalibFitPlots.C index 51d28a40f578d..445e726bd47ab 100644 --- a/Calibration/HcalCalibAlgos/macros/CalibFitPlots.C +++ b/Calibration/HcalCalibAlgos/macros/CalibFitPlots.C @@ -759,7 +759,7 @@ void FitHistExtended(const char* infile, const char* outfile,std::string prefix, if (total > 4) { sprintf (name, "%sOne", hist1->GetName()); TH1D* hist2 = (TH1D*)hist1->Clone(name); - results meanerr = fitOneGauss(hist2,false,debug); + results meanerr = fitOneGauss(hist2,true,debug); value = meanerr.mean; error = meanerr.errmean; width = meanerr.width; werror= meanerr.errwidth; double wbyv = width/value; @@ -1002,8 +1002,23 @@ void PlotHist(const char* infile, std::string prefix, std::string text, double xmin = hist->GetBinLowEdge(1); int nbin = hist->GetNbinsX(); double xmax = hist->GetBinLowEdge(nbin)+hist->GetBinWidth(nbin); + double mean(0), rms(0), total(0); + int kount(0); + for (int k=2; kGetBinContent(k); + double w = hist->GetBinError(k); + mean += (x*w); + rms += (x*x*w); + total += w; + ++kount; + } + mean /= total; + rms /= total; + double error = sqrt(rms-mean*mean)/total; line = new TLine(xmin,p0,xmax,p0); //etamin,1.0,etamax,1.0); - std::cout << xmin << ":" << xmax << ":" << p0 << std::endl; + std::cout << xmin << ":" << xmax << ":" << p0 << " Mean " << nbin + << ":" << kount << ":" << total <<":" << mean << ":" << rms + << ":" << error << std::endl; line->SetLineWidth(2); line->SetLineStyle(2); line->Draw("same"); diff --git a/Calibration/HcalCalibAlgos/macros/CalibMonitor.C b/Calibration/HcalCalibAlgos/macros/CalibMonitor.C index 8337df651e44d..3fa205f5dd422 100644 --- a/Calibration/HcalCalibAlgos/macros/CalibMonitor.C +++ b/Calibration/HcalCalibAlgos/macros/CalibMonitor.C @@ -798,6 +798,7 @@ void CalibMonitor::Loop() { std::vector kounts(kp1,0); std::vector kount50(20,0); for (Long64_t jentry=0; jentryGetEntry(jentry); nbytes += nb; @@ -1100,9 +1101,11 @@ bool CalibMonitor::goodTrack(double& eHcal, double &cuti, bool debug) { } correctEnergy(eHcal); select = ((t_qltyFlag) && (t_selectTk) && (t_hmaxNearP < cut) && - (t_eMipDR < 1.0)); + (t_eMipDR < 1.0) && (eHcal > 0.001)); if (debug) { - std::cout << " output " << eHcal << ":" << cut << ":" << select<< std::endl; + std::cout << " output " << select << " Based on " << t_qltyFlag << ":" + << t_selectTk << ":" << t_hmaxNearP << ":" << cut << ":" + << t_eMipDR << ":" << eHcal << std::endl; } return select; } @@ -1263,7 +1266,7 @@ void CalibMonitor::savePlot(const std::string& theName, bool append, bool all) { } void CalibMonitor::correctEnergy(double& eHcal) { - + bool debug(false); double pmom = (useGen_ && (t_gentrackP>0)) ? t_gentrackP : t_p; if ((corrPU_ < 0) && (pmom > 0)) { double ediff = (t_eHcal30-t_eHcal10); @@ -1286,7 +1289,13 @@ void CalibMonitor::correctEnergy(double& eHcal) { } ediff = (Etot3-Etot1); } - double fac = puFactor(-corrPU_,t_ieta,pmom,eHcal,ediff); + double fac = puFactor(-corrPU_,t_ieta,pmom,eHcal,ediff,false); + if (debug) { + double fac1 = puFactor(-corrPU_,t_ieta,pmom,eHcal,ediff,true); + double fac2 = puFactor(2,t_ieta,pmom,eHcal,ediff,true); + std::cout << "PU Factor for " << -corrPU_ << " = " << fac1 << "; for 2 = " + << fac2 << std::endl; + } eHcal *= fac; } else if (corrPU_ > 0) { eHcal = puFactorRho(corrPU_,t_ieta,t_rhoh,eHcal); @@ -2485,7 +2494,17 @@ void CalibPlotProperties::Loop() { if (plotHists_) { h_nvtx->Fill(t_nVtx); - if ((std::fabs(rat-1)<0.15) && (kp == kp50) && + bool bad(false); + for (unsigned int k=0; ksize(); ++k) { + unsigned int id = truncateId((*t_DetIds)[k],truncateFlag_,false); + double cfac = corrFactor_->getCorr(id); + if (cFactor_ != 0) + cfac *= cFactor_->getCorr(t_Run,(*t_DetIds)[k]); + double ener = cfac*(*t_HitEnergies)[k]; + if (corrPU_) correctEnergy(ener); + if (ener < 0.001) bad = true; + } + if ((!bad) && (std::fabs(rat-1)<0.15) && (kp == kp50) && ((std::abs(t_ieta) < 15) || (std::abs(t_ieta) > 17))) { float weight = (dataMC_ ? t_EventWeight : t_EventWeight*puweight(t_nVtx)); @@ -2493,7 +2512,7 @@ void CalibPlotProperties::Loop() { sel += weight; std::vector bv(7,0.0f), ev(7,0.0f); std::vector bnrec(7,0), enrec(7,0); - double eb(0), ee(0); + double eb(0), ee(0); for (unsigned int k=0; ksize(); ++k) { unsigned int id = truncateId((*t_DetIds)[k],truncateFlag_,false); double cfac = corrFactor_->getCorr(id); @@ -2559,7 +2578,7 @@ bool CalibPlotProperties::goodTrack(double& eHcal, double &cuti, bool debug) { } correctEnergy(eHcal); select = ((t_qltyFlag) && (t_selectTk) && (t_hmaxNearP < cut) && - (t_eMipDR < 100.0)); + (t_eMipDR < 100.0) && (eHcal > 0.001)); if (debug) { std::cout << " output " << eHcal << ":" << cut << ":" << select<< std::endl; } diff --git a/Calibration/HcalCalibAlgos/macros/CalibTree.C b/Calibration/HcalCalibAlgos/macros/CalibTree.C index 45e7600f7234a..6f03419354c05 100644 --- a/Calibration/HcalCalibAlgos/macros/CalibTree.C +++ b/Calibration/HcalCalibAlgos/macros/CalibTree.C @@ -128,6 +128,19 @@ class CalibTree { public : + struct myEntry { + myEntry (int k=0, double f0=0, double f1=0, + double f2=0) : kount(k), fact0(f0), fact1(f1), fact2(f2) {} + int kount; + double fact0, fact1, fact2; + }; + + struct energyCalor { + energyCalor (double e1=0, double e2=0, + double e3=0) : Etot(e1), Etot2(e2), ehcal(e3) {} + double Etot, Etot2, ehcal; + }; + CalibTree(const char *dupFileName, const char* rcorFileName, int truncateFlag, bool useMean, int runlo, int runhi, int phimin, int phimax, int zside, int nvxlo, int nvxhi, int sysmode, int rbx, int puCorr, @@ -152,19 +165,13 @@ public : bool useWeight, double fraction, bool debug); void fitPol0(TH1D* hist, bool debug); void highEtaFactors(int ietaMax, bool debug); + energyCalor energyHcal(double pmom, bool final); TChain *fChain; //!pointer to the analyzed TTree or TChain Int_t fCurrent;//!current Tree number in a TChain TH1D *h_pbyE, *h_cvg; TProfile *h_Ebyp_bfr, *h_Ebyp_aftr; - struct myEntry { - myEntry (int k=0, double f0=0, double f1=0, - double f2=0) : kount(k), fact0(f0), fact1(f1), fact2(f2) {} - int kount; - double fact0, fact1, fact2; - }; - private: // Declaration of leaf types @@ -660,112 +667,65 @@ Double_t CalibTree::Loop(int loop, TFile *fout, bool useweight, int nMin, double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p; if (goodTrack()) { ++ntkgood; - double Etot(0), Etot2(0); - for (unsigned int idet=0; idet<(*t_DetIds).size(); idet++) { - if (selectPhi((*t_DetIds)[idet])) { - unsigned int id = (*t_DetIds)[idet]; - double hitEn(0); - unsigned int detid = truncateId(id,truncateFlag_,false); - if (Cprev.find(detid) != Cprev.end()) - hitEn = Cprev[detid].first * (*t_HitEnergies)[idet]; - else - hitEn = (*t_HitEnergies)[idet]; - if (cFactor_) hitEn *= cFactor_->getCorr(t_Run,id); - Etot += hitEn; - Etot2 += ((*t_HitEnergies)[idet]); - } - } - // Now the outer cone - double Etot1(0), Etot3(0); - if (t_DetIds1 != 0 && t_DetIds3 != 0) { - for (unsigned int idet=0; idet<(*t_DetIds1).size(); idet++) { - if (selectPhi((*t_DetIds1)[idet])) { - unsigned int id = (*t_DetIds1)[idet]; - unsigned int detid = truncateId(id,truncateFlag_,false); - double hitEn(0); - if (Cprev.find(detid) != Cprev.end()) - hitEn = Cprev[detid].first * (*t_HitEnergies1)[idet]; - else - hitEn = (*t_HitEnergies1)[idet]; - if (cFactor_) hitEn *= cFactor_->getCorr(t_Run,id); - Etot1 += hitEn; - } + CalibTree::energyCalor en = energyHcal(pmom, true); + double evWt = (useweight) ? t_EventWeight : 1.0; + if (en.ehcal > 0.001) { + double pufac = (en.Etot > 0) ? (en.ehcal/en.Etot) : 1.0; + double ratio = en.ehcal/(pmom-t_eMipDR); + if (debug) std::cout << " Weights " << evWt << ":" << pufac + << " Energy " << en.Etot2 << ":" << en.Etot << ":" + << pmom << ":" << t_eMipDR << ":" << t_eHcal << ":" + << en.ehcal << " ratio " << ratio << std::endl; + if (loop==0) { + h_pbyE->Fill(ratio, evWt); + h_Ebyp_bfr->Fill(t_ieta, ratio, evWt); } - for (unsigned int idet=0; idet<(*t_DetIds3).size(); idet++) { - if (selectPhi((*t_DetIds3)[idet])) { - unsigned int id = (*t_DetIds3)[idet]; - unsigned int detid = truncateId(id,truncateFlag_,false); - double hitEn(0); - if (Cprev.find(detid) != Cprev.end()) - hitEn = Cprev[detid].first * (*t_HitEnergies3)[idet]; - else - hitEn = (*t_HitEnergies3)[idet]; - if (cFactor_) hitEn *= cFactor_->getCorr(t_Run,id); - Etot3 += hitEn; - } + if (last){ + h_Ebyp_aftr->Fill(t_ieta, ratio, evWt); } - } - eHcalDelta_ = Etot3-Etot1; - double evWt = (useweight) ? t_EventWeight : 1.0; - // PU correction only for loose isolation cut - double ehcal = ((puCorr_ == 0) ? Etot : - ((puCorr_ < 0) ? (Etot*puFactor(-puCorr_,t_ieta,pmom,Etot,eHcalDelta_)) : - puFactorRho(puCorr_,t_ieta,t_rhoh,Etot))); - double pufac = (Etot > 0) ? (ehcal/Etot) : 1.0; - double ratio = ehcal/(pmom-t_eMipDR); - if (debug) std::cout << " Weights " << evWt << ":" << pufac << " Energy " - << Etot2 << ":" << Etot << ":" << pmom << ":" - << t_eMipDR << ":" << t_eHcal << ":" << ehcal - << " ratio " << ratio << std::endl; - if (loop==0) { - h_pbyE->Fill(ratio, evWt); - h_Ebyp_bfr->Fill(t_ieta, ratio, evWt); - } - if (last){ - h_Ebyp_aftr->Fill(t_ieta, ratio, evWt); - } - bool l1c(true); - if (applyL1Cut != 0) l1c = ((t_mindR1 >= l1Cut) || - ((applyL1Cut == 1) && (t_DataType == 1))); - if ((rmin >=0 && ratio > rmin) && (rmax >= 0 && ratio < rmax) && l1c) { - for (unsigned int idet=0; idet<(*t_DetIds).size(); idet++) { - if (selectPhi((*t_DetIds)[idet])) { - unsigned int id = (*t_DetIds)[idet]; - unsigned int detid = truncateId(id,truncateFlag_,false); - double hitEn=0.0; - if (debug) { - std::cout << "idet " << idet << " detid/hitenergy : " - << std::hex << (*t_DetIds)[idet] << ":" << detid - << "/" << (*t_HitEnergies)[idet] << std::endl; - } - if (Cprev.find(detid) != Cprev.end()) - hitEn = Cprev[detid].first * (*t_HitEnergies)[idet]; - else - hitEn = (*t_HitEnergies)[idet]; - if (cFactor_) hitEn *= cFactor_->getCorr(t_Run,id); - double Wi = evWt * hitEn/Etot; - double Fac = (inverse) ? (ehcal/(pmom-t_eMipDR)) : - ((pmom-t_eMipDR)/ehcal); - double Fac2= Wi*Fac*Fac; - TH1D* hist(0); - std::map::const_iterator itr = histos.find(detid); - if (itr != histos.end()) hist = itr->second; - if (debug) { - std::cout << "Det Id " << std::hex << detid << std::dec - << " " << hist << std::endl; - } - if (hist != 0) hist->Fill(Fac, Wi);//////histola - Fac *= Wi; - if (SumW.find(detid) != SumW.end() ) { - Wi += SumW[detid].fact0; - Fac += SumW[detid].fact1; - Fac2+= SumW[detid].fact2; - int kount = SumW[detid].kount + 1; - SumW[detid] = myEntry(kount,Wi,Fac,Fac2); - nTrks[detid] += evWt; - } else { - SumW.insert(std::pair(detid,myEntry(1,Wi,Fac,Fac2))); - nTrks.insert(std::pair(detid, evWt)); + bool l1c(true); + if (applyL1Cut != 0) l1c = ((t_mindR1 >= l1Cut) || + ((applyL1Cut == 1) && (t_DataType == 1))); + if ((rmin >=0 && ratio > rmin) && (rmax >= 0 && ratio < rmax) && l1c) { + for (unsigned int idet=0; idet<(*t_DetIds).size(); idet++) { + if (selectPhi((*t_DetIds)[idet])) { + unsigned int id = (*t_DetIds)[idet]; + unsigned int detid = truncateId(id,truncateFlag_,false); + double hitEn=0.0; + if (debug) { + std::cout << "idet " << idet << " detid/hitenergy : " + << std::hex << (*t_DetIds)[idet] << ":" << detid + << "/" << (*t_HitEnergies)[idet] << std::endl; + } + if (Cprev.find(detid) != Cprev.end()) + hitEn = Cprev[detid].first * (*t_HitEnergies)[idet]; + else + hitEn = (*t_HitEnergies)[idet]; + if (cFactor_) hitEn *= cFactor_->getCorr(t_Run,id); + double Wi = evWt * hitEn/en.Etot; + double Fac = (inverse) ? (en.ehcal/(pmom-t_eMipDR)) : + ((pmom-t_eMipDR)/en.ehcal); + double Fac2= Wi*Fac*Fac; + TH1D* hist(0); + std::map::const_iterator itr = histos.find(detid); + if (itr != histos.end()) hist = itr->second; + if (debug) { + std::cout << "Det Id " << std::hex << detid << std::dec + << " " << hist << std::endl; + } + if (hist != 0) hist->Fill(Fac, Wi);//////histola + Fac *= Wi; + if (SumW.find(detid) != SumW.end() ) { + Wi += SumW[detid].fact0; + Fac += SumW[detid].fact1; + Fac2+= SumW[detid].fact2; + int kount = SumW[detid].kount + 1; + SumW[detid] = myEntry(kount,Wi,Fac,Fac2); + nTrks[detid] += evWt; + } else { + SumW.insert(std::pair(detid,myEntry(1,Wi,Fac,Fac2))); + nTrks.insert(std::pair(detid, evWt)); + } } } } @@ -1099,22 +1059,12 @@ void CalibTree::makeplots(double rmin, double rmax, int ietaMax, if (std::find(entries.begin(), entries.end(), jentry) != entries.end()) break; nb = fChain->GetEntry(jentry); nbytes += nb; if (goodTrack()) { - double Etot(0); - for (unsigned int idet=0; idet<(*t_DetIds).size(); idet++) { - double hitEn(0); - unsigned int id = (*t_DetIds)[idet]; - unsigned int detid = truncateId(id,truncateFlag_,false); - if (Cprev.find(detid) != Cprev.end()) - hitEn = Cprev[detid].first * (*t_HitEnergies)[idet]; - else - hitEn = (*t_HitEnergies)[idet]; - if (cFactor_) hitEn *= cFactor_->getCorr(t_Run,id); - Etot += hitEn; - } + double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p; + CalibTree::energyCalor en1 = energyHcal(pmom, false); + CalibTree::energyCalor en2 = energyHcal(pmom, true); double evWt = (useweight) ? t_EventWeight : 1.0; - double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p; - double ratioi = t_eHcal/(pmom-t_eMipDR); - double ratiof = Etot/(pmom-t_eMipDR); + double ratioi = en1.ehcal/(pmom-t_eMipDR); + double ratiof = en2.ehcal/(pmom-t_eMipDR); if (t_ieta >= -ietaMax && t_ieta <= ietaMax && t_ieta != 0) { if (ratioi>=rmin && ratioi<=rmax) { histos[0].first->Fill(ratioi,evWt); @@ -1215,3 +1165,63 @@ void CalibTree::highEtaFactors(int ietaMax, bool debug) { } } } + +CalibTree::energyCalor CalibTree::energyHcal(double pmom, bool final) { + + double etot = t_eHcal; + double etot2 = t_eHcal; + double ediff = (t_eHcal30 - t_eHcal10); + if (final) { + etot = etot2 = 0; + for (unsigned int idet=0; idet<(*t_DetIds).size(); idet++) { + if (selectPhi((*t_DetIds)[idet])) { + unsigned int id = (*t_DetIds)[idet]; + double hitEn(0); + unsigned int detid = truncateId(id,truncateFlag_,false); + if (Cprev.find(detid) != Cprev.end()) + hitEn = Cprev[detid].first * (*t_HitEnergies)[idet]; + else + hitEn = (*t_HitEnergies)[idet]; + if (cFactor_) hitEn *= cFactor_->getCorr(t_Run,id); + etot += hitEn; + etot2 += ((*t_HitEnergies)[idet]); + } + } + // Now the outer cone + double etot1(0), etot3(0); + if (t_DetIds1 != 0 && t_DetIds3 != 0) { + for (unsigned int idet=0; idet<(*t_DetIds1).size(); idet++) { + if (selectPhi((*t_DetIds1)[idet])) { + unsigned int id = (*t_DetIds1)[idet]; + unsigned int detid = truncateId(id,truncateFlag_,false); + double hitEn(0); + if (Cprev.find(detid) != Cprev.end()) + hitEn = Cprev[detid].first * (*t_HitEnergies1)[idet]; + else + hitEn = (*t_HitEnergies1)[idet]; + if (cFactor_) hitEn *= cFactor_->getCorr(t_Run,id); + etot1 += hitEn; + } + } + for (unsigned int idet=0; idet<(*t_DetIds3).size(); idet++) { + if (selectPhi((*t_DetIds3)[idet])) { + unsigned int id = (*t_DetIds3)[idet]; + unsigned int detid = truncateId(id,truncateFlag_,false); + double hitEn(0); + if (Cprev.find(detid) != Cprev.end()) + hitEn = Cprev[detid].first * (*t_HitEnergies3)[idet]; + else + hitEn = (*t_HitEnergies3)[idet]; + if (cFactor_) hitEn *= cFactor_->getCorr(t_Run,id); + etot3 += hitEn; + } + } + } + ediff = etot3-etot1; + } + // PU correction only for loose isolation cut + double ehcal = ((puCorr_ == 0) ? etot : + ((puCorr_ < 0) ? (etot*puFactor(-puCorr_,t_ieta,pmom,etot,ediff)) : + puFactorRho(puCorr_,t_ieta,t_rhoh,etot))); + return CalibTree::energyCalor(etot,etot2,ehcal); +}