CMS 3D CMS Logo

HGCalShowerShapeHelper.cc
Go to the documentation of this file.
2 
5 
7  : caloGeometryToken_{sumes.esConsumes()} {}
8 
9 void HGCalShowerShapeHelper::initPerEvent(const edm::EventSetup &iSetup, const std::vector<reco::PFRecHit> &pfRecHits) {
11  setPFRecHitPtrMap(pfRecHits);
12 }
13 
14 void HGCalShowerShapeHelper::initPerObject(const std::vector<std::pair<DetId, float> > &hitsAndFracs,
15  double minHitE,
16  double minHitET,
17  int minLayer,
18  int maxLayer,
19  DetId::Detector subDet) {
20  maxLayer = maxLayer <= 0 ? recHitTools_.lastLayerEE() : maxLayer;
21 
22  // Safety checks
23  nLayer_ = maxLayer - minLayer + 1;
24  assert(nLayer_ > 0);
25 
26  minHitE_ = minHitE;
27  minHitET_ = minHitET;
28  minHitET2_ = minHitET * minHitET;
30  maxLayer_ = maxLayer;
31  subDet_ = subDet;
32 
33  setFilteredHitsAndFractions(hitsAndFracs);
34 
36 }
37 
38 void HGCalShowerShapeHelper::setPFRecHitPtrMap(const std::vector<reco::PFRecHit> &recHits) {
39  pfRecHitPtrMap_.clear();
40 
41  for (const auto &recHit : recHits) {
42  pfRecHitPtrMap_[recHit.detId()] = &recHit;
43  }
44 }
45 
46 void HGCalShowerShapeHelper::setFilteredHitsAndFractions(const std::vector<std::pair<DetId, float> > &hitsAndFracs) {
47  hitsAndFracs_.clear();
48  hitEnergies_.clear();
49  hitEnergiesWithFracs_.clear();
50 
51  for (const auto &hnf : hitsAndFracs) {
52  DetId hitId = hnf.first;
53  float hitEfrac = hnf.second;
54 
55  int hitLayer = recHitTools_.getLayer(hitId);
56 
57  if (hitLayer > nLayer_) {
58  continue;
59  }
60 
61  if (hitId.det() != subDet_) {
62  continue;
63  }
64 
65  if (pfRecHitPtrMap_.find(hitId.rawId()) == pfRecHitPtrMap_.end()) {
66  continue;
67  }
68 
69  const reco::PFRecHit &recHit = *pfRecHitPtrMap_[hitId.rawId()];
70 
71  if (recHit.energy() < minHitE_) {
72  continue;
73  }
74 
75  if (recHit.pt2() < minHitET2_) {
76  continue;
77  }
78 
79  // Fill the vectors
80  hitsAndFracs_.push_back(hnf);
81  hitEnergies_.push_back(recHit.energy());
82  hitEnergiesWithFracs_.push_back(recHit.energy() * hitEfrac);
83  }
84 }
85 
87  layerEnergies_.clear();
88  layerEnergies_.resize(nLayer_);
89 
90  layerCentroids_.clear();
91  layerCentroids_.resize(nLayer_);
92 
93  centroid_.SetXYZ(0, 0, 0);
94 
95  int iHit = -1;
96  double totalW = 0.0;
97 
98  // Compute the centroid per layer
99  for (const auto &hnf : hitsAndFracs_) {
100  iHit++;
101 
102  DetId hitId = hnf.first;
103 
104  double weight = hitEnergies_[iHit];
105  totalW += weight;
106 
107  const auto &hitPos = recHitTools_.getPosition(hitId);
108  ROOT::Math::XYZVector hitXYZ(hitPos.x(), hitPos.y(), hitPos.z());
109 
110  centroid_ += weight * hitXYZ;
111 
112  int hitLayer = recHitTools_.getLayer(hitId) - 1;
113 
114  layerEnergies_[hitLayer] += weight;
115  layerCentroids_[hitLayer] += weight * hitXYZ;
116  }
117 
118  int iLayer = -1;
119 
120  for (auto &centroid : layerCentroids_) {
121  iLayer++;
122 
123  if (layerEnergies_[iLayer]) {
124  centroid /= layerEnergies_[iLayer];
125  }
126  }
127 
128  if (totalW) {
129  centroid_ /= totalW;
130  }
131 }
132 
135 }
136 
137 const double HGCalShowerShapeHelper::getRvar(double cylinderR, double energyNorm, bool useFractions, bool useCellSize) {
138  if (hitsAndFracs_.empty()) {
139  return 0.0;
140  }
141 
142  if (energyNorm <= 0.0) {
143  edm::LogWarning("HGCalShowerShapeHelper")
144  << "Encountered negative or zero energy for HGCal R-variable denominator: " << energyNorm << std::endl;
145  }
146 
147  double cylinderR2 = cylinderR * cylinderR;
148 
149  double rVar = 0.0;
150 
151  auto hitEnergyIter = useFractions ? hitEnergiesWithFracs_.begin() : hitEnergies_.begin();
152 
153  hitEnergyIter--;
154 
155  for (const auto &hnf : hitsAndFracs_) {
156  hitEnergyIter++;
157 
158  DetId hitId = hnf.first;
159 
160  int hitLayer = recHitTools_.getLayer(hitId) - 1;
161 
162  const auto &hitPos = recHitTools_.getPosition(hitId);
163  ROOT::Math::XYZVector hitXYZ(hitPos.x(), hitPos.y(), hitPos.z());
164 
165  auto distXYZ = hitXYZ - layerCentroids_[hitLayer];
166 
167  double r2 = distXYZ.x() * distXYZ.x() + distXYZ.y() * distXYZ.y();
168 
169  // Including the cell size seems to make the variable less sensitive to the HD/LD transition region
170  if (useCellSize) {
171  if (std::sqrt(r2) > cylinderR + getCellSize(hitId)) {
172  continue;
173  }
174  }
175 
176  else if (r2 > cylinderR2) {
177  continue;
178  }
179 
180  rVar += *hitEnergyIter;
181  }
182 
183  rVar /= energyNorm;
184 
185  return rVar;
186 }
187 
189  if (hitsAndFracs_.empty()) {
190  return ShowerWidths();
191  }
192 
193  double cylinderR2 = cylinderR * cylinderR;
194 
195  TMatrixD covMat(3, 3);
196 
197  double dxdx = 0.0;
198  double dydy = 0.0;
199  double dzdz = 0.0;
200 
201  double dxdy = 0.0;
202  double dydz = 0.0;
203  double dzdx = 0.0;
204 
205  double totalW = 0.0;
206 
207  auto hitEnergyIter = useFractions ? hitEnergiesWithFracs_.begin() : hitEnergies_.begin();
208 
209  int nHit = 0;
210  hitEnergyIter--;
211 
212  for (const auto &hnf : hitsAndFracs_) {
213  hitEnergyIter++;
214 
215  DetId hitId = hnf.first;
216 
217  const auto &hitPos = recHitTools_.getPosition(hitId);
218  ROOT::Math::XYZVector hitXYZ(hitPos.x(), hitPos.y(), hitPos.z());
219 
220  int hitLayer = recHitTools_.getLayer(hitId) - 1;
221 
222  ROOT::Math::XYZVector radXYZ = hitXYZ - layerCentroids_[hitLayer];
223 
224  double r2 = radXYZ.x() * radXYZ.x() + radXYZ.y() * radXYZ.y();
225 
226  if (r2 > cylinderR2) {
227  continue;
228  }
229 
230  ROOT::Math::XYZVector dXYZ = hitXYZ - centroid_;
231 
232  double weight = *hitEnergyIter;
233  totalW += weight;
234 
235  dxdx += weight * dXYZ.x() * dXYZ.x();
236  dydy += weight * dXYZ.y() * dXYZ.y();
237  dzdz += weight * dXYZ.z() * dXYZ.z();
238 
239  dxdy += weight * dXYZ.x() * dXYZ.y();
240  dydz += weight * dXYZ.y() * dXYZ.z();
241  dzdx += weight * dXYZ.z() * dXYZ.x();
242 
243  nHit++;
244  }
245 
246  if (!totalW || nHit < 2) {
247  return ShowerWidths();
248  }
249 
250  dxdx /= totalW;
251  dydy /= totalW;
252  dzdz /= totalW;
253 
254  dxdy /= totalW;
255  dydz /= totalW;
256  dzdx /= totalW;
257 
258  covMat(0, 0) = dxdx;
259  covMat(1, 1) = dydy;
260  covMat(2, 2) = dzdz;
261 
262  covMat(0, 1) = covMat(1, 0) = dxdy;
263  covMat(0, 2) = covMat(2, 0) = dzdx;
264  covMat(1, 2) = covMat(2, 1) = dydz;
265 
266  if (!covMat.Sum()) {
267  return ShowerWidths();
268  }
269 
270  // Get eigen values and vectors
271  TVectorD eigVals(3);
272  TMatrixD eigVecMat(3, 3);
273 
274  eigVecMat = covMat.EigenVectors(eigVals);
275 
276  ShowerWidths returnWidths;
277 
278  returnWidths.sigma2xx = dxdx;
279  returnWidths.sigma2yy = dydy;
280  returnWidths.sigma2zz = dzdz;
281 
282  returnWidths.sigma2xy = dxdy;
283  returnWidths.sigma2yz = dydz;
284  returnWidths.sigma2zx = dzdx;
285 
286  returnWidths.sigma2uu = eigVals(1);
287  returnWidths.sigma2vv = eigVals(2);
288  returnWidths.sigma2ww = eigVals(0);
289 
290  return returnWidths;
291 }
HGCalShowerShapeHelper::setFilteredHitsAndFractions
void setFilteredHitsAndFractions(const std::vector< std::pair< DetId, float > > &hitsAndFracs)
Definition: HGCalShowerShapeHelper.cc:46
HGCalShowerShapeHelper::nLayer_
int nLayer_
Definition: HGCalShowerShapeHelper.h:115
hgcal::RecHitTools::getSiThickIndex
int getSiThickIndex(const DetId &) const
Definition: RecHitTools.cc:202
HGCalShowerShapeHelper::ShowerWidths::sigma2yy
double sigma2yy
Definition: HGCalShowerShapeHelper.h:64
mps_merge.weight
weight
Definition: mps_merge.py:88
DetId::det
constexpr Detector det() const
get the detector field from this detid
Definition: DetId.h:46
HGCalShowerShapeHelper::ShowerWidths::sigma2uu
double sigma2uu
Definition: HGCalShowerShapeHelper.h:71
cms::cuda::assert
assert(be >=bs)
HGCalShowerShapeHelper::caloGeometryToken_
edm::ESGetToken< CaloGeometry, CaloGeometryRecord > caloGeometryToken_
Definition: HGCalShowerShapeHelper.h:118
HGCalShowerShapeHelper::initPerEvent
void initPerEvent(const edm::EventSetup &iSetup, const std::vector< reco::PFRecHit > &recHits)
Definition: HGCalShowerShapeHelper.cc:9
HGCalShowerShapeHelper::minHitET_
double minHitET_
Definition: HGCalShowerShapeHelper.h:111
HGCalShowerShapeHelper::centroid_
ROOT::Math::XYZVector centroid_
Definition: HGCalShowerShapeHelper.h:126
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:122
rpcPointValidation_cfi.recHit
recHit
Definition: rpcPointValidation_cfi.py:7
HGCalShowerShapeHelper::setPFRecHitPtrMap
void setPFRecHitPtrMap(const std::vector< reco::PFRecHit > &recHits)
Definition: HGCalShowerShapeHelper.cc:38
DetId
Definition: DetId.h:17
HGCalShowerShapeHelper::initPerObject
void initPerObject(const std::vector< std::pair< DetId, float > > &hitsAndFracs, double minHitE=0, double minHitET=0, int minLayer=1, int maxLayer=-1, DetId::Detector subDet=DetId::HGCalEE)
Definition: HGCalShowerShapeHelper.cc:14
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
HGCalShowerShapeHelper::ShowerWidths::sigma2vv
double sigma2vv
Definition: HGCalShowerShapeHelper.h:72
HGCalShowerShapeHelper::ShowerWidths::sigma2zz
double sigma2zz
Definition: HGCalShowerShapeHelper.h:65
HGCalShowerShapeHelper.h
HGCalShowerShapeHelper::ShowerWidths
Definition: HGCalShowerShapeHelper.h:62
HGCalShowerShapeHelper::hitEnergies_
std::vector< double > hitEnergies_
Definition: HGCalShowerShapeHelper.h:123
HGCalShowerShapeHelper::ShowerWidths::sigma2yz
double sigma2yz
Definition: HGCalShowerShapeHelper.h:68
HGCalShowerShapeHelper::hitsAndFracs_
std::vector< std::pair< DetId, float > > hitsAndFracs_
Definition: HGCalShowerShapeHelper.h:122
FastTrackerRecHitMaskProducer_cfi.recHits
recHits
Definition: FastTrackerRecHitMaskProducer_cfi.py:8
HGCalShowerShapeHelper::setLayerWiseInfo
void setLayerWiseInfo()
Definition: HGCalShowerShapeHelper.cc:86
HLT_FULL_cff.minLayer
minLayer
Definition: HLT_FULL_cff.py:53042
HGCalShowerShapeHelper::pfRecHitPtrMap_
std::unordered_map< uint32_t, const reco::PFRecHit * > pfRecHitPtrMap_
Definition: HGCalShowerShapeHelper.h:121
HGCalShowerShapeHelper::subDet_
DetId::Detector subDet_
Definition: HGCalShowerShapeHelper.h:116
HGCalShowerShapeHelper::recHitTools_
hgcal::RecHitTools recHitTools_
Definition: HGCalShowerShapeHelper.h:119
HGCalShowerShapeHelper::HGCalShowerShapeHelper
HGCalShowerShapeHelper(edm::ConsumesCollector &&sumes)
Definition: HGCalShowerShapeHelper.cc:6
hgcal::RecHitTools::lastLayerEE
unsigned int lastLayerEE(bool nose=false) const
Definition: RecHitTools.h:64
diffTwoXMLs.r2
r2
Definition: diffTwoXMLs.py:73
BeamSpotPI::dydz
Definition: BeamSpotPayloadInspectorHelper.h:37
HGCalShowerShapeHelper::ShowerWidths::sigma2xy
double sigma2xy
Definition: HGCalShowerShapeHelper.h:67
trackerHitRTTI::vector
Definition: trackerHitRTTI.h:21
DetId::Detector
Detector
Definition: DetId.h:24
edm::EventSetup
Definition: EventSetup.h:57
HGCalShowerShapeHelper::layerCentroids_
std::vector< ROOT::Math::XYZVector > layerCentroids_
Definition: HGCalShowerShapeHelper.h:128
HGCalShowerShapeHelper::minHitET2_
double minHitET2_
Definition: HGCalShowerShapeHelper.h:112
HGCalShowerShapeHelper::getRvar
const double getRvar(double cylinderR, double energyNorm, bool useFractions=true, bool useCellSize=true)
Definition: HGCalShowerShapeHelper.cc:137
HGCalShowerShapeHelper::ShowerWidths::sigma2ww
double sigma2ww
Definition: HGCalShowerShapeHelper.h:73
hgcal::RecHitTools::getLayer
unsigned int getLayer(DetId::Detector type, bool nose=false) const
Definition: RecHitTools.cc:294
edm::EventSetup::getData
bool getData(T &iHolder) const
Definition: EventSetup.h:120
DetId::rawId
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
HGCalShowerShapeHelper::layerEnergies_
std::vector< double > layerEnergies_
Definition: HGCalShowerShapeHelper.h:127
HGCalShowerShapeHelper::getPCAWidths
const ShowerWidths getPCAWidths(double cylinderR, bool useFractions=false)
Definition: HGCalShowerShapeHelper.cc:188
hgcal::RecHitTools::setGeometry
void setGeometry(CaloGeometry const &)
Definition: RecHitTools.cc:68
HGCalShowerShapeHelper::kHDWaferCellSize_
static const double kHDWaferCellSize_
Definition: HGCalShowerShapeHelper.h:58
HGCalShowerShapeHelper::hitEnergiesWithFracs_
std::vector< double > hitEnergiesWithFracs_
Definition: HGCalShowerShapeHelper.h:124
ROOT::Math::XYZVector
Transform3DPJ::Vector XYZVector
Definition: Transform3DPJ.cc:34
HGCalShowerShapeHelper::minHitE_
double minHitE_
Definition: HGCalShowerShapeHelper.h:110
HGCalShowerShapeHelper::getCellSize
const double getCellSize(DetId detId)
Definition: HGCalShowerShapeHelper.cc:133
reco::PFRecHit
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
Definition: PFRecHit.h:31
HGCalShowerShapeHelper::maxLayer_
int maxLayer_
Definition: HGCalShowerShapeHelper.h:114
HGCalShowerShapeHelper::ShowerWidths::sigma2zx
double sigma2zx
Definition: HGCalShowerShapeHelper.h:69
HGCalShowerShapeHelper::ShowerWidths::sigma2xx
double sigma2xx
Definition: HGCalShowerShapeHelper.h:63
hgcal::RecHitTools::getPosition
GlobalPoint getPosition(const DetId &id) const
Definition: RecHitTools.cc:126
HGCalShowerShapeHelper::kLDWaferCellSize_
static const double kLDWaferCellSize_
Definition: HGCalShowerShapeHelper.h:57
edm::ConsumesCollector
Definition: ConsumesCollector.h:45
weight
Definition: weight.py:1
HGCalShowerShapeHelper::minLayer_
int minLayer_
Definition: HGCalShowerShapeHelper.h:113