EgammaHFProducers/plugins/HFRecoEcalCandidateAlgo.cc

0001 /** \class HFRecoEcalCandidateProducers
0002  *
0003  *  \author Kevin Klapoetke (Minnesota)
0004  *
0005  * $Id:
0006  *
0007  */
0008
0009 #include "DataFormats/Math/interface/LorentzVector.h"
0010 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
0011 #include "HFRecoEcalCandidateAlgo.h"
0012 #include "HFEGammaSLCorrector.h"
0013
0014 using namespace std;
0015 using namespace reco;
0016
0017 HFRecoEcalCandidateAlgo::HFRecoEcalCandidateAlgo(bool correct,
0018                                                  double e9e25Cut,
0019                                                  double intercept2DCut,
0020                                                  double intercept2DSlope,
0021                                                  const std::vector<double>& e1e9Cut,
0022                                                  const std::vector<double>& eCOREe9Cut,
0023                                                  const std::vector<double>& eSeLCut,
0024                                                  const reco::HFValueStruct hfvv)
0025     : m_correct(correct),
0026       m_e9e25Cut(e9e25Cut),
0027       m_intercept2DCut(intercept2DCut),
0028       m_intercept2DSlope(intercept2DSlope),
0029       m_e1e9Cuthi(e1e9Cut[1]),
0030       m_eCOREe9Cuthi(eCOREe9Cut[1]),
0031       m_eSeLCuthi(eSeLCut[1]),
0032       m_e1e9Cutlo(e1e9Cut[0]),
0033       m_eCOREe9Cutlo(eCOREe9Cut[0]),
0034       m_eSeLCutlo(eSeLCut[0]),
0035       m_era(4),
0036       m_hfvv(hfvv) {}
0037
0038 RecoEcalCandidate HFRecoEcalCandidateAlgo::correctEPosition(const SuperCluster& original,
0039                                                             const HFEMClusterShape& shape,
0040                                                             int nvtx) const {
0041   double corEta = original.eta();
0042   //piece-wise log energy correction to eta
0043   double logel = log(shape.eLong3x3() / 100.0);
0044   double lowC = 0.00911;
0045   double hiC = 0.0193;
0046   double logSlope = 0.0146;
0047   double logIntercept = -0.00988;
0048   double sgn = (original.eta() > 0) ? (1.0) : (-1.0);
0049   if (logel <= 1.25)
0050     corEta += (lowC)*sgn;
0051   else if ((logel > 1.25) && (logel <= 2))
0052     corEta += (logIntercept + logSlope * logel) * sgn;
0053   else if (logel > 2)
0054     corEta += hiC * sgn;
0055   //poly-3 cell eta correction to eta (de)
0056   double p0 = 0.0125;
0057   double p1 = -0.0475;
0058   double p2 = 0.0732;
0059   double p3 = -0.0506;
0060   double de = sgn * (p0 + (p1 * (shape.CellEta())) + (p2 * (shape.CellEta() * shape.CellEta())) +
0061                      (p3 * (shape.CellEta() * shape.CellEta() * shape.CellEta())));
0062   corEta += de;
0063
0064   double corPhi = original.phi();
0065   //poly-3 cell phi correction to phi (dp)
0066   p0 = 0.0115;
0067   p1 = -0.0394;
0068   p2 = 0.0486;
0069   p3 = -0.0335;
0070   double dp = p0 + (p1 * (shape.CellPhi())) + (p2 * (shape.CellPhi() * shape.CellPhi())) +
0071               (p3 * (shape.CellPhi() * shape.CellPhi() * shape.CellPhi()));
0072   corPhi += dp;
0073
0074   double corEnergy = original.energy();
0075   //energy corrections
0076   //flat correction for using only long fibers
0077   double energyCorrect = 0.7397;
0078   corEnergy = corEnergy / energyCorrect;
0079   //corection based on ieta for pileup and general
0080   //find ieta
0081   int ieta = 0;
0082   double etabounds[30] = {
0083       2.846, 2.957, 3.132, 3.307, 3.482, 3.657, 3.833, 4.006, 4.184, 4.357, 4.532, 4.709, 4.882, 5.184};
0084   for (int kk = 0; kk <= 12; kk++) {
0085     if ((fabs(corEta) < etabounds[kk + 1]) && (fabs(corEta) > etabounds[kk])) {
0086       ieta = (corEta > 0) ? (kk + 29) : (-kk - 29);
0087     }
0088   }
0089   //check if ieta is in bounds, should be [-41, -29] U [29, 41]
0090   if (ieta != 0)
0091     corEnergy = (m_hfvv.PUSlope(ieta) * 1.0 * (nvtx - 1) + m_hfvv.PUIntercept(ieta) * 1.0) * corEnergy * 1.0 *
0092                 m_hfvv.EnCor(ieta);
0093
0094   //re-calculate the energy vector
0095   double corPx = corEnergy * cos(corPhi) / cosh(corEta);
0096   double corPy = corEnergy * sin(corPhi) / cosh(corEta);
0097   double corPz = corEnergy * tanh(corEta);
0098
0099   //make the candidate
0100   RecoEcalCandidate corCand(0, math::XYZTLorentzVector(corPx, corPy, corPz, corEnergy), math::XYZPoint(0, 0, 0));
0101
0102   return corCand;
0103 }
0104
0105 void HFRecoEcalCandidateAlgo::produce(const edm::Handle<SuperClusterCollection>& SuperClusters,
0106                                       const HFEMClusterShapeAssociationCollection& AssocShapes,
0107                                       RecoEcalCandidateCollection& RecoECand,
0108                                       int nvtx) const {
0109   //get super's and cluster shapes and associations
0110   for (unsigned int i = 0; i < SuperClusters->size(); ++i) {
0111     const SuperCluster& supClus = (*SuperClusters)[i];
0112     reco::SuperClusterRef theClusRef = edm::Ref<SuperClusterCollection>(SuperClusters, i);
0113     const HFEMClusterShapeRef clusShapeRef = AssocShapes.find(theClusRef)->val;
0114     const HFEMClusterShape& clusShape = *clusShapeRef;
0115
0116     // basic candidate
0117     double px = supClus.energy() * cos(supClus.phi()) / cosh(supClus.eta());
0118     double py = supClus.energy() * sin(supClus.phi()) / cosh(supClus.eta());
0119     double pz = supClus.energy() * tanh(supClus.eta());
0120     RecoEcalCandidate theCand(0, math::XYZTLorentzVector(px, py, pz, supClus.energy()), math::XYZPoint(0, 0, 0));
0121
0122     // correct it?
0123     if (m_correct)
0124       theCand = correctEPosition(supClus, clusShape, nvtx);
0125
0126     double e9e25 = clusShape.eLong3x3() / clusShape.eLong5x5();
0127     double e1e9 = clusShape.eLong1x1() / clusShape.eLong3x3();
0128     double eSeL = hf_egamma::eSeLCorrected(clusShape.eShort3x3(), clusShape.eLong3x3(), 4);
0129     double var2d = (clusShape.eCOREe9() - (eSeL * m_intercept2DSlope));
0130
0131     bool isAcceptable = true;
0132     isAcceptable = isAcceptable && (e9e25 > m_e9e25Cut);
0133     isAcceptable = isAcceptable && (var2d > m_intercept2DCut);
0134     isAcceptable = isAcceptable && ((e1e9 < m_e1e9Cuthi) && (e1e9 > m_e1e9Cutlo));
0135     isAcceptable = isAcceptable && ((clusShape.eCOREe9() < m_eCOREe9Cuthi) && (clusShape.eCOREe9() > m_eCOREe9Cutlo));
0136     isAcceptable = isAcceptable && ((clusShape.eSeL() < m_eSeLCuthi) && (clusShape.eSeL() > m_eSeLCutlo));
0137
0138     if (isAcceptable) {
0139       theCand.setSuperCluster(theClusRef);
0140       RecoECand.push_back(theCand);
0141     }
0142   }
0143 }