Back to home page

Project CMSSW displayed by LXR
 
 

    

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   // just calls reco versions
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_;  //this allows us to use the regression corrected H/E which is incorrect wrong
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   //check if we have specified an electron regression correction and
0108   //return the object unmodified if so
0109   if (!eleRegs_)
0110     return;
0111 
0112   const reco::SuperClusterRef& superClus = ele.superCluster();
0113 
0114   // skip HGCAL for now
0115   if (superClus->seed()->seed().det() == DetId::Forward)
0116     return;
0117 
0118   // do not apply corrections in case of missing info (slimmed MiniAOD electrons)
0119   bool rescaleDependentValues = superClus->clusters().isAvailable();
0120 
0121   //check if fbrem is filled as its needed for E/p combination so abort if its set to the default value
0122   //this will be the case for <5 (or current cuts) for miniAOD electrons
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   //bug here, it should include the ES, kept for backwards compat
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   //as the sample is trained flat in pt, the regression's only source of very high energy
0136   //electrons is in the high endcap and therefore it gives a very poor resolution estimate
0137   //to any electrons with this energy, regardless of their actual eta
0138   //hence this lovely hack
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   //check if we have specified an photon regression correction and
0164   //return the object unmodified if so
0165   if (!phoRegs_)
0166     return;
0167 
0168   const reco::SuperClusterRef& superClus = pho.superCluster();
0169 
0170   // skip HGCAL for now
0171   if (superClus->seed()->seed().det() == DetId::Forward)
0172     return;
0173 
0174   // do not apply corrections in case of missing info (happens for some slimmed MiniAOD photons)
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   //bug here, it should include the ES, kept for backwards compat
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   //see the electrons for explaination of this lovely feature
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   //the full5x5 is not regression corrected and thus is the correct one to use
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   //interestingly enough this differs from electrons where it uses cone based
0290   //naively Sam would have thought using cone based is even worse than tower based
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   //interestingly sigmaIEtaIPhi differs in defination here from
0298   //electron & sc definations of sigmaIEtaIPhi
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!