File indexing completed on 2023-03-17 11:17:31
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
0048
0049
0050 if (cfg_.ecalWeightsFromDB) {
0051 ecalRegBarrel_ = &es.getData(esGetTokens_.ecalRegBarrel);
0052 ecalRegEndcap_ = &es.getData(esGetTokens_.ecalRegEndcap);
0053 ecalRegErrorBarrel_ = &es.getData(esGetTokens_.ecalRegErrorBarrel);
0054 ecalRegErrorEndcap_ = &es.getData(esGetTokens_.ecalRegErrorEndcap);
0055 ecalRegressionInitialized_ = true;
0056 }
0057 if (cfg_.combinationWeightsFromDB) {
0058
0059 combinationReg_ = &es.getData(esGetTokens_.combinationReg);
0060 combinationRegressionInitialized_ = true;
0061 }
0062
0063
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
0146 float energyRelError = energyError / energy;
0147 float momentumRelError = momentumError / momentum;
0148 float errorRatio = energyRelError / momentumRelError;
0149
0150
0151 float eOverP = energy / momentum;
0152 float eOverPerror = eOverP * std::hypot(energyRelError, momentumRelError);
0153
0154
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
0171 float weight = 0.;
0172 if (eOverP > 0.025 &&
0173 fabs(momentum - energy) < 15. * sqrt(momentumError * momentumError +
0174 energyError * energyError))
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
0188
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 }