CMS 3D CMS Logo
ChargeIsolation.cc
Go to the documentation of this file.
1 #include "Calibration/IsolatedParticles/interface/CaloConstants.h"
2 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
3 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
4 #include "Calibration/IsolatedParticles/interface/FindDistCone.h"
5 #include "Calibration/IsolatedParticles/interface/MatrixECALDetIds.h"
6 #include "Calibration/IsolatedParticles/interface/MatrixHCALDetIds.h"
7 #include "FWCore/MessageLogger/interface/MessageLogger.h"
8 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
9 
10 namespace spr {
11 
12  double chargeIsolationEcal(unsigned int trkIndex,
13  std::vector<spr::propagatedTrackID>& vdetIds,
14  const CaloGeometry* geo,
15  const CaloTopology* caloTopology,
16  int ieta,
17  int iphi,
18  bool debug) {
19  const DetId coreDet = vdetIds[trkIndex].detIdECAL;
20  if (debug) {
21  if (coreDet.subdetId() == EcalBarrel)
22  edm::LogVerbatim("IsoTrack") << "DetId " << (EBDetId)(coreDet) << " Flag " << vdetIds[trkIndex].okECAL;
23  else
24  edm::LogVerbatim("IsoTrack") << "DetId " << (EEDetId)(coreDet) << " Flag " << vdetIds[trkIndex].okECAL;
25  }
26  double maxNearP = -1.0;
27  if (vdetIds[trkIndex].okECAL) {
28  std::vector<DetId> vdets = spr::matrixECALIds(coreDet, ieta, iphi, geo, caloTopology, debug);
29  if (debug)
30  edm::LogVerbatim("IsoTrack") << "chargeIsolationEcal:: eta/phi/dets " << ieta << " " << iphi << " "
31  << vdets.size();
32 
33  for (unsigned int indx = 0; indx < vdetIds.size(); ++indx) {
34  if (indx != trkIndex && vdetIds[indx].ok && vdetIds[indx].okECAL) {
35  const DetId anyCell = vdetIds[indx].detIdECAL;
36  if (!spr::chargeIsolation(anyCell, vdets)) {
37  const reco::Track* pTrack = &(*(vdetIds[indx].trkItr));
38  if (maxNearP < pTrack->p())
39  maxNearP = pTrack->p();
40  }
41  }
42  }
43  }
44  return maxNearP;
45  }
46 
47  double chargeIsolationGenEcal(unsigned int trkIndex,
48  std::vector<spr::propagatedGenParticleID>& vdetIds,
49  const CaloGeometry* geo,
50  const CaloTopology* caloTopology,
51  int ieta,
52  int iphi,
53  bool debug) {
54  const DetId coreDet = vdetIds[trkIndex].detIdECAL;
55  if (debug) {
56  if (coreDet.subdetId() == EcalBarrel)
57  edm::LogVerbatim("IsoTrack") << "DetId " << (EBDetId)(coreDet) << " Flag " << vdetIds[trkIndex].okECAL;
58  else
59  edm::LogVerbatim("IsoTrack") << "DetId " << (EEDetId)(coreDet) << " Flag " << vdetIds[trkIndex].okECAL;
60  }
61  double maxNearP = -1.0;
62  if (vdetIds[trkIndex].okECAL) {
63  std::vector<DetId> vdets = spr::matrixECALIds(coreDet, ieta, iphi, geo, caloTopology, debug);
64  if (debug)
65  edm::LogVerbatim("IsoTrack") << "chargeIsolationGenEcal:: eta/phi/dets " << ieta << " " << iphi << " "
66  << vdets.size();
67 
68  for (unsigned int indx = 0; indx < vdetIds.size(); ++indx) {
69  if (indx != trkIndex && vdetIds[indx].ok && vdetIds[indx].okECAL) {
70  const DetId anyCell = vdetIds[indx].detIdECAL;
71  if (!spr::chargeIsolation(anyCell, vdets)) {
72  const reco::GenParticle* pTrack = &(*(vdetIds[indx].trkItr));
73  if (maxNearP < pTrack->p())
74  maxNearP = pTrack->p();
75  }
76  }
77  }
78  }
79  return maxNearP;
80  }
81 
82  double chargeIsolationEcal(const DetId& coreDet,
83  reco::TrackCollection::const_iterator trkItr,
84  edm::Handle<reco::TrackCollection> trkCollection,
85  const CaloGeometry* geo,
86  const CaloTopology* caloTopology,
87  const MagneticField* bField,
88  int ieta,
89  int iphi,
90  const std::string& theTrackQuality,
91  bool debug) {
92  const CaloSubdetectorGeometry* barrelGeom = (geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel));
93  const CaloSubdetectorGeometry* endcapGeom = (geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap));
94 
95  std::vector<DetId> vdets = spr::matrixECALIds(coreDet, ieta, iphi, geo, caloTopology, debug);
96  if (debug)
97  edm::LogVerbatim("IsoTrack") << "chargeIsolation:: eta/phi/dets " << ieta << " " << iphi << " " << vdets.size();
98 
99  double maxNearP = -1.0;
100  reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
101 
102  // const DetId anyCell,
103  reco::TrackCollection::const_iterator trkItr2;
104  for (trkItr2 = trkCollection->begin(); trkItr2 != trkCollection->end(); ++trkItr2) {
105  const reco::Track* pTrack2 = &(*trkItr2);
106 
107  bool trkQuality = (trackQuality_ != reco::TrackBase::undefQuality) ? (pTrack2->quality(trackQuality_)) : true;
108  if ((trkItr2 != trkItr) && trkQuality) {
109  std::pair<math::XYZPoint, bool> info = spr::propagateECAL(pTrack2, bField);
110  const GlobalPoint point2(info.first.x(), info.first.y(), info.first.z());
111 
112  if (info.second) {
113  if (std::abs(point2.eta()) < spr::etaBEEcal) {
114  const DetId anyCell = barrelGeom->getClosestCell(point2);
115  if (!spr::chargeIsolation(anyCell, vdets)) {
116  if (debug)
117  edm::LogVerbatim("IsoTrack")
118  << "chargeIsolationEcal Cell " << (EBDetId)(anyCell) << " pt " << pTrack2->p();
119  if (maxNearP < pTrack2->p())
120  maxNearP = pTrack2->p();
121  }
122  } else {
123  if (endcapGeom) {
124  const DetId anyCell = endcapGeom->getClosestCell(point2);
125  if (!spr::chargeIsolation(anyCell, vdets)) {
126  if (debug)
127  edm::LogVerbatim("IsoTrack")
128  << "chargeIsolationEcal Cell " << (EEDetId)(anyCell) << " pt " << pTrack2->p();
129  if (maxNearP < pTrack2->p())
130  maxNearP = pTrack2->p();
131  }
132  }
133  }
134  } //info.second
135  }
136  }
137  return maxNearP;
138  }
139 
140  double chargeIsolationHcal(unsigned int trkIndex,
141  std::vector<spr::propagatedTrackID>& vdetIds,
142  const HcalTopology* topology,
143  int ieta,
144  int iphi,
145  bool debug) {
146  std::vector<DetId> dets(1, vdetIds[trkIndex].detIdHCAL);
147  if (debug) {
148  edm::LogVerbatim("IsoTrack") << "DetId " << (HcalDetId)(dets[0]) << " Flag " << vdetIds[trkIndex].okHCAL;
149  }
150 
151  double maxNearP = -1.0;
152  if (vdetIds[trkIndex].okHCAL) {
153  std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, false, debug);
154  if (debug)
155  edm::LogVerbatim("IsoTrack") << "chargeIsolationHcal:: eta/phi/dets " << ieta << " " << iphi << " "
156  << vdets.size();
157 
158  for (unsigned indx = 0; indx < vdetIds.size(); ++indx) {
159  if (indx != trkIndex && vdetIds[indx].ok && vdetIds[indx].okHCAL) {
160  const DetId anyCell = vdetIds[indx].detIdHCAL;
161  if (!spr::chargeIsolation(anyCell, vdets)) {
162  const reco::Track* pTrack = &(*(vdetIds[indx].trkItr));
163  if (debug)
164  edm::LogVerbatim("IsoTrack")
165  << "chargeIsolationHcal Cell " << (HcalDetId)(anyCell) << " pt " << pTrack->p();
166  if (maxNearP < pTrack->p())
167  maxNearP = pTrack->p();
168  }
169  }
170  }
171  }
172  return maxNearP;
173  }
174 
175  //===========================================================================================================
176 
177  double chargeIsolationHcal(reco::TrackCollection::const_iterator trkItr,
178  edm::Handle<reco::TrackCollection> trkCollection,
179  const DetId ClosestCell,
180  const HcalTopology* topology,
181  const CaloSubdetectorGeometry* gHB,
182  const MagneticField* bField,
183  int ieta,
184  int iphi,
185  const std::string& theTrackQuality,
186  bool debug) {
187  std::vector<DetId> dets(1, ClosestCell);
188  if (debug)
189  edm::LogVerbatim("IsoTrack") << (HcalDetId)ClosestCell;
190  std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, false, debug);
191 
192  if (debug) {
193  for (unsigned int i = 0; i < vdets.size(); i++) {
194  edm::LogVerbatim("IsoTrack") << "HcalDetId in " << 2 * ieta + 1 << "x" << 2 * iphi + 1 << " "
195  << static_cast<HcalDetId>(vdets[i]);
196  }
197  }
198 
199  double maxNearP = -1.0;
200  reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
201 
202  reco::TrackCollection::const_iterator trkItr2;
203  for (trkItr2 = trkCollection->begin(); trkItr2 != trkCollection->end(); ++trkItr2) {
204  const reco::Track* pTrack2 = &(*trkItr2);
205 
206  bool trkQuality = (trackQuality_ != reco::TrackBase::undefQuality) ? (pTrack2->quality(trackQuality_)) : true;
207  if ((trkItr2 != trkItr) && trkQuality) {
208  std::pair<math::XYZPoint, bool> info = spr::propagateHCAL(pTrack2, bField);
209  const GlobalPoint point2(info.first.x(), info.first.y(), info.first.z());
210 
211  if (debug)
212  edm::LogVerbatim("IsoTrack") << "Track2 (p,eta,phi) " << pTrack2->p() << " " << pTrack2->eta() << " "
213  << pTrack2->phi();
214  if (info.second) {
215  const DetId anyCell = gHB->getClosestCell(point2);
216  if (!spr::chargeIsolation(anyCell, vdets)) {
217  if (maxNearP < pTrack2->p())
218  maxNearP = pTrack2->p();
219  }
220 
221  if (debug)
222  edm::LogVerbatim("IsoTrack") << "maxNearP " << maxNearP << " thisCell " << static_cast<HcalDetId>(anyCell)
223  << " (" << info.first.x() << "," << info.first.y() << "," << info.first.z()
224  << ")";
225  }
226  }
227  }
228  return maxNearP;
229  }
230 
231  bool chargeIsolation(const DetId anyCell, std::vector<DetId>& vdets) {
232  bool isIsolated = true;
233  for (unsigned int i = 0; i < vdets.size(); i++) {
234  if (anyCell == vdets[i]) {
235  isIsolated = false;
236  break;
237  }
238  }
239  return isIsolated;
240  }
241 
242  double coneChargeIsolation(const edm::Event& iEvent,
243  const edm::EventSetup& iSetup,
244  reco::TrackCollection::const_iterator trkItr,
245  edm::Handle<reco::TrackCollection> trkCollection,
246  TrackDetectorAssociator& associator,
247  TrackAssociatorParameters& parameters_,
248  const std::string& theTrackQuality,
249  int& nNearTRKs,
250  int& nLayers_maxNearP,
251  int& trkQual_maxNearP,
252  double& maxNearP_goodTrk,
253  const GlobalPoint& hpoint1,
254  const GlobalVector& trackMom,
255  double dR) {
256  nNearTRKs = 0;
257  nLayers_maxNearP = 0;
258  trkQual_maxNearP = -1;
259  maxNearP_goodTrk = -999.0;
260  double maxNearP = -999.0;
261  reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
262 
263  // Iterate over tracks
264  reco::TrackCollection::const_iterator trkItr2;
265  for (trkItr2 = trkCollection->begin(); trkItr2 != trkCollection->end(); ++trkItr2) {
266  // Get track
267  const reco::Track* pTrack2 = &(*trkItr2);
268 
269  // Get track qual, nlayers, and hit pattern
270  bool trkQuality = (trackQuality_ != reco::TrackBase::undefQuality) ? (pTrack2->quality(trackQuality_)) : true;
271  if (trkQuality)
272  trkQual_maxNearP = 1;
273  const reco::HitPattern& hitp = pTrack2->hitPattern();
274  nLayers_maxNearP = hitp.trackerLayersWithMeasurement();
275 
276  // Skip if the neighboring track candidate is the iso-track
277  // candidate
278  if (trkItr2 != trkItr) {
279  // Get propagator
280  const FreeTrajectoryState fts2 =
281  associator.getFreeTrajectoryState(&iSetup.getData(parameters_.bFieldToken), *pTrack2);
282  TrackDetMatchInfo info2 = associator.associate(iEvent, iSetup, fts2, parameters_);
283 
284  // Make sure it reaches Hcal
285  if (info2.isGoodHcal) {
286  const GlobalPoint point2(info2.trkGlobPosAtHcal.x(), info2.trkGlobPosAtHcal.y(), info2.trkGlobPosAtHcal.z());
287 
288  int isConeChargedIso = spr::coneChargeIsolation(hpoint1, point2, trackMom, dR);
289 
290  if (isConeChargedIso == 0) {
291  nNearTRKs++;
292  if (maxNearP < pTrack2->p()) {
293  maxNearP = pTrack2->p();
294  if (trkQual_maxNearP > 0 && nLayers_maxNearP > 7 && maxNearP_goodTrk < pTrack2->p()) {
295  maxNearP_goodTrk = pTrack2->p();
296  }
297  }
298  }
299  }
300  }
301  } // Iterate over track loop
302 
303  return maxNearP;
304  }
305 
306  double chargeIsolationCone(unsigned int trkIndex,
307  std::vector<spr::propagatedTrackDirection>& trkDirs,
308  double dR,
309  int& nNearTRKs,
310  bool debug) {
311  double maxNearP = -1.0;
312  nNearTRKs = 0;
313  if (trkDirs[trkIndex].okHCAL) {
314  if (debug)
315  edm::LogVerbatim("IsoTrack") << "chargeIsolationCone with " << trkDirs.size() << " tracks ";
316  for (unsigned int indx = 0; indx < trkDirs.size(); ++indx) {
317  if (indx != trkIndex && trkDirs[indx].ok && trkDirs[indx].okHCAL) {
318  int isConeChargedIso = spr::coneChargeIsolation(
319  trkDirs[trkIndex].pointHCAL, trkDirs[indx].pointHCAL, trkDirs[trkIndex].directionHCAL, dR);
320  if (isConeChargedIso == 0) {
321  nNearTRKs++;
322  const reco::Track* pTrack = &(*(trkDirs[indx].trkItr));
323  if (maxNearP < pTrack->p())
324  maxNearP = pTrack->p();
325  }
326  }
327  }
328  }
329 
330  if (debug)
331  edm::LogVerbatim("IsoTrack") << "chargeIsolationCone Track " << trkDirs[trkIndex].okHCAL << " maxNearP "
332  << maxNearP;
333  return maxNearP;
334  }
335 
336  double chargeIsolationGenCone(unsigned int trkIndex,
337  std::vector<spr::propagatedGenParticleID>& trkDirs,
338  double dR,
339  int& nNearTRKs,
340  bool debug) {
341  double maxNearP = -1.0;
342  nNearTRKs = 0;
343  if (trkDirs[trkIndex].okHCAL) {
344  if (debug)
345  edm::LogVerbatim("IsoTrack") << "chargeIsolationCone with " << trkDirs.size() << " tracks ";
346  for (unsigned int indx = 0; indx < trkDirs.size(); ++indx) {
347  if (indx != trkIndex && trkDirs[indx].ok && trkDirs[indx].okHCAL) {
348  int isConeChargedIso = spr::coneChargeIsolation(
349  trkDirs[trkIndex].pointHCAL, trkDirs[indx].pointHCAL, trkDirs[trkIndex].directionHCAL, dR);
350  if (isConeChargedIso == 0) {
351  nNearTRKs++;
352  const reco::GenParticle* pTrack = &(*(trkDirs[indx].trkItr));
353  if (maxNearP < pTrack->p())
354  maxNearP = pTrack->p();
355  }
356  }
357  }
358  }
359 
360  if (debug)
361  edm::LogVerbatim("IsoTrack") << "chargeIsolationGenCone Track " << trkDirs[trkIndex].okHCAL << " maxNearP "
362  << maxNearP;
363  return maxNearP;
364  }
365 
366  std::pair<double, double> chargeIsolationCone(unsigned int trkIndex,
367  std::vector<spr::propagatedTrackDirection>& trkDirs,
368  double dR,
369  bool debug) {
370  double maxNearP = -1.0;
371  double sumP = 0;
372  if (trkDirs[trkIndex].okHCAL) {
373  if (debug)
374  edm::LogVerbatim("IsoTrack") << "chargeIsolationCone with " << trkDirs.size() << " tracks ";
375  for (unsigned int indx = 0; indx < trkDirs.size(); ++indx) {
376  if (indx != trkIndex && trkDirs[indx].ok && trkDirs[indx].okHCAL) {
377  int isConeChargedIso = spr::coneChargeIsolation(
378  trkDirs[trkIndex].pointHCAL, trkDirs[indx].pointHCAL, trkDirs[trkIndex].directionHCAL, dR);
379  if (isConeChargedIso == 0) {
380  const reco::Track* pTrack = &(*(trkDirs[indx].trkItr));
381  sumP += (pTrack->p());
382  if (maxNearP < pTrack->p())
383  maxNearP = pTrack->p();
384  }
385  }
386  }
387  }
388 
389  if (debug)
390  edm::LogVerbatim("IsoTrack") << "chargeIsolationCone Track " << trkDirs[trkIndex].okHCAL << " maxNearP "
391  << maxNearP << ":" << sumP;
392  return std::pair<double, double>(maxNearP, sumP);
393  }
394 
395  int coneChargeIsolation(const GlobalPoint& hpoint1,
396  const GlobalPoint& point2,
397  const GlobalVector& trackMom,
398  double dR) {
399  int isIsolated = 1;
400  if (spr::getDistInPlaneTrackDir(hpoint1, trackMom, point2) > dR)
401  isIsolated = 1;
402  else
403  isIsolated = 0;
404  return isIsolated;
405  }
406 
407  double chargeIsolation(const edm::Event& iEvent,
408  const edm::EventSetup& iSetup,
409  CaloNavigator<DetId>& theNavigator,
410  reco::TrackCollection::const_iterator trkItr,
411  edm::Handle<reco::TrackCollection> trkCollection,
412  const CaloSubdetectorGeometry* gEB,
413  const CaloSubdetectorGeometry* gEE,
414  TrackDetectorAssociator& associator,
415  TrackAssociatorParameters& parameters_,
416  int ieta,
417  int iphi,
418  const std::string& theTrackQuality,
419  bool debug) {
420  double maxNearP = -1.0;
421  reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
422 
423  // const DetId anyCell,
424  reco::TrackCollection::const_iterator trkItr2;
425  for (trkItr2 = trkCollection->begin(); trkItr2 != trkCollection->end(); ++trkItr2) {
426  const reco::Track* pTrack2 = &(*trkItr2);
427 
428  bool trkQuality = pTrack2->quality(trackQuality_);
429  if ((trkItr2 != trkItr) && trkQuality) {
430  const FreeTrajectoryState fts2 =
431  associator.getFreeTrajectoryState(&iSetup.getData(parameters_.bFieldToken), *pTrack2);
432  TrackDetMatchInfo info2 = associator.associate(iEvent, iSetup, fts2, parameters_);
433  const GlobalPoint point2(info2.trkGlobPosAtEcal.x(), info2.trkGlobPosAtEcal.y(), info2.trkGlobPosAtEcal.z());
434 
435  if (info2.isGoodEcal) {
436  if (std::abs(point2.eta()) < spr::etaBEEcal) {
437  const DetId anyCell = gEB->getClosestCell(point2);
438  if (debug)
439  edm::LogVerbatim("IsoTrack")
440  << "chargeIsolation:: EB cell " << (EBDetId)(anyCell) << " for pt " << pTrack2->p();
441  if (!spr::chargeIsolation(anyCell, theNavigator, ieta, iphi)) {
442  if (maxNearP < pTrack2->p())
443  maxNearP = pTrack2->p();
444  }
445  } else {
446  const DetId anyCell = gEE->getClosestCell(point2);
447  if (debug)
448  edm::LogVerbatim("IsoTrack")
449  << "chargeIsolation:: EE cell " << (EEDetId)(anyCell) << " for pt " << pTrack2->p();
450  if (!spr::chargeIsolation(anyCell, theNavigator, ieta, iphi)) {
451  if (maxNearP < pTrack2->p())
452  maxNearP = pTrack2->p();
453  }
454  }
455  } //info2.isGoodEcal
456  }
457  }
458  return maxNearP;
459  }
460 
461  //===========================================================================================================
462 
463  bool chargeIsolation(const DetId anyCell, CaloNavigator<DetId>& navigator, int ieta, int iphi) {
464  bool isIsolated = false;
465 
466  DetId thisDet;
467 
468  for (int dx = -ieta; dx < ieta + 1; ++dx) {
469  for (int dy = -iphi; dy < iphi + 1; ++dy) {
470  thisDet = navigator.offsetBy(dx, dy);
471  navigator.home();
472 
473  if (thisDet != DetId(0)) {
474  if (thisDet == anyCell) {
475  isIsolated = false;
476  return isIsolated;
477  }
478  }
479  }
480  }
481  return isIsolated;
482  }
483 
484  //===========================================================================================================
485 
486  double chargeIsolationEcal(const edm::Event& iEvent,
487  const edm::EventSetup& iSetup,
488  const DetId& coreDet,
489  reco::TrackCollection::const_iterator trkItr,
490  edm::Handle<reco::TrackCollection> trkCollection,
491  const CaloGeometry* geo,
492  const CaloTopology* caloTopology,
493  TrackDetectorAssociator& associator,
494  TrackAssociatorParameters& parameters_,
495  int ieta,
496  int iphi,
497  const std::string& theTrackQuality,
498  bool debug) {
499  const CaloSubdetectorGeometry* barrelGeom = (geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel));
500  const CaloSubdetectorGeometry* endcapGeom = (geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap));
501 
502  std::vector<DetId> vdets = spr::matrixECALIds(coreDet, ieta, iphi, geo, caloTopology, debug);
503  if (debug)
504  edm::LogVerbatim("IsoTrack") << "chargeIsolation:: eta/phi/dets " << ieta << " " << iphi << " " << vdets.size();
505 
506  double maxNearP = -1.0;
507  reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
508 
509  // const DetId anyCell,
510  reco::TrackCollection::const_iterator trkItr2;
511  for (trkItr2 = trkCollection->begin(); trkItr2 != trkCollection->end(); ++trkItr2) {
512  const reco::Track* pTrack2 = &(*trkItr2);
513 
514  bool trkQuality = pTrack2->quality(trackQuality_);
515  if ((trkItr2 != trkItr) && trkQuality) {
516  const FreeTrajectoryState fts2 =
517  associator.getFreeTrajectoryState(&iSetup.getData(parameters_.bFieldToken), *pTrack2);
518  TrackDetMatchInfo info2 = associator.associate(iEvent, iSetup, fts2, parameters_);
519  const GlobalPoint point2(info2.trkGlobPosAtEcal.x(), info2.trkGlobPosAtEcal.y(), info2.trkGlobPosAtEcal.z());
520 
521  if (info2.isGoodEcal) {
522  if (std::abs(point2.eta()) < spr::etaBEEcal) {
523  const DetId anyCell = barrelGeom->getClosestCell(point2);
524  if (debug)
525  edm::LogVerbatim("IsoTrack")
526  << "chargeIsolation:: EB cell " << (EBDetId)(anyCell) << " for pt " << pTrack2->p();
527  if (!spr::chargeIsolation(anyCell, vdets)) {
528  if (maxNearP < pTrack2->p())
529  maxNearP = pTrack2->p();
530  }
531  } else {
532  const DetId anyCell = endcapGeom->getClosestCell(point2);
533  if (debug)
534  edm::LogVerbatim("IsoTrack")
535  << "chargeIsolation:: EE cell " << (EEDetId)(anyCell) << " for pt " << pTrack2->p();
536  if (!spr::chargeIsolation(anyCell, vdets)) {
537  if (maxNearP < pTrack2->p())
538  maxNearP = pTrack2->p();
539  }
540  }
541  } //info2.isGoodEcal
542  }
543  }
544  return maxNearP;
545  }
546 
547  //===========================================================================================================
548 
549  double chargeIsolationHcal(const edm::Event& iEvent,
550  const edm::EventSetup& iSetup,
551  reco::TrackCollection::const_iterator trkItr,
552  edm::Handle<reco::TrackCollection> trkCollection,
553  const DetId ClosestCell,
554  const HcalTopology* topology,
555  const CaloSubdetectorGeometry* gHB,
556  TrackDetectorAssociator& associator,
557  TrackAssociatorParameters& parameters_,
558  int ieta,
559  int iphi,
560  const std::string& theTrackQuality,
561  bool debug) {
562  std::vector<DetId> dets(1, ClosestCell);
563 
564  if (debug)
565  edm::LogVerbatim("IsoTrack") << static_cast<HcalDetId>(ClosestCell);
566  std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, false, debug);
567 
568  if (debug) {
569  for (unsigned int i = 0; i < vdets.size(); i++) {
570  edm::LogVerbatim("IsoTrack") << "HcalDetId in " << 2 * ieta + 1 << "x" << 2 * iphi + 1 << " "
571  << static_cast<HcalDetId>(vdets[i]);
572  }
573  }
574 
575  double maxNearP = -1.0;
576  reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
577 
578  reco::TrackCollection::const_iterator trkItr2;
579  for (trkItr2 = trkCollection->begin(); trkItr2 != trkCollection->end(); ++trkItr2) {
580  const reco::Track* pTrack2 = &(*trkItr2);
581 
582  bool trkQuality = pTrack2->quality(trackQuality_);
583  if ((trkItr2 != trkItr) && trkQuality) {
584  const FreeTrajectoryState fts2 =
585  associator.getFreeTrajectoryState(&iSetup.getData(parameters_.bFieldToken), *pTrack2);
586  TrackDetMatchInfo info2 = associator.associate(iEvent, iSetup, fts2, parameters_);
587  const GlobalPoint point2(info2.trkGlobPosAtHcal.x(), info2.trkGlobPosAtHcal.y(), info2.trkGlobPosAtHcal.z());
588 
589  if (debug)
590  edm::LogVerbatim("IsoTrack") << "Track2 (p,eta,phi) " << pTrack2->p() << " " << pTrack2->eta() << " "
591  << pTrack2->phi();
592 
593  if (info2.isGoodHcal) {
594  const DetId anyCell = gHB->getClosestCell(point2);
595  if (debug)
596  edm::LogVerbatim("IsoTrack") << "chargeIsolation:: HCAL cell " << static_cast<HcalDetId>(anyCell)
597  << " for pt " << pTrack2->p();
598  if (!spr::chargeIsolation(anyCell, vdets)) {
599  if (maxNearP < pTrack2->p())
600  maxNearP = pTrack2->p();
601  }
602  if (debug) {
603  edm::LogVerbatim("IsoTrack") << "maxNearP " << maxNearP << " thisCell " << static_cast<HcalDetId>(anyCell)
604  << " (" << info2.trkGlobPosAtHcal.x() << "," << info2.trkGlobPosAtHcal.y()
605  << "," << info2.trkGlobPosAtHcal.z() << ")";
606  }
607  }
608  }
609  }
610  return maxNearP;
611  }
612 
613 } // namespace spr
reco::TrackBase::quality
bool quality(const TrackQuality) const
Track quality.
Definition: TrackBase.h:552
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
spr::coneChargeIsolation
double coneChargeIsolation(const edm::Event &iEvent, const edm::EventSetup &iSetup, reco::TrackCollection::const_iterator trkItr, edm::Handle< reco::TrackCollection > trkCollection, TrackDetectorAssociator &associator, TrackAssociatorParameters &parameters_, const std::string &theTrackQuality, int &nNearTRKs, int &nLayers_maxNearP, int &trkQual_maxNearP, double &maxNearP_goodTrk, const GlobalPoint &hpoint1, const GlobalVector &trackMom, double dR)
Definition: ChargeIsolation.cc:242
spr::etaBEEcal
static const double etaBEEcal
Definition: CaloConstants.h:12
info
static const TGPicture * info(bool iBackgroundIsBlack)
Definition: FWCollectionSummaryWidget.cc:153
spr::getDistInPlaneTrackDir
double getDistInPlaneTrackDir(const GlobalPoint &caloPoint, const GlobalVector &caloVector, const GlobalPoint &rechitPoint, bool debug=false)
Definition: FindDistCone.cc:9
mps_fire.i
i
Definition: mps_fire.py:429
edm::EventSetup::getData
T const & getData(const ESGetToken< T, R > &iToken) const noexcept(false)
Definition: EventSetup.h:119
spr::matrixECALIds
void matrixECALIds(const DetId &det, int ieta, int iphi, const CaloGeometry *geo, const CaloTopology *caloTopology, std::vector< DetId > &vdets, bool debug=false, bool igNoreTransition=true)
Definition: MatrixECALDetIds.cc:16
spr::chargeIsolation
bool chargeIsolation(const DetId anyCell, std::vector< DetId > &vdets)
Definition: ChargeIsolation.cc:231
reco::TrackBase::TrackQuality
TrackQuality
track quality
Definition: TrackBase.h:150
spr::propagateHCAL
std::pair< math::XYZPoint, bool > propagateHCAL(const reco::Track *, const MagneticField *, bool debug=false)
Definition: CaloPropagateTrack.cc:720
reco::TrackBase::p
double p() const
momentum vector magnitude
Definition: TrackBase.h:631
reco::HitPattern::trackerLayersWithMeasurement
int trackerLayersWithMeasurement() const
Definition: HitPattern.cc:534
spr::chargeIsolationGenCone
double chargeIsolationGenCone(unsigned int trkIndex, std::vector< spr::propagatedGenParticleID > &trackIDs, double dR, int &nNearTRKs, bool debug=false)
Definition: ChargeIsolation.cc:336
spr::chargeIsolationCone
double chargeIsolationCone(unsigned int trkIndex, std::vector< spr::propagatedTrackDirection > &trkDirs, double dR, int &nNearTRKs, bool debug=false)
Definition: ChargeIsolation.cc:306
spr::chargeIsolationGenEcal
double chargeIsolationGenEcal(unsigned int trkIndex, std::vector< spr::propagatedGenParticleID > &trackIDs, const CaloGeometry *geo, const CaloTopology *caloTopology, int ieta, int iphi, bool debug=false)
Definition: ChargeIsolation.cc:47
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
TrackDetMatchInfo::trkGlobPosAtHcal
math::XYZPoint trkGlobPosAtHcal
Definition: TrackDetMatchInfo.h:39
Calorimetry_cff.bField
bField
Definition: Calorimetry_cff.py:284
HLT_2023v11_cff.navigator
navigator
Definition: HLT_2023v11_cff.py:9207
spr::chargeIsolationEcal
double chargeIsolationEcal(unsigned int trkIndex, std::vector< spr::propagatedTrackID > &vdetIds, const CaloGeometry *geo, const CaloTopology *caloTopology, int ieta, int iphi, bool debug=false)
Definition: ChargeIsolation.cc:12
iEvent
int iEvent
Definition: GenABIO.cc:224
reco::LeafCandidate::p
double p() const final
magnitude of momentum vector
Definition: LeafCandidate.h:123
TrackAssociatorParameters::bFieldToken
edm::ESGetToken< MagneticField, IdealMagneticFieldRecord > bFieldToken
Definition: TrackAssociatorParameters.h:136
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
reco::TrackBase::phi
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:649
DetId::subdetId
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:48
spr::propagateECAL
std::pair< math::XYZPoint, bool > propagateECAL(const reco::Track *, const MagneticField *, bool debug=false)
Definition: CaloPropagateTrack.cc:660
CaloSubdetectorGeometry::getClosestCell
virtual DetId getClosestCell(const GlobalPoint &r) const
Definition: CaloSubdetectorGeometry.cc:44
reco::TrackBase::eta
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:652
DetId
Definition: DetId.h:17
reco::TrackBase::qualityByName
static TrackQuality qualityByName(const std::string &name)
Definition: TrackBase.cc:126
debug
#define debug
Definition: HDRShower.cc:19
reco::TrackBase::hitPattern
const HitPattern & hitPattern() const
Access the hit pattern, indicating in which Tracker layers the track has hits.
Definition: TrackBase.h:504
HGC3DClusterGenMatchSelector_cfi.dR
dR
Definition: HGC3DClusterGenMatchSelector_cfi.py:7
LCToCPAssociation_cfi.associator
associator
Definition: LCToCPAssociation_cfi.py:4
spr::matrixHCALIds
std::vector< DetId > matrixHCALIds(std::vector< DetId > &dets, const HcalTopology *topology, int ieta, int iphi, bool includeHO=false, bool debug=false)
Definition: MatrixHCALDetIds.cc:13
TrackDetMatchInfo::trkGlobPosAtEcal
math::XYZPoint trkGlobPosAtEcal
Track position at different parts of the calorimeter.
Definition: TrackDetMatchInfo.h:38
spr::chargeIsolationHcal
double chargeIsolationHcal(unsigned int trkIndex, std::vector< spr::propagatedTrackID > &vdetIds, const HcalTopology *topology, int ieta, int iphi, bool debug=false)
Definition: ChargeIsolation.cc:140
CaloGeometry::getSubdetectorGeometry
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
Definition: CaloGeometry.cc:34
AlCaHLTBitMon_ParallelJobs.p
def p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
Generated for CMSSW Reference Manual by doxygen 1.8.14