CMS 3D CMS Logo

HFRecoEcalCandidateAlgo.cc
Go to the documentation of this file.
1 
9 //#includes
13 #include "CLHEP/Vector/LorentzVector.h"
14 #include <iostream>
17 #include "HFEGammaSLCorrector.h"
18 
19 using namespace std;
20 using namespace reco;
21 
23  double e9e25Cut,
24  double intercept2DCut,
25  double intercept2DSlope,
26  const std::vector<double>& e1e9Cut,
27  const std::vector<double>& eCOREe9Cut,
28  const std::vector<double>& eSeLCut,
29  const reco::HFValueStruct hfvv)
30  : m_correct(correct),
31  m_e9e25Cut(e9e25Cut),
32  m_intercept2DCut(intercept2DCut),
33  m_intercept2DSlope(intercept2DSlope),
34  m_e1e9Cuthi(e1e9Cut[1]),
35  m_eCOREe9Cuthi(eCOREe9Cut[1]),
36  m_eSeLCuthi(eSeLCut[1]),
37  m_e1e9Cutlo(e1e9Cut[0]),
38  m_eCOREe9Cutlo(eCOREe9Cut[0]),
39  m_eSeLCutlo(eSeLCut[0]),
40  m_era(4),
41  m_hfvv(hfvv) {}
42 
44  const HFEMClusterShape& shape,
45  int nvtx) const {
46  double corEta = original.eta();
47  //piece-wise log energy correction to eta
48  double logel = log(shape.eLong3x3() / 100.0);
49  double lowC = 0.00911;
50  double hiC = 0.0193;
51  double logSlope = 0.0146;
52  double logIntercept = -0.00988;
53  double sgn = (original.eta() > 0) ? (1.0) : (-1.0);
54  if (logel <= 1.25)
55  corEta += (lowC)*sgn;
56  else if ((logel > 1.25) && (logel <= 2))
57  corEta += (logIntercept + logSlope * logel) * sgn;
58  else if (logel > 2)
59  corEta += hiC * sgn;
60  //poly-3 cell eta correction to eta (de)
61  double p0 = 0.0125;
62  double p1 = -0.0475;
63  double p2 = 0.0732;
64  double p3 = -0.0506;
65  double de = sgn * (p0 + (p1 * (shape.CellEta())) + (p2 * (shape.CellEta() * shape.CellEta())) +
66  (p3 * (shape.CellEta() * shape.CellEta() * shape.CellEta())));
67  corEta += de;
68 
69  double corPhi = original.phi();
70  //poly-3 cell phi correction to phi (dp)
71  p0 = 0.0115;
72  p1 = -0.0394;
73  p2 = 0.0486;
74  p3 = -0.0335;
75  double dp = p0 + (p1 * (shape.CellPhi())) + (p2 * (shape.CellPhi() * shape.CellPhi())) +
76  (p3 * (shape.CellPhi() * shape.CellPhi() * shape.CellPhi()));
77  corPhi += dp;
78 
79  double corEnergy = original.energy();
80  //energy corrections
81  //flat correction for using only long fibers
82  double energyCorrect = 0.7397;
83  corEnergy = corEnergy / energyCorrect;
84  //corection based on ieta for pileup and general
85  //find ieta
86  int ieta = 0;
87  double etabounds[30] = {
88  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};
89  for (int kk = 0; kk <= 12; kk++) {
90  if ((fabs(corEta) < etabounds[kk + 1]) && (fabs(corEta) > etabounds[kk])) {
91  ieta = (corEta > 0) ? (kk + 29) : (-kk - 29);
92  }
93  }
94  //check if ieta is in bounds, should be [-41, -29] U [29, 41]
95  if (ieta != 0)
96  corEnergy = (m_hfvv.PUSlope(ieta) * 1.0 * (nvtx - 1) + m_hfvv.PUIntercept(ieta) * 1.0) * corEnergy * 1.0 *
97  m_hfvv.EnCor(ieta);
98 
99  //re-calculate the energy vector
100  double corPx = corEnergy * cos(corPhi) / cosh(corEta);
101  double corPy = corEnergy * sin(corPhi) / cosh(corEta);
102  double corPz = corEnergy * tanh(corEta);
103 
104  //make the candidate
105  RecoEcalCandidate corCand(0, math::XYZTLorentzVector(corPx, corPy, corPz, corEnergy), math::XYZPoint(0, 0, 0));
106 
107  return corCand;
108 }
109 
111  const HFEMClusterShapeAssociationCollection& AssocShapes,
112  RecoEcalCandidateCollection& RecoECand,
113  int nvtx) const {
114  //get super's and cluster shapes and associations
115  for (unsigned int i = 0; i < SuperClusters->size(); ++i) {
116  const SuperCluster& supClus = (*SuperClusters)[i];
117  reco::SuperClusterRef theClusRef = edm::Ref<SuperClusterCollection>(SuperClusters, i);
118  const HFEMClusterShapeRef clusShapeRef = AssocShapes.find(theClusRef)->val;
119  const HFEMClusterShape& clusShape = *clusShapeRef;
120 
121  // basic candidate
122  double px = supClus.energy() * cos(supClus.phi()) / cosh(supClus.eta());
123  double py = supClus.energy() * sin(supClus.phi()) / cosh(supClus.eta());
124  double pz = supClus.energy() * tanh(supClus.eta());
125  RecoEcalCandidate theCand(0, math::XYZTLorentzVector(px, py, pz, supClus.energy()), math::XYZPoint(0, 0, 0));
126 
127  // correct it?
128  if (m_correct)
129  theCand = correctEPosition(supClus, clusShape, nvtx);
130 
131  double e9e25 = clusShape.eLong3x3() / clusShape.eLong5x5();
132  double e1e9 = clusShape.eLong1x1() / clusShape.eLong3x3();
133  double eSeL = hf_egamma::eSeLCorrected(clusShape.eShort3x3(), clusShape.eLong3x3(), 4);
134  double var2d = (clusShape.eCOREe9() - (eSeL * m_intercept2DSlope));
135 
136  bool isAcceptable = true;
137  isAcceptable = isAcceptable && (e9e25 > m_e9e25Cut);
138  isAcceptable = isAcceptable && (var2d > m_intercept2DCut);
139  isAcceptable = isAcceptable && ((e1e9 < m_e1e9Cuthi) && (e1e9 > m_e1e9Cutlo));
140  isAcceptable = isAcceptable && ((clusShape.eCOREe9() < m_eCOREe9Cuthi) && (clusShape.eCOREe9() > m_eCOREe9Cutlo));
141  isAcceptable = isAcceptable && ((clusShape.eSeL() < m_eSeLCuthi) && (clusShape.eSeL() > m_eSeLCutlo));
142 
143  if (isAcceptable) {
144  theCand.setSuperCluster(theClusRef);
145  RecoECand.push_back(theCand);
146  }
147  }
148 }
const reco::HFValueStruct m_hfvv
float sgn(float val)
Definition: FWPFMaths.cc:9
double eLong5x5() const
double eLong1x1() const
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
const_iterator find(const key_type &k) const
find element with specified reference key
double CellPhi() const
double PUIntercept(int ieta) const
void produce(const edm::Handle< reco::SuperClusterCollection > &SuperClusters, const reco::HFEMClusterShapeAssociationCollection &AssocShapes, reco::RecoEcalCandidateCollection &RecoECand, int nvtx) const
double eta() const
pseudorapidity of cluster centroid
Definition: CaloCluster.h:180
double EnCor(int ieta) const
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
reco::RecoEcalCandidate correctEPosition(const reco::SuperCluster &original, const reco::HFEMClusterShape &shape, int nvtx) const
double eShort3x3() const
double eLong3x3() const
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
double energy() const
cluster energy
Definition: CaloCluster.h:148
double p2[4]
Definition: TauolaWrapper.h:90
double PUSlope(int ieta) const
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
double eSeLCorrected(double es, double el, double pc, double px, double py)
std::vector< RecoEcalCandidate > RecoEcalCandidateCollection
collectin of RecoEcalCandidate objects
fixed size matrix
double p1[4]
Definition: TauolaWrapper.h:89
double phi() const
azimuthal angle of cluster centroid
Definition: CaloCluster.h:183
double CellEta() const
HFRecoEcalCandidateAlgo(bool correct, double e9e25Cut, double intercept2DCut, double intercept2DSlope, const std::vector< double > &e1e9Cut, const std::vector< double > &eCOREe9Cut, const std::vector< double > &eSeLCut, const reco::HFValueStruct hfvv)
double p3[4]
Definition: TauolaWrapper.h:91