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.);
60 faEtaBarrelEH = std::make_unique<TF1>(
"faEtaBarrelEH",
"[0]+[1]*exp(-x/[2])",1.,1000.);
64 fbEtaBarrelEH = std::make_unique<TF1>(
"fbEtaBarrelEH",
"[0]+[1]*exp(-x/[2])",1.,1000.);
68 faEtaBarrelH = std::make_unique<TF1>(
"faEtaBarrelH",
"[0]+[1]*x",1.,1000.);
71 fbEtaBarrelH = std::make_unique<TF1>(
"fbEtaBarrelH",
"[0]+[1]*exp(-x/[2])",1.,1000.);
76 faEndcap = std::make_unique<TF1>(
"faEndcap",
"[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",1.,1000.);
85 fbEndcap = std::make_unique<TF1>(
"fbEndcap",
"[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",1.,1000.);
94 fcEndcap = std::make_unique<TF1>(
"fcEndcap",
"[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",1.,1000.);
100 fcEndcap->SetParameter(5,0.0708707);
102 fcEndcap->SetParameter(7,-0.922411);
103 faEtaEndcapEH = std::make_unique<TF1>(
"faEtaEndcapEH",
"[0]+[1]*exp(-x/[2])",1.,1000.);
107 fbEtaEndcapEH = std::make_unique<TF1>(
"fbEtaEndcapEH",
"[0]+[1]*exp(-x/[2])",1.,1000.);
111 faEtaEndcapH = std::make_unique<TF1>(
"faEtaEndcapH",
"[0]+[1]*exp(-x/[2])+[3]*[3]*exp(-x*x/([4]*[4]))",1.,1000.);
117 fbEtaEndcapH = std::make_unique<TF1>(
"fbEtaEndcapH",
"[0]+[1]*exp(-x/[2])+[3]*[3]*exp(-x*x/([4]*[4]))",1.,1000.);
134 double etaCorrE = 1.;
135 double etaCorrH = 1.;
137 t =
min(999.9,
max(tt,e+h));
138 if ( t < 1. )
return;
141 if ( absEta < 1.48 ) {
148 if ( a < -0.25 || b < -0.25 ) {
155 t =
min(999.9,
max(tt, thresh+a*e+b*h));
158 if ( e > 0. && thresh > 0. ) {
165 if ( e > 0. && thresh > 0. )
167 if ( h > 0. && thresh > 0. ) {
179 if ( a < -0.25 || b < -0.25 ) {
186 t =
min(999.9,
max(tt, thresh+a*e+b*h));
189 double dEta =
std::abs( absEta - 1.5 );
190 double etaPow = dEta * dEta * dEta * dEta;
192 if ( e > 0. && thresh > 0. ) {
207 if ( e > 0. && thresh > 0. )
209 if ( h > 0. && thresh > 0. ) {
215 if ( e < 0. || h < 0. ) {
218 if ( e < 0. ) e = ee;
219 if ( h < 0. ) h = hh;
454 std::vector<double> &EclustersPS1,
455 std::vector<double> &EclustersPS2,
456 bool crackCorrection )
const {
457 double ePS1(std::accumulate(EclustersPS1.begin(), EclustersPS1.end(), 0.0));
458 double ePS2(std::accumulate(EclustersPS2.begin(), EclustersPS2.end(), 0.0));
459 return energyEm(clusterEcal, ePS1, ePS2, crackCorrection);
466 bool crackCorrection )
const {
467 double eEcal = clusterEcal.
energy();
474 double calibrated =
Ecorr(eEcal,ePS1,ePS2,eta,phi, crackCorrection);
480 std::vector<double> &EclustersPS1,
481 std::vector<double> &EclustersPS2,
482 double& ps1,
double& ps2,
483 bool crackCorrection)
const {
484 double ePS1(std::accumulate(EclustersPS1.begin(), EclustersPS1.end(), 0.0));
485 double ePS2(std::accumulate(EclustersPS2.begin(), EclustersPS2.end(), 0.0));
486 return energyEm(clusterEcal, ePS1, ePS2, ps1, ps2, crackCorrection);
489 double ePS1,
double ePS2,
490 double& ps1,
double& ps2,
491 bool crackCorrection)
const {
492 double eEcal = clusterEcal.
energy();
499 double calibrated =
Ecorr(eEcal,ePS1,ePS2,eta,phi,ps1,ps2,crackCorrection);
508 if(!out )
return out;
510 out<<
"PFEnergyCalibration -- "<<endl;
514 static std::map<std::string, PerformanceResult::ResultType> functType;
531 for(std::map<std::string,PerformanceResult::ResultType>::const_iterator
532 func = functType.begin();
533 func != functType.end();
536 cout <<
"Function: " <<
func->first << endl;
543 std::cout <<
"Default calibration functions : " << std::endl;
601 std::vector<double> fillcPhi() {
602 std::vector<double> retValue;
603 retValue.resize(18,0);
605 for(
unsigned i=1;
i<=17;++
i) retValue[
i]=retValue[0]-2*
i*pi/18;
611 const std::vector<double> cPhi = fillcPhi();
625 if(eta<0) phi +=delta_cPhi;
627 if (phi>=-
pi && phi<=
pi){
630 if (phi<cPhi[17] || phi>=cPhi[0]){
631 if (phi<0) phi+= 2*
pi;
632 m =
minimum(phi -cPhi[0],phi-cPhi[17]-2*
pi);
641 m=
minimum(phi-cPhi[i+1],phi-cPhi[i]);
650 std::cout<<
"Problem in dminphi"<<std::endl;
676 double result = (1+p1*TMath::Gaus(dminphi,p2,p3)+p4*TMath::Gaus(dminphi,p5,p6)+p7*TMath::Gaus(dminphi,p8,p9));
687 constexpr double a[] = {6.13349e-01, 5.08146e-01, 4.44480e-01, 3.3487e-01, 7.65627e-01};
688 constexpr double m[] = {-1.79514e-02, 4.44747e-01, 7.92824e-01, 1.14090e+00, 1.47464e+00};
689 constexpr double s[] = {7.92382e-03, 3.06028e-03, 3.36139e-03, 3.94521e-03, 8.63950e-04};
690 constexpr double sa[] = {1.27228e+01, 3.81517e-02, 1.63507e-01, -6.56480e-02, 1.87160e-01};
691 constexpr double ss[] = {5.48753e-02, -1.00223e-02, 2.22866e-03, 4.26288e-04, 2.67937e-03};
694 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]));
727 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);
758 constexpr double norm = (p1+p2*(2.6+1.656)/2);
779 constexpr double norm = (p4+p5*(2.6+1.656)/2);
781 double result = (1.0012+p0*TMath::Exp(-E/p3)+p1*TMath::Exp(-E/p2))*(p4+p5*eta)/norm;
797 constexpr double norm = (p1+p2*(2.6+1.656)/2);
799 double result = p0*(p1+p2*etaEcal)/norm;
817 bool crackCorrection )
const {
846 double result = E*(p0+p1*TMath::Exp(-E/p2))*(p3+p4*TMath::Gaus(eta,p6,p5)+p7*
eta)/norm;
858 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 ;
867 double result = E*(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
878 double gammaprime=
Gamma(etaEcal)/9
e-5;
890 outputPS2 = corrTotES - outputPS1;
896 outputPS1 = corrTotES - outputPS2;
900 outputPS1 = gammaprime*ePS1;
901 outputPS2 = gammaprime*
Alpha(etaEcal)*ePS2;
904 double E =
Beta(1.0155*eEcal+0.025*(ePS1+0.5976*ePS2)/9
e-5,etaEcal)*eEcal+outputPS1+outputPS2;
913 double corrfac=(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
916 double result = E*corrfac;
938 constexpr double norm = (p4+p5*(2.6+1.656)/2);
940 double result = eEcal*(p0+p1*TMath::Exp(-
TMath::Abs(eEcal-p3)/p2))*(p4+p5*eta)/norm;
971 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;
983 bool crackCorrection )
const {
995 if(eta <= endBarrel) result =
EcorrBarrel(eEcal,eta,phi,crackCorrection);
997 else if((eta < endPS) && ePS1==0 && ePS2==0) result =
EcorrPS_ePSNil(eEcal,eta);
998 else if(eta < endPS) result =
EcorrPS(eEcal,ePS1,ePS2,eta);
1002 else result = eEcal;
1004 if(result<eEcal) result=eEcal;
1023 if(eta <= endBarrel) result =
EcorrBarrel(eEcal,eta,phi,crackCorrection);
1025 else if((eta < endPS) && ePS1==0 && ePS2==0) result =
EcorrPS_ePSNil(eEcal,eta);
1026 else if(eta < endPS) result =
EcorrPS(eEcal,ePS1,ePS2,eta,ps1,ps2);
1030 else result = eEcal;
1032 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