CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
CaloPropagateTrack.cc
Go to the documentation of this file.
1 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
2 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
3 #include "DataFormats/GeometrySurface/interface/Plane.h"
4 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
5 #include "DataFormats/EcalDetId/interface/EBDetId.h"
6 #include "DataFormats/EcalDetId/interface/EEDetId.h"
7 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
8 
9 #include "FWCore/MessageLogger/interface/MessageLogger.h"
10 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
11 
12 #include "Calibration/IsolatedParticles/interface/CaloConstants.h"
13 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
14 
15 #include <iostream>
16 #include <sstream>
17 
18 namespace spr {
19 
20  std::vector<spr::propagatedTrackID> propagateCosmicCALO(edm::Handle<reco::TrackCollection>& trkCollection,
21  const CaloGeometry* geo,
22  const MagneticField* bField,
23  const std::string& theTrackQuality,
24  bool debug) {
25  const CaloSubdetectorGeometry* barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
26  const CaloSubdetectorGeometry* endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
27  const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
28  reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
29  std::vector<spr::propagatedTrackID> vdets;
30 
31  unsigned int indx;
32  reco::TrackCollection::const_iterator trkItr;
33  for (trkItr = trkCollection->begin(), indx = 0; trkItr != trkCollection->end(); ++trkItr, indx++) {
34  const reco::Track* pTrack = &(*trkItr);
35  spr::propagatedTrackID vdet;
36  vdet.trkItr = trkItr;
37  vdet.ok = (trackQuality_ != reco::TrackBase::undefQuality) ? (pTrack->quality(trackQuality_)) : true;
38  vdet.detIdECAL = DetId(0);
39  vdet.detIdHCAL = DetId(0);
40  vdet.detIdEHCAL = DetId(0);
41  if (debug)
42  edm::LogVerbatim("IsoTrack") << "Propagate track " << indx << " p " << trkItr->p() << " eta " << trkItr->eta()
43  << " phi " << trkItr->phi() << " Flag " << vdet.ok;
44  GlobalPoint vertex;
45  GlobalVector momentum;
46  int charge(pTrack->charge());
47  if (((pTrack->innerPosition()).Perp2()) < ((pTrack->outerPosition()).Perp2())) {
48  vertex = GlobalPoint(
49  ((pTrack->innerPosition()).X()), ((pTrack->innerPosition()).Y()), ((pTrack->innerPosition()).Z()));
50  momentum = GlobalVector(
51  ((pTrack->innerMomentum()).X()), ((pTrack->innerMomentum()).Y()), ((pTrack->innerMomentum()).Z()));
52  } else {
53  vertex = GlobalPoint(
54  ((pTrack->outerPosition()).X()), ((pTrack->outerPosition()).Y()), ((pTrack->outerPosition()).Z()));
55  momentum = GlobalVector(
56  ((pTrack->outerMomentum()).X()), ((pTrack->outerMomentum()).Y()), ((pTrack->outerMomentum()).Z()));
57  }
58  if (debug)
59  edm::LogVerbatim("IsoTrack") << "Track charge " << charge << " p " << momentum << " position " << vertex;
60  std::pair<math::XYZPoint, bool> info = spr::propagateECAL(vertex, momentum, charge, bField, debug);
61  if (debug)
62  edm::LogVerbatim("IsoTrack") << "Propagate to ECAL " << info.second << " at (" << info.first.x() << ", "
63  << info.first.y() << ", " << info.first.z() << ")";
64 
65  vdet.okECAL = info.second;
66  if (vdet.okECAL) {
67  const GlobalPoint point(info.first.x(), info.first.y(), info.first.z());
68  vdet.etaECAL = point.eta();
69  vdet.phiECAL = point.phi();
70  if (std::abs(point.eta()) < spr::etaBEEcal) {
71  vdet.detIdECAL = barrelGeom->getClosestCell(point);
72  } else {
73  if (endcapGeom)
74  vdet.detIdECAL = endcapGeom->getClosestCell(point);
75  else
76  vdet.okECAL = false;
77  }
78  vdet.detIdEHCAL = gHB->getClosestCell(point);
79  if (debug) {
80  std::ostringstream st1;
81  if (std::abs(point.eta()) < spr::etaBEEcal)
82  st1 << EBDetId(vdet.detIdECAL);
83  else
84  st1 << EEDetId(vdet.detIdECAL);
85  edm::LogVerbatim("IsoTrack") << "Point at ECAL (" << vdet.etaECAL << ", " << vdet.phiECAL << " " << st1.str()
86  << " " << HcalDetId(vdet.detIdEHCAL);
87  }
88  }
89  info = spr::propagateHCAL(vertex, momentum, charge, bField, debug);
90  if (debug)
91  edm::LogVerbatim("IsoTrack") << "Propagate to HCAL " << info.second << " at (" << info.first.x() << ", "
92  << info.first.y() << ", " << info.first.z() << ")";
93  vdet.okHCAL = info.second;
94  if (vdet.okHCAL) {
95  const GlobalPoint point(info.first.x(), info.first.y(), info.first.z());
96  vdet.etaHCAL = point.eta();
97  vdet.phiHCAL = point.phi();
98  vdet.detIdHCAL = gHB->getClosestCell(point);
99  }
100  if (debug) {
101  std::ostringstream st1;
102  if (vdet.detIdECAL.subdetId() == EcalBarrel)
103  st1 << (EBDetId)(vdet.detIdECAL);
104  else
105  st1 << (EEDetId)(vdet.detIdECAL);
106  edm::LogVerbatim("IsoTrack") << "Track [" << indx << "] Flag: " << vdet.ok << " ECAL (" << vdet.okECAL << ") "
107  << st1.str() << " HCAL (" << vdet.okHCAL << ") " << (HcalDetId)(vdet.detIdHCAL)
108  << " Or " << (HcalDetId)(vdet.detIdEHCAL);
109  }
110  vdets.push_back(vdet);
111  }
112 
113  if (debug) {
114  edm::LogVerbatim("IsoTrack") << "propagateCALO:: for " << vdets.size() << " tracks";
115  for (unsigned int i = 0; i < vdets.size(); ++i) {
116  std::ostringstream st1;
117  if (vdets[i].detIdECAL.subdetId() == EcalBarrel) {
118  st1 << (EBDetId)(vdets[i].detIdECAL);
119  } else {
120  st1 << (EEDetId)(vdets[i].detIdECAL);
121  }
122  edm::LogVerbatim("IsoTrack") << "Track [" << i << "] Flag: " << vdets[i].ok << " ECAL (" << vdets[i].okECAL
123  << ") " << st1.str() << " HCAL (" << vdets[i].okHCAL << ") "
124  << (HcalDetId)(vdets[i].detIdHCAL) << " Or " << (HcalDetId)(vdets[i].detIdEHCAL);
125  }
126  }
127  return vdets;
128  }
129 
130  std::vector<spr::propagatedTrackID> propagateCALO(edm::Handle<reco::TrackCollection>& trkCollection,
131  const CaloGeometry* geo,
132  const MagneticField* bField,
133  const std::string& theTrackQuality,
134  bool debug) {
135  std::vector<spr::propagatedTrackID> vdets;
136  spr::propagateCALO(trkCollection, geo, bField, theTrackQuality, vdets, debug);
137  return vdets;
138  }
139 
140  void propagateCALO(edm::Handle<reco::TrackCollection>& trkCollection,
141  const CaloGeometry* geo,
142  const MagneticField* bField,
143  const std::string& theTrackQuality,
144  std::vector<spr::propagatedTrackID>& vdets,
145  bool debug) {
146  const CaloSubdetectorGeometry* barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
147  const CaloSubdetectorGeometry* endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
148  const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
149  reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
150 
151  unsigned int indx;
152  reco::TrackCollection::const_iterator trkItr;
153  for (trkItr = trkCollection->begin(), indx = 0; trkItr != trkCollection->end(); ++trkItr, indx++) {
154  const reco::Track* pTrack = &(*trkItr);
155  spr::propagatedTrackID vdet;
156  vdet.trkItr = trkItr;
157  vdet.ok = (trackQuality_ != reco::TrackBase::undefQuality) ? (pTrack->quality(trackQuality_)) : true;
158  vdet.detIdECAL = DetId(0);
159  vdet.detIdHCAL = DetId(0);
160  vdet.detIdEHCAL = DetId(0);
161  if (debug)
162  edm::LogVerbatim("IsoTrack") << "Propagate track " << indx << " p " << trkItr->p() << " eta " << trkItr->eta()
163  << " phi " << trkItr->phi() << " Flag " << vdet.ok;
164  std::pair<math::XYZPoint, bool> info = spr::propagateECAL(pTrack, bField, debug);
165  vdet.okECAL = info.second;
166  if (vdet.okECAL) {
167  const GlobalPoint point(info.first.x(), info.first.y(), info.first.z());
168  vdet.etaECAL = point.eta();
169  vdet.phiECAL = point.phi();
170  if (std::abs(point.eta()) < spr::etaBEEcal) {
171  vdet.detIdECAL = barrelGeom->getClosestCell(point);
172  } else {
173  if (endcapGeom)
174  vdet.detIdECAL = endcapGeom->getClosestCell(point);
175  else
176  vdet.okECAL = false;
177  }
178  vdet.detIdEHCAL = gHB->getClosestCell(point);
179  }
180  info = spr::propagateHCAL(pTrack, bField, debug);
181  vdet.okHCAL = info.second;
182  if (vdet.okHCAL) {
183  const GlobalPoint point(info.first.x(), info.first.y(), info.first.z());
184  vdet.etaHCAL = point.eta();
185  vdet.phiHCAL = point.phi();
186  vdet.detIdHCAL = gHB->getClosestCell(point);
187  }
188  if (debug) {
189  std::ostringstream st1;
190  if (vdet.detIdECAL.subdetId() == EcalBarrel)
191  st1 << (EBDetId)(vdet.detIdECAL);
192  else
193  st1 << (EEDetId)(vdet.detIdECAL);
194  edm::LogVerbatim("IsoTrack") << "Track [" << indx << "] Flag: " << vdet.ok << " ECAL (" << vdet.okECAL << ") "
195  << st1.str() << " HCAL (" << vdet.okHCAL << ") " << (HcalDetId)(vdet.detIdHCAL)
196  << " Or " << (HcalDetId)(vdet.detIdEHCAL);
197  }
198  vdets.push_back(vdet);
199  }
200  if (debug) {
201  edm::LogVerbatim("IsoTrack") << "propagateCALO:: for " << vdets.size() << " tracks";
202  for (unsigned int i = 0; i < vdets.size(); ++i) {
203  std::ostringstream st1;
204  if (vdets[i].detIdECAL.subdetId() == EcalBarrel) {
205  st1 << (EBDetId)(vdets[i].detIdECAL);
206  } else {
207  st1 << (EEDetId)(vdets[i].detIdECAL);
208  }
209  edm::LogVerbatim("IsoTrack") << "Track [" << i << "] Flag: " << vdets[i].ok << " ECAL (" << vdets[i].okECAL
210  << ") " << st1.str() << " HCAL (" << vdets[i].okHCAL << ") "
211  << (HcalDetId)(vdets[i].detIdHCAL) << " Or " << (HcalDetId)(vdets[i].detIdEHCAL);
212  }
213  }
214  }
215 
216  void propagateCALO(edm::Handle<reco::TrackCollection>& trkCollection,
217  const CaloGeometry* geo,
218  const MagneticField* bField,
219  const std::string& theTrackQuality,
220  std::vector<spr::propagatedTrackDirection>& trkDir,
221  bool debug) {
222  const CaloSubdetectorGeometry* barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
223  const CaloSubdetectorGeometry* endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
224  const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
225  reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
226 
227  unsigned int indx;
228  reco::TrackCollection::const_iterator trkItr;
229  for (trkItr = trkCollection->begin(), indx = 0; trkItr != trkCollection->end(); ++trkItr, indx++) {
230  const reco::Track* pTrack = &(*trkItr);
231  spr::propagatedTrackDirection trkD;
232  trkD.trkItr = trkItr;
233  trkD.ok = (trackQuality_ != reco::TrackBase::undefQuality) ? (pTrack->quality(trackQuality_)) : true;
234  trkD.detIdECAL = DetId(0);
235  trkD.detIdHCAL = DetId(0);
236  trkD.detIdEHCAL = DetId(0);
237  if (debug)
238  edm::LogVerbatim("IsoTrack") << "Propagate track " << indx << " p " << trkItr->p() << " eta " << trkItr->eta()
239  << " phi " << trkItr->phi() << " Flag " << trkD.ok;
240  spr::propagatedTrack info = spr::propagateTrackToECAL(pTrack, bField, debug);
241  GlobalPoint point(info.point.x(), info.point.y(), info.point.z());
242  trkD.okECAL = info.ok;
243  trkD.pointECAL = point;
244  trkD.directionECAL = info.direction;
245  if (trkD.okECAL) {
246  if (std::abs(info.point.eta()) < spr::etaBEEcal) {
247  trkD.detIdECAL = barrelGeom->getClosestCell(point);
248  } else {
249  if (endcapGeom)
250  trkD.detIdECAL = endcapGeom->getClosestCell(point);
251  else
252  trkD.okECAL = false;
253  }
254  trkD.detIdEHCAL = gHB->getClosestCell(point);
255  }
256  info = spr::propagateTrackToHCAL(pTrack, bField, debug);
257  point = GlobalPoint(info.point.x(), info.point.y(), info.point.z());
258  trkD.okHCAL = info.ok;
259  trkD.pointHCAL = point;
260  trkD.directionHCAL = info.direction;
261  if (trkD.okHCAL) {
262  trkD.detIdHCAL = gHB->getClosestCell(point);
263  }
264  trkDir.push_back(trkD);
265  }
266  if (debug) {
267  edm::LogVerbatim("IsoTrack") << "propagateCALO:: for " << trkDir.size() << " tracks";
268  for (unsigned int i = 0; i < trkDir.size(); ++i) {
269  std::ostringstream st1, st2;
270  if (trkDir[i].okECAL) {
271  st1 << " point " << trkDir[i].pointECAL << " direction " << trkDir[i].directionECAL << " ";
272  if (trkDir[i].detIdECAL.subdetId() == EcalBarrel) {
273  st1 << (EBDetId)(trkDir[i].detIdECAL);
274  } else {
275  st1 << (EEDetId)(trkDir[i].detIdECAL);
276  }
277  }
278  if (trkDir[i].okHCAL) {
279  st2 << " point " << trkDir[i].pointHCAL << " direction " << trkDir[i].directionHCAL << " "
280  << (HcalDetId)(trkDir[i].detIdHCAL);
281  }
282  edm::LogVerbatim("IsoTrack") << "Track [" << i << "] Flag: " << trkDir[i].ok << " ECAL (" << trkDir[i].okECAL
283  << ")" << st1.str() << " HCAL (" << trkDir[i].okHCAL << ")" << st2.str() << " Or "
284  << (HcalDetId)(trkDir[i].detIdEHCAL);
285  }
286  }
287  }
288 
289  spr::propagatedTrackID propagateCALO(const reco::Track* pTrack,
290  const CaloGeometry* geo,
291  const MagneticField* bField,
292  bool debug) {
293  const CaloSubdetectorGeometry* barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
294  const CaloSubdetectorGeometry* endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
295  const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
296 
297  spr::propagatedTrackID vdet;
298  vdet.ok = true;
299  vdet.detIdECAL = DetId(0);
300  vdet.detIdHCAL = DetId(0);
301  vdet.detIdEHCAL = DetId(0);
302  if (debug)
303  edm::LogVerbatim("IsoTrack") << "Propagate track: p " << pTrack->p() << " eta " << pTrack->eta() << " phi "
304  << pTrack->phi() << " Flag " << vdet.ok;
305  std::pair<math::XYZPoint, bool> info = spr::propagateECAL(pTrack, bField, debug);
306  vdet.okECAL = info.second;
307  if (vdet.okECAL) {
308  const GlobalPoint point(info.first.x(), info.first.y(), info.first.z());
309  vdet.etaECAL = point.eta();
310  vdet.phiECAL = point.phi();
311  if (std::abs(point.eta()) < spr::etaBEEcal) {
312  vdet.detIdECAL = barrelGeom->getClosestCell(point);
313  } else {
314  if (endcapGeom)
315  vdet.detIdECAL = endcapGeom->getClosestCell(point);
316  else
317  vdet.okECAL = false;
318  }
319  vdet.detIdEHCAL = gHB->getClosestCell(point);
320  }
321  info = spr::propagateHCAL(pTrack, bField, debug);
322  vdet.okHCAL = info.second;
323  if (vdet.okHCAL) {
324  const GlobalPoint point(info.first.x(), info.first.y(), info.first.z());
325  vdet.etaHCAL = point.eta();
326  vdet.phiHCAL = point.phi();
327  vdet.detIdHCAL = gHB->getClosestCell(point);
328  }
329  if (debug) {
330  edm::LogVerbatim("IsoTrack") << "propagateCALO:: for 1 track";
331  std::ostringstream st1;
332  if (vdet.detIdECAL.subdetId() == EcalBarrel) {
333  st1 << (EBDetId)(vdet.detIdECAL);
334  } else {
335  st1 << (EEDetId)(vdet.detIdECAL);
336  }
337  edm::LogVerbatim("IsoTrack") << "Track [0] Flag: " << vdet.ok << " ECAL (" << vdet.okECAL << ") " << st1.str()
338  << " HCAL (" << vdet.okHCAL << ") " << (HcalDetId)(vdet.detIdHCAL) << " Or "
339  << (HcalDetId)(vdet.detIdEHCAL);
340  }
341  return vdet;
342  }
343 
344  std::vector<spr::propagatedGenTrackID> propagateCALO(const HepMC::GenEvent* genEvent,
345  const ParticleDataTable* pdt,
346  const CaloGeometry* geo,
347  const MagneticField* bField,
348  double etaMax,
349  bool debug) {
350  const CaloSubdetectorGeometry* barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
351  const CaloSubdetectorGeometry* endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
352  const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
353 
354  std::vector<spr::propagatedGenTrackID> trkDir;
355  unsigned int indx;
356  HepMC::GenEvent::particle_const_iterator p;
357  for (p = genEvent->particles_begin(), indx = 0; p != genEvent->particles_end(); ++p, ++indx) {
358  spr::propagatedGenTrackID trkD;
359  trkD.trkItr = p;
360  trkD.detIdECAL = DetId(0);
361  trkD.detIdHCAL = DetId(0);
362  trkD.detIdEHCAL = DetId(0);
363  trkD.pdgId = ((*p)->pdg_id());
364  trkD.charge = ((pdt->particle(trkD.pdgId))->ID().threeCharge()) / 3;
365  const GlobalVector momentum = GlobalVector((*p)->momentum().px(), (*p)->momentum().py(), (*p)->momentum().pz());
366  if (debug)
367  edm::LogVerbatim("IsoTrack") << "Propagate track " << indx << " pdg " << trkD.pdgId << " charge " << trkD.charge
368  << " p " << momentum;
369  // consider stable particles
370  if ((*p)->status() == 1 && std::abs((*p)->momentum().eta()) < etaMax) {
371  const GlobalPoint vertex = GlobalPoint(0.1 * (*p)->production_vertex()->position().x(),
372  0.1 * (*p)->production_vertex()->position().y(),
373  0.1 * (*p)->production_vertex()->position().z());
374  trkD.ok = true;
375  spr::propagatedTrack info = spr::propagateCalo(
376  vertex, momentum, trkD.charge, bField, spr::zFrontEE, spr::rFrontEB, spr::etaBEEcal, debug);
377  GlobalPoint point(info.point.x(), info.point.y(), info.point.z());
378  trkD.okECAL = info.ok;
379  trkD.pointECAL = point;
380  trkD.directionECAL = info.direction;
381  if (trkD.okECAL) {
382  if (std::abs(info.point.eta()) < spr::etaBEEcal) {
383  trkD.detIdECAL = barrelGeom->getClosestCell(point);
384  } else {
385  if (endcapGeom)
386  trkD.detIdECAL = endcapGeom->getClosestCell(point);
387  else
388  trkD.okECAL = false;
389  }
390  trkD.detIdEHCAL = gHB->getClosestCell(point);
391  }
392 
393  info = spr::propagateCalo(
394  vertex, momentum, trkD.charge, bField, spr::zFrontHE, spr::rFrontHB, spr::etaBEHcal, debug);
395  point = GlobalPoint(info.point.x(), info.point.y(), info.point.z());
396  trkD.okHCAL = info.ok;
397  trkD.pointHCAL = point;
398  trkD.directionHCAL = info.direction;
399  if (trkD.okHCAL) {
400  trkD.detIdHCAL = gHB->getClosestCell(point);
401  }
402  }
403  trkDir.push_back(trkD);
404  }
405  if (debug) {
406  edm::LogVerbatim("IsoTrack") << "propagateCALO:: for " << trkDir.size() << " tracks";
407  for (unsigned int i = 0; i < trkDir.size(); ++i) {
408  std::ostringstream st1, st2;
409  if (trkDir[i].okECAL) {
410  st1 << " point " << trkDir[i].pointECAL << " direction " << trkDir[i].directionECAL << " ";
411  if (trkDir[i].detIdECAL.subdetId() == EcalBarrel) {
412  st1 << (EBDetId)(trkDir[i].detIdECAL);
413  } else {
414  st1 << (EEDetId)(trkDir[i].detIdECAL);
415  }
416  }
417  st2 << " HCAL (" << trkDir[i].okHCAL << ")";
418  if (trkDir[i].okHCAL) {
419  st2 << " point " << trkDir[i].pointHCAL << " direction " << trkDir[i].directionHCAL << " "
420  << (HcalDetId)(trkDir[i].detIdHCAL);
421  }
422  edm::LogVerbatim("IsoTrack") << "Track [" << i << "] Flag: " << trkDir[i].ok << " ECAL (" << trkDir[i].okECAL
423  << ")" << st1.str() << st2.str() << " Or " << (HcalDetId)(trkDir[i].detIdEHCAL);
424  }
425  }
426  return trkDir;
427  }
428 
429  std::vector<spr::propagatedGenParticleID> propagateCALO(edm::Handle<reco::GenParticleCollection>& genParticles,
430  const ParticleDataTable* pdt,
431  const CaloGeometry* geo,
432  const MagneticField* bField,
433  double etaMax,
434  bool debug) {
435  const CaloSubdetectorGeometry* barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
436  const CaloSubdetectorGeometry* endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
437  const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
438 
439  std::vector<spr::propagatedGenParticleID> trkDir;
440  unsigned int indx;
441  reco::GenParticleCollection::const_iterator p;
442  for (p = genParticles->begin(), indx = 0; p != genParticles->end(); ++p, ++indx) {
443  spr::propagatedGenParticleID trkD;
444  trkD.trkItr = p;
445  trkD.detIdECAL = DetId(0);
446  trkD.detIdHCAL = DetId(0);
447  trkD.detIdEHCAL = DetId(0);
448  trkD.pdgId = (p->pdgId());
449  trkD.charge = p->charge();
450  const GlobalVector momentum = GlobalVector(p->momentum().x(), p->momentum().y(), p->momentum().z());
451  if (debug)
452  edm::LogVerbatim("IsoTrack") << "Propagate track " << indx << " pdg " << trkD.pdgId << " charge " << trkD.charge
453  << " p " << momentum;
454  // consider stable particles
455  if (p->status() == 1 && std::abs(momentum.eta()) < etaMax) {
456  const GlobalPoint vertex = GlobalPoint(p->vertex().x(), p->vertex().y(), p->vertex().z());
457  trkD.ok = true;
458  spr::propagatedTrack info = spr::propagateCalo(
459  vertex, momentum, trkD.charge, bField, spr::zFrontEE, spr::rFrontEB, spr::etaBEEcal, debug);
460  GlobalPoint point(info.point.x(), info.point.y(), info.point.z());
461  trkD.okECAL = info.ok;
462  trkD.pointECAL = point;
463  trkD.directionECAL = info.direction;
464  if (trkD.okECAL) {
465  if (std::abs(info.point.eta()) < spr::etaBEEcal) {
466  trkD.detIdECAL = barrelGeom->getClosestCell(point);
467  } else {
468  if (endcapGeom)
469  trkD.detIdECAL = endcapGeom->getClosestCell(point);
470  else
471  trkD.okECAL = false;
472  }
473  trkD.detIdEHCAL = gHB->getClosestCell(point);
474  }
475 
476  info = spr::propagateCalo(
477  vertex, momentum, trkD.charge, bField, spr::zFrontHE, spr::rFrontHB, spr::etaBEHcal, debug);
478  point = GlobalPoint(info.point.x(), info.point.y(), info.point.z());
479  trkD.okHCAL = info.ok;
480  trkD.pointHCAL = point;
481  trkD.directionHCAL = info.direction;
482  if (trkD.okHCAL) {
483  trkD.detIdHCAL = gHB->getClosestCell(point);
484  }
485  }
486  trkDir.push_back(trkD);
487  }
488  if (debug) {
489  edm::LogVerbatim("IsoTrack") << "propagateCALO:: for " << trkDir.size() << " tracks" << std::endl;
490  for (unsigned int i = 0; i < trkDir.size(); ++i) {
491  std::ostringstream st1, st2;
492  if (trkDir[i].okECAL) {
493  st1 << " point " << trkDir[i].pointECAL << " direction " << trkDir[i].directionECAL << " ";
494  if (trkDir[i].detIdECAL.subdetId() == EcalBarrel) {
495  st1 << (EBDetId)(trkDir[i].detIdECAL);
496  } else {
497  st1 << (EEDetId)(trkDir[i].detIdECAL);
498  }
499  }
500  st2 << " HCAL (" << trkDir[i].okHCAL << ")";
501  if (trkDir[i].okHCAL) {
502  st2 << " point " << trkDir[i].pointHCAL << " direction " << trkDir[i].directionHCAL << " "
503  << (HcalDetId)(trkDir[i].detIdHCAL);
504  }
505  edm::LogVerbatim("IsoTrack") << "Track [" << i << "] Flag: " << trkDir[i].ok << " ECAL (" << trkDir[i].okECAL
506  << ")" << st1.str() << st2.str() << " Or " << (HcalDetId)(trkDir[i].detIdEHCAL);
507  }
508  }
509  return trkDir;
510  }
511 
512  spr::propagatedTrackDirection propagateCALO(unsigned int thisTrk,
513  edm::Handle<edm::SimTrackContainer>& SimTk,
514  edm::Handle<edm::SimVertexContainer>& SimVtx,
515  const CaloGeometry* geo,
516  const MagneticField* bField,
517  bool debug) {
518  const CaloSubdetectorGeometry* barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
519  const CaloSubdetectorGeometry* endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
520  const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
521 
522  spr::trackAtOrigin trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
523  spr::propagatedTrackDirection trkD;
524  trkD.ok = trk.ok;
525  trkD.detIdECAL = DetId(0);
526  trkD.detIdHCAL = DetId(0);
527  trkD.detIdEHCAL = DetId(0);
528  if (debug)
529  edm::LogVerbatim("IsoTrack") << "Propagate track " << thisTrk << " charge " << trk.charge << " position "
530  << trk.position << " p " << trk.momentum << " Flag " << trkD.ok;
531  if (trkD.ok) {
532  spr::propagatedTrack info = spr::propagateCalo(
533  trk.position, trk.momentum, trk.charge, bField, spr::zFrontEE, spr::rFrontEB, spr::etaBEEcal, debug);
534  GlobalPoint point(info.point.x(), info.point.y(), info.point.z());
535  trkD.okECAL = info.ok;
536  trkD.pointECAL = point;
537  trkD.directionECAL = info.direction;
538  if (trkD.okECAL) {
539  if (std::abs(info.point.eta()) < spr::etaBEEcal) {
540  trkD.detIdECAL = barrelGeom->getClosestCell(point);
541  } else {
542  if (endcapGeom)
543  trkD.detIdECAL = endcapGeom->getClosestCell(point);
544  else
545  trkD.okECAL = false;
546  }
547  trkD.detIdEHCAL = gHB->getClosestCell(point);
548  }
549 
550  info = spr::propagateCalo(
551  trk.position, trk.momentum, trk.charge, bField, spr::zFrontHE, spr::rFrontHB, spr::etaBEHcal, debug);
552  point = GlobalPoint(info.point.x(), info.point.y(), info.point.z());
553  trkD.okHCAL = info.ok;
554  trkD.pointHCAL = point;
555  trkD.directionHCAL = info.direction;
556  if (trkD.okHCAL) {
557  trkD.detIdHCAL = gHB->getClosestCell(point);
558  }
559  }
560  if (debug) {
561  edm::LogVerbatim("IsoTrack") << "propagateCALO:: for track [" << thisTrk << "] Flag: " << trkD.ok << " ECAL ("
562  << trkD.okECAL << ") HCAL (" << trkD.okHCAL << ")";
563  std::ostringstream st1;
564  if (trkD.okECAL) {
565  st1 << "ECAL point " << trkD.pointECAL << " direction " << trkD.directionECAL << " ";
566  if (trkD.detIdECAL.subdetId() == EcalBarrel) {
567  st1 << (EBDetId)(trkD.detIdECAL);
568  } else {
569  st1 << (EEDetId)(trkD.detIdECAL);
570  }
571  }
572  if (trkD.okHCAL) {
573  st1 << " HCAL point " << trkD.pointHCAL << " direction " << trkD.directionHCAL << " "
574  << (HcalDetId)(trkD.detIdHCAL);
575  }
576  if (trkD.okECAL)
577  st1 << " Or " << (HcalDetId)(trkD.detIdEHCAL);
578  edm::LogVerbatim("IsoTrack") << st1.str();
579  }
580  return trkD;
581  }
582 
583  spr::propagatedTrackDirection propagateHCALBack(unsigned int thisTrk,
584  edm::Handle<edm::SimTrackContainer>& SimTk,
585  edm::Handle<edm::SimVertexContainer>& SimVtx,
586  const CaloGeometry* geo,
587  const MagneticField* bField,
588  bool debug) {
589  const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
590  spr::trackAtOrigin trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
591  spr::propagatedTrackDirection trkD;
592  trkD.ok = trk.ok;
593  trkD.detIdECAL = DetId(0);
594  trkD.detIdHCAL = DetId(0);
595  trkD.detIdEHCAL = DetId(0);
596  if (debug)
597  edm::LogVerbatim("IsoTrack") << "Propagate track " << thisTrk << " charge " << trk.charge << " position "
598  << trk.position << " p " << trk.momentum << " Flag " << trkD.ok;
599  if (trkD.ok) {
600  spr::propagatedTrack info = spr::propagateCalo(
601  trk.position, trk.momentum, trk.charge, bField, spr::zBackHE, spr::rBackHB, spr::etaBEHcal, debug);
602  const GlobalPoint point = GlobalPoint(info.point.x(), info.point.y(), info.point.z());
603  trkD.okHCAL = info.ok;
604  trkD.pointHCAL = point;
605  trkD.directionHCAL = info.direction;
606  if (trkD.okHCAL) {
607  trkD.detIdHCAL = gHB->getClosestCell(point);
608  }
609  }
610  if (debug) {
611  std::ostringstream st1;
612  if (trkD.okHCAL) {
613  st1 << " HCAL point " << trkD.pointHCAL << " direction " << trkD.directionHCAL << " "
614  << (HcalDetId)(trkD.detIdHCAL);
615  }
616  edm::LogVerbatim("IsoTrack") << "propagateCALO:: for track [" << thisTrk << "] Flag: " << trkD.ok << " ECAL ("
617  << trkD.okECAL << ") HCAL (" << trkD.okHCAL << ")" << st1.str();
618  }
619  return trkD;
620  }
621 
622  std::pair<bool, HcalDetId> propagateHCALBack(const reco::Track* track,
623  const CaloGeometry* geo,
624  const MagneticField* bField,
625  bool debug) {
626  const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
627  const GlobalPoint vertex(track->vx(), track->vy(), track->vz());
628  const GlobalVector momentum(track->px(), track->py(), track->pz());
629  int charge(track->charge());
630  spr::propagatedTrack info =
631  spr::propagateCalo(vertex, momentum, charge, bField, spr::zBackHE, spr::rBackHB, spr::etaBEHcal, debug);
632  if (info.ok) {
633  const GlobalPoint point = GlobalPoint(info.point.x(), info.point.y(), info.point.z());
634  return std::pair<bool, HcalDetId>(true, HcalDetId(gHB->getClosestCell(point)));
635  } else {
636  return std::pair<bool, HcalDetId>(false, HcalDetId());
637  }
638  }
639 
640  propagatedTrack propagateTrackToECAL(const reco::Track* track, const MagneticField* bfield, bool debug) {
641  const GlobalPoint vertex(track->vx(), track->vy(), track->vz());
642  const GlobalVector momentum(track->px(), track->py(), track->pz());
643  int charge(track->charge());
644  return spr::propagateCalo(vertex, momentum, charge, bfield, spr::zFrontEE, spr::rFrontEB, spr::etaBEEcal, debug);
645  }
646 
647  propagatedTrack propagateTrackToECAL(unsigned int thisTrk,
648  edm::Handle<edm::SimTrackContainer>& SimTk,
649  edm::Handle<edm::SimVertexContainer>& SimVtx,
650  const MagneticField* bfield,
651  bool debug) {
652  spr::trackAtOrigin trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
653  spr::propagatedTrack ptrk;
654  if (trk.ok)
655  ptrk = spr::propagateCalo(
656  trk.position, trk.momentum, trk.charge, bfield, spr::zFrontEE, spr::rFrontEB, spr::etaBEEcal, debug);
657  return ptrk;
658  }
659 
660  std::pair<math::XYZPoint, bool> propagateECAL(const reco::Track* track, const MagneticField* bfield, bool debug) {
661  const GlobalPoint vertex(track->vx(), track->vy(), track->vz());
662  const GlobalVector momentum(track->px(), track->py(), track->pz());
663  int charge(track->charge());
664  return spr::propagateECAL(vertex, momentum, charge, bfield, debug);
665  }
666 
667  std::pair<DetId, bool> propagateIdECAL(const HcalDetId& id,
668  const CaloGeometry* geo,
669  const MagneticField* bField,
670  bool debug) {
671  const HcalGeometry* gHB = static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel));
672  const GlobalPoint vertex(0, 0, 0);
673  const GlobalPoint hit(gHB->getPosition(id));
674  const GlobalVector momentum = GlobalVector(hit.x(), hit.y(), hit.z());
675  std::pair<math::XYZPoint, bool> info = propagateECAL(vertex, momentum, 0, bField, debug);
676  DetId eId(0);
677  if (info.second) {
678  const GlobalPoint point(info.first.x(), info.first.y(), info.first.z());
679  if (std::abs(point.eta()) < spr::etaBEEcal) {
680  const CaloSubdetectorGeometry* barrelGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
681  eId = barrelGeom->getClosestCell(point);
682  } else {
683  const CaloSubdetectorGeometry* endcapGeom = geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
684  if (endcapGeom)
685  eId = endcapGeom->getClosestCell(point);
686  else
687  info.second = false;
688  }
689  }
690  return std::pair<DetId, bool>(eId, info.second);
691  }
692 
693  std::pair<math::XYZPoint, bool> propagateECAL(
694  const GlobalPoint& vertex, const GlobalVector& momentum, int charge, const MagneticField* bfield, bool debug) {
695  spr::propagatedTrack track =
696  spr::propagateCalo(vertex, momentum, charge, bfield, spr::zFrontEE, spr::rFrontEB, spr::etaBEEcal, debug);
697  return std::pair<math::XYZPoint, bool>(track.point, track.ok);
698  }
699 
700  spr::propagatedTrack propagateTrackToHCAL(const reco::Track* track, const MagneticField* bfield, bool debug) {
701  const GlobalPoint vertex(track->vx(), track->vy(), track->vz());
702  const GlobalVector momentum(track->px(), track->py(), track->pz());
703  int charge(track->charge());
704  return spr::propagateCalo(vertex, momentum, charge, bfield, spr::zFrontHE, spr::rFrontHB, spr::etaBEHcal, debug);
705  }
706 
707  spr::propagatedTrack propagateTrackToHCAL(unsigned int thisTrk,
708  edm::Handle<edm::SimTrackContainer>& SimTk,
709  edm::Handle<edm::SimVertexContainer>& SimVtx,
710  const MagneticField* bfield,
711  bool debug) {
712  spr::trackAtOrigin trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
713  spr::propagatedTrack ptrk;
714  if (trk.ok)
715  ptrk = spr::propagateCalo(
716  trk.position, trk.momentum, trk.charge, bfield, spr::zFrontHE, spr::rFrontHB, spr::etaBEHcal, debug);
717  return ptrk;
718  }
719 
720  std::pair<math::XYZPoint, bool> propagateHCAL(const reco::Track* track, const MagneticField* bfield, bool debug) {
721  const GlobalPoint vertex(track->vx(), track->vy(), track->vz());
722  const GlobalVector momentum(track->px(), track->py(), track->pz());
723  int charge(track->charge());
724  return spr::propagateHCAL(vertex, momentum, charge, bfield, debug);
725  }
726 
727  std::pair<math::XYZPoint, bool> propagateHCAL(
728  const GlobalPoint& vertex, const GlobalVector& momentum, int charge, const MagneticField* bfield, bool debug) {
729  spr::propagatedTrack track =
730  spr::propagateCalo(vertex, momentum, charge, bfield, spr::zFrontHE, spr::rFrontHB, spr::etaBEHcal, debug);
731  return std::pair<math::XYZPoint, bool>(track.point, track.ok);
732  }
733 
734  std::pair<math::XYZPoint, bool> propagateTracker(const reco::Track* track, const MagneticField* bfield, bool debug) {
735  const GlobalPoint vertex(track->vx(), track->vy(), track->vz());
736  const GlobalVector momentum(track->px(), track->py(), track->pz());
737  int charge(track->charge());
738  spr::propagatedTrack track1 =
739  spr::propagateCalo(vertex, momentum, charge, bfield, spr::zBackTE, spr::rBackTB, spr::etaBETrak, debug);
740  return std::pair<math::XYZPoint, bool>(track1.point, track1.ok);
741  }
742 
743  std::pair<math::XYZPoint, double> propagateTrackerEnd(const reco::Track* track,
744  const MagneticField* bField,
745  bool debug) {
746  const GlobalPoint vertex(track->vx(), track->vy(), track->vz());
747  const GlobalVector momentum(track->px(), track->py(), track->pz());
748  int charge(track->charge());
749  float radius = track->outerPosition().Rho();
750  float zdist = track->outerPosition().Z();
751  if (debug)
752  edm::LogVerbatim("IsoTrack") << "propagateTrackerEnd:: Vertex " << vertex << " Momentum " << momentum
753  << " Charge " << charge << " Radius " << radius << " Z " << zdist;
754  FreeTrajectoryState fts(vertex, momentum, charge, bField);
755  Plane::PlanePointer endcap = Plane::build(Plane::PositionType(0, 0, zdist), Plane::RotationType());
756  Cylinder::CylinderPointer barrel =
757  Cylinder::build(Cylinder::PositionType(0, 0, 0), Cylinder::RotationType(), radius);
758 
759  AnalyticalPropagator myAP(bField, alongMomentum, 2 * M_PI);
760 
761  TrajectoryStateOnSurface tsose = myAP.propagate(fts, *endcap);
762  TrajectoryStateOnSurface tsosb = myAP.propagate(fts, *barrel);
763 
764  math::XYZPoint point(-999., -999., -999.);
765  bool ok = false;
766  GlobalVector direction(0, 0, 1);
767  if (tsosb.isValid() && std::abs(zdist) < spr::zFrontTE) {
768  point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
769  direction = tsosb.globalDirection();
770  ok = true;
771  } else if (tsose.isValid()) {
772  point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
773  direction = tsose.globalDirection();
774  ok = true;
775  }
776 
777  double length = -1;
778  if (ok) {
779  math::XYZPoint vDiff(point.x() - vertex.x(), point.y() - vertex.y(), point.z() - vertex.z());
780  double dphi = direction.phi() - momentum.phi();
781  double rdist = std::sqrt(vDiff.x() * vDiff.x() + vDiff.y() * vDiff.y());
782  double rat = 0.5 * dphi / std::sin(0.5 * dphi);
783  double dZ = vDiff.z();
784  double dS = rdist * rat; //dZ*momentum.z()/momentum.perp();
785  length = std::sqrt(dS * dS + dZ * dZ);
786  if (debug)
787  edm::LogVerbatim("IsoTrack") << "propagateTracker:: Barrel " << tsosb.isValid() << " Endcap " << tsose.isValid()
788  << " OverAll " << ok << " Point " << point << " RDist " << rdist << " dS " << dS
789  << " dS/pt " << rdist * rat / momentum.perp() << " zdist " << dZ << " dz/pz "
790  << dZ / momentum.z() << " Length " << length;
791  }
792 
793  return std::pair<math::XYZPoint, double>(point, length);
794  }
795 
796  spr::propagatedTrack propagateCalo(const GlobalPoint& tpVertex,
797  const GlobalVector& tpMomentum,
798  int tpCharge,
799  const MagneticField* bField,
800  float zdist,
801  float radius,
802  float corner,
803  bool debug) {
804  spr::propagatedTrack track;
805  if (debug)
806  edm::LogVerbatim("IsoTrack") << "propagateCalo:: Vertex " << tpVertex << " Momentum " << tpMomentum << " Charge "
807  << tpCharge << " Radius " << radius << " Z " << zdist << " Corner " << corner;
808  FreeTrajectoryState fts(tpVertex, tpMomentum, tpCharge, bField);
809 
810  Plane::PlanePointer lendcap = Plane::build(Plane::PositionType(0, 0, -zdist), Plane::RotationType());
811  Plane::PlanePointer rendcap = Plane::build(Plane::PositionType(0, 0, zdist), Plane::RotationType());
812 
813  Cylinder::CylinderPointer barrel =
814  Cylinder::build(Cylinder::PositionType(0, 0, 0), Cylinder::RotationType(), radius);
815 
816  AnalyticalPropagator myAP(bField, alongMomentum, 2 * M_PI);
817 
818  TrajectoryStateOnSurface tsose;
819  if (tpMomentum.eta() < 0) {
820  tsose = myAP.propagate(fts, *lendcap);
821  } else {
822  tsose = myAP.propagate(fts, *rendcap);
823  }
824 
825  TrajectoryStateOnSurface tsosb = myAP.propagate(fts, *barrel);
826 
827  track.ok = true;
828  if (tsose.isValid() && tsosb.isValid()) {
829  float absEta = std::abs(tsosb.globalPosition().eta());
830  if (absEta < corner) {
831  track.point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
832  track.direction = tsosb.globalDirection();
833  } else {
834  track.point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
835  track.direction = tsose.globalDirection();
836  }
837  } else if (tsose.isValid()) {
838  track.point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
839  track.direction = tsose.globalDirection();
840  } else if (tsosb.isValid()) {
841  track.point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
842  track.direction = tsosb.globalDirection();
843  } else {
844  track.point.SetXYZ(-999., -999., -999.);
845  track.direction = GlobalVector(0, 0, 1);
846  track.ok = false;
847  }
848  if (debug) {
849  edm::LogVerbatim("IsoTrack") << "propagateCalo:: Barrel " << tsosb.isValid() << " Endcap " << tsose.isValid()
850  << " OverAll " << track.ok << " Point " << track.point << " Direction "
851  << track.direction;
852  if (track.ok) {
853  math::XYZPoint vDiff(
854  track.point.x() - tpVertex.x(), track.point.y() - tpVertex.y(), track.point.z() - tpVertex.z());
855  double dphi = track.direction.phi() - tpMomentum.phi();
856  double rdist = std::sqrt(vDiff.x() * vDiff.x() + vDiff.y() * vDiff.y());
857  double pt = tpMomentum.perp();
858  double rat = 0.5 * dphi / std::sin(0.5 * dphi);
859  edm::LogVerbatim("IsoTrack") << "RDist " << rdist << " pt " << pt << " r/pt " << rdist * rat / pt << " zdist "
860  << vDiff.z() << " pz " << tpMomentum.z() << " z/pz " << vDiff.z() / tpMomentum.z()
861  << std::endl;
862  }
863  }
864  return track;
865  }
866 
867  spr::trackAtOrigin simTrackAtOrigin(unsigned int thisTrk,
868  edm::Handle<edm::SimTrackContainer>& SimTk,
869  edm::Handle<edm::SimVertexContainer>& SimVtx,
870  bool debug) {
871  spr::trackAtOrigin trk;
872 
873  edm::SimTrackContainer::const_iterator itr = SimTk->end();
874  for (edm::SimTrackContainer::const_iterator simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++) {
875  if (simTrkItr->trackId() == thisTrk) {
876  if (debug)
877  edm::LogVerbatim("IsoTrack") << "matched trackId (maximum occurance) " << thisTrk << " type "
878  << simTrkItr->type();
879  itr = simTrkItr;
880  break;
881  }
882  }
883 
884  if (itr != SimTk->end()) {
885  int vertIndex = itr->vertIndex();
886  if (vertIndex != -1 && vertIndex < (int)SimVtx->size()) {
887  edm::SimVertexContainer::const_iterator simVtxItr = SimVtx->begin();
888  for (int iv = 0; iv < vertIndex; iv++)
889  simVtxItr++;
890  const math::XYZTLorentzVectorD pos = simVtxItr->position();
891  const math::XYZTLorentzVectorD mom = itr->momentum();
892  trk.ok = true;
893  trk.charge = (int)(itr->charge());
894  trk.position = GlobalPoint(pos.x(), pos.y(), pos.z());
895  trk.momentum = GlobalVector(mom.x(), mom.y(), mom.z());
896  }
897  }
898  if (debug)
899  edm::LogVerbatim("IsoTrack") << "Track flag " << trk.ok << " Position " << trk.position << " Momentum "
900  << trk.momentum;
901  return trk;
902  }
903 
904  bool propagateHCAL(const reco::Track* track,
905  const CaloGeometry* geo,
906  const MagneticField* bField,
907  bool typeRZ,
908  const std::pair<double, double> rz,
909  bool debug) {
910  const GlobalPoint vertex(track->vx(), track->vy(), track->vz());
911  const GlobalVector momentum(track->px(), track->py(), track->pz());
912  int charge(track->charge());
913  if (debug)
914  edm::LogVerbatim("IsoTrack") << "Propagate track with charge " << charge << " position " << vertex << " p "
915  << momentum;
916  std::pair<HcalDetId, HcalDetId> ids = propagateHCAL(geo, bField, vertex, momentum, charge, typeRZ, rz, debug);
917  bool ok = ((ids.first != HcalDetId()) && (ids.first.ieta() == ids.second.ieta()) &&
918  (ids.first.iphi() == ids.second.iphi()));
919  return ok;
920  }
921 
922  bool propagateHCAL(unsigned int thisTrk,
923  edm::Handle<edm::SimTrackContainer>& SimTk,
924  edm::Handle<edm::SimVertexContainer>& SimVtx,
925  const CaloGeometry* geo,
926  const MagneticField* bField,
927  bool typeRZ,
928  const std::pair<double, double> rz,
929  bool debug) {
930  spr::trackAtOrigin trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
931  if (debug)
932  edm::LogVerbatim("IsoTrack") << "Propagate track " << thisTrk << " charge " << trk.charge << " position "
933  << trk.position << " p " << trk.momentum;
934  std::pair<HcalDetId, HcalDetId> ids =
935  propagateHCAL(geo, bField, trk.position, trk.momentum, trk.charge, typeRZ, rz, debug);
936  bool ok = ((ids.first != HcalDetId()) && (ids.first.ieta() == ids.second.ieta()) &&
937  (ids.first.iphi() == ids.second.iphi()));
938  return ok;
939  }
940 
941  std::pair<HcalDetId, HcalDetId> propagateHCAL(const CaloGeometry* geo,
942  const MagneticField* bField,
943  const GlobalPoint& vertex,
944  const GlobalVector& momentum,
945  int charge,
946  bool typeRZ,
947  const std::pair<double, double> rz,
948  bool debug) {
949  if (debug)
950  edm::LogVerbatim("IsoTrack") << "propagateCalo:: Vertex " << vertex << " Momentum " << momentum << " Charge "
951  << charge << " R/Z " << rz.first << " : " << rz.second << " Type " << typeRZ;
952  const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
953  FreeTrajectoryState fts(vertex, momentum, charge, bField);
954  AnalyticalPropagator myAP(bField, alongMomentum, 2 * M_PI);
955 
956  HcalDetId id1, id2;
957  for (int k = 0; k < 2; ++k) {
958  TrajectoryStateOnSurface tsos;
959  double rzv = (k == 0) ? rz.first : rz.second;
960  if (typeRZ) {
961  Cylinder::CylinderPointer barrel =
962  Cylinder::build(Cylinder::PositionType(0, 0, 0), Cylinder::RotationType(), rzv);
963  tsos = myAP.propagate(fts, *barrel);
964  } else {
965  Plane::PlanePointer endcap = Plane::build(Plane::PositionType(0, 0, rzv), Plane::RotationType());
966  tsos = myAP.propagate(fts, *endcap);
967  }
968 
969  if (tsos.isValid()) {
970  GlobalPoint point = tsos.globalPosition();
971  if (k == 0)
972  id1 = gHB->getClosestCell(point);
973  else
974  id2 = gHB->getClosestCell(point);
975  if (debug) {
976  std::ostringstream st1;
977  if (k == 0)
978  st1 << id1;
979  else
980  st1 << id2;
981  edm::LogVerbatim("IsoTrack") << "Iteration " << k << " Point " << point << " ID " << st1.str();
982  }
983  }
984  }
985  if (debug)
986  edm::LogVerbatim("IsoTrack") << "propagateCalo:: Front " << id1 << " Back " << id2;
987  return std::pair<HcalDetId, HcalDetId>(id1, id2);
988  }
989 } // namespace spr
reco::TrackBase::p
double p() const
momentum vector magnitude
Definition: TrackBase.h:631
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
CaloGeometry::getSubdetectorGeometry
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
Definition: CaloGeometry.cc:34
spr::etaBEEcal
static const double etaBEEcal
Definition: CaloConstants.h:12
spr::propagateCALO
std::vector< spr::propagatedTrackID > propagateCALO(edm::Handle< reco::TrackCollection > &trkCollection, const CaloGeometry *geo, const MagneticField *bField, const std::string &theTrackQuality, bool debug=false)
Definition: CaloPropagateTrack.cc:130
info
static const TGPicture * info(bool iBackgroundIsBlack)
Definition: FWCollectionSummaryWidget.cc:153
mps_fire.i
i
Definition: mps_fire.py:428
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:69
DiDispStaMuonMonitor_cfi.pt
tuple pt
Definition: DiDispStaMuonMonitor_cfi.py:39
gpuVertexFinder::iv
int32_t *__restrict__ iv
Definition: gpuClusterTracksDBSCAN.h:42
math::XYZTLorentzVectorD
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:14
spr::zBackTE
static const double zBackTE
Definition: CaloConstants.h:21
Reference_intrackfit_cff.barrel
barrel
Definition: Reference_intrackfit_cff.py:37
ParticleDataTable
HepPDT::ParticleDataTable ParticleDataTable
Definition: ParticleDataTable.h:8
spr::etaBEHcal
static const double etaBEHcal
Definition: CaloConstants.h:15
Reference_intrackfit_cff.endcap
endcap
Definition: Reference_intrackfit_cff.py:82
spr::propagateIdECAL
std::pair< DetId, bool > propagateIdECAL(const HcalDetId &id, const CaloGeometry *geo, const MagneticField *, bool debug=false)
Definition: CaloPropagateTrack.cc:667
spr::propagatedGenParticleID::trkItr
reco::GenParticleCollection::const_iterator trkItr
Definition: CaloPropagateTrack.h:87
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
spr::zFrontHE
static const double zFrontHE
Definition: CaloConstants.h:13
align::ID
uint32_t ID
Definition: Definitions.h:24
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
GlobalPoint
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
PV3DBase::y
T y() const
Definition: PV3DBase.h:60
X
#define X(str)
Definition: MuonsGrabber.cc:38
TrajectoryStateOnSurface::globalPosition
GlobalPoint globalPosition() const
Definition: TrajectoryStateOnSurface.h:65
reco::TrackBase::phi
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:649
spr::simTrackAtOrigin
spr::trackAtOrigin simTrackAtOrigin(unsigned int thisTrk, edm::Handle< edm::SimTrackContainer > &SimTk, edm::Handle< edm::SimVertexContainer > &SimVtx, bool debug=false)
Definition: CaloPropagateTrack.cc:867
spr::zBackHE
static const double zBackHE
Definition: CaloConstants.h:16
spr::propagateCalo
spr::propagatedTrack propagateCalo(const GlobalPoint &vertex, const GlobalVector &momentum, int charge, const MagneticField *, float zdist, float radius, float corner, bool debug=false)
Definition: CaloPropagateTrack.cc:796
HLT_FULL_cff.track
tuple track
Definition: HLT_FULL_cff.py:11953
reco::TrackBase::px
double px() const
x coordinate of momentum vector
Definition: TrackBase.h:640
genParticleCandidates2GenParticles_cfi.genParticles
tuple genParticles
Definition: genParticleCandidates2GenParticles_cfi.py:3
reco::Track::outerPosition
const math::XYZPoint & outerPosition() const
position of the outermost hit
Definition: Track.h:62
spr::rFrontHB
static const double rFrontHB
Definition: CaloConstants.h:14
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
spr::propagatedGenTrackID::trkItr
HepMC::GenEvent::particle_const_iterator trkItr
Definition: CaloPropagateTrack.h:74
reco::Track::innerPosition
const math::XYZPoint & innerPosition() const
position of the innermost hit
Definition: Track.h:56
spr::rBackHB
static const double rBackHB
Definition: CaloConstants.h:17
spr::propagatedTrackID::trkItr
reco::TrackCollection::const_iterator trkItr
Definition: CaloPropagateTrack.h:48
spr::zFrontEE
static const double zFrontEE
Definition: CaloConstants.h:9
reco::TrackBase::eta
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:652
spr::rBackTB
static const double rBackTB
Definition: CaloConstants.h:22
spr::propagateHCALBack
spr::propagatedTrackDirection propagateHCALBack(unsigned int thisTrk, edm::Handle< edm::SimTrackContainer > &SimTk, edm::Handle< edm::SimVertexContainer > &SimVtx, const CaloGeometry *geo, const MagneticField *bField, bool debug=false)
Definition: CaloPropagateTrack.cc:583
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
Plane::build
static PlanePointer build(Args &&...args)
Definition: Plane.h:33
PV3DBase::z
T z() const
Definition: PV3DBase.h:61
spr::rFrontEB
static const double rFrontEB
Definition: CaloConstants.h:10
if
if(conf_.getParameter< bool >("UseStripCablingDB"))
Definition: SiStripLAProfileBooker.cc:72
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
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
ecalTB2006H4_GenSimDigiReco_cfg.bField
tuple bField
Definition: ecalTB2006H4_GenSimDigiReco_cfg.py:236
MCTruth.genEvent
tuple genEvent
Definition: MCTruth.py:34
spr::propagateTracker
std::pair< math::XYZPoint, bool > propagateTracker(const reco::Track *, const MagneticField *, bool debug=false)
Definition: CaloPropagateTrack.cc:734
spr::propagateECAL
std::pair< math::XYZPoint, bool > propagateECAL(const reco::Track *, const MagneticField *, bool debug=false)
Definition: CaloPropagateTrack.cc:660
spr::propagateTrackToECAL
spr::propagatedTrack propagateTrackToECAL(const reco::Track *, const MagneticField *, bool debug=false)
Definition: CaloPropagateTrack.cc:640
HcalGeometry::getPosition
GlobalPoint getPosition(const DetId &id) const
Definition: HcalGeometry.cc:179
M_PI
#define M_PI
Definition: BXVectorInputProducer.cc:49
reco::TrackBase::pz
double pz() const
z coordinate of momentum vector
Definition: TrackBase.h:646
CaloSubdetectorGeometry::getClosestCell
virtual DetId getClosestCell(const GlobalPoint &r) const
Definition: CaloSubdetectorGeometry.cc:44
reco::TrackBase::vz
double vz() const
z coordinate of the reference point on track
Definition: TrackBase.h:661
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::Track::outerMomentum
const math::XYZVector & outerMomentum() const
momentum vector at the outermost hit position
Definition: Track.h:65
Cylinder::build
static CylinderPointer build(const PositionType &pos, const RotationType &rot, Scalar radius, Bounds *bounds=nullptr)
Definition: Cylinder.h:45
spr::zFrontTE
static const double zFrontTE
Definition: CaloConstants.h:20
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
spr::propagateCosmicCALO
std::vector< spr::propagatedTrackID > propagateCosmicCALO(edm::Handle< reco::TrackCollection > &trkCollection, const CaloGeometry *geo, const MagneticField *bField, const std::string &theTrackQuality, bool debug=false)
Definition: CaloPropagateTrack.cc:20
Propagator::propagate
TrajectoryStateOnSurface propagate(STA const &state, SUR const &surface) const
Definition: Propagator.h:50
reco::TrackBase::quality
bool quality(const TrackQuality) const
Track quality.
Definition: TrackBase.h:552
PV3DBase::eta
T eta() const
Definition: PV3DBase.h:73
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
reco::TrackBase::vy
double vy() const
y coordinate of the reference point on track
Definition: TrackBase.h:658
spr::propagateTrackerEnd
std::pair< math::XYZPoint, double > propagateTrackerEnd(const reco::Track *, const MagneticField *, bool debug=false)
Definition: CaloPropagateTrack.cc:743
reco::Track::innerMomentum
const math::XYZVector & innerMomentum() const
momentum vector at the innermost hit position
Definition: Track.h:59
spr::propagatedTrackDirection::trkItr
reco::TrackCollection::const_iterator trkItr
Definition: CaloPropagateTrack.h:61
reco::TrackBase::charge
int charge() const
track electric charge
Definition: TrackBase.h:596
Puppi_cff.etaMax
tuple etaMax
Definition: Puppi_cff.py:46
PV3DBase::x
T x() const
Definition: PV3DBase.h:59
spr::propagateTrackToHCAL
spr::propagatedTrack propagateTrackToHCAL(const reco::Track *, const MagneticField *, bool debug=false)
Definition: CaloPropagateTrack.cc:700
CosmicsPD_Skims.radius
tuple radius
Definition: CosmicsPD_Skims.py:135
spr::etaBETrak
static const double etaBETrak
Definition: CaloConstants.h:23
point
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
reco::TrackBase::py
double py() const
y coordinate of momentum vector
Definition: TrackBase.h:643
reco::TrackBase::vx
double vx() const
x coordinate of the reference point on track
Definition: TrackBase.h:655
GlobalVector
Global3DVector GlobalVector
Definition: GlobalVector.h:10
TrajectoryStateOnSurface::globalDirection
GlobalVector globalDirection() const
Definition: TrajectoryStateOnSurface.h:67
Generated for CMSSW Reference Manual by doxygen 1.8.5