EgammaClusterAlgos/src/EgammaSCEnergyCorrectionAlgo.cc

0001 //
0002 // Author: David Evans, Bristol
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   // A little bit of trivial info to be sure all is well
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   // Get the seed cluster
0032   const reco::CaloClusterPtr& seedC = cl.seed();
0033
0034   LogTrace("EgammaSCEnergyCorrectionAlgo") << "   Seed cluster energy... " << seedC->energy() << std::endl;
0035
0036   // Find the algorithm used to construct the basic clusters making up the supercluster
0037   LogTrace("EgammaSCEnergyCorrectionAlgo") << "   The seed cluster used algo " << theAlgo;
0038
0039   // Find the detector region of the supercluster
0040   // where is the seed cluster?
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   // Get number of crystals 2sigma above noise in seed basiccluster
0047   int nCryGT2Sigma = nCrystalsGT2Sigma(*seedC, rhc);
0048
0049   LogTrace("EgammaSCEnergyCorrectionAlgo") << "   nCryGT2Sigma " << nCryGT2Sigma << std::endl;
0050
0051   // Supercluster enery - seed basiccluster energy
0052   float bremsEnergy = cl.energy() - seedC->energy();
0053
0054   LogTrace("EgammaSCEnergyCorrectionAlgo") << "   bremsEnergy " << bremsEnergy << std::endl;
0055
0056   //Create the pointer ot class SuperClusterShapeAlgo
0057   //which calculates phiWidth and etaWidth
0058   SuperClusterShapeAlgo SCShape(&rhc, geometry);
0059
0060   double phiWidth = 0.;
0061   double etaWidth = 0.;
0062   //Calculate phiWidth & etaWidth for SuperClusters
0063   SCShape.Calculate_Covariances(cl);
0064   phiWidth = SCShape.phiWidth();
0065   etaWidth = SCShape.etaWidth();
0066
0067   // Calculate the new supercluster energy
0068   //as a function of number of crystals in the seed basiccluster for Endcap
0069   //or apply new Enegry SCale correction
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       //std::cout << "newEnergy="<<newEnergy<<std::endl;
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     //Apply f(nCry) correction on island algo and fixedMatrix algo
0092     newEnergy = seedC->energy() / fNCrystals(nCryGT2Sigma, theAlgo, theBase) + bremsEnergy;
0093   }
0094
0095   // Create a new supercluster with the corrected energy
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;  // extracted from fit to all endcap classes with Ebrem = 0.
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;  // extracted from fit to all endcap classes with Ebrem = 0.
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   //Rescale energy scale correction to take into account change in calibrated
0173   //RecHit definition introduced in CMSSW_1_5_0
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   // return number of crystals 2Sigma above
0189   // electronic noise
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     // need to get hit by DetID in order to get energy
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 // apply crack correction
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   }  // loop on BCs
0223
0224   reco::SuperCluster corrCl = cl;
0225   corrCl.setEnergy(cl.energy() * crackcor);
0226
0227   return corrCl;
0228 }
0229
0230 // apply local containment correction
0231 // Assume that the correction function provides correction for the seed Basic Cluster
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 }