CMS 3D CMS Logo

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