35 faBarrel = std::make_unique<TF1>(
"faBarrel",
"[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",0.,1000.);
41 faBarrel->SetParameter(5,0.000104857);
44 fbBarrel = std::make_unique<TF1>(
"fbBarrel",
"[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",0.,1000.);
53 fcBarrel = std::make_unique<TF1>(
"fcBarrel",
"[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",0.,1000.);
62 faEtaBarrelEH = std::make_unique<TF1>(
"faEtaBarrelEH",
"[0]+[1]*exp(-x/[2])",0.,1000.);
66 fbEtaBarrelEH = std::make_unique<TF1>(
"fbEtaBarrelEH",
"[0]+[1]*exp(-x/[2])",0.,1000.);
70 faEtaBarrelH = std::make_unique<TF1>(
"faEtaBarrelH",
"[0]+[1]*x",0.,1000.);
73 fbEtaBarrelH = std::make_unique<TF1>(
"fbEtaBarrelH",
"[0]+[1]*exp(-x/[2])",0.,1000.);
78 faEndcap = std::make_unique<TF1>(
"faEndcap",
"[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",0.,1000.);
87 fbEndcap = std::make_unique<TF1>(
"fbEndcap",
"[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",0.,1000.);
96 fcEndcap = std::make_unique<TF1>(
"fcEndcap",
"[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",0.,1000.);
104 fcEndcap->SetParameter(7,-0.0997326);
105 faEtaEndcapEH = std::make_unique<TF1>(
"faEtaEndcapEH",
"[0]+[1]*exp(-x/[2])",0.,1000.);
109 fbEtaEndcapEH = std::make_unique<TF1>(
"fbEtaEndcapEH",
"[0]+[1]*exp(-x/[2])",0.,1000.);
113 faEtaEndcapH = std::make_unique<TF1>(
"faEtaEndcapH",
"[0]+[1]*exp(-x/[2])+[3]*[3]*exp(-x*x/([4]*[4]))",0.,1000.);
119 fbEtaEndcapH = std::make_unique<TF1>(
"fbEtaEndcapH",
"[0]+[1]*exp(-x/[2])+[3]*[3]*exp(-x*x/([4]*[4]))",0.,1000.);
128 fcEtaBarrelH = std::make_unique<TF1>(
"fcEtaBarrelH",
"[3]*((x-[0])^[1])+[2]",0.,1000.);
134 fcEtaEndcapH = std::make_unique<TF1>(
"fcEtaEndcapH",
"[3]*((x-[0])^[1])+[2]",0.,1000.);
140 fdEtaEndcapH = std::make_unique<TF1>(
"fdEtaEndcapH",
"[3]*((x-[0])^[1])+[2]",0.,1000.);
146 fcEtaBarrelEH = std::make_unique<TF1>(
"fcEtaBarrelEH",
"[3]*((x-[0])^[1])+[2]",0.,1000.);
152 fcEtaEndcapEH = std::make_unique<TF1>(
"fcEtaEndcapEH",
"[3]*((x-[0])^[1])+[2]",0.,1000.);
158 fdEtaEndcapEH = std::make_unique<TF1>(
"fdEtaEndcapEH",
"[3]*((x-[0])^[1])+[2]",0.,1000.);
177 double etaCorrE = 1.;
178 double etaCorrH = 1.;
180 t =
min(999.9,
max(tt,e+h));
181 if ( t < 1. )
return;
184 if ( absEta < 1.48 ) {
191 if ( a < -0.25 || b < -0.25 ) {
198 t =
min(999.9,
max(tt, thresh+a*e+b*h));
201 if ( e > 0. && thresh > 0. ) {
210 if ( e > 0. && thresh > 0. )
212 if ( h > 0. && thresh > 0. ) {
224 if ( a < -0.25 || b < -0.25 ) {
231 t =
min(999.9,
max(tt, thresh+a*e+b*h));
235 double etaPow = dEta * dEta * dEta *
dEta;
238 if ( e > 0. && thresh > 0. ) {
260 if ( e > 0. && thresh > 0. )
262 if ( h > 0. && thresh > 0. ) {
268 if ( e < 0. || h < 0. ) {
271 if ( e < 0. ) e = ee;
272 if ( h < 0. ) h = hh;
582 std::vector<double> &EclustersPS1,
583 std::vector<double> &EclustersPS2,
584 bool crackCorrection )
const {
585 double ePS1(std::accumulate(EclustersPS1.begin(), EclustersPS1.end(), 0.0));
586 double ePS2(std::accumulate(EclustersPS2.begin(), EclustersPS2.end(), 0.0));
587 return energyEm(clusterEcal, ePS1, ePS2, crackCorrection);
594 bool crackCorrection )
const {
595 double eEcal = clusterEcal.
energy();
602 double calibrated =
Ecorr(eEcal,ePS1,ePS2,eta,phi, crackCorrection);
608 std::vector<double> &EclustersPS1,
609 std::vector<double> &EclustersPS2,
610 double& ps1,
double& ps2,
611 bool crackCorrection)
const {
612 double ePS1(std::accumulate(EclustersPS1.begin(), EclustersPS1.end(), 0.0));
613 double ePS2(std::accumulate(EclustersPS2.begin(), EclustersPS2.end(), 0.0));
614 return energyEm(clusterEcal, ePS1, ePS2, ps1, ps2, crackCorrection);
617 double ePS1,
double ePS2,
618 double& ps1,
double& ps2,
619 bool crackCorrection)
const {
620 double eEcal = clusterEcal.
energy();
627 double calibrated =
Ecorr(eEcal,ePS1,ePS2,eta,phi,ps1,ps2,crackCorrection);
636 if(!out )
return out;
638 out<<
"PFEnergyCalibration -- "<<endl;
642 static const std::map<std::string, PerformanceResult::ResultType> functType = {
667 for(std::map<std::string,PerformanceResult::ResultType>::const_iterator
668 func = functType.begin();
669 func != functType.end();
672 cout <<
"Function: " <<
func->first << endl;
679 std::cout <<
"Default calibration functions : " << std::endl;
738 std::vector<double> fillcPhi() {
739 std::vector<double> retValue;
740 retValue.resize(18,0);
742 for(
unsigned i=1;
i<=17;++
i) retValue[
i]=retValue[0]-2*
i*pi/18;
748 const std::vector<double> cPhi = fillcPhi();
762 if(eta<0) phi +=delta_cPhi;
764 if (phi>=-
pi && phi<=
pi){
767 if (phi<cPhi[17] || phi>=cPhi[0]){
768 if (phi<0) phi+= 2*
pi;
769 m =
minimum(phi -cPhi[0],phi-cPhi[17]-2*
pi);
778 m=
minimum(phi-cPhi[i+1],phi-cPhi[i]);
787 std::cout<<
"Problem in dminphi"<<std::endl;
813 double result = (1+p1*TMath::Gaus(dminphi,p2,p3)+p4*TMath::Gaus(dminphi,p5,p6)+p7*TMath::Gaus(dminphi,p8,p9));
824 constexpr double a[] = {6.13349e-01, 5.08146e-01, 4.44480e-01, 3.3487e-01, 7.65627e-01};
825 constexpr double m[] = {-1.79514e-02, 4.44747e-01, 7.92824e-01, 1.14090e+00, 1.47464e+00};
826 constexpr double s[] = {7.92382e-03, 3.06028e-03, 3.36139e-03, 3.94521e-03, 8.63950e-04};
827 constexpr double sa[] = {1.27228e+01, 3.81517e-02, 1.63507e-01, -6.56480e-02, 1.87160e-01};
828 constexpr double ss[] = {5.48753e-02, -1.00223e-02, 2.22866e-03, 4.26288e-04, 2.67937e-03};
831 for(
unsigned i=0;
i<=4;
i++) result+=a[
i]*TMath::Gaus(eta,m[
i],s[i])*(1+sa[
i]*
TMath::Sign(1.,eta-m[i])*TMath::Exp(-
TMath::Abs(eta-m[i])/ss[i]));
864 double result = (p0+1/(p1+p2*TMath::Power(E,p3))+p4*TMath::Exp(-E/p5)+p6*TMath::Exp(-E*E/(p7*p7)))*(1+p8*eta*eta);
895 constexpr double norm = (p1+p2*(2.6+1.656)/2);
916 constexpr double norm = (p4+p5*(2.6+1.656)/2);
918 double result = (1.0012+p0*TMath::Exp(-E/p3)+p1*TMath::Exp(-E/p2))*(p4+p5*eta)/norm;
934 constexpr double norm = (p1+p2*(2.6+1.656)/2);
936 double result = p0*(p1+p2*etaEcal)/norm;
954 bool crackCorrection )
const {
983 double result = E*(p0+p1*TMath::Exp(-E/p2))*(p3+p4*TMath::Gaus(eta,p6,p5)+p7*
eta)/norm;
995 double E =
Beta(1.0155*eEcal+0.025*(ePS1+0.5976*ePS2)/9
e-5,etaEcal)*eEcal+
Gamma(etaEcal)*(ePS1+
Alpha(etaEcal)*ePS2)/9
e-5 ;
1004 double result = E*(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
1015 double gammaprime=
Gamma(etaEcal)/9
e-5;
1023 else if(outputPS1 == 0 && outputPS2 == -1 &&
esEEInterCalib_ !=
nullptr){
1027 outputPS2 = corrTotES - outputPS1;
1029 else if(outputPS1 == -1 && outputPS2 == 0 &&
esEEInterCalib_ !=
nullptr){
1033 outputPS1 = corrTotES - outputPS2;
1037 outputPS1 = gammaprime*ePS1;
1038 outputPS2 = gammaprime*
Alpha(etaEcal)*ePS2;
1041 double E =
Beta(1.0155*eEcal+0.025*(ePS1+0.5976*ePS2)/9
e-5,etaEcal)*eEcal+outputPS1+outputPS2;
1050 double corrfac=(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
1053 double result = E*corrfac;
1075 constexpr double norm = (p4+p5*(2.6+1.656)/2);
1077 double result = eEcal*(p0+p1*TMath::Exp(-
TMath::Abs(eEcal-p3)/p2))*(p4+p5*eta)/norm;
1108 double result = E*(p0+p1*TMath::Exp(-(E-p3)/p2)+1/(p4+p5*TMath::Power(E,p6)))*(p7+p8*TMath::Gaus(eta,p9,p10)+p11*TMath::Gaus(eta,p12,p13))/norm;
1120 bool crackCorrection )
const {
1132 if(eta <= endBarrel) result =
EcorrBarrel(eEcal,eta,phi,crackCorrection);
1134 else if((eta < endPS) && ePS1==0 && ePS2==0) result =
EcorrPS_ePSNil(eEcal,eta);
1135 else if(eta < endPS) result =
EcorrPS(eEcal,ePS1,ePS2,eta);
1139 else result = eEcal;
1141 if(result<eEcal) result=eEcal;
1160 if(eta <= endBarrel) result =
EcorrBarrel(eEcal,eta,phi,crackCorrection);
1162 else if((eta < endPS) && ePS1==0 && ePS2==0) result =
EcorrPS_ePSNil(eEcal,eta);
1163 else if(eta < endPS) result =
EcorrPS(eEcal,ePS1,ePS2,eta,ps1,ps2);
1167 else result = eEcal;
1169 if(result<eEcal) result=eEcal;
std::unique_ptr< TF1 > faEtaBarrelH
std::unique_ptr< TF1 > fcEtaEndcapEH
const PerformancePayloadFromTFormula * pfCalibrations
double aEtaEndcapEH(double x) const
FWCore Framework interface EventSetupRecordImplementation h
Helper function to determine trigger accepts.
std::unique_ptr< TF1 > faEtaBarrelEH
double bBarrel(double x) const
Particle flow cluster, see clustering algorithm in PFClusterAlgo.
double cEtaEndcapEH(double x) const
double aEtaBarrelEH(double x) const
double aEndcap(double x) const
double bEndcap(double x) const
const ESEEIntercalibConstants * esEEInterCalib_
double bEtaEndcapH(double x) const
double Ecorr(double eEcal, double ePS1, double ePS2, double eta, double phi, bool crackCorrection=true) const
std::unique_ptr< TF1 > fbEtaBarrelEH
float getGammaLow2() const
double aEtaBarrelH(double x) const
double EcorrZoneBeforePS(double E, double eta) const
double cEtaBarrelH(double x) const
double Alpha(double eta) const
void energyEmHad(double t, double &e, double &h, double eta, double phi) const
double cEtaBarrelEH(double x) const
std::unique_ptr< TF1 > fcEtaBarrelEH
void calculatePositionREP()
computes posrep_ once and for all
double bEtaBarrelEH(double x) const
std::unique_ptr< TF1 > fbBarrel
std::unique_ptr< TF1 > faEtaEndcapEH
friend std::ostream & operator<<(std::ostream &out, const PFEnergyCalibration &calib)
double cEndcap(double x) const
const REPPoint & positionREP() const
cluster position: rho, eta, phi
double Beta(double E, double eta) const
std::unique_ptr< TF1 > faBarrel
Abs< T >::type abs(const T &t)
std::unique_ptr< TF1 > faEndcap
double dEtaEndcapEH(double x) const
std::unique_ptr< TF1 > fcEtaBarrelH
double energy() const
cluster energy
double bEtaEndcapEH(double x) const
double EcorrZoneAfterPS(double E, double eta) const
std::pair< int, edm::FunctionWithDict > OK
float getGammaLow3() const
bool insert(BinningVariables::BinningVariablesType, float)
float getGammaLow0() const
double bEtaBarrelH(double x) const
double dEtaEndcapH(double x) const
std::unique_ptr< TF1 > fbEtaBarrelH
std::unique_ptr< TF1 > fdEtaEndcapEH
std::unique_ptr< TF1 > faEtaEndcapH
double CorrBarrel(double E, double eta) const
double Gamma(double etaEcal) const
float getGammaLow1() const
double energyEm(const reco::PFCluster &clusterEcal, std::vector< double > &EclustersPS1, std::vector< double > &EclustersPS2, bool crackCorrection=true) const
double dCrackPhi(double phi, double eta) const
std::unique_ptr< TF1 > fcEtaEndcapH
std::unique_ptr< TF1 > fbEndcap
std::unique_ptr< TF1 > fcBarrel
std::unique_ptr< TF1 > fdEtaEndcapH
double CorrEta(double eta) const
void initializeCalibrationFunctions()
double aEtaEndcapH(double x) const
double aBarrel(double x) const
double CorrPhi(double phi, double eta) const
double EcorrPS_ePSNil(double eEcal, double eta) const
std::unique_ptr< TF1 > fbEtaEndcapEH
double cBarrel(double x) const
double cEtaEndcapH(double x) const
double EcorrBarrel(double E, double eta, double phi, bool crackCorrection=true) const
std::unique_ptr< TF1 > fcEndcap
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
std::unique_ptr< TF1 > fbEtaEndcapH
double EcorrPS(double eEcal, double ePS1, double ePS2, double etaEcal) const
double minimum(double a, double b) const