File indexing completed on 2023-03-17 11:17:55
0001 #include "RecoEgamma/EgammaTools/interface/EGEnergyCorrector.h"
0002 #include "CondFormats/GBRForest/interface/GBRForest.h"
0003 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0004 #include "DataFormats/Math/interface/deltaPhi.h"
0005 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0006 #include "DataFormats/VertexReco/interface/Vertex.h"
0007 #include "RecoEgamma/EgammaTools/interface/EcalClusterLocal.h"
0008 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0009
0010 using namespace reco;
0011
0012
0013 EGEnergyCorrector::EGEnergyCorrector(Initializer init) noexcept
0014 : fReadereb(std::move(init.readereb_)),
0015 fReaderebvariance(std::move(init.readerebvariance_)),
0016 fReaderee(std::move(init.readeree_)),
0017 fReadereevariance(std::move(init.readereevariance_)) {
0018 fVals.fill(0.0f);
0019 }
0020
0021
0022 std::pair<double, double> EGEnergyCorrector::CorrectedEnergyWithError(const Photon &p,
0023 const reco::VertexCollection &vtxcol,
0024 EcalClusterLazyTools &clustertools,
0025 CaloGeometry const &caloGeometry) {
0026 const SuperClusterRef s = p.superCluster();
0027 const CaloClusterPtr b = s->seed();
0028
0029
0030 CaloClusterPtr b2;
0031 double ebcmax = -99.;
0032 for (reco::CaloCluster_iterator bit = s->clustersBegin(); bit != s->clustersEnd(); ++bit) {
0033 const CaloClusterPtr bc = *bit;
0034 if (bc->energy() > ebcmax && bc != b) {
0035 b2 = bc;
0036 ebcmax = bc->energy();
0037 }
0038 }
0039
0040
0041 CaloClusterPtr bclast;
0042 double ebcmin = 1e6;
0043 for (reco::CaloCluster_iterator bit = s->clustersBegin(); bit != s->clustersEnd(); ++bit) {
0044 const CaloClusterPtr bc = *bit;
0045 if (bc->energy() < ebcmin && bc != b) {
0046 bclast = bc;
0047 ebcmin = bc->energy();
0048 }
0049 }
0050
0051
0052 CaloClusterPtr bclast2;
0053 ebcmin = 1e6;
0054 for (reco::CaloCluster_iterator bit = s->clustersBegin(); bit != s->clustersEnd(); ++bit) {
0055 const CaloClusterPtr bc = *bit;
0056 if (bc->energy() < ebcmin && bc != b && bc != bclast) {
0057 bclast2 = bc;
0058 ebcmin = bc->energy();
0059 }
0060 }
0061
0062 bool isbarrel = b->hitsAndFractions().at(0).first.subdetId() == EcalBarrel;
0063 bool hasbc2 = b2.isNonnull() && b2->energy() > 0.;
0064 bool hasbclast = bclast.isNonnull() && bclast->energy() > 0.;
0065 bool hasbclast2 = bclast2.isNonnull() && bclast2->energy() > 0.;
0066
0067 if (isbarrel) {
0068
0069 fVals[0] = s->rawEnergy();
0070 fVals[1] = p.r9();
0071 fVals[2] = s->eta();
0072 fVals[3] = s->phi();
0073 fVals[4] = p.e5x5() / s->rawEnergy();
0074 fVals[5] = p.hadronicOverEm();
0075 fVals[6] = s->etaWidth();
0076 fVals[7] = s->phiWidth();
0077
0078
0079 double bemax = clustertools.eMax(*b);
0080 double be2nd = clustertools.e2nd(*b);
0081 double betop = clustertools.eTop(*b);
0082 double bebottom = clustertools.eBottom(*b);
0083 double beleft = clustertools.eLeft(*b);
0084 double beright = clustertools.eRight(*b);
0085
0086 fVals[8] = b->eta() - s->eta();
0087 fVals[9] = reco::deltaPhi(b->phi(), s->phi());
0088 fVals[10] = b->energy() / s->rawEnergy();
0089 fVals[11] = clustertools.e3x3(*b) / b->energy();
0090 fVals[12] = clustertools.e5x5(*b) / b->energy();
0091 fVals[13] = sqrt(clustertools.localCovariances(*b)[0]);
0092 fVals[14] = sqrt(clustertools.localCovariances(*b)[2]);
0093 fVals[15] = clustertools.localCovariances(*b)[1];
0094 fVals[16] = bemax / b->energy();
0095 fVals[17] = log(be2nd / bemax);
0096 fVals[18] = log(betop / bemax);
0097 fVals[19] = log(bebottom / bemax);
0098 fVals[20] = log(beleft / bemax);
0099 fVals[21] = log(beright / bemax);
0100 fVals[22] = (betop - bebottom) / (betop + bebottom);
0101 fVals[23] = (beleft - beright) / (beleft + beright);
0102
0103 double bc2emax = hasbc2 ? clustertools.eMax(*b2) : 0.;
0104 double bc2e2nd = hasbc2 ? clustertools.e2nd(*b2) : 0.;
0105 double bc2etop = hasbc2 ? clustertools.eTop(*b2) : 0.;
0106 double bc2ebottom = hasbc2 ? clustertools.eBottom(*b2) : 0.;
0107 double bc2eleft = hasbc2 ? clustertools.eLeft(*b2) : 0.;
0108 double bc2eright = hasbc2 ? clustertools.eRight(*b2) : 0.;
0109
0110 fVals[24] = hasbc2 ? (b2->eta() - s->eta()) : 0.;
0111 fVals[25] = hasbc2 ? reco::deltaPhi(b2->phi(), s->phi()) : 0.;
0112 fVals[26] = hasbc2 ? b2->energy() / s->rawEnergy() : 0.;
0113 fVals[27] = hasbc2 ? clustertools.e3x3(*b2) / b2->energy() : 0.;
0114 fVals[28] = hasbc2 ? clustertools.e5x5(*b2) / b2->energy() : 0.;
0115 fVals[29] = hasbc2 ? sqrt(clustertools.localCovariances(*b2)[0]) : 0.;
0116 fVals[30] = hasbc2 ? sqrt(clustertools.localCovariances(*b2)[2]) : 0.;
0117 fVals[31] = hasbc2 ? clustertools.localCovariances(*b)[1] : 0.;
0118 fVals[32] = hasbc2 ? bc2emax / b2->energy() : 0.;
0119 fVals[33] = hasbc2 ? log(bc2e2nd / bc2emax) : 0.;
0120 fVals[34] = hasbc2 ? log(bc2etop / bc2emax) : 0.;
0121 fVals[35] = hasbc2 ? log(bc2ebottom / bc2emax) : 0.;
0122 fVals[36] = hasbc2 ? log(bc2eleft / bc2emax) : 0.;
0123 fVals[37] = hasbc2 ? log(bc2eright / bc2emax) : 0.;
0124 fVals[38] = hasbc2 ? (bc2etop - bc2ebottom) / (bc2etop + bc2ebottom) : 0.;
0125 fVals[39] = hasbc2 ? (bc2eleft - bc2eright) / (bc2eleft + bc2eright) : 0.;
0126
0127 fVals[40] = hasbclast ? (bclast->eta() - s->eta()) : 0.;
0128 fVals[41] = hasbclast ? reco::deltaPhi(bclast->phi(), s->phi()) : 0.;
0129 fVals[42] = hasbclast ? bclast->energy() / s->rawEnergy() : 0.;
0130 fVals[43] = hasbclast ? clustertools.e3x3(*bclast) / bclast->energy() : 0.;
0131 fVals[44] = hasbclast ? clustertools.e5x5(*bclast) / bclast->energy() : 0.;
0132 fVals[45] = hasbclast ? sqrt(clustertools.localCovariances(*bclast)[0]) : 0.;
0133 fVals[46] = hasbclast ? sqrt(clustertools.localCovariances(*bclast)[2]) : 0.;
0134 fVals[47] = hasbclast ? clustertools.localCovariances(*bclast)[1] : 0.;
0135
0136 fVals[48] = hasbclast2 ? (bclast2->eta() - s->eta()) : 0.;
0137 fVals[49] = hasbclast2 ? reco::deltaPhi(bclast2->phi(), s->phi()) : 0.;
0138 fVals[50] = hasbclast2 ? bclast2->energy() / s->rawEnergy() : 0.;
0139 fVals[51] = hasbclast2 ? clustertools.e3x3(*bclast2) / bclast2->energy() : 0.;
0140 fVals[52] = hasbclast2 ? clustertools.e5x5(*bclast2) / bclast2->energy() : 0.;
0141 fVals[53] = hasbclast2 ? sqrt(clustertools.localCovariances(*bclast2)[0]) : 0.;
0142 fVals[54] = hasbclast2 ? sqrt(clustertools.localCovariances(*bclast2)[2]) : 0.;
0143 fVals[55] = hasbclast2 ? clustertools.localCovariances(*bclast2)[1] : 0.;
0144
0145
0146
0147
0148 float betacry, bphicry, bthetatilt, bphitilt;
0149 int bieta, biphi;
0150 egammaTools::localEcalClusterCoordsEB(*b, caloGeometry, betacry, bphicry, bieta, biphi, bthetatilt, bphitilt);
0151
0152 fVals[56] = bieta;
0153 fVals[57] = biphi;
0154 fVals[58] = bieta % 5;
0155 fVals[59] = biphi % 2;
0156 fVals[60] = (std::abs(bieta) <= 25) * (bieta % 25) +
0157 (std::abs(bieta) > 25) *
0158 ((bieta - 25 * std::abs(bieta) / bieta) % 20);
0159 fVals[61] = biphi % 20;
0160 fVals[62] = betacry;
0161 fVals[63] = bphicry;
0162
0163
0164 float bc2etacry, bc2phicry, bc2thetatilt, bc2phitilt;
0165 int bc2ieta, bc2iphi;
0166 if (hasbc2)
0167 egammaTools::localEcalClusterCoordsEB(
0168 *b2, caloGeometry, bc2etacry, bc2phicry, bc2ieta, bc2iphi, bc2thetatilt, bc2phitilt);
0169
0170 fVals[64] = hasbc2 ? bc2ieta : 0.;
0171 fVals[65] = hasbc2 ? bc2iphi : 0.;
0172 fVals[66] = hasbc2 ? bc2ieta % 5 : 0.;
0173 fVals[67] = hasbc2 ? bc2iphi % 2 : 0.;
0174 fVals[68] = hasbc2 ? (std::abs(bc2ieta) <= 25) * (bc2ieta % 25) +
0175 (std::abs(bc2ieta) > 25) * ((bc2ieta - 25 * std::abs(bc2ieta) / bc2ieta) % 20)
0176 : 0.;
0177 fVals[69] = hasbc2 ? bc2iphi % 20 : 0.;
0178 fVals[70] = hasbc2 ? bc2etacry : 0.;
0179 fVals[71] = hasbc2 ? bc2phicry : 0.;
0180
0181 fVals[72] = vtxcol.size();
0182
0183 } else {
0184 fVals[0] = s->rawEnergy();
0185 fVals[1] = p.r9();
0186 fVals[2] = s->eta();
0187 fVals[3] = s->phi();
0188 fVals[4] = p.e5x5() / s->rawEnergy();
0189 fVals[5] = s->etaWidth();
0190 fVals[6] = s->phiWidth();
0191 fVals[7] = vtxcol.size();
0192 }
0193
0194 const double varscale = 1.253;
0195 double den;
0196 const GBRForest *reader;
0197 const GBRForest *readervar;
0198 if (isbarrel) {
0199 den = s->rawEnergy();
0200 reader = fReadereb.get();
0201 readervar = fReaderebvariance.get();
0202 } else {
0203 den = s->rawEnergy() + s->preshowerEnergy();
0204 reader = fReaderee.get();
0205 readervar = fReadereevariance.get();
0206 }
0207
0208 double ecor = reader->GetResponse(fVals.data()) * den;
0209 double ecorerr = readervar->GetResponse(fVals.data()) * den * varscale;
0210
0211
0212
0213 return {ecor, ecorerr};
0214 }
0215
0216
0217 std::pair<double, double> EGEnergyCorrector::CorrectedEnergyWithErrorV3(const Photon &p,
0218 const reco::VertexCollection &vtxcol,
0219 double rho,
0220 EcalClusterLazyTools &clustertools,
0221 CaloGeometry const &caloGeometry,
0222 bool applyRescale) {
0223 const SuperClusterRef s = p.superCluster();
0224 const CaloClusterPtr b = s->seed();
0225
0226 bool isbarrel = b->hitsAndFractions().at(0).first.subdetId() == EcalBarrel;
0227
0228
0229 fVals[0] = s->rawEnergy();
0230 fVals[1] = s->eta();
0231 fVals[2] = s->phi();
0232 fVals[3] = p.r9();
0233 fVals[4] = p.e5x5() / s->rawEnergy();
0234 fVals[5] = s->etaWidth();
0235 fVals[6] = s->phiWidth();
0236 fVals[7] = s->clustersSize();
0237 fVals[8] = p.hadTowOverEm();
0238 fVals[9] = rho;
0239 fVals[10] = vtxcol.size();
0240
0241
0242 double bemax = clustertools.eMax(*b);
0243 double be2nd = clustertools.e2nd(*b);
0244 double betop = clustertools.eTop(*b);
0245 double bebottom = clustertools.eBottom(*b);
0246 double beleft = clustertools.eLeft(*b);
0247 double beright = clustertools.eRight(*b);
0248
0249 double be2x5max = clustertools.e2x5Max(*b);
0250 double be2x5top = clustertools.e2x5Top(*b);
0251 double be2x5bottom = clustertools.e2x5Bottom(*b);
0252 double be2x5left = clustertools.e2x5Left(*b);
0253 double be2x5right = clustertools.e2x5Right(*b);
0254
0255 fVals[11] = b->eta() - s->eta();
0256 fVals[12] = reco::deltaPhi(b->phi(), s->phi());
0257 fVals[13] = b->energy() / s->rawEnergy();
0258 fVals[14] = clustertools.e3x3(*b) / b->energy();
0259 fVals[15] = clustertools.e5x5(*b) / b->energy();
0260 fVals[16] = sqrt(clustertools.localCovariances(*b)[0]);
0261 fVals[17] = sqrt(clustertools.localCovariances(*b)[2]);
0262 fVals[18] = clustertools.localCovariances(*b)[1];
0263 fVals[19] = bemax / b->energy();
0264 fVals[20] = be2nd / b->energy();
0265 fVals[21] = betop / b->energy();
0266 fVals[22] = bebottom / b->energy();
0267 fVals[23] = beleft / b->energy();
0268 fVals[24] = beright / b->energy();
0269 fVals[25] = be2x5max / b->energy();
0270 fVals[26] = be2x5top / b->energy();
0271 fVals[27] = be2x5bottom / b->energy();
0272 fVals[28] = be2x5left / b->energy();
0273 fVals[29] = be2x5right / b->energy();
0274
0275 if (isbarrel) {
0276
0277
0278
0279 float betacry, bphicry, bthetatilt, bphitilt;
0280 int bieta, biphi;
0281 egammaTools::localEcalClusterCoordsEB(*b, caloGeometry, betacry, bphicry, bieta, biphi, bthetatilt, bphitilt);
0282
0283 fVals[30] = bieta;
0284 fVals[31] = biphi;
0285 fVals[32] = bieta % 5;
0286 fVals[33] = biphi % 2;
0287 fVals[34] = (std::abs(bieta) <= 25) * (bieta % 25) +
0288 (std::abs(bieta) > 25) *
0289 ((bieta - 25 * std::abs(bieta) / bieta) % 20);
0290 fVals[35] = biphi % 20;
0291 fVals[36] = betacry;
0292 fVals[37] = bphicry;
0293
0294 } else {
0295
0296 fVals[30] = s->preshowerEnergy() / s->rawEnergy();
0297 }
0298
0299
0300
0301
0302
0303
0304 double den;
0305 const GBRForest *reader;
0306 const GBRForest *readervar;
0307 if (isbarrel) {
0308 den = s->rawEnergy();
0309 reader = fReadereb.get();
0310 readervar = fReaderebvariance.get();
0311 } else {
0312 den = s->rawEnergy() + s->preshowerEnergy();
0313 reader = fReaderee.get();
0314 readervar = fReadereevariance.get();
0315 }
0316
0317 double ecor = reader->GetResponse(fVals.data()) * den;
0318
0319
0320 if (applyRescale) {
0321 if (isbarrel) {
0322 fVals[3] = 1.0045 * p.r9() + 0.001;
0323 fVals[5] = 1.04302 * s->etaWidth() - 0.000618;
0324 fVals[6] = 1.00002 * s->phiWidth() - 0.000371;
0325 fVals[14] = fVals[3] * s->rawEnergy() / b->energy();
0326 if (fVals[15] <= 1.0)
0327 fVals[15] = std::min(1.0, 1.0022 * p.e5x5() / b->energy());
0328
0329 fVals[4] =
0330 fVals[15] * b->energy() / s->rawEnergy();
0331
0332 fVals[16] = 0.891832 * sqrt(clustertools.localCovariances(*b)[0]) + 0.0009133;
0333 fVals[17] = 0.993 * sqrt(clustertools.localCovariances(*b)[2]);
0334
0335 fVals[19] = 1.012 * bemax / b->energy();
0336 fVals[20] = 1.0 * be2nd / b->energy();
0337 fVals[21] = 0.94 * betop / b->energy();
0338 fVals[22] = 0.94 * bebottom / b->energy();
0339 fVals[23] = 0.94 * beleft / b->energy();
0340 fVals[24] = 0.94 * beright / b->energy();
0341 fVals[25] = 1.006 * be2x5max / b->energy();
0342 fVals[26] = 1.09 * be2x5top / b->energy();
0343 fVals[27] = 1.09 * be2x5bottom / b->energy();
0344 fVals[28] = 1.09 * be2x5left / b->energy();
0345 fVals[29] = 1.09 * be2x5right / b->energy();
0346
0347 } else {
0348 fVals[3] = 1.0086 * p.r9() - 0.0007;
0349 fVals[4] = std::min(1.0, 1.0022 * p.e5x5() / s->rawEnergy());
0350 fVals[5] = 0.903254 * s->etaWidth() + 0.001346;
0351 fVals[6] = 0.99992 * s->phiWidth() + 4.8e-07;
0352 fVals[13] =
0353 std::min(1.0, 1.0022 * b->energy() / s->rawEnergy());
0354
0355 fVals[14] = fVals[3] * s->rawEnergy() / b->energy();
0356
0357 fVals[16] = 0.9947 * sqrt(clustertools.localCovariances(*b)[0]) + 0.00003;
0358
0359 fVals[19] = 1.005 * bemax / b->energy();
0360 fVals[20] = 1.02 * be2nd / b->energy();
0361 fVals[21] = 0.96 * betop / b->energy();
0362 fVals[22] = 0.96 * bebottom / b->energy();
0363 fVals[23] = 0.96 * beleft / b->energy();
0364 fVals[24] = 0.96 * beright / b->energy();
0365 fVals[25] = 1.0075 * be2x5max / b->energy();
0366 fVals[26] = 1.13 * be2x5top / b->energy();
0367 fVals[27] = 1.13 * be2x5bottom / b->energy();
0368 fVals[28] = 1.13 * be2x5left / b->energy();
0369 fVals[29] = 1.13 * be2x5right / b->energy();
0370 }
0371 }
0372
0373 double ecorerr = readervar->GetResponse(fVals.data()) * den;
0374
0375
0376
0377 return {ecor, ecorerr};
0378 }