EgammaElectronAlgos/src/RegressionHelper.cc

0001 #include "RecoEgamma/EgammaElectronAlgos/interface/RegressionHelper.h"
0002 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
0003 #include "RecoEgamma/EgammaTools/interface/EcalRegressionData.h"
0004 #include "DataFormats/Math/interface/deltaR.h"
0005 #include "TVector2.h"
0006 #include "TFile.h"
0007
0008 RegressionHelper::ESGetTokens::ESGetTokens(Configuration const& cfg,
0009                                            bool useEcalReg,
0010                                            bool useCombinationReg,
0011                                            edm::ConsumesCollector& cc)
0012     : caloTopology{cc.esConsumes()}, caloGeometry{cc.esConsumes()} {
0013   if (useEcalReg && cfg.ecalWeightsFromDB) {
0014     ecalRegBarrel = cc.esConsumes<GBRForest, GBRWrapperRcd>(edm::ESInputTag("", cfg.ecalRegressionWeightLabels[0]));
0015     ecalRegEndcap = cc.esConsumes<GBRForest, GBRWrapperRcd>(edm::ESInputTag("", cfg.ecalRegressionWeightLabels[1]));
0016     ecalRegErrorBarrel =
0017         cc.esConsumes<GBRForest, GBRWrapperRcd>(edm::ESInputTag("", cfg.ecalRegressionWeightLabels[2]));
0018     ecalRegErrorEndcap =
0019         cc.esConsumes<GBRForest, GBRWrapperRcd>(edm::ESInputTag("", cfg.ecalRegressionWeightLabels[3]));
0020   }
0021   if (useCombinationReg && cfg.combinationWeightsFromDB) {
0022     combinationReg =
0023         cc.esConsumes<GBRForest, GBRWrapperRcd>(edm::ESInputTag("", cfg.combinationRegressionWeightLabels[0]));
0024   }
0025 }
0026
0027 RegressionHelper::RegressionHelper(const Configuration& config,
0028                                    bool useEcalReg,
0029                                    bool useCombinationReg,
0030                                    edm::ConsumesCollector& cc)
0031     : cfg_(config), esGetTokens_{cfg_, useEcalReg, useCombinationReg, cc} {}
0032
0033 void RegressionHelper::applyEcalRegression(reco::GsfElectron& ele,
0034                                            const reco::VertexCollection& vertices,
0035                                            const EcalRecHitCollection& rechitsEB,
0036                                            const EcalRecHitCollection& rechitsEE) const {
0037   double cor, err;
0038   getEcalRegression(*ele.superCluster(), vertices, rechitsEB, rechitsEE, cor, err);
0039   ele.setCorrectedEcalEnergy(cor * ele.superCluster()->correctedEnergy());
0040   ele.setCorrectedEcalEnergyError(err * ele.superCluster()->correctedEnergy());
0041 }
0042
0043 void RegressionHelper::checkSetup(const edm::EventSetup& es) {
0044   caloTopology_ = &es.getData(esGetTokens_.caloTopology);
0045   caloGeometry_ = &es.getData(esGetTokens_.caloGeometry);
0046
0047   // Ecal regression
0048
0049   // if at least one of the set of weights come from the DB
0050   if (cfg_.ecalWeightsFromDB) {
0051     ecalRegBarrel_ = &es.getData(esGetTokens_.ecalRegBarrel);            // ECAL barrel
0052     ecalRegEndcap_ = &es.getData(esGetTokens_.ecalRegEndcap);            // ECAL endcaps
0053     ecalRegErrorBarrel_ = &es.getData(esGetTokens_.ecalRegErrorBarrel);  // ECAL barrel error
0054     ecalRegErrorEndcap_ = &es.getData(esGetTokens_.ecalRegErrorEndcap);  // ECAL endcap error
0055     ecalRegressionInitialized_ = true;
0056   }
0057   if (cfg_.combinationWeightsFromDB) {
0058     // Combination
0059     combinationReg_ = &es.getData(esGetTokens_.combinationReg);
0060     combinationRegressionInitialized_ = true;
0061   }
0062
0063   // read weights from file - for debugging. Even if it is one single files, 4 files should b set in the vector
0064   if (!cfg_.ecalWeightsFromDB && !ecalRegressionInitialized_ && !cfg_.ecalRegressionWeightFiles.empty()) {
0065     TFile file0(edm::FileInPath(cfg_.ecalRegressionWeightFiles[0].c_str()).fullPath().c_str());
0066     ecalRegBarrel_ = (const GBRForest*)file0.Get(cfg_.ecalRegressionWeightLabels[0].c_str());
0067     file0.Close();
0068     TFile file1(edm::FileInPath(cfg_.ecalRegressionWeightFiles[1].c_str()).fullPath().c_str());
0069     ecalRegEndcap_ = (const GBRForest*)file1.Get(cfg_.ecalRegressionWeightLabels[1].c_str());
0070     file1.Close();
0071     TFile file2(edm::FileInPath(cfg_.ecalRegressionWeightFiles[2].c_str()).fullPath().c_str());
0072     ecalRegErrorBarrel_ = (const GBRForest*)file2.Get(cfg_.ecalRegressionWeightLabels[2].c_str());
0073     file2.Close();
0074     TFile file3(edm::FileInPath(cfg_.ecalRegressionWeightFiles[3].c_str()).fullPath().c_str());
0075     ecalRegErrorEndcap_ = (const GBRForest*)file3.Get(cfg_.ecalRegressionWeightLabels[3].c_str());
0076     ecalRegressionInitialized_ = true;
0077     file3.Close();
0078   }
0079
0080   if (!cfg_.combinationWeightsFromDB && !combinationRegressionInitialized_ &&
0081       !cfg_.combinationRegressionWeightFiles.empty()) {
0082     TFile file0(edm::FileInPath(cfg_.combinationRegressionWeightFiles[0].c_str()).fullPath().c_str());
0083     combinationReg_ = (const GBRForest*)file0.Get(cfg_.combinationRegressionWeightLabels[0].c_str());
0084     combinationRegressionInitialized_ = true;
0085     file0.Close();
0086   }
0087 }
0088
0089 void RegressionHelper::getEcalRegression(const reco::SuperCluster& sc,
0090                                          const reco::VertexCollection& vertices,
0091                                          const EcalRecHitCollection& rechitsEB,
0092                                          const EcalRecHitCollection& rechitsEE,
0093                                          double& energyFactor,
0094                                          double& errorFactor) const {
0095   energyFactor = -999.;
0096   errorFactor = -999.;
0097
0098   std::vector<float> rInputs;
0099   EcalRegressionData regData;
0100   regData.fill(sc, &rechitsEB, &rechitsEE, caloGeometry_, caloTopology_, &vertices);
0101   regData.fillVec(rInputs);
0102   if (sc.seed()->hitsAndFractions()[0].first.subdetId() == EcalBarrel) {
0103     energyFactor = ecalRegBarrel_->GetResponse(&rInputs[0]);
0104     errorFactor = ecalRegErrorBarrel_->GetResponse(&rInputs[0]);
0105   } else if (sc.seed()->hitsAndFractions()[0].first.subdetId() == EcalEndcap) {
0106     energyFactor = ecalRegEndcap_->GetResponse(&rInputs[0]);
0107     errorFactor = ecalRegErrorEndcap_->GetResponse(&rInputs[0]);
0108   } else {
0109     throw cms::Exception("RegressionHelper::calculateRegressedEnergy")
0110         << "Supercluster seed is either EB nor EE!" << std::endl;
0111   }
0112 }
0113
0114 void RegressionHelper::applyCombinationRegression(reco::GsfElectron& ele) const {
0115   float energy = ele.correctedEcalEnergy();
0116   float energyError = ele.correctedEcalEnergyError();
0117   float momentum = ele.trackMomentumAtVtx().R();
0118   float momentumError = ele.trackMomentumError();
0119   int elClass = -1;
0120
0121   switch (ele.classification()) {
0122     case reco::GsfElectron::GOLDEN:
0123       elClass = 0;
0124       break;
0125     case reco::GsfElectron::BIGBREM:
0126       elClass = 1;
0127       break;
0128     case reco::GsfElectron::BADTRACK:
0129       elClass = 2;
0130       break;
0131     case reco::GsfElectron::SHOWERING:
0132       elClass = 3;
0133       break;
0134     case reco::GsfElectron::GAP:
0135       elClass = 4;
0136       break;
0137     default:
0138       elClass = -1;
0139   }
0140
0141   bool isEcalDriven = ele.ecalDriven();
0142   bool isTrackerDriven = ele.trackerDrivenSeed();
0143   bool isEB = ele.isEB();
0144
0145   // compute relative errors and ratio of errors
0146   float energyRelError = energyError / energy;
0147   float momentumRelError = momentumError / momentum;
0148   float errorRatio = energyRelError / momentumRelError;
0149
0150   // calculate E/p and corresponding error
0151   float eOverP = energy / momentum;
0152   float eOverPerror = eOverP * std::hypot(energyRelError, momentumRelError);
0153
0154   // fill input variables
0155   std::vector<float> regressionInputs;
0156   regressionInputs.resize(11, 0.);
0157
0158   regressionInputs[0] = energy;
0159   regressionInputs[1] = energyRelError;
0160   regressionInputs[2] = momentum;
0161   regressionInputs[3] = momentumRelError;
0162   regressionInputs[4] = errorRatio;
0163   regressionInputs[5] = eOverP;
0164   regressionInputs[6] = eOverPerror;
0165   regressionInputs[7] = static_cast<float>(isEcalDriven);
0166   regressionInputs[8] = static_cast<float>(isTrackerDriven);
0167   regressionInputs[9] = static_cast<float>(elClass);
0168   regressionInputs[10] = static_cast<float>(isEB);
0169
0170   // retrieve combination weight
0171   float weight = 0.;
0172   if (eOverP > 0.025 &&
0173       fabs(momentum - energy) < 15. * sqrt(momentumError * momentumError +
0174                                            energyError * energyError))  // protection against crazy track measurement
0175   {
0176     weight = combinationReg_->GetResponse(regressionInputs.data());
0177     if (weight > 1.)
0178       weight = 1.;
0179     else if (weight < 0.)
0180       weight = 0.;
0181   }
0182
0183   float combinedMomentum = weight * momentum + (1. - weight) * energy;
0184   float combinedMomentumError =
0185       sqrt(weight * weight * momentumError * momentumError + (1. - weight) * (1. - weight) * energyError * energyError);
0186
0187   // FIXME : pure tracker electrons have track momentum error of 999.
0188   // If the combination try to combine such electrons then the original combined momentum is kept
0189   if (momentumError != 999. || weight == 0.) {
0190     math::XYZTLorentzVector oldMomentum = ele.p4();
0191     math::XYZTLorentzVector newMomentum(oldMomentum.x() * combinedMomentum / oldMomentum.t(),
0192                                         oldMomentum.y() * combinedMomentum / oldMomentum.t(),
0193                                         oldMomentum.z() * combinedMomentum / oldMomentum.t(),
0194                                         combinedMomentum);
0195
0196     ele.setP4(reco::GsfElectron::P4_COMBINATION, newMomentum, combinedMomentumError, true);
0197     ele.setMass(0.0);
0198   }
0199 }