EgammaTools/src/EGEnergyCorrector.cc

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();  //seed  basic cluster
0028
0029   //highest energy basic cluster excluding seed basic cluster
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   //lowest energy basic cluster excluding seed (for pileup mitigation)
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   //2nd lowest energy basic cluster excluding seed (for pileup mitigation)
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     //basic supercluster variables
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     //seed basic cluster variables
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]);  //sigietaieta
0092     fVals[14] = sqrt(clustertools.localCovariances(*b)[2]);  //sigiphiiphi
0093     fVals[15] = clustertools.localCovariances(*b)[1];        //sigietaiphi
0094     fVals[16] = bemax / b->energy();                         //crystal energy ratio gap variables
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     //local coordinates and crystal indices
0146
0147     //seed cluster
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;      //crystal ieta
0153     fVals[57] = biphi;      //crystal iphi
0154     fVals[58] = bieta % 5;  //submodule boundary eta symmetry
0155     fVals[59] = biphi % 2;  //submodule boundary phi symmetry
0156     fVals[60] = (std::abs(bieta) <= 25) * (bieta % 25) +
0157                 (std::abs(bieta) > 25) *
0158                     ((bieta - 25 * std::abs(bieta) / bieta) % 20);  //module boundary eta approximate symmetry
0159     fVals[61] = biphi % 20;                                         //module boundary phi symmetry
0160     fVals[62] = betacry;  //local coordinates with respect to closest crystal center at nominal shower depth
0161     fVals[63] = bphicry;
0162
0163     //2nd cluster (meaningful gap corrections for converted photons)
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   //printf("ecor = %5f, ecorerr = %5f\n",ecor,ecorerr);
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();  //seed  basic cluster
0225
0226   bool isbarrel = b->hitsAndFractions().at(0).first.subdetId() == EcalBarrel;
0227
0228   //basic supercluster variables
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   //seed basic cluster variables
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]);  //sigietaieta
0261   fVals[17] = sqrt(clustertools.localCovariances(*b)[2]);  //sigiphiiphi
0262   fVals[18] = clustertools.localCovariances(*b)[1];        //sigietaiphi
0263   fVals[19] = bemax / b->energy();                         //crystal energy ratio gap variables
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();  //crystal energy ratio gap variables
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     //local coordinates and crystal indices (barrel only)
0277
0278     //seed cluster
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;      //crystal ieta
0284     fVals[31] = biphi;      //crystal iphi
0285     fVals[32] = bieta % 5;  //submodule boundary eta symmetry
0286     fVals[33] = biphi % 2;  //submodule boundary phi symmetry
0287     fVals[34] = (std::abs(bieta) <= 25) * (bieta % 25) +
0288                 (std::abs(bieta) > 25) *
0289                     ((bieta - 25 * std::abs(bieta) / bieta) % 20);  //module boundary eta approximate symmetry
0290     fVals[35] = biphi % 20;                                         //module boundary phi symmetry
0291     fVals[36] = betacry;  //local coordinates with respect to closest crystal center at nominal shower depth
0292     fVals[37] = bphicry;
0293
0294   } else {
0295     //preshower energy ratio (endcap only)
0296     fVals[30] = s->preshowerEnergy() / s->rawEnergy();
0297   }
0298
0299   //   if (isbarrel) {
0300   //     for (int i=0; i<38; ++i) printf("%i: %5f\n",i,fVals[i]);
0301   //   }
0302   //   else for (int i=0; i<31; ++i) printf("%i: %5f\n",i,fVals[i]);
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   //apply shower shape rescaling - for Monte Carlo only, and only for calculation of energy uncertainty
0320   if (applyRescale) {
0321     if (isbarrel) {
0322       fVals[3] = 1.0045 * p.r9() + 0.001;                   //r9
0323       fVals[5] = 1.04302 * s->etaWidth() - 0.000618;        //etawidth
0324       fVals[6] = 1.00002 * s->phiWidth() - 0.000371;        //phiwidth
0325       fVals[14] = fVals[3] * s->rawEnergy() / b->energy();  //compute consistent e3x3/eseed after r9 rescaling
0326       if (fVals[15] <= 1.0)  // rescale e5x5/eseed only if value is <=1.0, don't allow scaled values to exceed 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();  // compute consistent e5x5()/rawEnergy() after e5x5/eseed resacling
0331
0332       fVals[16] = 0.891832 * sqrt(clustertools.localCovariances(*b)[0]) + 0.0009133;  //sigietaieta
0333       fVals[17] = 0.993 * sqrt(clustertools.localCovariances(*b)[2]);                 //sigiphiiphi
0334
0335       fVals[19] = 1.012 * bemax / b->energy();  //crystal energy ratio gap variables
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();  //crystal energy ratio gap variables
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;                           //r9
0349       fVals[4] = std::min(1.0, 1.0022 * p.e5x5() / s->rawEnergy());  //e5x5/rawenergy
0350       fVals[5] = 0.903254 * s->etaWidth() + 0.001346;                //etawidth
0351       fVals[6] = 0.99992 * s->phiWidth() + 4.8e-07;                  //phiwidth
0352       fVals[13] =
0353           std::min(1.0, 1.0022 * b->energy() / s->rawEnergy());  //eseed/rawenergy (practically equivalent to e5x5)
0354
0355       fVals[14] = fVals[3] * s->rawEnergy() / b->energy();  //compute consistent e3x3/eseed after r9 rescaling
0356
0357       fVals[16] = 0.9947 * sqrt(clustertools.localCovariances(*b)[0]) + 0.00003;  //sigietaieta
0358
0359       fVals[19] = 1.005 * bemax / b->energy();  //crystal energy ratio gap variables
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();  //crystal energy ratio gap variables
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   //printf("ecor = %5f, ecorerr = %5f\n",ecor,ecorerr);
0376
0377   return {ecor, ecorerr};
0378 }