Skip to content

Commit

Permalink
Merge pull request #28960 from guitargeek/PFEnergyCalibration_1
Browse files Browse the repository at this point in the history 
Slightly refactor PF endcap cluster calibration
  • Loading branch information
@cmsbuild
cmsbuild committed Mar 9, 2020
2 parents 761c1de + 6b13a95 commit 48dccbc
Show file tree
Hide file tree
Showing 5 changed files with 110 additions and 210 deletions.
57 changes: 10 additions & 47 deletions RecoEcal/EgammaClusterAlgos/src/PFECALSuperClusterAlgo.cc
View file
Expand Up @@ -316,62 +316,25 @@ void PFECALSuperClusterAlgo::buildSuperCluster(CalibClusterPtr& seed, CalibClust
// need the vector of raw pointers for a PF width class
std::vector<const reco::PFCluster*> 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<double> 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<reco::PFCluster>());
const auto clustops = std::equal_range(EEtoPS_->begin(), EEtoPS_->end(), ee_key_val, sortByKey);
std::vector<reco::PFCluster const*> psClusterPointers;
for (auto i_ps = clustops.first; i_ps != clustops.second; ++i_ps) {
edm::Ptr<reco::PFCluster> 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.)
Expand Down
53 changes: 26 additions & 27 deletions RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h
View file
Expand Up @@ -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 <TF1.h>

// -*- C++ -*-
//
Expand Down Expand Up @@ -42,34 +43,24 @@ class PFEnergyCalibration {
public:
PFEnergyCalibration();

~PFEnergyCalibration();

// ecal calibration for photons
double energyEm(const reco::PFCluster& clusterEcal,
std::vector<double>& EclustersPS1,
std::vector<double>& 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<double>& EclustersPS1,
std::vector<double>& 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<reco::PFCluster const*> 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; }

Expand All @@ -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<TF1> faBarrel;
Expand Down Expand Up @@ -110,7 +109,6 @@ class PFEnergyCalibration {
std::unique_ptr<TF1> fcEtaEndcapH;
std::unique_ptr<TF1> fdEtaEndcapH;

private:
double minimum(double a, double b) const;
double dCrackPhi(double phi, double eta) const;
double CorrPhi(double phi, double eta) const;
Expand Down Expand Up @@ -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
124 changes: 58 additions & 66 deletions RecoParticleFlow/PFClusterTools/src/PFEnergyCalibration.cc
View file
Expand Up @@ -7,23 +7,13 @@
#include <TMath.h>
#include <cmath>
#include <vector>
#include <TF1.h>
#include <map>
#include <algorithm>
#include <numeric>

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<TF1>(
"faBarrel", "[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))", 0., 1000.);
Expand Down Expand Up @@ -161,13 +151,52 @@ void PFEnergyCalibration::initializeCalibrationFunctions() {
fdEtaEndcapEH->SetParameter(3, 1.0);
}

PFEnergyCalibration::CalibratedEndcapPFClusterEnergies PFEnergyCalibration::calibrateEndcapClusterEnergies(
reco::PFCluster const& eeCluster,
std::vector<reco::PFCluster const*> 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);
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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) {
Expand Down Expand Up @@ -485,58 +514,21 @@ double PFEnergyCalibration::dEtaEndcapEH(double x) const {
return fdEtaEndcapEH->Eval(x);
}
}
//
double PFEnergyCalibration::energyEm(const reco::PFCluster& clusterEcal,
std::vector<double>& EclustersPS1,
std::vector<double>& 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 = "<<eEcal<<" eta = "<<eta<<" phi = "<<phi<<std::endl;
return calibrated;
return Ecorr(clusterEcal.energy(), ePS1, ePS2, clusterEcal.eta(), clusterEcal.phi(), crackCorrection);
}

double PFEnergyCalibration::energyEm(const reco::PFCluster& clusterEcal,
std::vector<double>& EclustersPS1,
std::vector<double>& EclustersPS2,
double& ps1,
double& ps2,
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, ps1, ps2, crackCorrection);
}
double PFEnergyCalibration::energyEm(const reco::PFCluster& clusterEcal,
double ePS1,
double ePS2,
double& ps1,
double& ps2,
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, ps1, ps2, crackCorrection);
// if(eEcal!=0 && calibrated==0) std::cout<<"Eecal = "<<eEcal<<" eta = "<<eta<<" phi = "<<phi<<std::endl;
return calibrated;
return Ecorr(clusterEcal.energy(), ePS1, ePS2, clusterEcal.eta(), clusterEcal.phi(), ps1, ps2, crackCorrection);
}

std::ostream& operator<<(std::ostream& out, const PFEnergyCalibration& calib) {
Expand Down Expand Up @@ -622,7 +614,7 @@ std::ostream& operator<<(std::ostream& out, const PFEnergyCalibration& calib) {

//useful to compute the signed distance to the closest crack in the barrel
double PFEnergyCalibration::minimum(double a, double b) const {
if (TMath::Abs(b) < TMath::Abs(a))
if (std::abs(b) < std::abs(a))
a = b;
return a;
}
Expand Down Expand Up @@ -719,7 +711,7 @@ double PFEnergyCalibration::CorrEta(double eta) const {

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]));
(1 + sa[i] * TMath::Sign(1., eta - m[i]) * TMath::Exp(-std::abs(eta - m[i]) / ss[i]));

return result;
}
Expand Down Expand Up @@ -937,7 +929,7 @@ double PFEnergyCalibration::EcorrPS_ePSNil(double eEcal, double eta) const {
//so that <feta()> = 1
constexpr double norm = (p4 + p5 * (2.6 + 1.656) / 2);

double result = eEcal * (p0 + p1 * TMath::Exp(-TMath::Abs(eEcal - p3) / p2)) * (p4 + p5 * eta) / norm;
double result = eEcal * (p0 + p1 * TMath::Exp(-std::abs(eEcal - p3) / p2)) * (p4 + p5 * eta) / norm;

return result;
}
Expand Down Expand Up @@ -981,7 +973,7 @@ double PFEnergyCalibration::Ecorr(

double result = 0;

eta = TMath::Abs(eta);
eta = std::abs(eta);

if (eEcal > 0) {
if (eta <= endBarrel)
Expand Down Expand Up @@ -1021,7 +1013,7 @@ double PFEnergyCalibration::Ecorr(double eEcal,

double result = 0;

eta = TMath::Abs(eta);
eta = std::abs(eta);

if (eEcal > 0) {
if (eta <= endBarrel)
Expand Down

0 comments on commit 48dccbc

Please sign in to comment.