EgammaClusterAlgos/src/SCEnergyCorrectorSemiParm.cc

0001 #include "RecoEcal/EgammaClusterAlgos/interface/SCEnergyCorrectorSemiParm.h"
0002
0003 #include "FWCore/Utilities/interface/isFinite.h"
0004 #include "FWCore/Utilities/interface/Transition.h"
0005 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0006 #include "DataFormats/VertexReco/interface/Vertex.h"
0007 #include "DataFormats/Math/interface/deltaPhi.h"
0008 #include "RecoEcal/EgammaCoreTools/interface/EcalTools.h"
0009 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
0010 #include "RecoEgamma/EgammaTools/interface/EgammaHGCALIDParamDefaults.h"
0011
0012 #include <vdt/vdtMath.h>
0013
0014 using namespace reco;
0015
0016 namespace {
0017   //bool is if a valid dr was found, float is the dr
0018   std::pair<bool, float> getMaxDRNonSeedCluster(const reco::SuperCluster& sc) {
0019     float maxDR2 = 0.;
0020     const edm::Ptr<reco::CaloCluster>& seedClus = sc.seed();
0021
0022     for (const auto& clus : sc.clusters()) {
0023       if (clus == seedClus) {
0024         continue;
0025       }
0026
0027       // find cluster with max dR
0028       const double dr2 = reco::deltaR2(*clus, *seedClus);
0029       if (dr2 > maxDR2) {
0030         maxDR2 = dr2;
0031       }
0032     }
0033     return {sc.clustersSize() != 1, sc.clustersSize() != 1 ? std::sqrt(maxDR2) : 999.};
0034   }
0035   template <typename T>
0036   int countRecHits(const T& recHitHandle, float threshold) {
0037     int count = 0;
0038     if (recHitHandle.isValid()) {
0039       for (const auto& recHit : *recHitHandle) {
0040         if (recHit.energy() > threshold) {
0041           count++;
0042         }
0043       }
0044     }
0045     return count;
0046   }
0047 }  // namespace
0048
0049 SCEnergyCorrectorSemiParm::SCEnergyCorrectorSemiParm()
0050     : caloTopo_(nullptr),
0051       caloGeom_(nullptr),
0052       isHLT_(false),
0053       isPhaseII_(false),
0054       regTrainedWithPS_(true),
0055       applySigmaIetaIphiBug_(false),
0056       nHitsAboveThresholdEB_(0),
0057       nHitsAboveThresholdEE_(0),
0058       nHitsAboveThresholdHG_(0),
0059       hitsEnergyThreshold_(-1.),
0060       hgcalCylinderR_(0.) {}
0061
0062 SCEnergyCorrectorSemiParm::SCEnergyCorrectorSemiParm(const edm::ParameterSet& iConfig, edm::ConsumesCollector cc)
0063     : SCEnergyCorrectorSemiParm() {
0064   setTokens(iConfig, cc);
0065 }
0066
0067 void SCEnergyCorrectorSemiParm::fillPSetDescription(edm::ParameterSetDescription& desc) {
0068   desc.add<bool>("isHLT", false);
0069   desc.add<bool>("isPhaseII", false);
0070   desc.add<bool>("regTrainedWithPS", true);
0071   desc.add<bool>("applySigmaIetaIphiBug", false);
0072   desc.add<edm::InputTag>("ecalRecHitsEE", edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0073   desc.add<edm::InputTag>("ecalRecHitsEB", edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0074   desc.add<std::string>("regressionKeyEB", "pfscecal_EBCorrection_offline_v2");
0075   desc.add<std::string>("regressionKeyEE", "pfscecal_EECorrection_offline_v2");
0076   desc.add<std::string>("uncertaintyKeyEB", "pfscecal_EBUncertainty_offline_v2");
0077   desc.add<std::string>("uncertaintyKeyEE", "pfscecal_EEUncertainty_offline_v2");
0078   desc.add<double>("regressionMinEB", 0.2);
0079   desc.add<double>("regressionMaxEB", 2.0);
0080   desc.add<double>("regressionMinEE", 0.2);
0081   desc.add<double>("regressionMaxEE", 2.0);
0082   desc.add<double>("uncertaintyMinEB", 0.0002);
0083   desc.add<double>("uncertaintyMaxEB", 0.5);
0084   desc.add<double>("uncertaintyMinEE", 0.0002);
0085   desc.add<double>("uncertaintyMaxEE", 0.5);
0086   desc.add<edm::InputTag>("vertexCollection", edm::InputTag("offlinePrimaryVertices"));
0087   desc.add<double>("eRecHitThreshold", 1.);
0088   desc.add<edm::InputTag>("hgcalRecHits", edm::InputTag());
0089   desc.add<double>("hgcalCylinderR", EgammaHGCALIDParamDefaults::kRCylinder);
0090 }
0091
0092 edm::ParameterSetDescription SCEnergyCorrectorSemiParm::makePSetDescription() {
0093   edm::ParameterSetDescription desc;
0094   fillPSetDescription(desc);
0095   return desc;
0096 }
0097
0098 void SCEnergyCorrectorSemiParm::setEventSetup(const edm::EventSetup& es) {
0099   caloTopo_ = &es.getData(caloTopoToken_);
0100   caloGeom_ = &es.getData(caloGeomToken_);
0101
0102   regParamBarrel_.setForests(es);
0103   regParamEndcap_.setForests(es);
0104
0105   if (isPhaseII_) {
0106     hgcalShowerShapes_.initPerSetup(es);
0107   }
0108 }
0109
0110 void SCEnergyCorrectorSemiParm::setEvent(const edm::Event& event) {
0111   event.getByToken(tokenEBRecHits_, recHitsEB_);
0112   if (!isPhaseII_) {
0113     event.getByToken(tokenEERecHits_, recHitsEE_);
0114   } else {
0115     event.getByToken(tokenHgcalRecHits_, recHitsHgcal_);
0116     hgcalShowerShapes_.initPerEvent(*recHitsHgcal_);
0117   }
0118   if (isHLT_ || isPhaseII_) {
0119     //note countRecHits checks the validity of the handle and returns 0
0120     //if invalid so its okay to call on all rec-hit collections here
0121     nHitsAboveThresholdEB_ = countRecHits(recHitsEB_, hitsEnergyThreshold_);
0122     nHitsAboveThresholdEE_ = countRecHits(recHitsEE_, hitsEnergyThreshold_);
0123     nHitsAboveThresholdHG_ = countRecHits(recHitsHgcal_, hitsEnergyThreshold_);
0124   }
0125   if (!isHLT_) {
0126     event.getByToken(tokenVertices_, vertices_);
0127   }
0128 }
0129
0130 std::pair<double, double> SCEnergyCorrectorSemiParm::getCorrections(const reco::SuperCluster& sc) const {
0131   std::pair<double, double> corrEnergyAndRes = {-1, -1};
0132
0133   const auto regData = getRegData(sc);
0134   if (regData.empty()) {
0135     //supercluster has no valid regression, return default values
0136     return corrEnergyAndRes;
0137   }
0138   DetId seedId = sc.seed()->seed();
0139   const auto& regParam = getRegParam(seedId);
0140
0141   double mean = regParam.mean(regData);
0142   double sigma = regParam.sigma(regData);
0143
0144   double energyCorr = mean * sc.rawEnergy();
0145   if (isHLT_ && sc.seed()->seed().det() == DetId::Ecal && seedId.subdetId() == EcalEndcap && !regTrainedWithPS_) {
0146     energyCorr += sc.preshowerEnergy();
0147   }
0148   double resolutionEst = sigma * energyCorr;
0149
0150   corrEnergyAndRes.first = energyCorr;
0151   corrEnergyAndRes.second = resolutionEst;
0152
0153   return corrEnergyAndRes;
0154 }
0155
0156 void SCEnergyCorrectorSemiParm::modifyObject(reco::SuperCluster& sc) const {
0157   std::pair<double, double> cor = getCorrections(sc);
0158   if (cor.first < 0)
0159     return;
0160   sc.setEnergy(cor.first);
0161   sc.setCorrectedEnergy(cor.first);
0162   if (cor.second >= 0) {
0163     sc.setCorrectedEnergyUncertainty(cor.second);
0164   }
0165 }
0166
0167 std::vector<float> SCEnergyCorrectorSemiParm::getRegData(const reco::SuperCluster& sc) const {
0168   switch (sc.seed()->seed().det()) {
0169     case DetId::Ecal:
0170       if (isPhaseII_ && sc.seed()->seed().subdetId() == EcalEndcap) {
0171         throw cms::Exception("ConfigError") << " Error in SCEnergyCorrectorSemiParm: "
0172                                             << " running over events with EcalEndcap clusters while enabling "
0173                                                "isPhaseII, please set isPhaseII = False in regression config";
0174       }
0175       return isHLT_ ? getRegDataECALHLTV1(sc) : getRegDataECALV1(sc);
0176     case DetId::HGCalEE:
0177       if (!isPhaseII_) {
0178         throw cms::Exception("ConfigError") << " Error in SCEnergyCorrectorSemiParm: "
0179                                             << " running over PhaseII events without enabling isPhaseII, please set "
0180                                                "isPhaseII = True in regression config";
0181       }
0182       return isHLT_ ? getRegDataHGCALHLTV1(sc) : getRegDataHGCALV1(sc);
0183     default:
0184       return std::vector<float>();
0185   }
0186 }
0187
0188 void SCEnergyCorrectorSemiParm::RegParam::setForests(const edm::EventSetup& setup) {
0189   meanForest_ = &setup.getData(meanForestToken_);
0190   sigmaForest_ = &setup.getData(sigmaForestToken_);
0191 }
0192
0193 double SCEnergyCorrectorSemiParm::RegParam::mean(const std::vector<float>& data) const {
0194   return meanForest_ ? meanOutTrans_(meanForest_->GetResponse(data.data())) : -1;
0195 }
0196
0197 double SCEnergyCorrectorSemiParm::RegParam::sigma(const std::vector<float>& data) const {
0198   return sigmaForest_ ? sigmaOutTrans_(sigmaForest_->GetResponse(data.data())) : -1;
0199 }
0200
0201 std::vector<float> SCEnergyCorrectorSemiParm::getRegDataECALV1(const reco::SuperCluster& sc) const {
0202   std::vector<float> eval(30, 0.);
0203
0204   const reco::CaloCluster& seedCluster = *(sc.seed());
0205   const bool iseb = seedCluster.hitsAndFractions()[0].first.subdetId() == EcalBarrel;
0206   const EcalRecHitCollection* recHits = iseb ? recHitsEB_.product() : recHitsEE_.product();
0207
0208   const double raw_energy = sc.rawEnergy();
0209   const int numberOfClusters = sc.clusters().size();
0210
0211   const auto& localCovariances = EcalClusterTools::localCovariances(seedCluster, recHits, caloTopo_);
0212
0213   const float eLeft = EcalClusterTools::eLeft(seedCluster, recHits, caloTopo_);
0214   const float eRight = EcalClusterTools::eRight(seedCluster, recHits, caloTopo_);
0215   const float eTop = EcalClusterTools::eTop(seedCluster, recHits, caloTopo_);
0216   const float eBottom = EcalClusterTools::eBottom(seedCluster, recHits, caloTopo_);
0217
0218   float sigmaIetaIeta = sqrt(localCovariances[0]);
0219   float sigmaIetaIphi = std::numeric_limits<float>::max();
0220   float sigmaIphiIphi = std::numeric_limits<float>::max();
0221
0222   if (!edm::isNotFinite(localCovariances[2]))
0223     sigmaIphiIphi = sqrt(localCovariances[2]);
0224
0225   // extra shower shapes
0226   const float see_by_spp = sigmaIetaIeta * (applySigmaIetaIphiBug_ ? std::numeric_limits<float>::max() : sigmaIphiIphi);
0227   if (see_by_spp > 0) {
0228     sigmaIetaIphi = localCovariances[1] / see_by_spp;
0229   } else if (localCovariances[1] > 0) {
0230     sigmaIetaIphi = 1.f;
0231   } else {
0232     sigmaIetaIphi = -1.f;
0233   }
0234
0235   // calculate sub-cluster variables
0236   std::vector<float> clusterRawEnergy;
0237   clusterRawEnergy.resize(std::max(3, numberOfClusters), 0);
0238   std::vector<float> clusterDEtaToSeed;
0239   clusterDEtaToSeed.resize(std::max(3, numberOfClusters), 0);
0240   std::vector<float> clusterDPhiToSeed;
0241   clusterDPhiToSeed.resize(std::max(3, numberOfClusters), 0);
0242   float clusterMaxDR = 999.;
0243   float clusterMaxDRDPhi = 999.;
0244   float clusterMaxDRDEta = 999.;
0245   float clusterMaxDRRawEnergy = 0.;
0246
0247   size_t iclus = 0;
0248   float maxDR = 0;
0249   edm::Ptr<reco::CaloCluster> pclus;
0250   const edm::Ptr<reco::CaloCluster>& theseed = sc.seed();
0251   // loop over all clusters that aren't the seed
0252   auto clusend = sc.clustersEnd();
0253   for (auto clus = sc.clustersBegin(); clus != clusend; ++clus) {
0254     pclus = *clus;
0255
0256     if (theseed == pclus)
0257       continue;
0258     clusterRawEnergy[iclus] = pclus->energy();
0259     clusterDPhiToSeed[iclus] = reco::deltaPhi(pclus->phi(), theseed->phi());
0260     clusterDEtaToSeed[iclus] = pclus->eta() - theseed->eta();
0261
0262     // find cluster with max dR
0263     const auto the_dr = reco::deltaR(*pclus, *theseed);
0264     if (the_dr > maxDR) {
0265       maxDR = the_dr;
0266       clusterMaxDR = maxDR;
0267       clusterMaxDRDPhi = clusterDPhiToSeed[iclus];
0268       clusterMaxDRDEta = clusterDEtaToSeed[iclus];
0269       clusterMaxDRRawEnergy = clusterRawEnergy[iclus];
0270     }
0271     ++iclus;
0272   }
0273
0274   eval[0] = vertices_->size();
0275   eval[1] = raw_energy;
0276   eval[2] = sc.etaWidth();
0277   eval[3] = sc.phiWidth();
0278   eval[4] = EcalClusterTools::e3x3(seedCluster, recHits, caloTopo_) / raw_energy;
0279   eval[5] = seedCluster.energy() / raw_energy;
0280   eval[6] = EcalClusterTools::eMax(seedCluster, recHits) / raw_energy;
0281   eval[7] = EcalClusterTools::e2nd(seedCluster, recHits) / raw_energy;
0282   eval[8] = (eLeft + eRight != 0.f ? (eLeft - eRight) / (eLeft + eRight) : 0.f);
0283   eval[9] = (eTop + eBottom != 0.f ? (eTop - eBottom) / (eTop + eBottom) : 0.f);
0284   eval[10] = sigmaIetaIeta;
0285   eval[11] = sigmaIetaIphi;
0286   eval[12] = sigmaIphiIphi;
0287   eval[13] = std::max(0, numberOfClusters - 1);
0288   eval[14] = clusterMaxDR;
0289   eval[15] = clusterMaxDRDPhi;
0290   eval[16] = clusterMaxDRDEta;
0291   eval[17] = clusterMaxDRRawEnergy / raw_energy;
0292   eval[18] = clusterRawEnergy[0] / raw_energy;
0293   eval[19] = clusterRawEnergy[1] / raw_energy;
0294   eval[20] = clusterRawEnergy[2] / raw_energy;
0295   eval[21] = clusterDPhiToSeed[0];
0296   eval[22] = clusterDPhiToSeed[1];
0297   eval[23] = clusterDPhiToSeed[2];
0298   eval[24] = clusterDEtaToSeed[0];
0299   eval[25] = clusterDEtaToSeed[1];
0300   eval[26] = clusterDEtaToSeed[2];
0301   if (iseb) {
0302     EBDetId ebseedid(seedCluster.seed());
0303     eval[27] = ebseedid.ieta();
0304     eval[28] = ebseedid.iphi();
0305   } else {
0306     EEDetId eeseedid(seedCluster.seed());
0307     eval[27] = eeseedid.ix();
0308     eval[28] = eeseedid.iy();
0309     //seed cluster eta is only needed for the 106X Ultra Legacy regressions
0310     //and was not used in the 74X regression however as its just an extra varaible
0311     //at the end, its harmless to add for the 74X regression
0312     eval[29] = seedCluster.eta();
0313   }
0314   return eval;
0315 }
0316
0317 std::vector<float> SCEnergyCorrectorSemiParm::getRegDataECALHLTV1(const reco::SuperCluster& sc) const {
0318   std::vector<float> eval(7, 0.);
0319   auto maxDRNonSeedClus = getMaxDRNonSeedCluster(sc);
0320   const float clusterMaxDR = maxDRNonSeedClus.first ? maxDRNonSeedClus.second : 999.;
0321
0322   const reco::CaloCluster& seedCluster = *(sc.seed());
0323   const bool iseb = seedCluster.hitsAndFractions()[0].first.subdetId() == EcalBarrel;
0324   const EcalRecHitCollection* recHits = iseb ? recHitsEB_.product() : recHitsEE_.product();
0325
0326   eval[0] = nHitsAboveThresholdEB_ + nHitsAboveThresholdEE_;
0327   eval[1] = sc.eta();
0328   eval[2] = sc.phiWidth();
0329   eval[3] = EcalClusterTools::e3x3(seedCluster, recHits, caloTopo_) / sc.rawEnergy();
0330   eval[4] = std::max(0, static_cast<int>(sc.clusters().size()) - 1);
0331   eval[5] = clusterMaxDR;
0332   eval[6] = sc.rawEnergy();
0333
0334   return eval;
0335 }
0336
0337 std::vector<float> SCEnergyCorrectorSemiParm::getRegDataHGCALV1(const reco::SuperCluster& sc) const {
0338   std::vector<float> eval(17, 0.);
0339
0340   auto ssCalc = hgcalShowerShapes_.createCalc(sc);
0341   auto pcaWidths = ssCalc.getPCAWidths(hgcalCylinderR_);
0342   auto energyHighestHits = ssCalc.getEnergyHighestHits(2);
0343
0344   auto maxDRNonSeedClus = getMaxDRNonSeedCluster(sc);
0345   const float clusterMaxDR = maxDRNonSeedClus.first ? maxDRNonSeedClus.second : 999.;
0346
0347   eval[0] = sc.rawEnergy();
0348   eval[1] = sc.eta();
0349   eval[2] = sc.etaWidth();
0350   eval[3] = sc.phiWidth();
0351   eval[4] = sc.clusters().size();
0352   eval[5] = sc.hitsAndFractions().size();
0353   eval[6] = clusterMaxDR;
0354   eval[7] = sc.eta() - sc.seed()->eta();
0355   eval[8] = reco::deltaPhi(sc.phi(), sc.seed()->phi());
0356   eval[9] = energyHighestHits[0] / sc.rawEnergy();
0357   eval[10] = energyHighestHits[1] / sc.rawEnergy();
0358   eval[11] = std::sqrt(pcaWidths.sigma2uu);
0359   eval[12] = std::sqrt(pcaWidths.sigma2vv);
0360   eval[13] = std::sqrt(pcaWidths.sigma2ww);
0361   eval[14] = ssCalc.getRvar(hgcalCylinderR_, sc.rawEnergy());
0362   eval[15] = sc.seed()->energy() / sc.rawEnergy();
0363   eval[16] = nHitsAboveThresholdEB_ + nHitsAboveThresholdHG_;
0364
0365   return eval;
0366 }
0367
0368 std::vector<float> SCEnergyCorrectorSemiParm::getRegDataHGCALHLTV1(const reco::SuperCluster& sc) const {
0369   std::vector<float> eval(7, 0.);
0370   const float clusterMaxDR = getMaxDRNonSeedCluster(sc).second;
0371
0372   auto ssCalc = hgcalShowerShapes_.createCalc(sc);
0373
0374   eval[0] = sc.rawEnergy();
0375   eval[1] = sc.eta();
0376   eval[2] = sc.phiWidth();
0377   eval[3] = std::max(0, static_cast<int>(sc.clusters().size()) - 1);
0378   eval[4] = ssCalc.getRvar(hgcalCylinderR_);
0379   eval[5] = clusterMaxDR;
0380   eval[6] = nHitsAboveThresholdEB_ + nHitsAboveThresholdHG_;
0381
0382   return eval;
0383 }