CMS 3D CMS Logo

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