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;
192 if ( e > 0. && thresh > 0. ) {
197 etaPow = dEta * dEta;
200 etaPow = dEta * dEta * dEta * dEta;
214 if ( e > 0. && thresh > 0. )
216 if ( h > 0. && thresh > 0. ) {
222 if ( e < 0. || h < 0. ) {
225 if ( e < 0. ) e = ee;
226 if ( h < 0. ) h = hh;
461 std::vector<double> &EclustersPS1,
462 std::vector<double> &EclustersPS2,
463 bool crackCorrection )
const {
464 double ePS1(std::accumulate(EclustersPS1.begin(), EclustersPS1.end(), 0.0));
465 double ePS2(std::accumulate(EclustersPS2.begin(), EclustersPS2.end(), 0.0));
466 return energyEm(clusterEcal, ePS1, ePS2, crackCorrection);
473 bool crackCorrection )
const {
474 double eEcal = clusterEcal.
energy();
481 double calibrated =
Ecorr(eEcal,ePS1,ePS2,eta,phi, crackCorrection);
487 std::vector<double> &EclustersPS1,
488 std::vector<double> &EclustersPS2,
489 double& ps1,
double& ps2,
490 bool crackCorrection)
const {
491 double ePS1(std::accumulate(EclustersPS1.begin(), EclustersPS1.end(), 0.0));
492 double ePS2(std::accumulate(EclustersPS2.begin(), EclustersPS2.end(), 0.0));
493 return energyEm(clusterEcal, ePS1, ePS2, ps1, ps2, crackCorrection);
496 double ePS1,
double ePS2,
497 double& ps1,
double& ps2,
498 bool crackCorrection)
const {
499 double eEcal = clusterEcal.
energy();
506 double calibrated =
Ecorr(eEcal,ePS1,ePS2,eta,phi,ps1,ps2,crackCorrection);
515 if(!out )
return out;
517 out<<
"PFEnergyCalibration -- "<<endl;
521 static std::map<std::string, PerformanceResult::ResultType> functType;
538 for(std::map<std::string,PerformanceResult::ResultType>::const_iterator
539 func = functType.begin();
540 func != functType.end();
543 cout <<
"Function: " <<
func->first << endl;
550 std::cout <<
"Default calibration functions : " << std::endl;
608 std::vector<double> fillcPhi() {
609 std::vector<double> retValue;
610 retValue.resize(18,0);
612 for(
unsigned i=1;
i<=17;++
i) retValue[
i]=retValue[0]-2*
i*pi/18;
618 const std::vector<double> cPhi = fillcPhi();
632 if(eta<0) phi +=delta_cPhi;
634 if (phi>=-
pi && phi<=
pi){
637 if (phi<cPhi[17] || phi>=cPhi[0]){
638 if (phi<0) phi+= 2*
pi;
639 m =
minimum(phi -cPhi[0],phi-cPhi[17]-2*
pi);
648 m=
minimum(phi-cPhi[i+1],phi-cPhi[i]);
657 std::cout<<
"Problem in dminphi"<<std::endl;
683 double result = (1+p1*TMath::Gaus(dminphi,p2,p3)+p4*TMath::Gaus(dminphi,p5,p6)+p7*TMath::Gaus(dminphi,p8,p9));
694 constexpr double a[] = {6.13349e-01, 5.08146e-01, 4.44480e-01, 3.3487e-01, 7.65627e-01};
695 constexpr double m[] = {-1.79514e-02, 4.44747e-01, 7.92824e-01, 1.14090e+00, 1.47464e+00};
696 constexpr double s[] = {7.92382e-03, 3.06028e-03, 3.36139e-03, 3.94521e-03, 8.63950e-04};
697 constexpr double sa[] = {1.27228e+01, 3.81517e-02, 1.63507e-01, -6.56480e-02, 1.87160e-01};
698 constexpr double ss[] = {5.48753e-02, -1.00223e-02, 2.22866e-03, 4.26288e-04, 2.67937e-03};
701 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]));
734 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);
765 constexpr double norm = (p1+p2*(2.6+1.656)/2);
786 constexpr double norm = (p4+p5*(2.6+1.656)/2);
788 double result = (1.0012+p0*TMath::Exp(-E/p3)+p1*TMath::Exp(-E/p2))*(p4+p5*eta)/norm;
804 constexpr double norm = (p1+p2*(2.6+1.656)/2);
806 double result = p0*(p1+p2*etaEcal)/norm;
824 bool crackCorrection )
const {
853 double result = E*(p0+p1*TMath::Exp(-E/p2))*(p3+p4*TMath::Gaus(eta,p6,p5)+p7*
eta)/norm;
865 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 ;
874 double result = E*(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
885 double gammaprime=
Gamma(etaEcal)/9
e-5;
893 else if(outputPS1 == 0 && outputPS2 == -1 &&
esEEInterCalib_ !=
nullptr){
897 outputPS2 = corrTotES - outputPS1;
899 else if(outputPS1 == -1 && outputPS2 == 0 &&
esEEInterCalib_ !=
nullptr){
903 outputPS1 = corrTotES - outputPS2;
907 outputPS1 = gammaprime*ePS1;
908 outputPS2 = gammaprime*
Alpha(etaEcal)*ePS2;
911 double E =
Beta(1.0155*eEcal+0.025*(ePS1+0.5976*ePS2)/9
e-5,etaEcal)*eEcal+outputPS1+outputPS2;
920 double corrfac=(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
923 double result = E*corrfac;
945 constexpr double norm = (p4+p5*(2.6+1.656)/2);
947 double result = eEcal*(p0+p1*TMath::Exp(-
TMath::Abs(eEcal-p3)/p2))*(p4+p5*eta)/norm;
978 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;
990 bool crackCorrection )
const {
1002 if(eta <= endBarrel) result =
EcorrBarrel(eEcal,eta,phi,crackCorrection);
1004 else if((eta < endPS) && ePS1==0 && ePS2==0) result =
EcorrPS_ePSNil(eEcal,eta);
1005 else if(eta < endPS) result =
EcorrPS(eEcal,ePS1,ePS2,eta);
1009 else result = eEcal;
1011 if(result<eEcal) result=eEcal;
1030 if(eta <= endBarrel) result =
EcorrBarrel(eEcal,eta,phi,crackCorrection);
1032 else if((eta < endPS) && ePS1==0 && ePS2==0) result =
EcorrPS_ePSNil(eEcal,eta);
1033 else if(eta < endPS) result =
EcorrPS(eEcal,ePS1,ePS2,eta,ps1,ps2);
1037 else result = eEcal;
1039 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