33 faBarrel = std::make_unique<TF1>(
"faBarrel",
"[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",1.,1000.);
38 faBarrel->SetParameter(4,-0.00759275);
39 faBarrel->SetParameter(5,0.00373563);
42 fbBarrel = std::make_unique<TF1>(
"fbBarrel",
"[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",1.,1000.);
51 fcBarrel = std::make_unique<TF1>(
"fcBarrel",
"[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",1.,1000.);
59 faEtaBarrelEH = std::make_unique<TF1>(
"faEtaBarrelEH",
"[0]+[1]*exp(-x/[2])",1.,1000.);
63 fbEtaBarrelEH = std::make_unique<TF1>(
"fbEtaBarrelEH",
"[0]+[1]*exp(-x/[2])",1.,1000.);
67 faEtaBarrelH = std::make_unique<TF1>(
"faEtaBarrelH",
"[0]+[1]*x",1.,1000.);
70 fbEtaBarrelH = std::make_unique<TF1>(
"fbEtaBarrelH",
"[0]+[1]*exp(-x/[2])",1.,1000.);
75 faEndcap = std::make_unique<TF1>(
"faEndcap",
"[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",1.,1000.);
84 fbEndcap = std::make_unique<TF1>(
"fbEndcap",
"[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",1.,1000.);
93 fcEndcap = std::make_unique<TF1>(
"fcEndcap",
"[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",1.,1000.);
101 fcEndcap->SetParameter(7,-0.901421);
102 faEtaEndcapEH = std::make_unique<TF1>(
"faEtaEndcapEH",
"[0]+[1]*exp(-x/[2])",1.,1000.);
106 fbEtaEndcapEH = std::make_unique<TF1>(
"fbEtaEndcapEH",
"[0]+[1]*exp(-x/[2])",1.,1000.);
110 faEtaEndcapH = std::make_unique<TF1>(
"faEtaEndcapH",
"[0]+[1]*exp(-x/[2])+[3]*[3]*exp(-x*x/([4]*[4]))",1.,1000.);
116 fbEtaEndcapH = std::make_unique<TF1>(
"fbEtaEndcapH",
"[0]+[1]*exp(-x/[2])+[3]*[3]*exp(-x*x/([4]*[4]))",1.,1000.);
133 double etaCorrE = 1.;
134 double etaCorrH = 1.;
136 t =
min(999.9,
max(tt,e+h));
137 if ( t < 1. )
return;
140 if ( absEta < 1.48 ) {
147 if ( a < -0.25 || b < -0.25 ) {
154 t =
min(999.9,
max(tt, thresh+a*e+b*h));
157 if ( e > 0. && thresh > 0. ) {
164 if ( e > 0. && thresh > 0. )
166 if ( h > 0. && thresh > 0. ) {
178 if ( a < -0.25 || b < -0.25 ) {
185 t =
min(999.9,
max(tt, thresh+a*e+b*h));
188 double dEta =
std::abs( absEta - 1.5 );
189 double etaPow = dEta * dEta * dEta * dEta;
191 if ( e > 0. && thresh > 0. ) {
206 if ( e > 0. && thresh > 0. )
208 if ( h > 0. && thresh > 0. ) {
214 if ( e < 0. || h < 0. ) {
217 if ( e < 0. ) e = ee;
218 if ( h < 0. ) h = hh;
453 std::vector<double> &EclustersPS1,
454 std::vector<double> &EclustersPS2,
455 bool crackCorrection )
const {
456 double ePS1(std::accumulate(EclustersPS1.begin(), EclustersPS1.end(), 0.0));
457 double ePS2(std::accumulate(EclustersPS2.begin(), EclustersPS2.end(), 0.0));
458 return energyEm(clusterEcal, ePS1, ePS2, crackCorrection);
465 bool crackCorrection )
const {
466 double eEcal = clusterEcal.
energy();
473 double calibrated =
Ecorr(eEcal,ePS1,ePS2,eta,phi, crackCorrection);
479 std::vector<double> &EclustersPS1,
480 std::vector<double> &EclustersPS2,
481 double& ps1,
double& ps2,
482 bool crackCorrection)
const {
483 double ePS1(std::accumulate(EclustersPS1.begin(), EclustersPS1.end(), 0.0));
484 double ePS2(std::accumulate(EclustersPS2.begin(), EclustersPS2.end(), 0.0));
485 return energyEm(clusterEcal, ePS1, ePS2, ps1, ps2, crackCorrection);
488 double ePS1,
double ePS2,
489 double& ps1,
double& ps2,
490 bool crackCorrection)
const {
491 double eEcal = clusterEcal.
energy();
498 double calibrated =
Ecorr(eEcal,ePS1,ePS2,eta,phi,ps1,ps2,crackCorrection);
507 if(!out )
return out;
509 out<<
"PFEnergyCalibration -- "<<endl;
513 static std::map<std::string, PerformanceResult::ResultType> functType;
530 for(std::map<std::string,PerformanceResult::ResultType>::const_iterator
531 func = functType.begin();
532 func != functType.end();
535 cout <<
"Function: " <<
func->first << endl;
542 std::cout <<
"Default calibration functions : " << std::endl;
600 std::vector<double> fillcPhi() {
601 std::vector<double> retValue;
602 retValue.resize(18,0);
604 for(
unsigned i=1;
i<=17;++
i) retValue[
i]=retValue[0]-2*
i*pi/18;
610 const std::vector<double> cPhi = fillcPhi();
624 if(eta<0) phi +=delta_cPhi;
626 if (phi>=-
pi && phi<=
pi){
629 if (phi<cPhi[17] || phi>=cPhi[0]){
630 if (phi<0) phi+= 2*
pi;
631 m =
minimum(phi -cPhi[0],phi-cPhi[17]-2*
pi);
640 m=
minimum(phi-cPhi[i+1],phi-cPhi[i]);
649 std::cout<<
"Problem in dminphi"<<std::endl;
675 double result = (1+p1*TMath::Gaus(dminphi,p2,p3)+p4*TMath::Gaus(dminphi,p5,p6)+p7*TMath::Gaus(dminphi,p8,p9));
686 constexpr double a[] = {6.13349e-01, 5.08146e-01, 4.44480e-01, 3.3487e-01, 7.65627e-01};
687 constexpr double m[] = {-1.79514e-02, 4.44747e-01, 7.92824e-01, 1.14090e+00, 1.47464e+00};
688 constexpr double s[] = {7.92382e-03, 3.06028e-03, 3.36139e-03, 3.94521e-03, 8.63950e-04};
689 constexpr double sa[] = {1.27228e+01, 3.81517e-02, 1.63507e-01, -6.56480e-02, 1.87160e-01};
690 constexpr double ss[] = {5.48753e-02, -1.00223e-02, 2.22866e-03, 4.26288e-04, 2.67937e-03};
693 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]));
726 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);
757 constexpr double norm = (p1+p2*(2.6+1.656)/2);
778 constexpr double norm = (p4+p5*(2.6+1.656)/2);
780 double result = (1.0012+p0*TMath::Exp(-E/p3)+p1*TMath::Exp(-E/p2))*(p4+p5*eta)/norm;
796 constexpr double norm = (p1+p2*(2.6+1.656)/2);
798 double result = p0*(p1+p2*etaEcal)/norm;
816 bool crackCorrection )
const {
845 double result = E*(p0+p1*TMath::Exp(-E/p2))*(p3+p4*TMath::Gaus(eta,p6,p5)+p7*
eta)/norm;
857 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 ;
866 double result = E*(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
877 double gammaprime=
Gamma(etaEcal)/9
e-5;
889 outputPS2 = corrTotES - outputPS1;
895 outputPS1 = corrTotES - outputPS2;
899 outputPS1 = gammaprime*ePS1;
900 outputPS2 = gammaprime*
Alpha(etaEcal)*ePS2;
903 double E =
Beta(1.0155*eEcal+0.025*(ePS1+0.5976*ePS2)/9
e-5,etaEcal)*eEcal+outputPS1+outputPS2;
912 double corrfac=(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
915 double result = E*corrfac;
937 constexpr double norm = (p4+p5*(2.6+1.656)/2);
939 double result = eEcal*(p0+p1*TMath::Exp(-
TMath::Abs(eEcal-p3)/p2))*(p4+p5*eta)/norm;
970 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;
982 bool crackCorrection )
const {
994 if(eta <= endBarrel) result =
EcorrBarrel(eEcal,eta,phi,crackCorrection);
996 else if((eta < endPS) && ePS1==0 && ePS2==0) result =
EcorrPS_ePSNil(eEcal,eta);
997 else if(eta < endPS) result =
EcorrPS(eEcal,ePS1,ePS2,eta);
1001 else result = eEcal;
1003 if(result<eEcal) result=eEcal;
1022 if(eta <= endBarrel) result =
EcorrBarrel(eEcal,eta,phi,crackCorrection);
1024 else if((eta < endPS) && ePS1==0 && ePS2==0) result =
EcorrPS_ePSNil(eEcal,eta);
1025 else if(eta < endPS) result =
EcorrPS(eEcal,ePS1,ePS2,eta,ps1,ps2);
1029 else result = eEcal;
1031 if(result<eEcal) result=eEcal;
std::unique_ptr< TF1 > faEtaBarrelH
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 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 Alpha(double eta) const
void energyEmHad(double t, double &e, double &h, double eta, double phi) const
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 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
std::unique_ptr< TF1 > fbEtaBarrelH
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 > fbEndcap
std::unique_ptr< TF1 > fcBarrel
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 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