Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Slightly refactor PF endcap cluster calibration #28960

Merged
merged 4 commits into from Mar 9, 2020
Merged

Slightly refactor PF endcap cluster calibration #28960

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
f04135f
don't store vectors just to accumulate them
guitargeek Feb 12, 2020
9fb176b
factored out calibrateEndcapClusterEnergies function
guitargeek Feb 13, 2020
a932054
address review comments
guitargeek Mar 3, 2020
6b13a95
make one of the energyEm functions private
guitargeek Mar 5, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
57 changes: 10 additions & 47 deletions RecoEcal/EgammaClusterAlgos/src/PFECALSuperClusterAlgo.cc
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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