diff --git a/RecoEcal/EgammaClusterAlgos/src/PFECALSuperClusterAlgo.cc b/RecoEcal/EgammaClusterAlgos/src/PFECALSuperClusterAlgo.cc index cf2e01546a220..f9519556cc1c5 100644 --- a/RecoEcal/EgammaClusterAlgos/src/PFECALSuperClusterAlgo.cc +++ b/RecoEcal/EgammaClusterAlgos/src/PFECALSuperClusterAlgo.cc @@ -316,62 +316,25 @@ void PFECALSuperClusterAlgo::buildSuperCluster(CalibClusterPtr& seed, CalibClust // need the vector of raw pointers for a PF width class std::vector bare_ptrs; // calculate necessary parameters and build the SC - double posX(0), posY(0), posZ(0), corrSCEnergy(0), corrPS1Energy(0), corrPS2Energy(0), ePS1(0), ePS2(0), - energyweight(0), energyweighttot(0); - std::vector ps1_energies, ps2_energies; - int condP1(1), condP2(1); + double posX(0), posY(0), posZ(0), corrSCEnergy(0), corrPS1Energy(0), corrPS2Energy(0), energyweight(0), + energyweighttot(0); for (auto& clus : clustered) { - ePS1 = ePS2 = 0; + double ePS1 = 0.0; + double ePS2 = 0.0; energyweight = clus->energy_nocalib(); bare_ptrs.push_back(clus->the_ptr().get()); // update EE calibrated super cluster energies if (isEE) { - ePS1 = ePS2 = 0; - condP1 = condP2 = 1; - ps1_energies.clear(); - ps2_energies.clear(); auto ee_key_val = std::make_pair(clus->the_ptr().key(), edm::Ptr()); const auto clustops = std::equal_range(EEtoPS_->begin(), EEtoPS_->end(), ee_key_val, sortByKey); + std::vector psClusterPointers; for (auto i_ps = clustops.first; i_ps != clustops.second; ++i_ps) { - edm::Ptr psclus(i_ps->second); - - auto const& recH_Frac = psclus->recHitFractions(); - - switch (psclus->layer()) { - case PFLayer::PS1: - ps1_energies.push_back(psclus->energy()); - for (auto const& recH : recH_Frac) { - ESDetId strip1 = recH.recHitRef()->detId(); - if (strip1 != ESDetId(0)) { - ESChannelStatusMap::const_iterator status_p1 = channelStatus_->getMap().find(strip1); - // getStatusCode() == 1 => dead channel - //apply correction if all recHits in dead region - if (status_p1->getStatusCode() == 0) - condP1 = 0; //active - } - } - break; - case PFLayer::PS2: - ps2_energies.push_back(psclus->energy()); - for (auto const& recH : recH_Frac) { - ESDetId strip2 = recH.recHitRef()->detId(); - if (strip2 != ESDetId(0)) { - ESChannelStatusMap::const_iterator status_p2 = channelStatus_->getMap().find(strip2); - if (status_p2->getStatusCode() == 0) - condP2 = 0; - } - } - break; - default: - break; - } + psClusterPointers.push_back(i_ps->second.get()); } - if (condP1 == 1) - ePS1 = -1.; - if (condP2 == 1) - ePS2 = -1.; - _pfEnergyCalibration->energyEm( - *(clus->the_ptr()), ps1_energies, ps2_energies, ePS1, ePS2, applyCrackCorrections_); + auto calibratedEnergies = _pfEnergyCalibration->calibrateEndcapClusterEnergies( + *(clus->the_ptr()), psClusterPointers, *channelStatus_, applyCrackCorrections_); + ePS1 = calibratedEnergies.ps1Energy; + ePS2 = calibratedEnergies.ps2Energy; } if (ePS1 == -1.) diff --git a/RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h b/RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h index 981702f7febc8..585faf1a2a5b9 100644 --- a/RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h +++ b/RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h @@ -5,8 +5,9 @@ #include "CondFormats/PhysicsToolsObjects/interface/PerformancePayloadFromTFormula.h" #include "CondFormats/ESObjects/interface/ESEEIntercalibConstants.h" +#include "CondFormats/ESObjects/interface/ESChannelStatus.h" -class TF1; +#include // -*- C++ -*- // @@ -42,34 +43,24 @@ class PFEnergyCalibration { public: PFEnergyCalibration(); - ~PFEnergyCalibration(); - // ecal calibration for photons - double energyEm(const reco::PFCluster& clusterEcal, - std::vector& EclustersPS1, - std::vector& EclustersPS2, - bool crackCorrection = true) const; double energyEm(const reco::PFCluster& clusterEcal, double ePS1, double ePS2, bool crackCorrection = true) const; - double energyEm(const reco::PFCluster& clusterEcal, - std::vector& EclustersPS1, - std::vector& EclustersPS2, - double& ps1, - double& ps2, - bool crackCorrection = true) const; - double energyEm(const reco::PFCluster& clusterEcal, - double ePS1, - double ePS2, - double& ps1, - double& ps2, - bool crackCorrection = true) const; + struct CalibratedEndcapPFClusterEnergies { + double clusterEnergy = 0.; + double ps1Energy = 0.; + double ps2Energy = 0.; + }; + + CalibratedEndcapPFClusterEnergies calibrateEndcapClusterEnergies( + reco::PFCluster const& eeCluster, + std::vector const& psClusterPointers, + ESChannelStatus const& channelStatus, + bool applyCrackCorrections) const; // ECAL+HCAL (abc) calibration, with E and eta dependent coefficients, for hadrons void energyEmHad(double t, double& e, double& h, double eta, double phi) const; - // Initialize default E- and eta-dependent coefficient functional form - void initializeCalibrationFunctions(); - // Set the run-dependent calibration functions from the global tag void setCalibrationFunctions(const PerformancePayloadFromTFormula* thePFCal) { pfCalibrations = thePFCal; } @@ -79,10 +70,18 @@ class PFEnergyCalibration { friend std::ostream& operator<<(std::ostream& out, const PFEnergyCalibration& calib); -protected: +private: + // ecal calibration for photons + double energyEm(const reco::PFCluster& clusterEcal, + double ePS1, + double ePS2, + double& ps1, + double& ps2, + bool crackCorrection = true) const; + // Calibration functions from global tag - const PerformancePayloadFromTFormula* pfCalibrations; - const ESEEIntercalibConstants* esEEInterCalib_; + const PerformancePayloadFromTFormula* pfCalibrations = nullptr; + const ESEEIntercalibConstants* esEEInterCalib_ = nullptr; // Barrel calibration (eta 0.00 -> 1.48) std::unique_ptr faBarrel; @@ -110,7 +109,6 @@ class PFEnergyCalibration { std::unique_ptr fcEtaEndcapH; std::unique_ptr fdEtaEndcapH; -private: double minimum(double a, double b) const; double dCrackPhi(double phi, double eta) const; double CorrPhi(double phi, double eta) const; @@ -159,7 +157,8 @@ class PFEnergyCalibration { double dEtaEndcapH(double x) const; // Threshold correction (offset) - double threshE, threshH; + const double threshE = 3.5; + const double threshH = 2.5; }; #endif diff --git a/RecoParticleFlow/PFClusterTools/src/PFEnergyCalibration.cc b/RecoParticleFlow/PFClusterTools/src/PFEnergyCalibration.cc index 2800f42034fe0..e654a0fd5bca1 100644 --- a/RecoParticleFlow/PFClusterTools/src/PFEnergyCalibration.cc +++ b/RecoParticleFlow/PFClusterTools/src/PFEnergyCalibration.cc @@ -7,23 +7,13 @@ #include #include #include -#include #include #include #include using namespace std; -PFEnergyCalibration::PFEnergyCalibration() : pfCalibrations(nullptr), esEEInterCalib_(nullptr) { - initializeCalibrationFunctions(); -} - -PFEnergyCalibration::~PFEnergyCalibration() {} - -void PFEnergyCalibration::initializeCalibrationFunctions() { - threshE = 3.5; - threshH = 2.5; - +PFEnergyCalibration::PFEnergyCalibration() { //calibChrisClean.C calibration parameters bhumika Nov, 2018 faBarrel = std::make_unique( "faBarrel", "[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))", 0., 1000.); @@ -161,13 +151,52 @@ void PFEnergyCalibration::initializeCalibrationFunctions() { fdEtaEndcapEH->SetParameter(3, 1.0); } +PFEnergyCalibration::CalibratedEndcapPFClusterEnergies PFEnergyCalibration::calibrateEndcapClusterEnergies( + reco::PFCluster const& eeCluster, + std::vector const& psClusterPointers, + ESChannelStatus const& channelStatus, + bool applyCrackCorrections) const { + double ps1_energy_sum = 0.; + double ps2_energy_sum = 0.; + bool condP1 = true; + bool condP2 = true; + + for (auto const& psclus : psClusterPointers) { + bool cond = true; + for (auto const& recH : psclus->recHitFractions()) { + auto strip = recH.recHitRef()->detId(); + if (strip != ESDetId(0)) { + //getStatusCode() == 0 => active channel + // apply correction if all recHits are dead + if (channelStatus.getMap().find(strip)->getStatusCode() == 0) { + cond = false; + break; + } + } + } + + if (psclus->layer() == PFLayer::PS1) { + ps1_energy_sum += psclus->energy(); + condP1 &= cond; + } else if (psclus->layer() == PFLayer::PS2) { + ps2_energy_sum += psclus->energy(); + condP2 &= cond; + } + } + + double ePS1 = condP1 ? -1. : 0.; + double ePS2 = condP2 ? -1. : 0.; + + double cluscalibe = energyEm(eeCluster, ps1_energy_sum, ps2_energy_sum, ePS1, ePS2, applyCrackCorrections); + + return {cluscalibe, ePS1, ePS2}; +} + void PFEnergyCalibration::energyEmHad(double t, double& e, double& h, double eta, double phi) const { // Use calorimetric energy as true energy for neutral particles - double tt = t; - double ee = e; - double hh = h; - double a = 1.; - double b = 1.; + const double tt = t; + const double ee = e; + const double hh = h; double etaCorrE = 1.; double etaCorrH = 1.; auto absEta = std::abs(eta); @@ -178,8 +207,8 @@ void PFEnergyCalibration::energyEmHad(double t, double& e, double& h, double eta // Barrel calibration if (absEta < 1.48) { // The energy correction - a = e > 0. ? aBarrel(t) : 1.; - b = e > 0. ? bBarrel(t) : cBarrel(t); + double a = e > 0. ? aBarrel(t) : 1.; + double b = e > 0. ? bBarrel(t) : cBarrel(t); double thresh = e > 0. ? threshE : threshH; // Protection against negative calibration @@ -209,8 +238,8 @@ void PFEnergyCalibration::energyEmHad(double t, double& e, double& h, double eta // Endcap calibration } else { // The energy correction - a = e > 0. ? aEndcap(t) : 1.; - b = e > 0. ? bEndcap(t) : cEndcap(t); + double a = e > 0. ? aEndcap(t) : 1.; + double b = e > 0. ? bEndcap(t) : cEndcap(t); double thresh = e > 0. ? threshE : threshH; // Protection against negative calibration @@ -224,8 +253,8 @@ void PFEnergyCalibration::energyEmHad(double t, double& e, double& h, double eta t = min(999.9, max(tt, thresh + a * e + b * h)); // The angular correction - double dEta = std::abs(absEta - 1.5); - double etaPow = dEta * dEta * dEta * dEta; + const double dEta = std::abs(absEta - 1.5); + const double etaPow = dEta * dEta * dEta * dEta; if (e > 0. && thresh > 0.) { if (absEta < 2.5) { @@ -485,58 +514,21 @@ double PFEnergyCalibration::dEtaEndcapEH(double x) const { return fdEtaEndcapEH->Eval(x); } } -// -double PFEnergyCalibration::energyEm(const reco::PFCluster& clusterEcal, - std::vector& EclustersPS1, - std::vector& EclustersPS2, - bool crackCorrection) const { - double ePS1(std::accumulate(EclustersPS1.begin(), EclustersPS1.end(), 0.0)); - double ePS2(std::accumulate(EclustersPS2.begin(), EclustersPS2.end(), 0.0)); - return energyEm(clusterEcal, ePS1, ePS2, crackCorrection); -} double PFEnergyCalibration::energyEm(const reco::PFCluster& clusterEcal, double ePS1, double ePS2, bool crackCorrection) const { - double eEcal = clusterEcal.energy(); - //temporaty ugly fix - reco::PFCluster myPFCluster = clusterEcal; - myPFCluster.calculatePositionREP(); - double eta = myPFCluster.positionREP().eta(); - double phi = myPFCluster.positionREP().phi(); - - double calibrated = Ecorr(eEcal, ePS1, ePS2, eta, phi, crackCorrection); - // if(eEcal!=0 && calibrated==0) std::cout<<"Eecal = "<