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 "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
8 
9 #include <iostream>
10 
11 //#define EDM_ML_DEBUG
12 
13 namespace spr {
14 
15  double chargeIsolationEcal(unsigned int trkIndex,
16  std::vector<spr::propagatedTrackID>& vdetIds,
17  const CaloGeometry* geo,
18  const CaloTopology* caloTopology,
19  int ieta,
20  int iphi,
21  bool debug) {
22  const DetId coreDet = vdetIds[trkIndex].detIdECAL;
23 #ifdef EDM_ML_DEBUG
24  if (debug) {
25  if (coreDet.subdetId() == EcalBarrel)
26  std::cout << "DetId " << (EBDetId)(coreDet) << " Flag " << vdetIds[trkIndex].okECAL << std::endl;
27  else
28  std::cout << "DetId " << (EEDetId)(coreDet) << " Flag " << vdetIds[trkIndex].okECAL << std::endl;
29  }
30 #endif
31  double maxNearP = -1.0;
32  if (vdetIds[trkIndex].okECAL) {
33  std::vector<DetId> vdets = spr::matrixECALIds(coreDet, ieta, iphi, geo, caloTopology, debug);
34 #ifdef EDM_ML_DEBUG
35  if (debug)
36  std::cout << "chargeIsolationEcal:: eta/phi/dets " << ieta << " " << iphi << " " << vdets.size() << std::endl;
37 #endif
38 
39  for (unsigned int indx = 0; indx < vdetIds.size(); ++indx) {
40  if (indx != trkIndex && vdetIds[indx].ok && vdetIds[indx].okECAL) {
41  const DetId anyCell = vdetIds[indx].detIdECAL;
42  if (!spr::chargeIsolation(anyCell, vdets)) {
43  const reco::Track* pTrack = &(*(vdetIds[indx].trkItr));
44  if (maxNearP < pTrack->p())
45  maxNearP = pTrack->p();
46  }
47  }
48  }
49  }
50  return maxNearP;
51  }
52 
53  double chargeIsolationEcal(const DetId& coreDet,
54  reco::TrackCollection::const_iterator trkItr,
55  edm::Handle<reco::TrackCollection> trkCollection,
56  const CaloGeometry* geo,
57  const CaloTopology* caloTopology,
58  const MagneticField* bField,
59  int ieta,
60  int iphi,
61  const std::string& theTrackQuality,
62  bool debug) {
63  const CaloSubdetectorGeometry* barrelGeom = (geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel));
64  const CaloSubdetectorGeometry* endcapGeom = (geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap));
65 
66  std::vector<DetId> vdets = spr::matrixECALIds(coreDet, ieta, iphi, geo, caloTopology, debug);
67 #ifdef EDM_ML_DEBUG
68  if (debug)
69  std::cout << "chargeIsolation:: eta/phi/dets " << ieta << " " << iphi << " " << vdets.size() << std::endl;
70 #endif
71  double maxNearP = -1.0;
72  reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
73 
74  // const DetId anyCell,
75  reco::TrackCollection::const_iterator trkItr2;
76  for (trkItr2 = trkCollection->begin(); trkItr2 != trkCollection->end(); ++trkItr2) {
77  const reco::Track* pTrack2 = &(*trkItr2);
78 
79  bool trkQuality = (trackQuality_ != reco::TrackBase::undefQuality) ? (pTrack2->quality(trackQuality_)) : true;
80  if ((trkItr2 != trkItr) && trkQuality) {
81  std::pair<math::XYZPoint, bool> info = spr::propagateECAL(pTrack2, bField);
82  const GlobalPoint point2(info.first.x(), info.first.y(), info.first.z());
83 
84  if (info.second) {
85  if (std::abs(point2.eta()) < spr::etaBEEcal) {
86  const DetId anyCell = barrelGeom->getClosestCell(point2);
87  if (!spr::chargeIsolation(anyCell, vdets)) {
88 #ifdef EDM_ML_DEBUG
89  if (debug)
90  std::cout << "chargeIsolationEcal Cell " << (EBDetId)(anyCell) << " pt " << pTrack2->p() << std::endl;
91 #endif
92  if (maxNearP < pTrack2->p())
93  maxNearP = pTrack2->p();
94  }
95  } else {
96  if (endcapGeom) {
97  const DetId anyCell = endcapGeom->getClosestCell(point2);
98  if (!spr::chargeIsolation(anyCell, vdets)) {
99 #ifdef EDM_ML_DEBUG
100  if (debug)
101  std::cout << "chargeIsolationEcal Cell " << (EEDetId)(anyCell) << " pt " << pTrack2->p() << std::endl;
102 #endif
103  if (maxNearP < pTrack2->p())
104  maxNearP = pTrack2->p();
105  }
106  }
107  }
108  } //info.second
109  }
110  }
111  return maxNearP;
112  }
113 
114  double chargeIsolationHcal(unsigned int trkIndex,
115  std::vector<spr::propagatedTrackID>& vdetIds,
116  const HcalTopology* topology,
117  int ieta,
118  int iphi,
119  bool debug) {
120  std::vector<DetId> dets(1, vdetIds[trkIndex].detIdHCAL);
121 #ifdef EDM_ML_DEBUG
122  if (debug) {
123  std::cout << "DetId " << (HcalDetId)(dets[0]) << " Flag " << vdetIds[trkIndex].okHCAL << std::endl;
124  }
125 #endif
126  double maxNearP = -1.0;
127  if (vdetIds[trkIndex].okHCAL) {
128  std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, false, debug);
129 #ifdef EDM_ML_DEBUG
130  if (debug)
131  std::cout << "chargeIsolationHcal:: eta/phi/dets " << ieta << " " << iphi << " " << vdets.size() << std::endl;
132 #endif
133  for (unsigned indx = 0; indx < vdetIds.size(); ++indx) {
134  if (indx != trkIndex && vdetIds[indx].ok && vdetIds[indx].okHCAL) {
135  const DetId anyCell = vdetIds[indx].detIdHCAL;
136  if (!spr::chargeIsolation(anyCell, vdets)) {
137  const reco::Track* pTrack = &(*(vdetIds[indx].trkItr));
138 #ifdef EDM_ML_DEBUG
139  if (debug)
140  std::cout << "chargeIsolationHcal Cell " << (HcalDetId)(anyCell) << " pt " << pTrack->p() << std::endl;
141 #endif
142  if (maxNearP < pTrack->p())
143  maxNearP = pTrack->p();
144  }
145  }
146  }
147  }
148  return maxNearP;
149  }
150 
151  //===========================================================================================================
152 
153  double chargeIsolationHcal(reco::TrackCollection::const_iterator trkItr,
154  edm::Handle<reco::TrackCollection> trkCollection,
155  const DetId ClosestCell,
156  const HcalTopology* topology,
157  const CaloSubdetectorGeometry* gHB,
158  const MagneticField* bField,
159  int ieta,
160  int iphi,
161  const std::string& theTrackQuality,
162  bool debug) {
163  std::vector<DetId> dets(1, ClosestCell);
164 #ifdef EDM_ML_DEBUG
165  if (debug)
166  std::cout << (HcalDetId)ClosestCell << std::endl;
167 #endif
168  std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, false, debug);
169 
170 #ifdef EDM_ML_DEBUG
171  if (debug) {
172  for (unsigned int i = 0; i < vdets.size(); i++) {
173  std::cout << "HcalDetId in " << 2 * ieta + 1 << "x" << 2 * iphi + 1 << " " << (HcalDetId)vdets[i] << std::endl;
174  }
175  }
176 #endif
177  double maxNearP = -1.0;
178  reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
179 
180  reco::TrackCollection::const_iterator trkItr2;
181  for (trkItr2 = trkCollection->begin(); trkItr2 != trkCollection->end(); ++trkItr2) {
182  const reco::Track* pTrack2 = &(*trkItr2);
183 
184  bool trkQuality = (trackQuality_ != reco::TrackBase::undefQuality) ? (pTrack2->quality(trackQuality_)) : true;
185  if ((trkItr2 != trkItr) && trkQuality) {
186  std::pair<math::XYZPoint, bool> info = spr::propagateHCAL(pTrack2, bField);
187  const GlobalPoint point2(info.first.x(), info.first.y(), info.first.z());
188 
189 #ifdef EDM_ML_DEBUG
190  if (debug) {
191  std::cout << "Track2 (p,eta,phi) " << pTrack2->p() << " " << pTrack2->eta() << " " << pTrack2->phi()
192  << std::endl;
193  }
194 #endif
195  if (info.second) {
196  const DetId anyCell = gHB->getClosestCell(point2);
197  if (!spr::chargeIsolation(anyCell, vdets)) {
198  if (maxNearP < pTrack2->p())
199  maxNearP = pTrack2->p();
200  }
201 #ifdef EDM_ML_DEBUG
202  if (debug) {
203  std::cout << "maxNearP " << maxNearP << " thisCell " << (HcalDetId)anyCell << " (" << info.first.x() << ","
204  << info.first.y() << "," << info.first.z() << ")" << std::endl;
205  }
206 #endif
207  }
208  }
209  }
210  return maxNearP;
211  }
212 
213  bool chargeIsolation(const DetId anyCell, std::vector<DetId>& vdets) {
214  bool isIsolated = true;
215  for (unsigned int i = 0; i < vdets.size(); i++) {
216  if (anyCell == vdets[i]) {
217  isIsolated = false;
218  break;
219  }
220  }
221  return isIsolated;
222  }
223 
224  double coneChargeIsolation(const edm::Event& iEvent,
225  const edm::EventSetup& iSetup,
226  reco::TrackCollection::const_iterator trkItr,
227  edm::Handle<reco::TrackCollection> trkCollection,
228  TrackDetectorAssociator& associator,
229  TrackAssociatorParameters& parameters_,
230  const std::string& theTrackQuality,
231  int& nNearTRKs,
232  int& nLayers_maxNearP,
233  int& trkQual_maxNearP,
234  double& maxNearP_goodTrk,
235  const GlobalPoint& hpoint1,
236  const GlobalVector& trackMom,
237  double dR) {
238  nNearTRKs = 0;
239  nLayers_maxNearP = 0;
240  trkQual_maxNearP = -1;
241  maxNearP_goodTrk = -999.0;
242  double maxNearP = -999.0;
243  reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
244 
245  // Iterate over tracks
246  reco::TrackCollection::const_iterator trkItr2;
247  for (trkItr2 = trkCollection->begin(); trkItr2 != trkCollection->end(); ++trkItr2) {
248  // Get track
249  const reco::Track* pTrack2 = &(*trkItr2);
250 
251  // Get track qual, nlayers, and hit pattern
252  bool trkQuality = (trackQuality_ != reco::TrackBase::undefQuality) ? (pTrack2->quality(trackQuality_)) : true;
253  if (trkQuality)
254  trkQual_maxNearP = 1;
255  const reco::HitPattern& hitp = pTrack2->hitPattern();
256  nLayers_maxNearP = hitp.trackerLayersWithMeasurement();
257 
258  // Skip if the neighboring track candidate is the iso-track
259  // candidate
260  if (trkItr2 != trkItr) {
261  // Get propagator
262  const FreeTrajectoryState fts2 = associator.getFreeTrajectoryState(iSetup, *pTrack2);
263  TrackDetMatchInfo info2 = associator.associate(iEvent, iSetup, fts2, parameters_);
264 
265  // Make sure it reaches Hcal
266  if (info2.isGoodHcal) {
267  const GlobalPoint point2(info2.trkGlobPosAtHcal.x(), info2.trkGlobPosAtHcal.y(), info2.trkGlobPosAtHcal.z());
268 
269  int isConeChargedIso = spr::coneChargeIsolation(hpoint1, point2, trackMom, dR);
270 
271  if (isConeChargedIso == 0) {
272  nNearTRKs++;
273  if (maxNearP < pTrack2->p()) {
274  maxNearP = pTrack2->p();
275  if (trkQual_maxNearP > 0 && nLayers_maxNearP > 7 && maxNearP_goodTrk < pTrack2->p()) {
276  maxNearP_goodTrk = pTrack2->p();
277  }
278  }
279  }
280  }
281  }
282  } // Iterate over track loop
283 
284  return maxNearP;
285  }
286 
287  double chargeIsolationCone(unsigned int trkIndex,
288  std::vector<spr::propagatedTrackDirection>& trkDirs,
289  double dR,
290  int& nNearTRKs,
291  bool debug) {
292  double maxNearP = -1.0;
293  nNearTRKs = 0;
294  if (trkDirs[trkIndex].okHCAL) {
295 #ifdef EDM_ML_DEBUG
296  if (debug)
297  std::cout << "chargeIsolationCone with " << trkDirs.size() << " tracks " << std::endl;
298 #endif
299  for (unsigned int indx = 0; indx < trkDirs.size(); ++indx) {
300  if (indx != trkIndex && trkDirs[indx].ok && trkDirs[indx].okHCAL) {
301  int isConeChargedIso = spr::coneChargeIsolation(
302  trkDirs[trkIndex].pointHCAL, trkDirs[indx].pointHCAL, trkDirs[trkIndex].directionHCAL, dR);
303  if (isConeChargedIso == 0) {
304  nNearTRKs++;
305  const reco::Track* pTrack = &(*(trkDirs[indx].trkItr));
306  if (maxNearP < pTrack->p())
307  maxNearP = pTrack->p();
308  }
309  }
310  }
311  }
312 #ifdef EDM_ML_DEBUG
313  if (debug)
314  std::cout << "chargeIsolationCone Track " << trkDirs[trkIndex].okHCAL << " maxNearP " << maxNearP << std::endl;
315 #endif
316  return maxNearP;
317  }
318 
319  std::pair<double, double> chargeIsolationCone(unsigned int trkIndex,
320  std::vector<spr::propagatedTrackDirection>& trkDirs,
321  double dR,
322  bool debug) {
323  double maxNearP = -1.0;
324  double sumP = 0;
325  if (trkDirs[trkIndex].okHCAL) {
326 #ifdef EDM_ML_DEBUG
327  if (debug)
328  std::cout << "chargeIsolationCone with " << trkDirs.size() << " tracks " << std::endl;
329 #endif
330  for (unsigned int indx = 0; indx < trkDirs.size(); ++indx) {
331  if (indx != trkIndex && trkDirs[indx].ok && trkDirs[indx].okHCAL) {
332  int isConeChargedIso = spr::coneChargeIsolation(
333  trkDirs[trkIndex].pointHCAL, trkDirs[indx].pointHCAL, trkDirs[trkIndex].directionHCAL, dR);
334  if (isConeChargedIso == 0) {
335  const reco::Track* pTrack = &(*(trkDirs[indx].trkItr));
336  sumP += (pTrack->p());
337  if (maxNearP < pTrack->p())
338  maxNearP = pTrack->p();
339  }
340  }
341  }
342  }
343 #ifdef EDM_ML_DEBUG
344  if (debug)
345  std::cout << "chargeIsolationCone Track " << trkDirs[trkIndex].okHCAL << " maxNearP " << maxNearP << ":" << sumP
346  << std::endl;
347 #endif
348  return std::pair<double, double>(maxNearP, sumP);
349  }
350 
351  int coneChargeIsolation(const GlobalPoint& hpoint1,
352  const GlobalPoint& point2,
353  const GlobalVector& trackMom,
354  double dR) {
355  int isIsolated = 1;
356  if (spr::getDistInPlaneTrackDir(hpoint1, trackMom, point2) > dR)
357  isIsolated = 1;
358  else
359  isIsolated = 0;
360  return isIsolated;
361  }
362 
363 } // namespace spr
Vector3DBase
Definition: Vector3DBase.h:8
mps_fire.i
i
Definition: mps_fire.py:355
FreeTrajectoryState.h
spr
Definition: CaloConstants.h:6
reco::TrackBase::p
double p() const
momentum vector magnitude
Definition: TrackBase.h:605
EBDetId
Definition: EBDetId.h:17
CaloConstants.h
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
reco::TrackBase::undefQuality
Definition: TrackBase.h:151
gather_cfg.cout
cout
Definition: gather_cfg.py:144
HcalTopology
Definition: HcalTopology.h:26
TrackDetMatchInfo::trkGlobPosAtHcal
math::XYZPoint trkGlobPosAtHcal
Definition: TrackDetMatchInfo.h:40
CaloGeometry::getSubdetectorGeometry
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
Definition: CaloGeometry.cc:34
reco::TrackBase::TrackQuality
TrackQuality
track quality
Definition: TrackBase.h:150
info
static const TGPicture * info(bool iBackgroundIsBlack)
Definition: FWCollectionSummaryWidget.cc:152
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
convertSQLiteXML.ok
bool ok
Definition: convertSQLiteXML.py:98
edm::Handle< reco::TrackCollection >
spr::getDistInPlaneTrackDir
double getDistInPlaneTrackDir(const GlobalPoint &caloPoint, const GlobalVector &caloVector, const GlobalPoint &rechitPoint, bool debug=false)
Definition: FindDistCone.cc:12
TrackDetMatchInfo::isGoodHcal
bool isGoodHcal
Definition: TrackDetMatchInfo.h:48
CaloTopology
Definition: CaloTopology.h:19
EcalBarrel
Definition: EcalSubdetector.h:10
LEDCalibrationChannels.iphi
iphi
Definition: LEDCalibrationChannels.py:64
DetId
Definition: DetId.h:17
CaloGeometry
Definition: CaloGeometry.h:21
reco::HitPattern
Definition: HitPattern.h:147
debug
#define debug
Definition: HDRShower.cc:19
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:15
reco::HitPattern::trackerLayersWithMeasurement
int trackerLayersWithMeasurement() const
Definition: HitPattern.cc:513
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:114
ecaldqm::topology
const CaloTopology * topology(nullptr)
ChargeIsolation.h
reco::Track
Definition: Track.h:27
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:224
MatrixECALDetIds.h
spr::etaBEEcal
static const double etaBEEcal
Definition: CaloConstants.h:12
ctfWithMaterialTrackMCMatch_cfi.associator
associator
Definition: ctfWithMaterialTrackMCMatch_cfi.py:7
Point3DBase< float, GlobalTag >
spr::chargeIsolationCone
double chargeIsolationCone(unsigned int trkIndex, std::vector< spr::propagatedTrackDirection > &trkDirs, double dR, int &nNearTRKs, bool debug=false)
Definition: ChargeIsolation.cc:287
reco::TrackBase::phi
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:620
DetId::subdetId
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector'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:699
EEDetId
Definition: EEDetId.h:14
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
EcalEndcap
Definition: EcalSubdetector.h:10
MatrixHCALDetIds.h
LEDCalibrationChannels.ieta
ieta
Definition: LEDCalibrationChannels.py:63
spr::chargeIsolation
bool chargeIsolation(const DetId anyCell, std::vector< DetId > &vdets)
Definition: ChargeIsolation.cc:213
TrackDetectorAssociator
Definition: TrackDetectorAssociator.h:49
reco::TrackBase::eta
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:623
HcalDetId
Definition: HcalDetId.h:12
iEvent
int iEvent
Definition: GenABIO.cc:224
edm::EventSetup
Definition: EventSetup.h:57
CaloSubdetectorGeometry::getClosestCell
virtual DetId getClosestCell(const GlobalPoint &r) const
Definition: CaloSubdetectorGeometry.cc:44
DetId::Ecal
Definition: DetId.h:27
reco::TrackBase::qualityByName
static TrackQuality qualityByName(const std::string &name)
Definition: TrackBase.cc:126
reco::TrackBase::hitPattern
const HitPattern & hitPattern() const
Access the hit pattern, indicating in which Tracker layers the track has hits.
Definition: TrackBase.h:489
TrackDetMatchInfo
Definition: TrackDetMatchInfo.h:14
FreeTrajectoryState
Definition: FreeTrajectoryState.h:27
Calorimetry_cff.bField
bField
Definition: Calorimetry_cff.py:292
CaloPropagateTrack.h
FindDistCone.h
CaloSubdetectorGeometry
Definition: CaloSubdetectorGeometry.h:22
spr::propagateHCAL
std::pair< math::XYZPoint, bool > propagateHCAL(const reco::Track *, const MagneticField *, bool debug=false)
Definition: CaloPropagateTrack.cc:759
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
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:15
HGC3DClusterGenMatchSelector_cfi.dR
dR
Definition: HGC3DClusterGenMatchSelector_cfi.py:7
edm::Event
Definition: Event.h:73
MagneticField
Definition: MagneticField.h:19
reco::TrackBase::quality
bool quality(const TrackQuality) const
Track quality.
Definition: TrackBase.h:537
TrackAssociatorParameters
Definition: TrackAssociatorParameters.h:34
Generated for CMSSW Reference Manual by doxygen 1.8.16