File indexing completed on 2024-04-06 12:25:17
0001
0002
0003
0004 #include "RecoHI/HiEgammaAlgos/interface/HiEgammaSCEnergyCorrectionAlgo.h"
0005 #include "RecoEcal/EgammaCoreTools/interface/SuperClusterShapeAlgo.h"
0006 #include <iostream>
0007 #include <string>
0008 #include <vector>
0009
0010 HiEgammaSCEnergyCorrectionAlgo::HiEgammaSCEnergyCorrectionAlgo(float noise,
0011 const edm::ParameterSet& pSet,
0012 HiEgammaSCEnergyCorrectionAlgo::VerbosityLevel verbosity)
0013 : sigmaElectronicNoise_(noise),
0014 verbosity_(verbosity),
0015
0016
0017 p_fEta_(pSet.getParameter<std::vector<double> >("fEtaVect")),
0018 p_fBremTh_(pSet.getParameter<std::vector<double> >("fBremThVect")),
0019 p_fBrem_(pSet.getParameter<std::vector<double> >("fBremVect")),
0020 p_fEtEta_(pSet.getParameter<std::vector<double> >("fEtEtaVect")),
0021
0022
0023 minR9Barrel_(pSet.getParameter<double>("minR9Barrel")),
0024 minR9Endcap_(pSet.getParameter<double>("minR9Endcap")),
0025
0026
0027 maxR9_(pSet.getParameter<double>("maxR9")) {}
0028
0029 reco::SuperCluster HiEgammaSCEnergyCorrectionAlgo::applyCorrection(const reco::SuperCluster& cl,
0030 const EcalRecHitCollection& rhc,
0031 reco::CaloCluster::AlgoId algoId,
0032 const CaloSubdetectorGeometry& geometry,
0033 const CaloTopology& topology) const {
0034
0035 if (verbosity_ <= pINFO) {
0036 std::cout << " HiEgammaSCEnergyCorrectionAlgo::applyCorrection" << std::endl;
0037 std::cout << " SC has energy " << cl.energy() << std::endl;
0038 std::cout << " Will correct now.... " << std::endl;
0039 }
0040
0041
0042 const reco::CaloClusterPtr& seedC = cl.seed();
0043
0044 if (verbosity_ <= pINFO) {
0045 std::cout << " Seed cluster energy... " << seedC->energy() << std::endl;
0046 }
0047
0048
0049 reco::CaloClusterPtrVector clusters_v;
0050
0051 if (verbosity_ <= pINFO)
0052 std::cout << " Constituent cluster energies... ";
0053
0054 for (reco::CaloCluster_iterator cluster = cl.clustersBegin(); cluster != cl.clustersEnd(); cluster++) {
0055 clusters_v.push_back(*cluster);
0056 if (verbosity_ <= pINFO)
0057 std::cout << (*cluster)->energy() << ", ";
0058 }
0059 if (verbosity_ <= pINFO)
0060 std::cout << std::endl;
0061
0062
0063 if (verbosity_ <= pINFO) {
0064 std::cout << " The seed cluster used algo " << algoId;
0065 }
0066
0067
0068
0069 std::vector<std::pair<DetId, float> > const& seedHits = seedC->hitsAndFractions();
0070 EcalSubdetector theBase = EcalSubdetector(seedHits.at(0).first.subdetId());
0071 if (verbosity_ <= pINFO) {
0072 std::cout << " seed cluster location == " << theBase << std::endl;
0073 }
0074
0075
0076 int nCryGT2Sigma = nCrystalsGT2Sigma(*seedC, rhc);
0077 if (verbosity_ <= pINFO) {
0078 std::cout << " nCryGT2Sigma " << nCryGT2Sigma << std::endl;
0079 }
0080
0081
0082 float bremsEnergy = cl.energy() - seedC->energy();
0083 if (verbosity_ <= pINFO) {
0084 std::cout << " bremsEnergy " << bremsEnergy << std::endl;
0085 }
0086
0087
0088
0089 SuperClusterShapeAlgo SCShape(&rhc, &geometry);
0090
0091 double phiWidth = 0.;
0092 double etaWidth = 0.;
0093
0094
0095 SCShape.Calculate_Covariances(cl);
0096 phiWidth = SCShape.phiWidth();
0097 etaWidth = SCShape.etaWidth();
0098
0099
0100 float e3x3 = EcalClusterTools::e3x3(*(cl.seed()), &rhc, &topology);
0101 float e5x5 = EcalClusterTools::e5x5(*(cl.seed()), &rhc, &topology);
0102 float r9 = e3x3 / cl.rawEnergy();
0103
0104
0105
0106
0107 float newEnergy = 0;
0108
0109
0110 if ((r9 < minR9Barrel_ && theBase == EcalBarrel) || (r9 < minR9Endcap_ && theBase == EcalEndcap)) {
0111
0112 newEnergy = (cl.rawEnergy()) / fEta(cl.eta(), algoId, theBase) / fBrem(phiWidth / etaWidth, algoId, theBase) /
0113 fEtEta(cl.energy() / cosh(cl.eta()), cl.eta(), algoId, theBase);
0114
0115 } else {
0116 if (r9 < maxR9_) {
0117
0118 newEnergy = e5x5 / fEta(cl.eta(), algoId, theBase);
0119 } else {
0120
0121 newEnergy = cl.rawEnergy();
0122 }
0123 }
0124
0125 if (newEnergy > 2 * cl.rawEnergy())
0126 newEnergy = cl.rawEnergy();
0127
0128
0129 if (verbosity_ <= pINFO) {
0130 std::cout << " UNCORRECTED SC has energy... " << cl.energy() << std::endl;
0131 std::cout << " CORRECTED SC has energy... " << newEnergy << std::endl;
0132 std::cout << " Size..." << cl.size() << std::endl;
0133 std::cout << " Seed nCryGT2Sigma Size..." << nCryGT2Sigma << std::endl;
0134 }
0135
0136 reco::SuperCluster corrCl(newEnergy,
0137 math::XYZPoint(cl.position().X(), cl.position().Y(), cl.position().Z()),
0138 cl.seed(),
0139 clusters_v,
0140 cl.preshowerEnergy());
0141
0142
0143 corrCl.setFlags(cl.flags());
0144 corrCl.setPhiWidth(phiWidth);
0145 corrCl.setEtaWidth(etaWidth);
0146
0147 return corrCl;
0148 }
0149
0150 float HiEgammaSCEnergyCorrectionAlgo::fEtEta(float et,
0151 float eta,
0152 reco::CaloCluster::AlgoId algoId,
0153 EcalSubdetector theBase) const {
0154 int offset = 0;
0155 float factor;
0156 if ((theBase == EcalBarrel) && (algoId == reco::CaloCluster::island)) {
0157 offset = 0;
0158 } else if ((theBase == EcalEndcap) && (algoId == reco::CaloCluster::island)) {
0159 offset = 7;
0160 }
0161
0162
0163 factor = (p_fEtEta_[0 + offset] + p_fEtEta_[1 + offset] * sqrt(et));
0164
0165 factor *= (p_fEtEta_[2 + offset] + p_fEtEta_[3 + offset] * fabs(eta) + p_fEtEta_[4 + offset] * eta * eta +
0166 p_fEtEta_[5 + offset] * eta * eta * fabs(eta) + +p_fEtEta_[6 + offset] * eta * eta * eta * eta);
0167
0168
0169 if (factor < 0.66f)
0170 factor = 0.66f;
0171 if (factor > 1.5f)
0172 factor = 1.5f;
0173
0174 return factor;
0175 }
0176
0177 float HiEgammaSCEnergyCorrectionAlgo::fEta(float eta, reco::CaloCluster::AlgoId algoId, EcalSubdetector theBase) const {
0178 int offset = 0;
0179 float factor;
0180 if ((theBase == EcalBarrel) && (algoId == reco::CaloCluster::island)) {
0181 offset = 0;
0182 } else if ((theBase == EcalEndcap) && (algoId == reco::CaloCluster::island)) {
0183 offset = 3;
0184 }
0185
0186 factor = (p_fEta_[0 + offset] + p_fEta_[1 + offset] * fabs(eta) + p_fEta_[2 + offset] * eta * eta);
0187
0188
0189 if (factor < 0.66f)
0190 factor = 0.66f;
0191 if (factor > 1.5f)
0192 factor = 1.5f;
0193
0194 return factor;
0195 }
0196
0197 float HiEgammaSCEnergyCorrectionAlgo::fBrem(float brem,
0198 reco::CaloCluster::AlgoId algoId,
0199 EcalSubdetector theBase) const {
0200 int det = 0;
0201 int offset = 0;
0202 float factor;
0203 if ((theBase == EcalBarrel) && (algoId == reco::CaloCluster::island)) {
0204 det = 0;
0205 offset = 0;
0206 } else if ((theBase == EcalEndcap) && (algoId == reco::CaloCluster::island)) {
0207 det = 1;
0208 offset = 6;
0209 }
0210
0211 if (brem < p_fBremTh_[det]) {
0212 factor = p_fBrem_[0 + offset] + p_fBrem_[1 + offset] * brem + p_fBrem_[2 + offset] * brem * brem;
0213 } else {
0214 factor = p_fBrem_[3 + offset] + p_fBrem_[4 + offset] * brem + p_fBrem_[5 + offset] * brem * brem;
0215 };
0216
0217
0218 if (factor < 0.66f)
0219 factor = 0.66f;
0220 if (factor > 1.5f)
0221 factor = 1.5f;
0222
0223 return factor;
0224 }
0225
0226
0227
0228 float HiEgammaSCEnergyCorrectionAlgo::fNCrystals(int nCry,
0229 reco::CaloCluster::AlgoId algoId,
0230 EcalSubdetector theBase) const {
0231 float x = (float)nCry;
0232 float result = 1.f;
0233
0234 if ((theBase == EcalBarrel) && (algoId == reco::CaloCluster::island)) {
0235 float const p0 = 0.682554f;
0236 float const p1 = 0.0253013f;
0237 float const p2 = -0.0007907f;
0238 float const p3 = 1.166e-5f;
0239 float const p4 = -6.7387e-8f;
0240 if (x < 10.f)
0241 x = 10.f;
0242 if (x < 40.f)
0243 result = p0 + x * (p1 + x * (p2 + x * (p3 + x * p4)));
0244 else
0245 result = 1.f;
0246 }
0247
0248 else if ((theBase == EcalEndcap) && (algoId == reco::CaloCluster::island)) {
0249 float const p0 = 0.712185f;
0250 float const p1 = 0.0273609f;
0251 float const p2 = -0.00103818f;
0252 float const p3 = 2.01828e-05f;
0253 float const p4 = -1.71438e-07f;
0254 if (x < 10.f)
0255 x = 10.f;
0256 if (x < 40.f)
0257 result = p0 + x * (p1 + x * (p2 + x * (p3 + x * p4)));
0258 else
0259 result = 1.f;
0260 }
0261
0262 else {
0263 if (verbosity_ <= pINFO) {
0264 std::cout << "trying to correct unknown cluster!!!" << std::endl;
0265 }
0266 }
0267
0268 if (result > 1.5f)
0269 result = 1.5f;
0270 if (result < 0.5f)
0271 result = 0.5f;
0272
0273 return result;
0274 }
0275
0276 int HiEgammaSCEnergyCorrectionAlgo::nCrystalsGT2Sigma(reco::BasicCluster const& seed,
0277 EcalRecHitCollection const& rhc) const {
0278
0279
0280
0281 std::vector<std::pair<DetId, float> > const& hits = seed.hitsAndFractions();
0282
0283 if (verbosity_ <= pINFO) {
0284 std::cout << " HiEgammaSCEnergyCorrectionAlgo::nCrystalsGT2Sigma" << std::endl;
0285 std::cout << " Will calculate number of crystals above 2sigma noise" << std::endl;
0286 std::cout << " sigmaElectronicNoise = " << sigmaElectronicNoise_ << std::endl;
0287 std::cout << " There are " << hits.size() << " recHits" << std::endl;
0288 }
0289
0290 int nCry = 0;
0291 for (std::vector<std::pair<DetId, float> >::const_iterator hit = hits.begin(); hit != hits.end(); ++hit) {
0292
0293 EcalRecHitCollection::const_iterator aHit = rhc.find((*hit).first);
0294
0295 if (aHit->energy() > 2.f * sigmaElectronicNoise_)
0296 ++nCry;
0297 }
0298
0299 if (verbosity_ <= pINFO) {
0300 std::cout << " " << nCry << " of these above 2sigma noise" << std::endl;
0301 }
0302
0303 return nCry;
0304 }