HiEgammaAlgos/src/HiEgammaSCEnergyCorrectionAlgo.cc

0001 //
0002 // Author: David Evans, Bristol
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       // Parameters
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       // Min R9
0023       minR9Barrel_(pSet.getParameter<double>("minR9Barrel")),
0024       minR9Endcap_(pSet.getParameter<double>("minR9Endcap")),
0025
0026       // Max R9
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   // Print out a little bit of trivial info to be sure all is well
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   // Get the seed cluster
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   // Get the constituent clusters
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   // Find the algorithm used to construct the basic clusters making up the supercluster
0063   if (verbosity_ <= pINFO) {
0064     std::cout << "   The seed cluster used algo " << algoId;
0065   }
0066
0067   // Find the detector region of the supercluster
0068   // where is the seed cluster?
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   // Get number of crystals 2sigma above noise in seed basiccluster
0076   int nCryGT2Sigma = nCrystalsGT2Sigma(*seedC, rhc);
0077   if (verbosity_ <= pINFO) {
0078     std::cout << "   nCryGT2Sigma " << nCryGT2Sigma << std::endl;
0079   }
0080
0081   // Supercluster enery - seed basiccluster energy
0082   float bremsEnergy = cl.energy() - seedC->energy();
0083   if (verbosity_ <= pINFO) {
0084     std::cout << "   bremsEnergy " << bremsEnergy << std::endl;
0085   }
0086
0087   // Create a SuperClusterShapeAlgo
0088   // which calculates phiWidth and etaWidth
0089   SuperClusterShapeAlgo SCShape(&rhc, &geometry);
0090
0091   double phiWidth = 0.;
0092   double etaWidth = 0.;
0093
0094   // Calculate phiWidth & etaWidth for SuperClusters
0095   SCShape.Calculate_Covariances(cl);
0096   phiWidth = SCShape.phiWidth();
0097   etaWidth = SCShape.etaWidth();
0098
0099   // Calculate r9 and 5x5 energy
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   // Calculate the new supercluster energy
0105   // as a function of number of crystals in the seed basiccluster for Endcap
0106   // lslslsor apply new Enegry SCale correction
0107   float newEnergy = 0;
0108
0109   // if r9 > maxR9_ -> uncaptured brem.
0110   if ((r9 < minR9Barrel_ && theBase == EcalBarrel) || (r9 < minR9Endcap_ && theBase == EcalEndcap)) {
0111     // if r9 is greater than threshold, then use the SC raw energy
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       // use 5x5 energy if r9 < threshold
0118       newEnergy = e5x5 / fEta(cl.eta(), algoId, theBase);
0119     } else {
0120       // it comes from a uncaptured brem, doesn't correct
0121       newEnergy = cl.rawEnergy();
0122     }
0123   }
0124
0125   if (newEnergy > 2 * cl.rawEnergy())
0126     newEnergy = cl.rawEnergy();  // avoid very large correction due to the uncaptured brem
0127
0128   // Create a new supercluster with the corrected energy
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   //set the flags, although we should implement a ctor in SuperCluster
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   // Et dependent correction
0163   factor = (p_fEtEta_[0 + offset] + p_fEtEta_[1 + offset] * sqrt(et));
0164   // eta dependent correction
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   // Constraint correction factor
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   // Constraint correction factor
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   // Constraint correction factor
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 //   char *var ="rawEnergy/cosh(genMatchedEta)/(1.01606-0.0162668*abs(eta))/genMatchedPt/(1.022-0.02812*phiWidth/etaWidth+0.001637*phiWidth*phiWidth/etaWidth/etaWidth)/((0.682554+0.0253013*scSize-(0.0007907)*scSize*scSize+(1.166e-5)*scSize*scSize*scSize-(6.7387e-8)*scSize*scSize*scSize*scSize)*(scSize<40)+(scSize>=40))/((1.016-0.009877*((clustersSize<=4)*clustersSize+(clustersSize>4)*4)))";
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   // return number of crystals 2Sigma above
0279   // electronic noise
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     // need to get hit by DetID in order to get energy
0293     EcalRecHitCollection::const_iterator aHit = rhc.find((*hit).first);
0294     // better the hit to exists....
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 }