CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros 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 "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
10 #include "Geometry/EcalAlgo/interface/EcalEndcapGeometry.h"
11 
12 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
13 
14 #include <iostream>
15 
16 namespace spr{
17 
18  std::vector<spr::propagatedTrackID> propagateCALO(edm::Handle<reco::TrackCollection>& trkCollection, const CaloGeometry* geo, const MagneticField* bField, std::string & theTrackQuality, bool debug) {
19 
20  std::vector<spr::propagatedTrackID> vdets;
21  spr::propagateCALO(trkCollection,geo,bField,theTrackQuality, vdets, debug);
22  return vdets;
23  }
24 
25  void propagateCALO(edm::Handle<reco::TrackCollection>& trkCollection, const CaloGeometry* geo, const MagneticField* bField, std::string & theTrackQuality, std::vector<spr::propagatedTrackID>& vdets, bool debug) {
26 
27  const EcalBarrelGeometry *barrelGeom = (dynamic_cast< const EcalBarrelGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel)));
28  const EcalEndcapGeometry *endcapGeom = (dynamic_cast< const EcalEndcapGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap)));
29  const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);
30  reco::TrackBase::TrackQuality trackQuality_=reco::TrackBase::qualityByName(theTrackQuality);
31 
32  unsigned int indx;
33  reco::TrackCollection::const_iterator trkItr;
34  for (trkItr = trkCollection->begin(),indx=0; trkItr != trkCollection->end(); ++trkItr,indx++) {
35  const reco::Track* pTrack = &(*trkItr);
36  spr::propagatedTrackID vdet;
37  vdet.trkItr = trkItr;
38  vdet.ok = (pTrack->quality(trackQuality_));
39  vdet.detIdECAL = DetId(0);
40  vdet.detIdHCAL = DetId(0);
41  vdet.detIdEHCAL= DetId(0);
42  if (debug) std::cout << "Propagate track " << indx << " p " << trkItr->p() << " eta " << trkItr->eta() << " phi " << trkItr->phi() << " Flag " << vdet.ok << std::endl;
43 
44  std::pair<math::XYZPoint,bool> info = spr::propagateECAL (pTrack, bField, debug);
45  vdet.okECAL = info.second;
46  if (vdet.okECAL) {
47  const GlobalPoint point(info.first.x(),info.first.y(),info.first.z());
48  vdet.etaECAL = point.eta();
49  vdet.phiECAL = point.phi();
50  if (std::abs(point.eta())<1.479) {
51  vdet.detIdECAL = barrelGeom->getClosestCell(point);
52  } else {
53  vdet.detIdECAL = endcapGeom->getClosestCell(point);
54  }
55  vdet.detIdEHCAL = gHB->getClosestCell(point);
56  }
57  info = spr::propagateHCAL (pTrack, bField, debug);
58  vdet.okHCAL = info.second;
59  if (vdet.okHCAL) {
60  const GlobalPoint point(info.first.x(),info.first.y(),info.first.z());
61  vdet.etaHCAL = point.eta();
62  vdet.phiHCAL = point.phi();
63  vdet.detIdHCAL = gHB->getClosestCell(point);
64  }
65 
66  vdets.push_back(vdet);
67  }
68 
69  if (debug) {
70  std::cout << "propagateCALO:: for " << vdets.size() << " tracks" << std::endl;
71  for (unsigned int i=0; i<vdets.size(); ++i) {
72  std::cout << "Track [" << i << "] Flag: " << vdets[i].ok << " ECAL (" << vdets[i].okECAL << ") ";
73  if (vdets[i].detIdECAL.subdetId() == EcalBarrel) {
74  std::cout << (EBDetId)(vdets[i].detIdECAL);
75  } else {
76  std::cout << (EEDetId)(vdets[i].detIdECAL);
77  }
78  std::cout << " HCAL (" << vdets[i].okHCAL << ") " << (HcalDetId)(vdets[i].detIdHCAL) << " Or " << (HcalDetId)(vdets[i].detIdEHCAL) << std::endl;
79  }
80  }
81  }
82 
83  void propagateCALO(edm::Handle<reco::TrackCollection>& trkCollection, const CaloGeometry* geo, const MagneticField* bField, std::string & theTrackQuality, std::vector<spr::propagatedTrackDirection>& trkDir, bool debug) {
84 
85  const EcalBarrelGeometry *barrelGeom = (dynamic_cast< const EcalBarrelGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel)));
86  const EcalEndcapGeometry *endcapGeom = (dynamic_cast< const EcalEndcapGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap)));
87  const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);
88  reco::TrackBase::TrackQuality trackQuality_=reco::TrackBase::qualityByName(theTrackQuality);
89 
90  unsigned int indx;
91  reco::TrackCollection::const_iterator trkItr;
92  for (trkItr = trkCollection->begin(),indx=0; trkItr != trkCollection->end(); ++trkItr,indx++) {
93  const reco::Track* pTrack = &(*trkItr);
94  spr::propagatedTrackDirection trkD;
95  trkD.trkItr = trkItr;
96  trkD.ok = (pTrack->quality(trackQuality_));
97  trkD.detIdECAL = DetId(0);
98  trkD.detIdHCAL = DetId(0);
99  trkD.detIdEHCAL= DetId(0);
100  if (debug) std::cout << "Propagate track " << indx << " p " << trkItr->p() << " eta " << trkItr->eta() << " phi " << trkItr->phi() << " Flag " << trkD.ok << std::endl;
101 
102  spr::propagatedTrack info = spr::propagateTrackToECAL (pTrack, bField, debug);
103  GlobalPoint point(info.point.x(),info.point.y(),info.point.z());
104  trkD.okECAL = info.ok;
105  trkD.pointECAL = point;
106  trkD.directionECAL = info.direction;
107  if (trkD.okECAL) {
108  if (std::abs(info.point.eta())<1.479) {
109  trkD.detIdECAL = barrelGeom->getClosestCell(point);
110  } else {
111  trkD.detIdECAL = endcapGeom->getClosestCell(point);
112  }
113  trkD.detIdEHCAL = gHB->getClosestCell(point);
114  }
115  info = spr::propagateTrackToHCAL (pTrack, bField, debug);
116  point = GlobalPoint(info.point.x(),info.point.y(),info.point.z());
117  trkD.okHCAL = info.ok;
118  trkD.pointHCAL = point;
119  trkD.directionHCAL = info.direction;
120  if (trkD.okHCAL) {
121  trkD.detIdHCAL = gHB->getClosestCell(point);
122  }
123  trkDir.push_back(trkD);
124  }
125 
126  if (debug) {
127  std::cout << "propagateCALO:: for " << trkDir.size() << " tracks" << std::endl;
128  for (unsigned int i=0; i<trkDir.size(); ++i) {
129  std::cout << "Track [" << i << "] Flag: " << trkDir[i].ok << " ECAL (" << trkDir[i].okECAL << ")";
130  if (trkDir[i].okECAL) {
131  std::cout << " point " << trkDir[i].pointECAL << " direction "
132  << trkDir[i].directionECAL << " ";
133  if (trkDir[i].detIdECAL.subdetId() == EcalBarrel) {
134  std::cout << (EBDetId)(trkDir[i].detIdECAL);
135  } else {
136  std::cout << (EEDetId)(trkDir[i].detIdECAL);
137  }
138  }
139  std::cout << " HCAL (" << trkDir[i].okHCAL << ")";
140  if (trkDir[i].okHCAL) {
141  std::cout << " point " << trkDir[i].pointHCAL << " direction "
142  << trkDir[i].directionHCAL << " "
143  << (HcalDetId)(trkDir[i].detIdHCAL);
144  }
145  std::cout << " Or " << (HcalDetId)(trkDir[i].detIdEHCAL) << std::endl;
146  }
147  }
148  }
149 
150  std::vector<spr::propagatedGenTrackID> propagateCALO(const HepMC::GenEvent * genEvent, edm::ESHandle<ParticleDataTable>& pdt, const CaloGeometry* geo, const MagneticField* bField, double etaMax, bool debug) {
151 
152  const EcalBarrelGeometry *barrelGeom = (dynamic_cast< const EcalBarrelGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel)));
153  const EcalEndcapGeometry *endcapGeom = (dynamic_cast< const EcalEndcapGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap)));
154  const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);
155 
156  std::vector<spr::propagatedGenTrackID> trkDir;
157  unsigned int indx;
158  HepMC::GenEvent::particle_const_iterator p;
159  for (p=genEvent->particles_begin(),indx=0; p != genEvent->particles_end(); ++p,++indx) {
160  spr::propagatedGenTrackID trkD;
161  trkD.trkItr = p;
162  trkD.detIdECAL = DetId(0);
163  trkD.detIdHCAL = DetId(0);
164  trkD.detIdEHCAL= DetId(0);
165  trkD.pdgId = ((*p)->pdg_id());
166  trkD.charge = ((pdt->particle(trkD.pdgId))->ID().threeCharge())/3;
167  GlobalVector momentum = GlobalVector((*p)->momentum().px(), (*p)->momentum().py(), (*p)->momentum().pz());
168  if (debug) std::cout << "Propagate track " << indx << " pdg " << trkD.pdgId << " charge " << trkD.charge << " p " << momentum << std::endl;
169 
170  // consider stable particles
171  if ( (*p)->status()==1 && std::abs((*p)->momentum().eta()) < etaMax ) {
172  GlobalPoint vertex = GlobalPoint(0.1*(*p)->production_vertex()->position().x(),
173  0.1*(*p)->production_vertex()->position().y(),
174  0.1*(*p)->production_vertex()->position().z());
175  trkD.ok = true;
176  spr::propagatedTrack info = spr::propagateCalo (vertex, momentum, trkD.charge, bField, 319.2, 129.4, 1.479, debug);
177  GlobalPoint point(info.point.x(),info.point.y(),info.point.z());
178  trkD.okECAL = info.ok;
179  trkD.pointECAL = point;
180  trkD.directionECAL = info.direction;
181  if (trkD.okECAL) {
182  if (std::abs(info.point.eta())<1.479) {
183  trkD.detIdECAL = barrelGeom->getClosestCell(point);
184  } else {
185  trkD.detIdECAL = endcapGeom->getClosestCell(point);
186  }
187  trkD.detIdEHCAL = gHB->getClosestCell(point);
188  }
189 
190  info = spr::propagateCalo (vertex, momentum, trkD.charge, bField, 402.7, 180.7, 1.392, debug);
191  point = GlobalPoint(info.point.x(),info.point.y(),info.point.z());
192  trkD.okHCAL = info.ok;
193  trkD.pointHCAL = point;
194  trkD.directionHCAL = info.direction;
195  if (trkD.okHCAL) {
196  trkD.detIdHCAL = gHB->getClosestCell(point);
197  }
198  }
199  trkDir.push_back(trkD);
200  }
201 
202  if (debug) {
203  std::cout << "propagateCALO:: for " << trkDir.size() << " tracks" << std::endl;
204  for (unsigned int i=0; i<trkDir.size(); ++i) {
205  if (trkDir[i].okECAL) std::cout << "Track [" << i << "] Flag: " << trkDir[i].ok << " ECAL (" << trkDir[i].okECAL << ")";
206  if (trkDir[i].okECAL) {
207  std::cout << " point " << trkDir[i].pointECAL << " direction "
208  << trkDir[i].directionECAL << " ";
209  if (trkDir[i].detIdECAL.subdetId() == EcalBarrel) {
210  std::cout << (EBDetId)(trkDir[i].detIdECAL);
211  } else {
212  std::cout << (EEDetId)(trkDir[i].detIdECAL);
213  }
214  }
215  if (trkDir[i].okECAL) std::cout << " HCAL (" << trkDir[i].okHCAL << ")";
216  if (trkDir[i].okHCAL) {
217  std::cout << " point " << trkDir[i].pointHCAL << " direction "
218  << trkDir[i].directionHCAL << " "
219  << (HcalDetId)(trkDir[i].detIdHCAL);
220  }
221  if (trkDir[i].okECAL) std::cout << " Or " << (HcalDetId)(trkDir[i].detIdEHCAL) << std::endl;
222  }
223  }
224  return trkDir;
225  }
226 
227  std::vector<spr::propagatedGenParticleID> propagateCALO(edm::Handle<reco::GenParticleCollection>& genParticles, edm::ESHandle<ParticleDataTable>& pdt, const CaloGeometry* geo, const MagneticField* bField, double etaMax, bool debug) {
228 
229  const EcalBarrelGeometry *barrelGeom = (dynamic_cast< const EcalBarrelGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel)));
230  const EcalEndcapGeometry *endcapGeom = (dynamic_cast< const EcalEndcapGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap)));
231  const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);
232 
233  std::vector<spr::propagatedGenParticleID> trkDir;
234  unsigned int indx;
235  reco::GenParticleCollection::const_iterator p;
236  for (p=genParticles->begin(),indx=0; p != genParticles->end(); ++p,++indx) {
237  spr::propagatedGenParticleID trkD;
238  trkD.trkItr = p;
239  trkD.detIdECAL = DetId(0);
240  trkD.detIdHCAL = DetId(0);
241  trkD.detIdEHCAL= DetId(0);
242  trkD.pdgId = (p->pdgId());
243  trkD.charge = p->charge();
244  GlobalVector momentum = GlobalVector(p->momentum().x(), p->momentum().y(), p->momentum().z());
245  if (debug) std::cout << "Propagate track " << indx << " pdg " << trkD.pdgId << " charge " << trkD.charge << " p " << momentum << std::endl;
246 
247  // consider stable particles
248  if ( p->status()==1 && std::abs(momentum.eta()) < etaMax ) {
249  GlobalPoint vertex = GlobalPoint(p->vertex().x(), p->vertex().y(), p->vertex().z());
250  trkD.ok = true;
251  spr::propagatedTrack info = spr::propagateCalo (vertex, momentum, trkD.charge, bField, 319.2, 129.4, 1.479, debug);
252  GlobalPoint point(info.point.x(),info.point.y(),info.point.z());
253  trkD.okECAL = info.ok;
254  trkD.pointECAL = point;
255  trkD.directionECAL = info.direction;
256  if (trkD.okECAL) {
257  if (std::abs(info.point.eta())<1.479) {
258  trkD.detIdECAL = barrelGeom->getClosestCell(point);
259  } else {
260  trkD.detIdECAL = endcapGeom->getClosestCell(point);
261  }
262  trkD.detIdEHCAL = gHB->getClosestCell(point);
263  }
264 
265  info = spr::propagateCalo (vertex, momentum, trkD.charge, bField, 402.7, 180.7, 1.392, debug);
266  point = GlobalPoint(info.point.x(),info.point.y(),info.point.z());
267  trkD.okHCAL = info.ok;
268  trkD.pointHCAL = point;
269  trkD.directionHCAL = info.direction;
270  if (trkD.okHCAL) {
271  trkD.detIdHCAL = gHB->getClosestCell(point);
272  }
273  }
274  trkDir.push_back(trkD);
275  }
276 
277  if (debug) {
278  std::cout << "propagateCALO:: for " << trkDir.size() << " tracks" << std::endl;
279  for (unsigned int i=0; i<trkDir.size(); ++i) {
280  if (trkDir[i].okECAL) std::cout << "Track [" << i << "] Flag: " << trkDir[i].ok << " ECAL (" << trkDir[i].okECAL << ")";
281  if (trkDir[i].okECAL) {
282  std::cout << " point " << trkDir[i].pointECAL << " direction "
283  << trkDir[i].directionECAL << " ";
284  if (trkDir[i].detIdECAL.subdetId() == EcalBarrel) {
285  std::cout << (EBDetId)(trkDir[i].detIdECAL);
286  } else {
287  std::cout << (EEDetId)(trkDir[i].detIdECAL);
288  }
289  }
290  if (trkDir[i].okECAL) std::cout << " HCAL (" << trkDir[i].okHCAL << ")";
291  if (trkDir[i].okHCAL) {
292  std::cout << " point " << trkDir[i].pointHCAL << " direction "
293  << trkDir[i].directionHCAL << " "
294  << (HcalDetId)(trkDir[i].detIdHCAL);
295  }
296  if (trkDir[i].okECAL) std::cout << " Or " << (HcalDetId)(trkDir[i].detIdEHCAL) << std::endl;
297  }
298  }
299  return trkDir;
300  }
301 
302  spr::propagatedTrackDirection propagateCALO(unsigned int thisTrk, edm::Handle<edm::SimTrackContainer>& SimTk, edm::Handle<edm::SimVertexContainer>& SimVtx, const CaloGeometry* geo, const MagneticField* bField, bool debug) {
303 
304  const EcalBarrelGeometry *barrelGeom = (dynamic_cast< const EcalBarrelGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel)));
305  const EcalEndcapGeometry *endcapGeom = (dynamic_cast< const EcalEndcapGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap)));
306  const CaloSubdetectorGeometry* gHB = geo->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);
307 
308  spr::trackAtOrigin trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
309  spr::propagatedTrackDirection trkD;
310  trkD.ok = trk.ok;
311  trkD.detIdECAL = DetId(0);
312  trkD.detIdHCAL = DetId(0);
313  trkD.detIdEHCAL= DetId(0);
314  if (debug) std::cout << "Propagate track " << thisTrk << " charge " << trk.charge << " position " << trk.position << " p " << trk.momentum << " Flag " << trkD.ok << std::endl;
315 
316  if (trkD.ok) {
317  spr::propagatedTrack info = spr::propagateCalo (trk.position, trk.momentum, trk.charge, bField, 319.2, 129.4, 1.479, debug);
318  GlobalPoint point(info.point.x(),info.point.y(),info.point.z());
319  trkD.okECAL = info.ok;
320  trkD.pointECAL = point;
321  trkD.directionECAL = info.direction;
322  if (trkD.okECAL) {
323  if (std::abs(info.point.eta())<1.479) {
324  trkD.detIdECAL = barrelGeom->getClosestCell(point);
325  } else {
326  trkD.detIdECAL = endcapGeom->getClosestCell(point);
327  }
328  trkD.detIdEHCAL = gHB->getClosestCell(point);
329  }
330 
331  info = spr::propagateCalo (trk.position, trk.momentum, trk.charge, bField, 402.7, 180.7, 1.392, debug);
332  point = GlobalPoint(info.point.x(),info.point.y(),info.point.z());
333  trkD.okHCAL = info.ok;
334  trkD.pointHCAL = point;
335  trkD.directionHCAL = info.direction;
336  if (trkD.okHCAL) {
337  trkD.detIdHCAL = gHB->getClosestCell(point);
338  }
339  }
340 
341  if (debug) {
342  std::cout << "propagateCALO:: for track [" << thisTrk << "] Flag: " << trkD.ok << " ECAL (" << trkD.okECAL << ") HCAL (" << trkD.okHCAL << ")" << std::endl;
343  if (trkD.okECAL) {
344  std::cout << "ECAL point " << trkD.pointECAL << " direction "
345  << trkD.directionECAL << " ";
346  if (trkD.detIdECAL.subdetId() == EcalBarrel) {
347  std::cout << (EBDetId)(trkD.detIdECAL);
348  } else {
349  std::cout << (EEDetId)(trkD.detIdECAL);
350  }
351  }
352  if (trkD.okHCAL) {
353  std::cout << " HCAL point " << trkD.pointHCAL << " direction "
354  << trkD.directionHCAL << " " << (HcalDetId)(trkD.detIdHCAL);
355  }
356  if (trkD.okECAL) std::cout << " Or " << (HcalDetId)(trkD.detIdEHCAL);
357  std::cout << std::endl;
358  }
359 
360  return trkD;
361  }
362 
363  propagatedTrack propagateTrackToECAL(const reco::Track *track, const MagneticField* bfield, bool debug) {
364  GlobalPoint vertex (track->vx(), track->vy(), track->vz());
365  GlobalVector momentum (track->px(), track->py(), track->pz());
366  int charge (track->charge());
367  return spr::propagateCalo (vertex, momentum, charge, bfield, 319.2, 129.4, 1.479, debug);
368  }
369 
370  propagatedTrack propagateTrackToECAL(unsigned int thisTrk, edm::Handle<edm::SimTrackContainer>& SimTk, edm::Handle<edm::SimVertexContainer>& SimVtx, const MagneticField* bfield, bool debug) {
371 
372  spr::trackAtOrigin trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
373  spr::propagatedTrack ptrk;
374  if (trk.ok)
375  ptrk = spr::propagateCalo (trk.position, trk.momentum, trk.charge, bfield, 319.2, 129.4, 1.479, debug);
376  return ptrk;
377  }
378 
379  std::pair<math::XYZPoint,bool> propagateECAL(const reco::Track *track, const MagneticField* bfield, bool debug) {
380  GlobalPoint vertex (track->vx(), track->vy(), track->vz());
381  GlobalVector momentum (track->px(), track->py(), track->pz());
382  int charge (track->charge());
383  return spr::propagateECAL (vertex, momentum, charge, bfield, debug);
384  }
385 
386  std::pair<math::XYZPoint,bool> propagateECAL(const GlobalPoint& vertex, const GlobalVector& momentum, int charge, const MagneticField* bfield, bool debug) {
387  spr::propagatedTrack track = spr::propagateCalo (vertex, momentum, charge, bfield, 319.2, 129.4, 1.479, debug);
388  return std::pair<math::XYZPoint,bool>(track.point,track.ok);
389  }
390 
391  spr::propagatedTrack propagateTrackToHCAL(const reco::Track *track, const MagneticField* bfield, bool debug) {
392  GlobalPoint vertex (track->vx(), track->vy(), track->vz());
393  GlobalVector momentum (track->px(), track->py(), track->pz());
394  int charge (track->charge());
395  return spr::propagateCalo (vertex, momentum, charge, bfield, 402.7, 180.7, 1.392, debug);
396  }
397 
398  spr::propagatedTrack propagateTrackToHCAL(unsigned int thisTrk, edm::Handle<edm::SimTrackContainer>& SimTk, edm::Handle<edm::SimVertexContainer>& SimVtx, const MagneticField* bfield, bool debug) {
399  spr::trackAtOrigin trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
400  spr::propagatedTrack ptrk;
401  if (trk.ok)
402  ptrk = spr::propagateCalo (trk.position, trk.momentum, trk.charge, bfield, 402.7, 180.7, 1.392, debug);
403  return ptrk;
404  }
405 
406  std::pair<math::XYZPoint,bool> propagateHCAL(const reco::Track *track, const MagneticField* bfield, bool debug) {
407  GlobalPoint vertex (track->vx(), track->vy(), track->vz());
408  GlobalVector momentum (track->px(), track->py(), track->pz());
409  int charge (track->charge());
410  return spr::propagateHCAL (vertex, momentum, charge, bfield, debug);
411  }
412 
413  std::pair<math::XYZPoint,bool> propagateHCAL(const GlobalPoint& vertex, const GlobalVector& momentum, int charge, const MagneticField* bfield, bool debug) {
414  spr::propagatedTrack track = spr::propagateCalo (vertex, momentum, charge, bfield, 402.7, 180.7, 1.392, debug);
415  return std::pair<math::XYZPoint,bool>(track.point,track.ok);
416  }
417 
418  std::pair<math::XYZPoint,bool> propagateTracker(const reco::Track *track, const MagneticField* bfield, bool debug) {
419  GlobalPoint vertex (track->vx(), track->vy(), track->vz());
420  GlobalVector momentum (track->px(), track->py(), track->pz());
421  int charge (track->charge());
422  spr::propagatedTrack track1 = spr::propagateCalo (vertex, momentum, charge, bfield, 290.0, 109.0, 1.705, debug);
423  return std::pair<math::XYZPoint,bool>(track1.point,track1.ok);
424  }
425 
426  std::pair<math::XYZPoint,double> propagateTrackerEnd(const reco::Track *track, const MagneticField* bField, bool debug) {
427 
428  GlobalPoint vertex (track->vx(), track->vy(), track->vz());
429  GlobalVector momentum (track->px(), track->py(), track->pz());
430  int charge (track->charge());
431  float radius = track->outerPosition().Rho();
432  float zdist = track->outerPosition().Z();
433  if (debug) std::cout << "propagateTrackerEnd:: Vertex " << vertex << " Momentum " << momentum << " Charge " << charge << " Radius " << radius << " Z " << zdist << std::endl;
434 
435  FreeTrajectoryState fts (vertex, momentum, charge, bField);
436  Plane::PlanePointer endcap = Plane::build(Plane::PositionType (0, 0, zdist), Plane::RotationType());
437  Cylinder::CylinderPointer barrel = Cylinder::build(Cylinder::PositionType (0, 0, 0), Cylinder::RotationType (), radius);
438 
439  AnalyticalPropagator myAP (bField, alongMomentum, 2*M_PI);
440 
441  TrajectoryStateOnSurface tsose = myAP.propagate(fts, *endcap);
442  TrajectoryStateOnSurface tsosb = myAP.propagate(fts, *barrel);
443 
444  math::XYZPoint point(-999.,-999.,-999.);
445  bool ok=false;
446  GlobalVector direction(0,0,1);
447  if (tsosb.isValid() && std::abs(zdist) < 110) {
448  point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
449  direction = tsosb.globalDirection();
450  ok = true;
451  } else if (tsose.isValid()) {
452  point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
453  direction = tsose.globalDirection();
454  ok = true;
455  }
456 
457  double length = -1;
458  if (ok) {
459  math::XYZPoint vDiff(point.x()-vertex.x(), point.y()-vertex.y(), point.z()-vertex.z());
460  double dphi = direction.phi()-momentum.phi();
461  double rdist = std::sqrt(vDiff.x()*vDiff.x()+vDiff.y()*vDiff.y());
462  double rat = 0.5*dphi/std::sin(0.5*dphi);
463  double dZ = vDiff.z();
464  double dS = rdist*rat; //dZ*momentum.z()/momentum.perp();
465  length = std::sqrt(dS*dS+dZ*dZ);
466  if (debug)
467  std::cout << "propagateTracker:: Barrel " << tsosb.isValid() << " Endcap " << tsose.isValid() << " OverAll " << ok << " Point " << point << " RDist " << rdist << " dS " << dS << " dS/pt " << rdist*rat/momentum.perp() << " zdist " << dZ << " dz/pz " << dZ/momentum.z() << " Length " << length << std::endl;
468  }
469 
470  return std::pair<math::XYZPoint,double>(point,length);
471  }
472 
473  spr::propagatedTrack propagateCalo(const GlobalPoint& tpVertex, const GlobalVector& tpMomentum, int tpCharge, const MagneticField* bField, float zdist, float radius, float corner, bool debug) {
474 
475  spr::propagatedTrack track;
476  if (debug) std::cout << "propagateCalo:: Vertex " << tpVertex << " Momentum " << tpMomentum << " Charge " << tpCharge << " Radius " << radius << " Z " << zdist << " Corner " << corner << std::endl;
477  FreeTrajectoryState fts (tpVertex, tpMomentum, tpCharge, bField);
478 
479  Plane::PlanePointer lendcap = Plane::build(Plane::PositionType (0, 0, -zdist), Plane::RotationType());
480  Plane::PlanePointer rendcap = Plane::build(Plane::PositionType (0, 0, zdist), Plane::RotationType());
481 
482  Cylinder::CylinderPointer barrel = Cylinder::build(Cylinder::PositionType (0, 0, 0), Cylinder::RotationType (), radius);
483 
484  AnalyticalPropagator myAP (bField, alongMomentum, 2*M_PI);
485 
486  TrajectoryStateOnSurface tsose;
487  if (tpMomentum.eta() < 0) {
488  tsose = myAP.propagate(fts, *lendcap);
489  } else {
490  tsose = myAP.propagate(fts, *rendcap);
491  }
492 
493  TrajectoryStateOnSurface tsosb = myAP.propagate(fts, *barrel);
494 
495  track.ok=true;
496  if (tsose.isValid() && tsosb.isValid()) {
497  float absEta = std::abs(tsosb.globalPosition().eta());
498  if (absEta < corner) {
499  track.point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
500  track.direction = tsosb.globalDirection();
501  } else {
502  track.point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
503  track.direction = tsose.globalDirection();
504  }
505  } else if (tsose.isValid()) {
506  track.point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
507  track.direction = tsose.globalDirection();
508  } else if (tsosb.isValid()) {
509  track.point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
510  track.direction = tsosb.globalDirection();
511  } else {
512  track.point.SetXYZ(-999., -999., -999.);
513  track.direction = GlobalVector(0,0,1);
514  track.ok = false;
515  }
516  if (debug) {
517  std::cout << "propagateCalo:: Barrel " << tsosb.isValid() << " Endcap " << tsose.isValid() << " OverAll " << track.ok << " Point " << track.point << " Direction " << track.direction << std::endl;
518  if (track.ok) {
519  math::XYZPoint vDiff(track.point.x()-tpVertex.x(), track.point.y()-tpVertex.y(), track.point.z()-tpVertex.z());
520  double dphi = track.direction.phi()-tpMomentum.phi();
521  double rdist = std::sqrt(vDiff.x()*vDiff.x()+vDiff.y()*vDiff.y());
522  double pt = tpMomentum.perp();
523  double rat = 0.5*dphi/std::sin(0.5*dphi);
524  std::cout << "RDist " << rdist << " pt " << pt << " r/pt " << rdist*rat/pt << " zdist " << vDiff.z() << " pz " << tpMomentum.z() << " z/pz " << vDiff.z()/tpMomentum.z() << std::endl;
525  }
526  }
527  return track;
528  }
529 
530  spr::trackAtOrigin simTrackAtOrigin(unsigned int thisTrk, edm::Handle<edm::SimTrackContainer>& SimTk, edm::Handle<edm::SimVertexContainer>& SimVtx, bool debug) {
531 
532  spr::trackAtOrigin trk;
533 
534  edm::SimTrackContainer::const_iterator itr = SimTk->end();
535  for (edm::SimTrackContainer::const_iterator simTrkItr = SimTk->begin(); simTrkItr!= SimTk->end(); simTrkItr++) {
536  if ( simTrkItr->trackId() == thisTrk ) {
537  if (debug) std::cout << "matched trackId (maximum occurance) " << thisTrk << " type " << simTrkItr->type() << std::endl;
538  itr = simTrkItr;
539  break;
540  }
541  }
542 
543  if (itr != SimTk->end()) {
544  int vertIndex = itr->vertIndex();
545  if (vertIndex != -1 && vertIndex < (int)SimVtx->size()) {
546  edm::SimVertexContainer::const_iterator simVtxItr= SimVtx->begin();
547  for (int iv=0; iv<vertIndex; iv++) simVtxItr++;
548  const math::XYZTLorentzVectorD pos = simVtxItr->position();
549  const math::XYZTLorentzVectorD mom = itr->momentum();
550  trk.ok = true;
551  trk.charge = (int)(itr->charge());
552  trk.position = GlobalPoint(pos.x(), pos.y(), pos.z());
553  trk.momentum = GlobalVector(mom.x(), mom.y(), mom.z());
554  }
555  }
556  if (debug) std::cout << "Track flag " << trk.ok << " Position " << trk.position << " Momentum " << trk.momentum << std::endl;;
557  return trk;
558  }
559 
560 }
EcalEndcapGeometry::getClosestCell
virtual DetId getClosestCell(const GlobalPoint &r) const
Definition: EcalEndcapGeometry.cc:230
CaloGeometry::getSubdetectorGeometry
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
Definition: CaloGeometry.cc:43
i
int i
Definition: DBlmapReader.cc:9
spr::propagateCALO
std::vector< spr::propagatedTrackID > propagateCALO(edm::Handle< reco::TrackCollection > &trkCollection, const CaloGeometry *geo, const MagneticField *bField, std::string &theTrackQuality, bool debug=false)
Definition: CaloPropagateTrack.cc:18
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:71
math::XYZTLorentzVectorD
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:15
AnalyticalPropagator::propagate
TrajectoryStateOnSurface propagate(const FreeTrajectoryState &fts, const Plane &plane) const
propagation to plane
Definition: AnalyticalPropagator.h:48
Reference_intrackfit_cff.barrel
barrel
Definition: Reference_intrackfit_cff.py:37
EcalBarrelGeometry::getClosestCell
virtual DetId getClosestCell(const GlobalPoint &r) const
Definition: EcalBarrelGeometry.cc:66
Reference_intrackfit_cff.endcap
endcap
Definition: Reference_intrackfit_cff.py:82
spr::propagatedGenParticleID::trkItr
reco::GenParticleCollection::const_iterator trkItr
Definition: CaloPropagateTrack.h:80
reco::TrackBase::TrackQuality
TrackQuality
track quality
Definition: TrackBase.h:95
spr::propagateHCAL
std::pair< math::XYZPoint, bool > propagateHCAL(const reco::Track *, const MagneticField *, bool debug=false)
Definition: CaloPropagateTrack.cc:406
align::ID
uint32_t ID
Definition: Definitions.h:26
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:68
GlobalPoint
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
PV3DBase::y
T y() const
Definition: PV3DBase.h:62
abs
#define abs(x)
Definition: mlp_lapack.h:159
TrajectoryStateOnSurface::globalPosition
GlobalPoint globalPosition() const
Definition: TrajectoryStateOnSurface.h:143
Plane::build
static PlanePointer build(const PositionType &pos, const RotationType &rot, MediumProperties *mp=0)
Definition: Plane.h:25
spr::simTrackAtOrigin
spr::trackAtOrigin simTrackAtOrigin(unsigned int thisTrk, edm::Handle< edm::SimTrackContainer > &SimTk, edm::Handle< edm::SimVertexContainer > &SimVtx, bool debug=false)
Definition: CaloPropagateTrack.cc:530
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:473
reco::TrackBase::px
double px() const
x coordinate of momentum vector
Definition: TrackBase.h:133
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:47
DeDxDiscriminatorTools::charge
double charge(const std::vector< uint8_t > &Ampls)
Definition: DeDxDiscriminatorTools.cc:43
spr::propagatedGenTrackID::trkItr
HepMC::GenEvent::particle_const_iterator trkItr
Definition: CaloPropagateTrack.h:67
spr::propagatedTrackID::trkItr
reco::TrackCollection::const_iterator trkItr
Definition: CaloPropagateTrack.h:45
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
PV3DBase::z
T z() const
Definition: PV3DBase.h:63
ecalTB2006H4_GenSimDigiReco_cfg.bField
tuple bField
Definition: ecalTB2006H4_GenSimDigiReco_cfg.py:239
MCTruth.genEvent
tuple genEvent
Definition: MCTruth.py:33
spr::propagateTracker
std::pair< math::XYZPoint, bool > propagateTracker(const reco::Track *, const MagneticField *, bool debug=false)
Definition: CaloPropagateTrack.cc:418
spr::propagateECAL
std::pair< math::XYZPoint, bool > propagateECAL(const reco::Track *, const MagneticField *, bool debug=false)
Definition: CaloPropagateTrack.cc:379
spr::propagateTrackToECAL
spr::propagatedTrack propagateTrackToECAL(const reco::Track *, const MagneticField *, bool debug=false)
Definition: CaloPropagateTrack.cc:363
DetId::subdetId
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:39
reco::TrackBase::pz
double pz() const
z coordinate of momentum vector
Definition: TrackBase.h:137
CaloSubdetectorGeometry::getClosestCell
virtual DetId getClosestCell(const GlobalPoint &r) const
Definition: CaloSubdetectorGeometry.cc:55
reco::TrackBase::vz
double vz() const
z coordinate of the reference point on track
Definition: TrackBase.h:147
DetId
Definition: DetId.h:20
reco::TrackBase::qualityByName
static TrackQuality qualityByName(const std::string &name)
Definition: TrackBase.cc:46
M_PI
#define M_PI
Definition: BFit3D.cc:3
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:13
Cylinder::build
static CylinderPointer build(const PositionType &pos, const RotationType &rot, Scalar radius, MediumProperties *mp=0)
Definition: Cylinder.h:29
reco::TrackBase::quality
bool quality(const TrackQuality) const
Track quality.
Definition: TrackBase.h:377
PV3DBase::eta
T eta() const
Definition: PV3DBase.h:75
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
reco::TrackBase::vy
double vy() const
y coordinate of the reference point on track
Definition: TrackBase.h:145
jptDQMConfig_cff.etaMax
float etaMax
Definition: jptDQMConfig_cff.py:9
spr::propagateTrackerEnd
std::pair< math::XYZPoint, double > propagateTrackerEnd(const reco::Track *, const MagneticField *, bool debug=false)
Definition: CaloPropagateTrack.cc:426
spr::propagatedTrackDirection::trkItr
reco::TrackCollection::const_iterator trkItr
Definition: CaloPropagateTrack.h:54
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
reco::TrackBase::charge
int charge() const
track electric charge
Definition: TrackBase.h:113
debug
#define debug
Definition: MEtoEDMFormat.h:34
PV3DBase::x
T x() const
Definition: PV3DBase.h:61
spr::propagateTrackToHCAL
spr::propagatedTrack propagateTrackToHCAL(const reco::Track *, const MagneticField *, bool debug=false)
Definition: CaloPropagateTrack.cc:391
CosmicsPD_Skims.radius
tuple radius
Definition: CosmicsPD_Skims.py:135
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:135
reco::TrackBase::vx
double vx() const
x coordinate of the reference point on track
Definition: TrackBase.h:143
GlobalVector
Global3DVector GlobalVector
Definition: GlobalVector.h:10
TrajectoryStateOnSurface::globalDirection
GlobalVector globalDirection() const
Definition: TrajectoryStateOnSurface.h:149
Generated for CMSSW Reference Manual by doxygen 1.8.5