File indexing completed on 2024-04-06 12:24:38
0001
0002
0003
0004
0005 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0006 #include "RecoEcal/EgammaClusterAlgos/interface/EgammaSCEnergyCorrectionAlgo.h"
0007 #include "RecoEcal/EgammaCoreTools/interface/SuperClusterShapeAlgo.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009
0010 #include <string>
0011 #include <vector>
0012
0013 EgammaSCEnergyCorrectionAlgo::EgammaSCEnergyCorrectionAlgo(float noise) : sigmaElectronicNoise_{noise} {}
0014
0015 reco::SuperCluster EgammaSCEnergyCorrectionAlgo::applyCorrection(const reco::SuperCluster& cl,
0016 const EcalRecHitCollection& rhc,
0017 reco::CaloCluster::AlgoId theAlgo,
0018 const CaloSubdetectorGeometry* geometry,
0019 EcalClusterFunctionBaseClass* energyCorrectionFunction,
0020 std::string energyCorrectorName_,
0021 const int modeEB_,
0022 const int modeEE_) {
0023
0024
0025 {
0026 LogTrace("EgammaSCEnergyCorrectionAlgo") << "::applyCorrection" << std::endl;
0027 LogTrace() << " SC has energy " << cl.energy() << std::endl;
0028 LogTrace() << " Will correct now.... " << std::endl;
0029 }
0030
0031
0032 const reco::CaloClusterPtr& seedC = cl.seed();
0033
0034 LogTrace("EgammaSCEnergyCorrectionAlgo") << " Seed cluster energy... " << seedC->energy() << std::endl;
0035
0036
0037 LogTrace("EgammaSCEnergyCorrectionAlgo") << " The seed cluster used algo " << theAlgo;
0038
0039
0040
0041 std::vector<std::pair<DetId, float> > seedHits = seedC->hitsAndFractions();
0042 EcalSubdetector theBase = EcalSubdetector(seedHits.at(0).first.subdetId());
0043
0044 LogTrace("EgammaSCEnergyCorrectionAlgo") << " seed cluster location == " << theBase << std::endl;
0045
0046
0047 int nCryGT2Sigma = nCrystalsGT2Sigma(*seedC, rhc);
0048
0049 LogTrace("EgammaSCEnergyCorrectionAlgo") << " nCryGT2Sigma " << nCryGT2Sigma << std::endl;
0050
0051
0052 float bremsEnergy = cl.energy() - seedC->energy();
0053
0054 LogTrace("EgammaSCEnergyCorrectionAlgo") << " bremsEnergy " << bremsEnergy << std::endl;
0055
0056
0057
0058 SuperClusterShapeAlgo SCShape(&rhc, geometry);
0059
0060 double phiWidth = 0.;
0061 double etaWidth = 0.;
0062
0063 SCShape.Calculate_Covariances(cl);
0064 phiWidth = SCShape.phiWidth();
0065 etaWidth = SCShape.etaWidth();
0066
0067
0068
0069
0070 float newEnergy = 0;
0071
0072 reco::SuperCluster tmp = cl;
0073 tmp.setPhiWidth(phiWidth);
0074 tmp.setEtaWidth(etaWidth);
0075
0076 if (theAlgo == reco::CaloCluster::hybrid || theAlgo == reco::CaloCluster::dynamicHybrid) {
0077 if (energyCorrectorName_ == "EcalClusterEnergyCorrection")
0078 newEnergy = tmp.rawEnergy() + energyCorrectionFunction->getValue(tmp, modeEB_);
0079 if (energyCorrectorName_ == "EcalClusterEnergyCorrectionObjectSpecific") {
0080
0081 newEnergy = energyCorrectionFunction->getValue(tmp, modeEB_);
0082 }
0083
0084 } else if (theAlgo == reco::CaloCluster::multi5x5) {
0085 if (energyCorrectorName_ == "EcalClusterEnergyCorrection")
0086 newEnergy = tmp.rawEnergy() + tmp.preshowerEnergy() + energyCorrectionFunction->getValue(tmp, modeEE_);
0087 if (energyCorrectorName_ == "EcalClusterEnergyCorrectionObjectSpecific")
0088 newEnergy = energyCorrectionFunction->getValue(tmp, modeEE_);
0089
0090 } else {
0091
0092 newEnergy = seedC->energy() / fNCrystals(nCryGT2Sigma, theAlgo, theBase) + bremsEnergy;
0093 }
0094
0095
0096
0097 LogTrace("EgammaSCEnergyCorrectionAlgo") << " UNCORRECTED SC has energy... " << cl.energy() << std::endl;
0098 LogTrace("EgammaSCEnergyCorrectionAlgo") << " CORRECTED SC has energy... " << newEnergy << std::endl;
0099
0100 reco::SuperCluster corrCl = cl;
0101
0102 corrCl.setEnergy(newEnergy);
0103 corrCl.setPhiWidth(phiWidth);
0104 corrCl.setEtaWidth(etaWidth);
0105
0106 return corrCl;
0107 }
0108
0109 float EgammaSCEnergyCorrectionAlgo::fNCrystals(int nCry,
0110 reco::CaloCluster::AlgoId theAlgo,
0111 EcalSubdetector theBase) const {
0112 float p0 = 0, p1 = 0, p2 = 0, p3 = 0, p4 = 0;
0113 float x = nCry;
0114 float result = 1.f;
0115
0116 if ((theBase == EcalBarrel) && (theAlgo == reco::CaloCluster::hybrid)) {
0117 if (nCry <= 10) {
0118 p0 = 6.32879e-01f;
0119 p1 = 1.14893e-01f;
0120 p2 = -2.45705e-02f;
0121 p3 = 2.53074e-03f;
0122 p4 = -9.29654e-05f;
0123 } else if (nCry > 10 && nCry <= 30) {
0124 p0 = 6.93196e-01f;
0125 p1 = 4.44034e-02f;
0126 p2 = -2.82229e-03f;
0127 p3 = 8.19495e-05f;
0128 p4 = -8.96645e-07f;
0129 } else {
0130 p0 = 5.65474e+00f;
0131 p1 = -6.31640e-01f;
0132 p2 = 3.14218e-02f;
0133 p3 = -6.84256e-04f;
0134 p4 = 5.50659e-06f;
0135 }
0136 if (x > 40.f)
0137 x = 40.f;
0138 }
0139
0140 else if ((theBase == EcalEndcap) && (theAlgo == reco::CaloCluster::hybrid)) {
0141 LogTrace("EgammaSCEnergyCorrectionAlgo") << "ERROR! HybridEFRYsc called" << std::endl;
0142
0143 return 1.f;
0144 }
0145
0146 else if ((theBase == EcalBarrel) && (theAlgo == reco::CaloCluster::island)) {
0147 p0 = 4.69976e-01f;
0148 p1 = 1.45900e-01f;
0149 p2 = -1.61359e-02f;
0150 p3 = 7.99423e-04f;
0151 p4 = -1.47873e-05f;
0152 if (x > 16.f)
0153 x = 16.f;
0154 }
0155
0156 else if ((theBase == EcalEndcap) && (theAlgo == reco::CaloCluster::island)) {
0157 p0 = 4.69976e-01f;
0158 p1 = 1.45900e-01f;
0159 p2 = -1.61359e-02f;
0160 p3 = 7.99423e-04f;
0161 p4 = -1.47873e-05f;
0162 if (x > 16.f)
0163 x = 16.f;
0164 }
0165
0166 else {
0167 LogTrace("EgammaSCEnergyCorrectionAlgo") << "trying to correct unknown cluster!!!" << std::endl;
0168 return 1.f;
0169 }
0170 result = p0 + x * (p1 + x * (p2 + x * (p3 + x * p4)));
0171
0172
0173
0174 float const ebfact = 1.f / 0.965f;
0175 float const eefact = 1.f / 0.975f;
0176
0177 if (theBase == EcalBarrel) {
0178 result *= ebfact;
0179 } else {
0180 result *= eefact;
0181 }
0182
0183 return result;
0184 }
0185
0186 int EgammaSCEnergyCorrectionAlgo::nCrystalsGT2Sigma(reco::BasicCluster const& seed,
0187 EcalRecHitCollection const& rhc) const {
0188
0189
0190
0191 std::vector<std::pair<DetId, float> > const& hits = seed.hitsAndFractions();
0192
0193 LogTrace("EgammaSCEnergyCorrectionAlgo") << " EgammaSCEnergyCorrectionAlgo::nCrystalsGT2Sigma" << std::endl;
0194 LogTrace("EgammaSCEnergyCorrectionAlgo") << " Will calculate number of crystals above 2sigma noise" << std::endl;
0195 LogTrace("EgammaSCEnergyCorrectionAlgo") << " sigmaElectronicNoise = " << sigmaElectronicNoise_ << std::endl;
0196 LogTrace("EgammaSCEnergyCorrectionAlgo") << " There are " << hits.size() << " recHits" << std::endl;
0197
0198 int nCry = 0;
0199 std::vector<std::pair<DetId, float> >::const_iterator hit;
0200 EcalRecHitCollection::const_iterator aHit;
0201 for (hit = hits.begin(); hit != hits.end(); hit++) {
0202
0203 aHit = rhc.find((*hit).first);
0204 if ((*aHit).energy() > 2.f * sigmaElectronicNoise_)
0205 nCry++;
0206 }
0207
0208 LogTrace("EgammaSCEnergyCorrectionAlgo") << " " << nCry << " of these above 2sigma noise" << std::endl;
0209
0210 return nCry;
0211 }
0212
0213
0214
0215 reco::SuperCluster EgammaSCEnergyCorrectionAlgo::applyCrackCorrection(
0216 const reco::SuperCluster& cl, EcalClusterFunctionBaseClass* crackCorrectionFunction) {
0217 double crackcor = 1.;
0218
0219 for (reco::CaloCluster_iterator cIt = cl.clustersBegin(); cIt != cl.clustersEnd(); ++cIt) {
0220 const reco::CaloClusterPtr cc = *cIt;
0221 crackcor *= ((cl.rawEnergy() + cc->energy() * (crackCorrectionFunction->getValue(*cc) - 1.)) / cl.rawEnergy());
0222 }
0223
0224 reco::SuperCluster corrCl = cl;
0225 corrCl.setEnergy(cl.energy() * crackcor);
0226
0227 return corrCl;
0228 }
0229
0230
0231
0232
0233 reco::SuperCluster EgammaSCEnergyCorrectionAlgo::applyLocalContCorrection(
0234 reco::SuperCluster const& cl, BasicClusterFunction localContCorrectionFunction) {
0235 const EcalRecHitCollection dummy;
0236
0237 const reco::CaloClusterPtr& seedBC = cl.seed();
0238 float seedBCene = seedBC->energy();
0239 float correctedSeedBCene = localContCorrectionFunction(*seedBC, dummy) * seedBCene;
0240
0241 reco::SuperCluster correctedSC = cl;
0242 correctedSC.setEnergy(cl.energy() - seedBCene + correctedSeedBCene);
0243
0244 return correctedSC;
0245 }