File indexing completed on 2024-07-16 22:52:33
0001 #include "CommonTools/CandAlgos/interface/ModifyObjectValueBase.h"
0002
0003 #include "FWCore/Utilities/interface/InputTag.h"
0004 #include "FWCore/Utilities/interface/EDGetToken.h"
0005 #include "FWCore/Framework/interface/ESHandle.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007
0008 #include "DataFormats/Common/interface/ValueMap.h"
0009 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0010 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0011 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0012 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0013
0014 #include "RecoEgamma/EgammaTools/interface/EgammaRegressionContainer.h"
0015 #include "RecoEgamma/EgammaTools/interface/EpCombinationTool.h"
0016 #include "RecoEgamma/EgammaTools/interface/EcalClusterLocal.h"
0017 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0018 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0019
0020 #include <vdt/vdtMath.h>
0021
0022 class EGRegressionModifierV3 : public ModifyObjectValueBase {
0023 public:
0024 struct EleRegs {
0025 EleRegs(const edm::ParameterSet& iConfig, edm::ConsumesCollector& cc);
0026 void setEventContent(const edm::EventSetup& iSetup);
0027 EgammaRegressionContainer ecalOnlyMean;
0028 EgammaRegressionContainer ecalOnlySigma;
0029 EpCombinationTool epComb;
0030 };
0031
0032 struct PhoRegs {
0033 PhoRegs(const edm::ParameterSet& iConfig, edm::ConsumesCollector& cc);
0034 void setEventContent(const edm::EventSetup& iSetup);
0035 EgammaRegressionContainer ecalOnlyMean;
0036 EgammaRegressionContainer ecalOnlySigma;
0037 };
0038
0039 EGRegressionModifierV3(const edm::ParameterSet& conf, edm::ConsumesCollector& cc);
0040 ~EGRegressionModifierV3() override;
0041
0042 void setEvent(const edm::Event&) final;
0043 void setEventContent(const edm::EventSetup&) final;
0044
0045 void modifyObject(reco::GsfElectron&) const final;
0046 void modifyObject(reco::Photon&) const final;
0047
0048
0049 void modifyObject(pat::Electron&) const final;
0050 void modifyObject(pat::Photon&) const final;
0051
0052 private:
0053 std::array<float, 32> getRegData(const reco::GsfElectron& ele) const;
0054 std::array<float, 32> getRegData(const reco::Photon& pho) const;
0055 void getSeedCrysCoord(const reco::CaloCluster& clus, int& iEtaOrX, int& iPhiOrY) const;
0056
0057 std::unique_ptr<EleRegs> eleRegs_;
0058 std::unique_ptr<PhoRegs> phoRegs_;
0059
0060 float rhoValue_;
0061 edm::EDGetTokenT<double> rhoToken_;
0062
0063 bool useClosestToCentreSeedCrysDef_;
0064 bool useBuggedHOverE_;
0065 float maxRawEnergyForLowPtEBSigma_;
0066 float maxRawEnergyForLowPtEESigma_;
0067 edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;
0068 edm::ESHandle<CaloGeometry> caloGeomHandle_;
0069 };
0070
0071 DEFINE_EDM_PLUGIN(ModifyObjectValueFactory, EGRegressionModifierV3, "EGRegressionModifierV3");
0072
0073 EGRegressionModifierV3::EGRegressionModifierV3(const edm::ParameterSet& conf, edm::ConsumesCollector& cc)
0074 : ModifyObjectValueBase(conf),
0075 rhoValue_(0.),
0076 rhoToken_(cc.consumes(conf.getParameter<edm::InputTag>("rhoTag"))),
0077 useClosestToCentreSeedCrysDef_(conf.getParameter<bool>("useClosestToCentreSeedCrysDef")),
0078 useBuggedHOverE_(conf.getParameter<bool>("useBuggedHOverE")),
0079 maxRawEnergyForLowPtEBSigma_(conf.getParameter<double>("maxRawEnergyForLowPtEBSigma")),
0080 maxRawEnergyForLowPtEESigma_(conf.getParameter<double>("maxRawEnergyForLowPtEESigma")) {
0081 if (conf.exists("eleRegs")) {
0082 eleRegs_ = std::make_unique<EleRegs>(conf.getParameterSet("eleRegs"), cc);
0083 }
0084 if (conf.exists("phoRegs")) {
0085 phoRegs_ = std::make_unique<PhoRegs>(conf.getParameterSet("phoRegs"), cc);
0086 }
0087 if (useClosestToCentreSeedCrysDef_) {
0088 caloGeomToken_ = cc.esConsumes();
0089 }
0090 }
0091
0092 EGRegressionModifierV3::~EGRegressionModifierV3() {}
0093
0094 void EGRegressionModifierV3::setEvent(const edm::Event& evt) { rhoValue_ = evt.get(rhoToken_); }
0095
0096 void EGRegressionModifierV3::setEventContent(const edm::EventSetup& iSetup) {
0097 if (eleRegs_)
0098 eleRegs_->setEventContent(iSetup);
0099 if (phoRegs_)
0100 phoRegs_->setEventContent(iSetup);
0101 if (useClosestToCentreSeedCrysDef_) {
0102 caloGeomHandle_ = iSetup.getHandle(caloGeomToken_);
0103 }
0104 }
0105
0106 void EGRegressionModifierV3::modifyObject(reco::GsfElectron& ele) const {
0107
0108
0109 if (!eleRegs_)
0110 return;
0111
0112 const reco::SuperClusterRef& superClus = ele.superCluster();
0113
0114
0115 if (superClus->seed()->seed().det() == DetId::Forward)
0116 return;
0117
0118
0119 bool rescaleDependentValues = superClus->clusters().isAvailable();
0120
0121
0122
0123 if (ele.fbrem() == reco::GsfElectron::ClassificationVariables::kDefaultValue)
0124 return;
0125
0126 auto regData = getRegData(ele);
0127 const float rawEnergy = superClus->rawEnergy();
0128 const float rawESEnergy = superClus->preshowerEnergy();
0129
0130 const float rawEt = rawEnergy * superClus->position().rho() / superClus->position().r();
0131 const bool isSaturated = ele.nSaturatedXtals() != 0;
0132
0133 const float ecalMean = eleRegs_->ecalOnlyMean(rawEt, ele.isEB(), isSaturated, regData.data());
0134 const float ecalMeanCorr = ecalMean > 0. ? ecalMean : 1.0;
0135
0136
0137
0138
0139 if (ele.isEB() && maxRawEnergyForLowPtEBSigma_ >= 0 && eleRegs_->ecalOnlySigma.useLowEtBin(rawEt, isSaturated)) {
0140 regData[0] = std::min(regData[0], maxRawEnergyForLowPtEBSigma_);
0141 }
0142 if (!ele.isEB() && maxRawEnergyForLowPtEESigma_ >= 0 && eleRegs_->ecalOnlySigma.useLowEtBin(rawEt, isSaturated)) {
0143 regData[0] = std::min(regData[0], maxRawEnergyForLowPtEESigma_);
0144 }
0145 const float ecalSigma = eleRegs_->ecalOnlySigma(rawEt, ele.isEB(), isSaturated, regData.data());
0146
0147 const float corrEnergy = (rawEnergy + rawESEnergy) * ecalMeanCorr;
0148 const float corrEnergyErr = corrEnergy * ecalSigma;
0149
0150 ele.setCorrectedEcalEnergy(corrEnergy, rescaleDependentValues);
0151 ele.setCorrectedEcalEnergyError(corrEnergyErr);
0152
0153 std::pair<float, float> combEnergyAndErr = eleRegs_->epComb.combine(ele);
0154 const math::XYZTLorentzVector newP4 = ele.p4() * combEnergyAndErr.first / ele.p4().t();
0155 ele.correctMomentum(newP4, ele.trackMomentumError(), combEnergyAndErr.second);
0156 }
0157
0158 void EGRegressionModifierV3::modifyObject(pat::Electron& ele) const {
0159 modifyObject(static_cast<reco::GsfElectron&>(ele));
0160 }
0161
0162 void EGRegressionModifierV3::modifyObject(reco::Photon& pho) const {
0163
0164
0165 if (!phoRegs_)
0166 return;
0167
0168 const reco::SuperClusterRef& superClus = pho.superCluster();
0169
0170
0171 if (superClus->seed()->seed().det() == DetId::Forward)
0172 return;
0173
0174
0175 if (!superClus->clusters().isAvailable())
0176 return;
0177
0178 auto regData = getRegData(pho);
0179
0180 const float rawEnergy = superClus->rawEnergy();
0181 const float rawESEnergy = superClus->preshowerEnergy();
0182
0183 const float rawEt = rawEnergy * superClus->position().rho() / superClus->position().r();
0184 const bool isSaturated = pho.nSaturatedXtals();
0185 const float ecalMean = phoRegs_->ecalOnlyMean(rawEt, pho.isEB(), isSaturated, regData.data());
0186 const float ecalMeanCorr = ecalMean > 0. ? ecalMean : 1.0;
0187
0188
0189 if (pho.isEB() && maxRawEnergyForLowPtEBSigma_ >= 0 && phoRegs_->ecalOnlySigma.useLowEtBin(rawEt, isSaturated)) {
0190 regData[0] = std::min(regData[0], maxRawEnergyForLowPtEBSigma_);
0191 }
0192 if (!pho.isEB() && maxRawEnergyForLowPtEESigma_ >= 0 && phoRegs_->ecalOnlySigma.useLowEtBin(rawEt, isSaturated)) {
0193 regData[0] = std::min(regData[0], maxRawEnergyForLowPtEESigma_);
0194 }
0195 const float ecalSigma = phoRegs_->ecalOnlySigma(rawEt, pho.isEB(), isSaturated, regData.data());
0196
0197 const double corrEnergy = (rawEnergy + rawESEnergy) * ecalMeanCorr;
0198 const double corrEnergyErr = corrEnergy * ecalSigma;
0199
0200 pho.setCorrectedEnergy(reco::Photon::P4type::regression2, corrEnergy, corrEnergyErr, true);
0201 }
0202
0203 void EGRegressionModifierV3::modifyObject(pat::Photon& pho) const { modifyObject(static_cast<reco::Photon&>(pho)); }
0204
0205 std::array<float, 32> EGRegressionModifierV3::getRegData(const reco::GsfElectron& ele) const {
0206 std::array<float, 32> data;
0207
0208 const reco::SuperClusterRef& superClus = ele.superCluster();
0209 const edm::Ptr<reco::CaloCluster>& seedClus = superClus->seed();
0210
0211 const bool isEB = ele.isEB();
0212 const double rawEnergy = superClus->rawEnergy();
0213 const double rawESEnergy = superClus->preshowerEnergy();
0214 const int numberOfClusters = superClus->clusters().size();
0215 const auto& ssFull5x5 = ele.full5x5_showerShape();
0216
0217 float e5x5Inverse = ssFull5x5.e5x5 != 0. ? vdt::fast_inv(ssFull5x5.e5x5) : 0.;
0218
0219 data[0] = rawEnergy;
0220 data[1] = superClus->etaWidth();
0221 data[2] = superClus->phiWidth();
0222 data[3] = superClus->seed()->energy() / rawEnergy;
0223 data[4] = ssFull5x5.e5x5 / rawEnergy;
0224
0225 data[5] = useBuggedHOverE_ ? ele.hcalOverEcalBc() : ele.full5x5_hcalOverEcalBc();
0226 data[6] = rhoValue_;
0227 data[7] = seedClus->eta() - superClus->position().Eta();
0228 data[8] = reco::deltaPhi(seedClus->phi(), superClus->position().Phi());
0229 data[9] = ssFull5x5.r9;
0230 data[10] = ssFull5x5.sigmaIetaIeta;
0231 data[11] = ssFull5x5.sigmaIetaIphi;
0232 data[12] = ssFull5x5.sigmaIphiIphi;
0233 data[13] = ssFull5x5.eMax * e5x5Inverse;
0234 data[14] = ssFull5x5.e2nd * e5x5Inverse;
0235 data[15] = ssFull5x5.eTop * e5x5Inverse;
0236 data[16] = ssFull5x5.eBottom * e5x5Inverse;
0237 data[17] = ssFull5x5.eLeft * e5x5Inverse;
0238 data[18] = ssFull5x5.eRight * e5x5Inverse;
0239 data[19] = ssFull5x5.e2x5Max * e5x5Inverse;
0240 data[20] = ssFull5x5.e2x5Left * e5x5Inverse;
0241 data[21] = ssFull5x5.e2x5Right * e5x5Inverse;
0242 data[22] = ssFull5x5.e2x5Top * e5x5Inverse;
0243 data[23] = ssFull5x5.e2x5Bottom * e5x5Inverse;
0244 data[24] = ele.nSaturatedXtals();
0245 data[25] = std::max(0, numberOfClusters);
0246
0247 if (isEB) {
0248 int iEta, iPhi;
0249 getSeedCrysCoord(*seedClus, iEta, iPhi);
0250 int signIEta = iEta > 0 ? +1 : -1;
0251 data[26] = iEta;
0252 data[27] = iPhi;
0253 data[28] = (iEta - signIEta) % 5;
0254 data[29] = (iPhi - 1) % 2;
0255 const int iEtaCorr = iEta - (iEta > 0 ? +1 : -1);
0256 const int iEtaCorr26 = iEta - (iEta > 0 ? +26 : -26);
0257 data[30] = std::abs(iEta) <= 25 ? iEtaCorr % 20 : iEtaCorr26 % 20;
0258 data[31] = (iPhi - 1) % 20;
0259 } else {
0260 int iX, iY;
0261 getSeedCrysCoord(*seedClus, iX, iY);
0262 data[26] = iX;
0263 data[27] = iY;
0264 data[28] = rawESEnergy / rawEnergy;
0265 }
0266
0267 return data;
0268 }
0269
0270 std::array<float, 32> EGRegressionModifierV3::getRegData(const reco::Photon& pho) const {
0271 std::array<float, 32> data;
0272
0273 const reco::SuperClusterRef& superClus = pho.superCluster();
0274 const edm::Ptr<reco::CaloCluster>& seedClus = superClus->seed();
0275
0276 const bool isEB = pho.isEB();
0277 const double rawEnergy = superClus->rawEnergy();
0278 const double rawESEnergy = superClus->preshowerEnergy();
0279 const int numberOfClusters = superClus->clusters().size();
0280 const auto& ssFull5x5 = pho.full5x5_showerShapeVariables();
0281
0282 float e5x5Inverse = ssFull5x5.e5x5 != 0. ? vdt::fast_inv(ssFull5x5.e5x5) : 0.;
0283
0284 data[0] = rawEnergy;
0285 data[1] = superClus->etaWidth();
0286 data[2] = superClus->phiWidth();
0287 data[3] = superClus->seed()->energy() / rawEnergy;
0288 data[4] = ssFull5x5.e5x5 / rawEnergy;
0289
0290
0291 data[5] = pho.hadronicOverEm();
0292 data[6] = rhoValue_;
0293 data[7] = seedClus->eta() - superClus->position().Eta();
0294 data[8] = reco::deltaPhi(seedClus->phi(), superClus->position().Phi());
0295 data[9] = pho.full5x5_r9();
0296 data[10] = ssFull5x5.sigmaIetaIeta;
0297
0298
0299 data[11] = ssFull5x5.sigmaIetaIphi;
0300 data[12] = ssFull5x5.sigmaIphiIphi;
0301 data[13] = ssFull5x5.maxEnergyXtal * e5x5Inverse;
0302 data[14] = ssFull5x5.e2nd * e5x5Inverse;
0303 data[15] = ssFull5x5.eTop * e5x5Inverse;
0304 data[16] = ssFull5x5.eBottom * e5x5Inverse;
0305 data[17] = ssFull5x5.eLeft * e5x5Inverse;
0306 data[18] = ssFull5x5.eRight * e5x5Inverse;
0307 data[19] = ssFull5x5.e2x5Max * e5x5Inverse;
0308 data[20] = ssFull5x5.e2x5Left * e5x5Inverse;
0309 data[21] = ssFull5x5.e2x5Right * e5x5Inverse;
0310 data[22] = ssFull5x5.e2x5Top * e5x5Inverse;
0311 data[23] = ssFull5x5.e2x5Bottom * e5x5Inverse;
0312 data[24] = pho.nSaturatedXtals();
0313 data[25] = std::max(0, numberOfClusters);
0314
0315 if (isEB) {
0316 int iEta, iPhi;
0317 getSeedCrysCoord(*seedClus, iEta, iPhi);
0318 data[26] = iEta;
0319 data[27] = iPhi;
0320 int signIEta = iEta > 0 ? +1 : -1;
0321 data[28] = (iEta - signIEta) % 5;
0322 data[29] = (iPhi - 1) % 2;
0323 const int iEtaCorr = iEta - (iEta > 0 ? +1 : -1);
0324 const int iEtaCorr26 = iEta - (iEta > 0 ? +26 : -26);
0325 data[30] = std::abs(iEta) <= 25 ? iEtaCorr % 20 : iEtaCorr26 % 20;
0326 data[31] = (iPhi - 1) % 20;
0327 } else {
0328 int iX, iY;
0329 getSeedCrysCoord(*seedClus, iX, iY);
0330 data[26] = iX;
0331 data[27] = iY;
0332 data[28] = rawESEnergy / rawEnergy;
0333 }
0334
0335 return data;
0336 }
0337
0338 void EGRegressionModifierV3::getSeedCrysCoord(const reco::CaloCluster& clus, int& iEtaOrX, int& iPhiOrY) const {
0339 iEtaOrX = 0;
0340 iPhiOrY = 0;
0341
0342 const bool isEB = clus.seed().subdetId() == EcalBarrel;
0343
0344 if (useClosestToCentreSeedCrysDef_) {
0345 float dummy;
0346 if (isEB) {
0347 egammaTools::localEcalClusterCoordsEB(clus, *caloGeomHandle_, dummy, dummy, iEtaOrX, iPhiOrY, dummy, dummy);
0348 } else {
0349 egammaTools::localEcalClusterCoordsEE(clus, *caloGeomHandle_, dummy, dummy, iEtaOrX, iPhiOrY, dummy, dummy);
0350 }
0351 } else {
0352 if (isEB) {
0353 const EBDetId ebId(clus.seed());
0354 iEtaOrX = ebId.ieta();
0355 iPhiOrY = ebId.iphi();
0356 } else {
0357 const EEDetId eeId(clus.seed());
0358 iEtaOrX = eeId.ix();
0359 iPhiOrY = eeId.iy();
0360 }
0361 }
0362 }
0363
0364 EGRegressionModifierV3::EleRegs::EleRegs(const edm::ParameterSet& iConfig, edm::ConsumesCollector& cc)
0365 : ecalOnlyMean(iConfig.getParameterSet("ecalOnlyMean"), cc),
0366 ecalOnlySigma(iConfig.getParameterSet("ecalOnlySigma"), cc),
0367 epComb(iConfig.getParameterSet("epComb"), std::move(cc)) {}
0368
0369 void EGRegressionModifierV3::EleRegs::setEventContent(const edm::EventSetup& iSetup) {
0370 ecalOnlyMean.setEventContent(iSetup);
0371 ecalOnlySigma.setEventContent(iSetup);
0372 epComb.setEventContent(iSetup);
0373 }
0374
0375 EGRegressionModifierV3::PhoRegs::PhoRegs(const edm::ParameterSet& iConfig, edm::ConsumesCollector& cc)
0376 : ecalOnlyMean(iConfig.getParameterSet("ecalOnlyMean"), cc),
0377 ecalOnlySigma(iConfig.getParameterSet("ecalOnlySigma"), cc) {}
0378
0379 void EGRegressionModifierV3::PhoRegs::setEventContent(const edm::EventSetup& iSetup) {
0380 ecalOnlyMean.setEventContent(iSetup);
0381 ecalOnlySigma.setEventContent(iSetup);
0382 }