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.
8 
11 
13 
14 #include <iostream>
15 
16 namespace spr{
17 
18  std::vector<spr::propagatedTrackID> propagateCosmicCALO(edm::Handle<reco::TrackCollection>& trkCollection, const CaloGeometry* geo, const MagneticField* bField, std::string & theTrackQuality, bool debug) {
19 
20  const EcalBarrelGeometry *barrelGeom = (dynamic_cast< const EcalBarrelGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel)));
21  const EcalEndcapGeometry *endcapGeom = (dynamic_cast< const EcalEndcapGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap)));
24  std::vector<spr::propagatedTrackID> vdets;
25 
26  unsigned int indx;
27  reco::TrackCollection::const_iterator trkItr;
28  for (trkItr = trkCollection->begin(),indx=0; trkItr != trkCollection->end(); ++trkItr,indx++) {
29  const reco::Track* pTrack = &(*trkItr);
31  vdet.trkItr = trkItr;
32  vdet.ok = (pTrack->quality(trackQuality_));
33  vdet.detIdECAL = DetId(0);
34  vdet.detIdHCAL = DetId(0);
35  vdet.detIdEHCAL= DetId(0);
36  if (debug) std::cout << "Propagate track " << indx << " p " << trkItr->p() << " eta " << trkItr->eta() << " phi " << trkItr->phi() << " Flag " << vdet.ok << std::endl;
37 
38  GlobalPoint vertex;
39  GlobalVector momentum;
40  int charge (pTrack->charge());
41  if (((pTrack->innerPosition()).Perp2()) <
42  ((pTrack->outerPosition()).Perp2())) {
43  vertex = GlobalPoint(((pTrack->innerPosition()).X()),
44  ((pTrack->innerPosition()).Y()),
45  ((pTrack->innerPosition()).Z()));
46  momentum = GlobalVector(((pTrack->innerMomentum()).X()),
47  ((pTrack->innerMomentum()).Y()),
48  ((pTrack->innerMomentum()).Z()));
49  } else {
50  vertex = GlobalPoint(((pTrack->outerPosition()).X()),
51  ((pTrack->outerPosition()).Y()),
52  ((pTrack->outerPosition()).Z()));
53  momentum = GlobalVector(((pTrack->outerMomentum()).X()),
54  ((pTrack->outerMomentum()).Y()),
55  ((pTrack->outerMomentum()).Z()));
56  }
57  if (debug) std::cout << "Track charge " << charge << " p " << momentum << " position " << vertex << std::endl;
58 
59  std::pair<math::XYZPoint,bool> info =
60  spr::propagateECAL (vertex, momentum, charge, bField, debug);
61 
62  vdet.okECAL = info.second;
63  if (vdet.okECAL) {
64  const GlobalPoint point(info.first.x(),info.first.y(),info.first.z());
65  vdet.etaECAL = point.eta();
66  vdet.phiECAL = point.phi();
67  if (std::abs(point.eta())<1.479) {
68  vdet.detIdECAL = barrelGeom->getClosestCell(point);
69  } else {
70  vdet.detIdECAL = endcapGeom->getClosestCell(point);
71  }
72  vdet.detIdEHCAL = gHB->getClosestCell(point);
73  }
74  info = spr::propagateHCAL (vertex, momentum, charge, bField, debug);
75  vdet.okHCAL = info.second;
76  if (vdet.okHCAL) {
77  const GlobalPoint point(info.first.x(),info.first.y(),info.first.z());
78  vdet.etaHCAL = point.eta();
79  vdet.phiHCAL = point.phi();
80  vdet.detIdHCAL = gHB->getClosestCell(point);
81  }
82  if (debug) {
83  std::cout << "Track [" << indx << "] Flag: " << vdet.ok << " ECAL ("
84  << vdet.okECAL << ") ";
85  if (vdet.detIdECAL.subdetId() == EcalBarrel)
86  std::cout << (EBDetId)(vdet.detIdECAL);
87  else
88  std::cout << (EEDetId)(vdet.detIdECAL);
89  std::cout << " HCAL (" << vdet.okHCAL << ") " << (HcalDetId)(vdet.detIdHCAL) << " Or " << (HcalDetId)(vdet.detIdEHCAL) << std::endl;
90  }
91  vdets.push_back(vdet);
92  }
93 
94  if (debug) {
95  std::cout << "propagateCALO:: for " << vdets.size() << " tracks" << std::endl;
96  for (unsigned int i=0; i<vdets.size(); ++i) {
97  std::cout << "Track [" << i << "] Flag: " << vdets[i].ok << " ECAL (" << vdets[i].okECAL << ") ";
98  if (vdets[i].detIdECAL.subdetId() == EcalBarrel) {
99  std::cout << (EBDetId)(vdets[i].detIdECAL);
100  } else {
101  std::cout << (EEDetId)(vdets[i].detIdECAL);
102  }
103  std::cout << " HCAL (" << vdets[i].okHCAL << ") " << (HcalDetId)(vdets[i].detIdHCAL) << " Or " << (HcalDetId)(vdets[i].detIdEHCAL) << std::endl;
104  }
105  }
106  return vdets;
107  }
108 
109  std::vector<spr::propagatedTrackID> propagateCALO(edm::Handle<reco::TrackCollection>& trkCollection, const CaloGeometry* geo, const MagneticField* bField, std::string & theTrackQuality, bool debug) {
110 
111  std::vector<spr::propagatedTrackID> vdets;
112  spr::propagateCALO(trkCollection,geo,bField,theTrackQuality, vdets, debug);
113  return vdets;
114  }
115 
116  void propagateCALO(edm::Handle<reco::TrackCollection>& trkCollection, const CaloGeometry* geo, const MagneticField* bField, std::string & theTrackQuality, std::vector<spr::propagatedTrackID>& vdets, bool debug) {
117 
118  const EcalBarrelGeometry *barrelGeom = (dynamic_cast< const EcalBarrelGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel)));
119  const EcalEndcapGeometry *endcapGeom = (dynamic_cast< const EcalEndcapGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap)));
121  reco::TrackBase::TrackQuality trackQuality_=reco::TrackBase::qualityByName(theTrackQuality);
122 
123  unsigned int indx;
124  reco::TrackCollection::const_iterator trkItr;
125  for (trkItr = trkCollection->begin(),indx=0; trkItr != trkCollection->end(); ++trkItr,indx++) {
126  const reco::Track* pTrack = &(*trkItr);
128  vdet.trkItr = trkItr;
129  vdet.ok = (pTrack->quality(trackQuality_));
130  vdet.detIdECAL = DetId(0);
131  vdet.detIdHCAL = DetId(0);
132  vdet.detIdEHCAL= DetId(0);
133  if (debug) std::cout << "Propagate track " << indx << " p " << trkItr->p() << " eta " << trkItr->eta() << " phi " << trkItr->phi() << " Flag " << vdet.ok << std::endl;
134 
135  std::pair<math::XYZPoint,bool> info = spr::propagateECAL (pTrack, bField, debug);
136  vdet.okECAL = info.second;
137  if (vdet.okECAL) {
138  const GlobalPoint point(info.first.x(),info.first.y(),info.first.z());
139  vdet.etaECAL = point.eta();
140  vdet.phiECAL = point.phi();
141  if (std::abs(point.eta())<1.479) {
142  vdet.detIdECAL = barrelGeom->getClosestCell(point);
143  } else {
144  vdet.detIdECAL = endcapGeom->getClosestCell(point);
145  }
146  vdet.detIdEHCAL = gHB->getClosestCell(point);
147  }
148  info = spr::propagateHCAL (pTrack, bField, debug);
149  vdet.okHCAL = info.second;
150  if (vdet.okHCAL) {
151  const GlobalPoint point(info.first.x(),info.first.y(),info.first.z());
152  vdet.etaHCAL = point.eta();
153  vdet.phiHCAL = point.phi();
154  vdet.detIdHCAL = gHB->getClosestCell(point);
155  }
156  if (debug) {
157  std::cout << "Track [" << indx << "] Flag: " << vdet.ok << " ECAL ("
158  << vdet.okECAL << ") ";
159  if (vdet.detIdECAL.subdetId() == EcalBarrel)
160  std::cout << (EBDetId)(vdet.detIdECAL);
161  else
162  std::cout << (EEDetId)(vdet.detIdECAL);
163  std::cout << " HCAL (" << vdet.okHCAL << ") " << (HcalDetId)(vdet.detIdHCAL) << " Or " << (HcalDetId)(vdet.detIdEHCAL) << std::endl;
164  }
165  vdets.push_back(vdet);
166  }
167 
168  if (debug) {
169  std::cout << "propagateCALO:: for " << vdets.size() << " tracks" << std::endl;
170  for (unsigned int i=0; i<vdets.size(); ++i) {
171  std::cout << "Track [" << i << "] Flag: " << vdets[i].ok << " ECAL (" << vdets[i].okECAL << ") ";
172  if (vdets[i].detIdECAL.subdetId() == EcalBarrel) {
173  std::cout << (EBDetId)(vdets[i].detIdECAL);
174  } else {
175  std::cout << (EEDetId)(vdets[i].detIdECAL);
176  }
177  std::cout << " HCAL (" << vdets[i].okHCAL << ") " << (HcalDetId)(vdets[i].detIdHCAL) << " Or " << (HcalDetId)(vdets[i].detIdEHCAL) << std::endl;
178  }
179  }
180  }
181 
182  void propagateCALO(edm::Handle<reco::TrackCollection>& trkCollection, const CaloGeometry* geo, const MagneticField* bField, std::string & theTrackQuality, std::vector<spr::propagatedTrackDirection>& trkDir, bool debug) {
183 
184  const EcalBarrelGeometry *barrelGeom = (dynamic_cast< const EcalBarrelGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel)));
185  const EcalEndcapGeometry *endcapGeom = (dynamic_cast< const EcalEndcapGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap)));
187  reco::TrackBase::TrackQuality trackQuality_=reco::TrackBase::qualityByName(theTrackQuality);
188 
189  unsigned int indx;
190  reco::TrackCollection::const_iterator trkItr;
191  for (trkItr = trkCollection->begin(),indx=0; trkItr != trkCollection->end(); ++trkItr,indx++) {
192  const reco::Track* pTrack = &(*trkItr);
194  trkD.trkItr = trkItr;
195  trkD.ok = (pTrack->quality(trackQuality_));
196  trkD.detIdECAL = DetId(0);
197  trkD.detIdHCAL = DetId(0);
198  trkD.detIdEHCAL= DetId(0);
199  if (debug) std::cout << "Propagate track " << indx << " p " << trkItr->p() << " eta " << trkItr->eta() << " phi " << trkItr->phi() << " Flag " << trkD.ok << std::endl;
200 
201  spr::propagatedTrack info = spr::propagateTrackToECAL (pTrack, bField, debug);
202  GlobalPoint point(info.point.x(),info.point.y(),info.point.z());
203  trkD.okECAL = info.ok;
204  trkD.pointECAL = point;
205  trkD.directionECAL = info.direction;
206  if (trkD.okECAL) {
207  if (std::abs(info.point.eta())<1.479) {
208  trkD.detIdECAL = barrelGeom->getClosestCell(point);
209  } else {
210  trkD.detIdECAL = endcapGeom->getClosestCell(point);
211  }
212  trkD.detIdEHCAL = gHB->getClosestCell(point);
213  }
214  info = spr::propagateTrackToHCAL (pTrack, bField, debug);
215  point = GlobalPoint(info.point.x(),info.point.y(),info.point.z());
216  trkD.okHCAL = info.ok;
217  trkD.pointHCAL = point;
218  trkD.directionHCAL = info.direction;
219  if (trkD.okHCAL) {
220  trkD.detIdHCAL = gHB->getClosestCell(point);
221  }
222  trkDir.push_back(trkD);
223  }
224 
225  if (debug) {
226  std::cout << "propagateCALO:: for " << trkDir.size() << " tracks" << std::endl;
227  for (unsigned int i=0; i<trkDir.size(); ++i) {
228  std::cout << "Track [" << i << "] Flag: " << trkDir[i].ok << " ECAL (" << trkDir[i].okECAL << ")";
229  if (trkDir[i].okECAL) {
230  std::cout << " point " << trkDir[i].pointECAL << " direction "
231  << trkDir[i].directionECAL << " ";
232  if (trkDir[i].detIdECAL.subdetId() == EcalBarrel) {
233  std::cout << (EBDetId)(trkDir[i].detIdECAL);
234  } else {
235  std::cout << (EEDetId)(trkDir[i].detIdECAL);
236  }
237  }
238  std::cout << " HCAL (" << trkDir[i].okHCAL << ")";
239  if (trkDir[i].okHCAL) {
240  std::cout << " point " << trkDir[i].pointHCAL << " direction "
241  << trkDir[i].directionHCAL << " "
242  << (HcalDetId)(trkDir[i].detIdHCAL);
243  }
244  std::cout << " Or " << (HcalDetId)(trkDir[i].detIdEHCAL) << std::endl;
245  }
246  }
247  }
248 
250 
251  const EcalBarrelGeometry *barrelGeom = (dynamic_cast< const EcalBarrelGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel)));
252  const EcalEndcapGeometry *endcapGeom = (dynamic_cast< const EcalEndcapGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap)));
254 
256  vdet.ok = true;
257  vdet.detIdECAL = DetId(0);
258  vdet.detIdHCAL = DetId(0);
259  vdet.detIdEHCAL= DetId(0);
260  if (debug) std::cout << "Propagate track: p " << pTrack->p() << " eta " << pTrack->eta() << " phi " << pTrack->phi() << " Flag " << vdet.ok << std::endl;
261 
262  std::pair<math::XYZPoint,bool> info = spr::propagateECAL (pTrack, bField, debug);
263  vdet.okECAL = info.second;
264  if (vdet.okECAL) {
265  const GlobalPoint point(info.first.x(),info.first.y(),info.first.z());
266  vdet.etaECAL = point.eta();
267  vdet.phiECAL = point.phi();
268  if (std::abs(point.eta())<1.479) {
269  vdet.detIdECAL = barrelGeom->getClosestCell(point);
270  } else {
271  vdet.detIdECAL = endcapGeom->getClosestCell(point);
272  }
273  vdet.detIdEHCAL = gHB->getClosestCell(point);
274  }
275  info = spr::propagateHCAL (pTrack, bField, debug);
276  vdet.okHCAL = info.second;
277  if (vdet.okHCAL) {
278  const GlobalPoint point(info.first.x(),info.first.y(),info.first.z());
279  vdet.etaHCAL = point.eta();
280  vdet.phiHCAL = point.phi();
281  vdet.detIdHCAL = gHB->getClosestCell(point);
282  }
283 
284  if (debug) {
285  std::cout << "propagateCALO:: for 1 track" << std::endl;
286  std::cout << "Track [0] Flag: " << vdet.ok << " ECAL (" << vdet.okECAL << ") ";
287  if (vdet.detIdECAL.subdetId() == EcalBarrel) {
288  std::cout << (EBDetId)(vdet.detIdECAL);
289  } else {
290  std::cout << (EEDetId)(vdet.detIdECAL);
291  }
292  std::cout << " HCAL (" << vdet.okHCAL << ") " << (HcalDetId)(vdet.detIdHCAL) << " Or " << (HcalDetId)(vdet.detIdEHCAL) << std::endl;
293  }
294  return vdet;
295  }
296 
297  std::vector<spr::propagatedGenTrackID> propagateCALO(const HepMC::GenEvent * genEvent, edm::ESHandle<ParticleDataTable>& pdt, const CaloGeometry* geo, const MagneticField* bField, double etaMax, bool debug) {
298 
299  const EcalBarrelGeometry *barrelGeom = (dynamic_cast< const EcalBarrelGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel)));
300  const EcalEndcapGeometry *endcapGeom = (dynamic_cast< const EcalEndcapGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap)));
302 
303  std::vector<spr::propagatedGenTrackID> trkDir;
304  unsigned int indx;
305  HepMC::GenEvent::particle_const_iterator p;
306  for (p=genEvent->particles_begin(),indx=0; p != genEvent->particles_end(); ++p,++indx) {
308  trkD.trkItr = p;
309  trkD.detIdECAL = DetId(0);
310  trkD.detIdHCAL = DetId(0);
311  trkD.detIdEHCAL= DetId(0);
312  trkD.pdgId = ((*p)->pdg_id());
313  trkD.charge = ((pdt->particle(trkD.pdgId))->ID().threeCharge())/3;
314  GlobalVector momentum = GlobalVector((*p)->momentum().px(), (*p)->momentum().py(), (*p)->momentum().pz());
315  if (debug) std::cout << "Propagate track " << indx << " pdg " << trkD.pdgId << " charge " << trkD.charge << " p " << momentum << std::endl;
316 
317  // consider stable particles
318  if ( (*p)->status()==1 && std::abs((*p)->momentum().eta()) < etaMax ) {
319  GlobalPoint vertex = GlobalPoint(0.1*(*p)->production_vertex()->position().x(),
320  0.1*(*p)->production_vertex()->position().y(),
321  0.1*(*p)->production_vertex()->position().z());
322  trkD.ok = true;
323  spr::propagatedTrack info = spr::propagateCalo (vertex, momentum, trkD.charge, bField, 319.2, 129.4, 1.479, debug);
324  GlobalPoint point(info.point.x(),info.point.y(),info.point.z());
325  trkD.okECAL = info.ok;
326  trkD.pointECAL = point;
327  trkD.directionECAL = info.direction;
328  if (trkD.okECAL) {
329  if (std::abs(info.point.eta())<1.479) {
330  trkD.detIdECAL = barrelGeom->getClosestCell(point);
331  } else {
332  trkD.detIdECAL = endcapGeom->getClosestCell(point);
333  }
334  trkD.detIdEHCAL = gHB->getClosestCell(point);
335  }
336 
337  info = spr::propagateCalo (vertex, momentum, trkD.charge, bField, 402.7, 180.7, 1.392, debug);
338  point = GlobalPoint(info.point.x(),info.point.y(),info.point.z());
339  trkD.okHCAL = info.ok;
340  trkD.pointHCAL = point;
341  trkD.directionHCAL = info.direction;
342  if (trkD.okHCAL) {
343  trkD.detIdHCAL = gHB->getClosestCell(point);
344  }
345  }
346  trkDir.push_back(trkD);
347  }
348 
349  if (debug) {
350  std::cout << "propagateCALO:: for " << trkDir.size() << " tracks" << std::endl;
351  for (unsigned int i=0; i<trkDir.size(); ++i) {
352  if (trkDir[i].okECAL) std::cout << "Track [" << i << "] Flag: " << trkDir[i].ok << " ECAL (" << trkDir[i].okECAL << ")";
353  if (trkDir[i].okECAL) {
354  std::cout << " point " << trkDir[i].pointECAL << " direction "
355  << trkDir[i].directionECAL << " ";
356  if (trkDir[i].detIdECAL.subdetId() == EcalBarrel) {
357  std::cout << (EBDetId)(trkDir[i].detIdECAL);
358  } else {
359  std::cout << (EEDetId)(trkDir[i].detIdECAL);
360  }
361  }
362  if (trkDir[i].okECAL) std::cout << " HCAL (" << trkDir[i].okHCAL << ")";
363  if (trkDir[i].okHCAL) {
364  std::cout << " point " << trkDir[i].pointHCAL << " direction "
365  << trkDir[i].directionHCAL << " "
366  << (HcalDetId)(trkDir[i].detIdHCAL);
367  }
368  if (trkDir[i].okECAL) std::cout << " Or " << (HcalDetId)(trkDir[i].detIdEHCAL) << std::endl;
369  }
370  }
371  return trkDir;
372  }
373 
374  std::vector<spr::propagatedGenParticleID> propagateCALO(edm::Handle<reco::GenParticleCollection>& genParticles, edm::ESHandle<ParticleDataTable>& pdt, const CaloGeometry* geo, const MagneticField* bField, double etaMax, bool debug) {
375 
376  const EcalBarrelGeometry *barrelGeom = (dynamic_cast< const EcalBarrelGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel)));
377  const EcalEndcapGeometry *endcapGeom = (dynamic_cast< const EcalEndcapGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap)));
379 
380  std::vector<spr::propagatedGenParticleID> trkDir;
381  unsigned int indx;
382  reco::GenParticleCollection::const_iterator p;
383  for (p=genParticles->begin(),indx=0; p != genParticles->end(); ++p,++indx) {
385  trkD.trkItr = p;
386  trkD.detIdECAL = DetId(0);
387  trkD.detIdHCAL = DetId(0);
388  trkD.detIdEHCAL= DetId(0);
389  trkD.pdgId = (p->pdgId());
390  trkD.charge = p->charge();
391  GlobalVector momentum = GlobalVector(p->momentum().x(), p->momentum().y(), p->momentum().z());
392  if (debug) std::cout << "Propagate track " << indx << " pdg " << trkD.pdgId << " charge " << trkD.charge << " p " << momentum << std::endl;
393 
394  // consider stable particles
395  if ( p->status()==1 && std::abs(momentum.eta()) < etaMax ) {
396  GlobalPoint vertex = GlobalPoint(p->vertex().x(), p->vertex().y(), p->vertex().z());
397  trkD.ok = true;
398  spr::propagatedTrack info = spr::propagateCalo (vertex, momentum, trkD.charge, bField, 319.2, 129.4, 1.479, debug);
399  GlobalPoint point(info.point.x(),info.point.y(),info.point.z());
400  trkD.okECAL = info.ok;
401  trkD.pointECAL = point;
402  trkD.directionECAL = info.direction;
403  if (trkD.okECAL) {
404  if (std::abs(info.point.eta())<1.479) {
405  trkD.detIdECAL = barrelGeom->getClosestCell(point);
406  } else {
407  trkD.detIdECAL = endcapGeom->getClosestCell(point);
408  }
409  trkD.detIdEHCAL = gHB->getClosestCell(point);
410  }
411 
412  info = spr::propagateCalo (vertex, momentum, trkD.charge, bField, 402.7, 180.7, 1.392, debug);
413  point = GlobalPoint(info.point.x(),info.point.y(),info.point.z());
414  trkD.okHCAL = info.ok;
415  trkD.pointHCAL = point;
416  trkD.directionHCAL = info.direction;
417  if (trkD.okHCAL) {
418  trkD.detIdHCAL = gHB->getClosestCell(point);
419  }
420  }
421  trkDir.push_back(trkD);
422  }
423 
424  if (debug) {
425  std::cout << "propagateCALO:: for " << trkDir.size() << " tracks" << std::endl;
426  for (unsigned int i=0; i<trkDir.size(); ++i) {
427  if (trkDir[i].okECAL) std::cout << "Track [" << i << "] Flag: " << trkDir[i].ok << " ECAL (" << trkDir[i].okECAL << ")";
428  if (trkDir[i].okECAL) {
429  std::cout << " point " << trkDir[i].pointECAL << " direction "
430  << trkDir[i].directionECAL << " ";
431  if (trkDir[i].detIdECAL.subdetId() == EcalBarrel) {
432  std::cout << (EBDetId)(trkDir[i].detIdECAL);
433  } else {
434  std::cout << (EEDetId)(trkDir[i].detIdECAL);
435  }
436  }
437  if (trkDir[i].okECAL) std::cout << " HCAL (" << trkDir[i].okHCAL << ")";
438  if (trkDir[i].okHCAL) {
439  std::cout << " point " << trkDir[i].pointHCAL << " direction "
440  << trkDir[i].directionHCAL << " "
441  << (HcalDetId)(trkDir[i].detIdHCAL);
442  }
443  if (trkDir[i].okECAL) std::cout << " Or " << (HcalDetId)(trkDir[i].detIdEHCAL) << std::endl;
444  }
445  }
446  return trkDir;
447  }
448 
450 
451  const EcalBarrelGeometry *barrelGeom = (dynamic_cast< const EcalBarrelGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel)));
452  const EcalEndcapGeometry *endcapGeom = (dynamic_cast< const EcalEndcapGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap)));
454 
455  spr::trackAtOrigin trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
457  trkD.ok = trk.ok;
458  trkD.detIdECAL = DetId(0);
459  trkD.detIdHCAL = DetId(0);
460  trkD.detIdEHCAL= DetId(0);
461  if (debug) std::cout << "Propagate track " << thisTrk << " charge " << trk.charge << " position " << trk.position << " p " << trk.momentum << " Flag " << trkD.ok << std::endl;
462 
463  if (trkD.ok) {
464  spr::propagatedTrack info = spr::propagateCalo (trk.position, trk.momentum, trk.charge, bField, 319.2, 129.4, 1.479, debug);
465  GlobalPoint point(info.point.x(),info.point.y(),info.point.z());
466  trkD.okECAL = info.ok;
467  trkD.pointECAL = point;
468  trkD.directionECAL = info.direction;
469  if (trkD.okECAL) {
470  if (std::abs(info.point.eta())<1.479) {
471  trkD.detIdECAL = barrelGeom->getClosestCell(point);
472  } else {
473  trkD.detIdECAL = endcapGeom->getClosestCell(point);
474  }
475  trkD.detIdEHCAL = gHB->getClosestCell(point);
476  }
477 
478  info = spr::propagateCalo (trk.position, trk.momentum, trk.charge, bField, 402.7, 180.7, 1.392, debug);
479  point = GlobalPoint(info.point.x(),info.point.y(),info.point.z());
480  trkD.okHCAL = info.ok;
481  trkD.pointHCAL = point;
482  trkD.directionHCAL = info.direction;
483  if (trkD.okHCAL) {
484  trkD.detIdHCAL = gHB->getClosestCell(point);
485  }
486  }
487 
488  if (debug) {
489  std::cout << "propagateCALO:: for track [" << thisTrk << "] Flag: " << trkD.ok << " ECAL (" << trkD.okECAL << ") HCAL (" << trkD.okHCAL << ")" << std::endl;
490  if (trkD.okECAL) {
491  std::cout << "ECAL point " << trkD.pointECAL << " direction "
492  << trkD.directionECAL << " ";
493  if (trkD.detIdECAL.subdetId() == EcalBarrel) {
494  std::cout << (EBDetId)(trkD.detIdECAL);
495  } else {
496  std::cout << (EEDetId)(trkD.detIdECAL);
497  }
498  }
499  if (trkD.okHCAL) {
500  std::cout << " HCAL point " << trkD.pointHCAL << " direction "
501  << trkD.directionHCAL << " " << (HcalDetId)(trkD.detIdHCAL);
502  }
503  if (trkD.okECAL) std::cout << " Or " << (HcalDetId)(trkD.detIdEHCAL);
504  std::cout << std::endl;
505  }
506 
507  return trkD;
508  }
509 
511  GlobalPoint vertex (track->vx(), track->vy(), track->vz());
512  GlobalVector momentum (track->px(), track->py(), track->pz());
513  int charge (track->charge());
514  return spr::propagateCalo (vertex, momentum, charge, bfield, 319.2, 129.4, 1.479, debug);
515  }
516 
518 
519  spr::trackAtOrigin trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
521  if (trk.ok)
522  ptrk = spr::propagateCalo (trk.position, trk.momentum, trk.charge, bfield, 319.2, 129.4, 1.479, debug);
523  return ptrk;
524  }
525 
526  std::pair<math::XYZPoint,bool> propagateECAL(const reco::Track *track, const MagneticField* bfield, bool debug) {
527  GlobalPoint vertex (track->vx(), track->vy(), track->vz());
528  GlobalVector momentum (track->px(), track->py(), track->pz());
529  int charge (track->charge());
530  return spr::propagateECAL (vertex, momentum, charge, bfield, debug);
531  }
532 
533  std::pair<math::XYZPoint,bool> propagateECAL(const GlobalPoint& vertex, const GlobalVector& momentum, int charge, const MagneticField* bfield, bool debug) {
534  spr::propagatedTrack track = spr::propagateCalo (vertex, momentum, charge, bfield, 319.2, 129.4, 1.479, debug);
535  return std::pair<math::XYZPoint,bool>(track.point,track.ok);
536  }
537 
539  GlobalPoint vertex (track->vx(), track->vy(), track->vz());
540  GlobalVector momentum (track->px(), track->py(), track->pz());
541  int charge (track->charge());
542  return spr::propagateCalo (vertex, momentum, charge, bfield, 402.7, 180.7, 1.392, debug);
543  }
544 
546  spr::trackAtOrigin trk = spr::simTrackAtOrigin(thisTrk, SimTk, SimVtx, debug);
548  if (trk.ok)
549  ptrk = spr::propagateCalo (trk.position, trk.momentum, trk.charge, bfield, 402.7, 180.7, 1.392, debug);
550  return ptrk;
551  }
552 
553  std::pair<math::XYZPoint,bool> propagateHCAL(const reco::Track *track, const MagneticField* bfield, bool debug) {
554  GlobalPoint vertex (track->vx(), track->vy(), track->vz());
555  GlobalVector momentum (track->px(), track->py(), track->pz());
556  int charge (track->charge());
557  return spr::propagateHCAL (vertex, momentum, charge, bfield, debug);
558  }
559 
560  std::pair<math::XYZPoint,bool> propagateHCAL(const GlobalPoint& vertex, const GlobalVector& momentum, int charge, const MagneticField* bfield, bool debug) {
561  spr::propagatedTrack track = spr::propagateCalo (vertex, momentum, charge, bfield, 402.7, 180.7, 1.392, debug);
562  return std::pair<math::XYZPoint,bool>(track.point,track.ok);
563  }
564 
565  std::pair<math::XYZPoint,bool> propagateTracker(const reco::Track *track, const MagneticField* bfield, bool debug) {
566  GlobalPoint vertex (track->vx(), track->vy(), track->vz());
567  GlobalVector momentum (track->px(), track->py(), track->pz());
568  int charge (track->charge());
569  spr::propagatedTrack track1 = spr::propagateCalo (vertex, momentum, charge, bfield, 290.0, 109.0, 1.705, debug);
570  return std::pair<math::XYZPoint,bool>(track1.point,track1.ok);
571  }
572 
573  std::pair<math::XYZPoint,double> propagateTrackerEnd(const reco::Track *track, const MagneticField* bField, bool debug) {
574 
575  GlobalPoint vertex (track->vx(), track->vy(), track->vz());
576  GlobalVector momentum (track->px(), track->py(), track->pz());
577  int charge (track->charge());
578  float radius = track->outerPosition().Rho();
579  float zdist = track->outerPosition().Z();
580  if (debug) std::cout << "propagateTrackerEnd:: Vertex " << vertex << " Momentum " << momentum << " Charge " << charge << " Radius " << radius << " Z " << zdist << std::endl;
581 
582  FreeTrajectoryState fts (vertex, momentum, charge, bField);
585 
586  AnalyticalPropagator myAP (bField, alongMomentum, 2*M_PI);
587 
588  TrajectoryStateOnSurface tsose = myAP.propagate(fts, *endcap);
589  TrajectoryStateOnSurface tsosb = myAP.propagate(fts, *barrel);
590 
591  math::XYZPoint point(-999.,-999.,-999.);
592  bool ok=false;
593  GlobalVector direction(0,0,1);
594  if (tsosb.isValid() && std::abs(zdist) < 110) {
595  point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
596  direction = tsosb.globalDirection();
597  ok = true;
598  } else if (tsose.isValid()) {
599  point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
600  direction = tsose.globalDirection();
601  ok = true;
602  }
603 
604  double length = -1;
605  if (ok) {
606  math::XYZPoint vDiff(point.x()-vertex.x(), point.y()-vertex.y(), point.z()-vertex.z());
607  double dphi = direction.phi()-momentum.phi();
608  double rdist = std::sqrt(vDiff.x()*vDiff.x()+vDiff.y()*vDiff.y());
609  double rat = 0.5*dphi/std::sin(0.5*dphi);
610  double dZ = vDiff.z();
611  double dS = rdist*rat; //dZ*momentum.z()/momentum.perp();
612  length = std::sqrt(dS*dS+dZ*dZ);
613  if (debug)
614  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;
615  }
616 
617  return std::pair<math::XYZPoint,double>(point,length);
618  }
619 
620  spr::propagatedTrack propagateCalo(const GlobalPoint& tpVertex, const GlobalVector& tpMomentum, int tpCharge, const MagneticField* bField, float zdist, float radius, float corner, bool debug) {
621 
622  spr::propagatedTrack track;
623  if (debug) std::cout << "propagateCalo:: Vertex " << tpVertex << " Momentum " << tpMomentum << " Charge " << tpCharge << " Radius " << radius << " Z " << zdist << " Corner " << corner << std::endl;
624  FreeTrajectoryState fts (tpVertex, tpMomentum, tpCharge, bField);
625 
628 
630 
631  AnalyticalPropagator myAP (bField, alongMomentum, 2*M_PI);
632 
634  if (tpMomentum.eta() < 0) {
635  tsose = myAP.propagate(fts, *lendcap);
636  } else {
637  tsose = myAP.propagate(fts, *rendcap);
638  }
639 
640  TrajectoryStateOnSurface tsosb = myAP.propagate(fts, *barrel);
641 
642  track.ok=true;
643  if (tsose.isValid() && tsosb.isValid()) {
644  float absEta = std::abs(tsosb.globalPosition().eta());
645  if (absEta < corner) {
646  track.point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
647  track.direction = tsosb.globalDirection();
648  } else {
649  track.point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
650  track.direction = tsose.globalDirection();
651  }
652  } else if (tsose.isValid()) {
653  track.point.SetXYZ(tsose.globalPosition().x(), tsose.globalPosition().y(), tsose.globalPosition().z());
654  track.direction = tsose.globalDirection();
655  } else if (tsosb.isValid()) {
656  track.point.SetXYZ(tsosb.globalPosition().x(), tsosb.globalPosition().y(), tsosb.globalPosition().z());
657  track.direction = tsosb.globalDirection();
658  } else {
659  track.point.SetXYZ(-999., -999., -999.);
660  track.direction = GlobalVector(0,0,1);
661  track.ok = false;
662  }
663  if (debug) {
664  std::cout << "propagateCalo:: Barrel " << tsosb.isValid() << " Endcap " << tsose.isValid() << " OverAll " << track.ok << " Point " << track.point << " Direction " << track.direction << std::endl;
665  if (track.ok) {
666  math::XYZPoint vDiff(track.point.x()-tpVertex.x(), track.point.y()-tpVertex.y(), track.point.z()-tpVertex.z());
667  double dphi = track.direction.phi()-tpMomentum.phi();
668  double rdist = std::sqrt(vDiff.x()*vDiff.x()+vDiff.y()*vDiff.y());
669  double pt = tpMomentum.perp();
670  double rat = 0.5*dphi/std::sin(0.5*dphi);
671  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;
672  }
673  }
674  return track;
675  }
676 
678 
679  spr::trackAtOrigin trk;
680 
681  edm::SimTrackContainer::const_iterator itr = SimTk->end();
682  for (edm::SimTrackContainer::const_iterator simTrkItr = SimTk->begin(); simTrkItr!= SimTk->end(); simTrkItr++) {
683  if ( simTrkItr->trackId() == thisTrk ) {
684  if (debug) std::cout << "matched trackId (maximum occurance) " << thisTrk << " type " << simTrkItr->type() << std::endl;
685  itr = simTrkItr;
686  break;
687  }
688  }
689 
690  if (itr != SimTk->end()) {
691  int vertIndex = itr->vertIndex();
692  if (vertIndex != -1 && vertIndex < (int)SimVtx->size()) {
693  edm::SimVertexContainer::const_iterator simVtxItr= SimVtx->begin();
694  for (int iv=0; iv<vertIndex; iv++) simVtxItr++;
695  const math::XYZTLorentzVectorD pos = simVtxItr->position();
696  const math::XYZTLorentzVectorD mom = itr->momentum();
697  trk.ok = true;
698  trk.charge = (int)(itr->charge());
699  trk.position = GlobalPoint(pos.x(), pos.y(), pos.z());
700  trk.momentum = GlobalVector(mom.x(), mom.y(), mom.z());
701  }
702  }
703  if (debug) std::cout << "Track flag " << trk.ok << " Position " << trk.position << " Momentum " << trk.momentum << std::endl;;
704  return trk;
705  }
706 
707 }
const double Z[kNumberCalorimeter]
virtual DetId getClosestCell(const GlobalPoint &r) const
double p() const
momentum vector magnitude
Definition: TrackBase.h:578
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
Definition: CaloGeometry.cc:43
virtual FreeTrajectoryState propagate(const FreeTrajectoryState &ftsStart, const GlobalPoint &pDest) const final
Definition: Propagator.h:119
int i
Definition: DBlmapReader.cc:9
static const TGPicture * info(bool iBackgroundIsBlack)
std::vector< spr::propagatedTrackID > propagateCALO(edm::Handle< reco::TrackCollection > &trkCollection, const CaloGeometry *geo, const MagneticField *bField, std::string &theTrackQuality, bool debug=false)
T perp() const
Definition: PV3DBase.h:72
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:14
virtual DetId getClosestCell(const GlobalPoint &r) const
reco::GenParticleCollection::const_iterator trkItr
TrackQuality
track quality
Definition: TrackBase.h:149
std::pair< math::XYZPoint, bool > propagateHCAL(const reco::Track *, const MagneticField *, bool debug=false)
uint32_t ID
Definition: Definitions.h:26
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
T y() const
Definition: PV3DBase.h:63
#define X(str)
Definition: MuonsGrabber.cc:48
GlobalPoint globalPosition() const
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:608
spr::trackAtOrigin simTrackAtOrigin(unsigned int thisTrk, edm::Handle< edm::SimTrackContainer > &SimTk, edm::Handle< edm::SimVertexContainer > &SimVtx, bool debug=false)
spr::propagatedTrack propagateCalo(const GlobalPoint &vertex, const GlobalVector &momentum, int charge, const MagneticField *, float zdist, float radius, float corner, bool debug=false)
double px() const
x coordinate of momentum vector
Definition: TrackBase.h:590
const math::XYZPoint & outerPosition() const
position of the outermost hit
Definition: Track.h:65
HepMC::GenEvent::particle_const_iterator trkItr
const math::XYZPoint & innerPosition() const
position of the innermost hit
Definition: Track.h:55
static CylinderPointer build(const PositionType &pos, const RotationType &rot, Scalar radius, Bounds *bounds=0)
Definition: Cylinder.h:51
reco::TrackCollection::const_iterator trkItr
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:614
T sqrt(T t)
Definition: SSEVec.h:48
T phi() const
Definition: Phi.h:41
T z() const
Definition: PV3DBase.h:64
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
tuple genEvent
Definition: MCTruth.py:33
std::pair< math::XYZPoint, bool > propagateTracker(const reco::Track *, const MagneticField *, bool debug=false)
std::pair< math::XYZPoint, bool > propagateECAL(const reco::Track *, const MagneticField *, bool debug=false)
spr::propagatedTrack propagateTrackToECAL(const reco::Track *, const MagneticField *, bool debug=false)
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:37
#define M_PI
double pz() const
z coordinate of momentum vector
Definition: TrackBase.h:602
virtual DetId getClosestCell(const GlobalPoint &r) const
double vz() const
z coordinate of the reference point on track
Definition: TrackBase.h:632
Definition: DetId.h:18
static TrackQuality qualityByName(const std::string &name)
Definition: TrackBase.cc:123
#define debug
Definition: HDRShower.cc:19
const math::XYZVector & outerMomentum() const
momentum vector at the outermost hit position
Definition: Track.h:70
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
bool quality(const TrackQuality) const
Track quality.
Definition: TrackBase.h:473
T eta() const
Definition: PV3DBase.h:76
double vy() const
y coordinate of the reference point on track
Definition: TrackBase.h:626
std::pair< math::XYZPoint, double > propagateTrackerEnd(const reco::Track *, const MagneticField *, bool debug=false)
const math::XYZVector & innerMomentum() const
momentum vector at the innermost hit position
Definition: Track.h:60
std::vector< spr::propagatedTrackID > propagateCosmicCALO(edm::Handle< reco::TrackCollection > &trkCollection, const CaloGeometry *geo, const MagneticField *bField, std::string &theTrackQuality, bool debug=false)
reco::TrackCollection::const_iterator trkItr
tuple cout
Definition: gather_cfg.py:121
int charge() const
track electric charge
Definition: TrackBase.h:530
T x() const
Definition: PV3DBase.h:62
spr::propagatedTrack propagateTrackToHCAL(const reco::Track *, const MagneticField *, bool debug=false)
*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
double py() const
y coordinate of momentum vector
Definition: TrackBase.h:596
double vx() const
x coordinate of the reference point on track
Definition: TrackBase.h:620
Global3DVector GlobalVector
Definition: GlobalVector.h:10
GlobalVector globalDirection() const