CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
TrackDetectorAssociator.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: TrackAssociator
4 // Class: TrackDetectorAssociator
5 //
6 /*
7 
8  Description: <one line class summary>
9 
10  Implementation:
11  <Notes on implementation>
12 */
13 //
14 // Original Author: Dmytro Kovalskyi
15 // Created: Fri Apr 21 10:59:41 PDT 2006
16 //
17 //
18 
23 
27 
42 
43 // calorimeter and muon infos
52 
54 #include <stack>
55 #include <set>
56 
58 #include "Math/VectorUtil.h"
59 #include <algorithm>
60 
67 
68 using namespace reco;
69 
71  ivProp_ = nullptr;
72  useDefaultPropagator_ = false;
73 }
74 
76 
78  ivProp_ = ptr;
79  cachedTrajectory_.setPropagator(ivProp_);
80 }
81 
82 void TrackDetectorAssociator::useDefaultPropagator() { useDefaultPropagator_ = true; }
83 
85  // access the calorimeter geometry
86  theCaloGeometry_ = &iSetup.getData(parameters.theCaloGeometryToken);
87 
88  // get the tracking Geometry
89  theTrackingGeometry_ = &iSetup.getData(parameters.theTrackingGeometryToken);
90 
91  if (useDefaultPropagator_ && (!defProp_ || theMagneticFieldWatcher_.check(iSetup))) {
92  // setup propagator
93  const MagneticField* bField = &iSetup.getData(parameters.bFieldToken);
94 
95  auto prop = std::make_unique<SteppingHelixPropagator>(bField, anyDirection);
96  prop->setMaterialMode(false);
97  prop->applyRadX0Correction(true);
98  // prop->setDebug(true); // tmp
99  defProp_ = std::move(prop);
100  setPropagator(defProp_.get());
101  }
102 
103  ecalDetIdAssociator_ = &iSetup.getData(parameters.ecalDetIdAssociatorToken);
104  hcalDetIdAssociator_ = &iSetup.getData(parameters.hcalDetIdAssociatorToken);
105  hoDetIdAssociator_ = &iSetup.getData(parameters.hoDetIdAssociatorToken);
106  caloDetIdAssociator_ = &iSetup.getData(parameters.caloDetIdAssociatorToken);
107  muonDetIdAssociator_ = &iSetup.getData(parameters.muonDetIdAssociatorToken);
108  preshowerDetIdAssociator_ = &iSetup.getData(parameters.preshowerDetIdAssociatorToken);
109 }
110 
112  const edm::EventSetup& iSetup,
113  const FreeTrajectoryState& fts,
115  return associate(iEvent, iSetup, parameters, &fts);
116 }
117 
119  const edm::EventSetup& iSetup,
121  const FreeTrajectoryState* innerState,
122  const FreeTrajectoryState* outerState) {
124  if (!parameters.useEcal && !parameters.useCalo && !parameters.useHcal && !parameters.useHO && !parameters.useMuon &&
125  !parameters.usePreshower)
126  throw cms::Exception("ConfigurationError")
127  << "Configuration error! No subdetector was selected for the track association.";
128 
129  SteppingHelixStateInfo trackOrigin(*innerState);
130  info.stateAtIP = *innerState;
131  cachedTrajectory_.setStateAtIP(trackOrigin);
132 
133  init(iSetup, parameters);
134  // get track trajectory
135  // ECAL points (EB+EE)
136  // If the phi angle between a track entrance and exit points is more
137  // than 2 crystals, it is possible that the track will cross 3 crystals
138  // and therefore one has to check at least 3 points along the track
139  // trajectory inside ECAL. In order to have a chance to cross 4 crystalls
140  // in the barrel, a track should have P_t as low as 3 GeV or smaller
141  // If it's necessary, number of points along trajectory can be increased
142 
143  info.setCaloGeometry(theCaloGeometry_);
144 
145  cachedTrajectory_.reset_trajectory();
146  // estimate propagation outer boundaries based on
147  // requested sub-detector information. For now limit
148  // propagation region only if muon matching is not
149  // requested.
150  double HOmaxR = hoDetIdAssociator_->volume().maxR();
151  double HOmaxZ = hoDetIdAssociator_->volume().maxZ();
152  double minR = ecalDetIdAssociator_->volume().minR();
153  double minZ = preshowerDetIdAssociator_->volume().minZ();
154  cachedTrajectory_.setMaxHORadius(HOmaxR);
155  cachedTrajectory_.setMaxHOLength(HOmaxZ * 2.);
156  cachedTrajectory_.setMinDetectorRadius(minR);
157  cachedTrajectory_.setMinDetectorLength(minZ * 2.);
158 
159  if (parameters.useMuon) {
160  double maxR = muonDetIdAssociator_->volume().maxR();
161  double maxZ = muonDetIdAssociator_->volume().maxZ();
162  cachedTrajectory_.setMaxDetectorRadius(maxR);
163  cachedTrajectory_.setMaxDetectorLength(maxZ * 2.);
164  } else {
165  cachedTrajectory_.setMaxDetectorRadius(HOmaxR);
166  cachedTrajectory_.setMaxDetectorLength(HOmaxZ * 2.);
167  }
168 
169  // If track extras exist and outerState is before HO maximum, then use outerState
170  if (outerState) {
171  if (outerState->position().perp() < HOmaxR && std::abs(outerState->position().z()) < HOmaxZ) {
172  LogTrace("TrackAssociator") << "Using outerState as trackOrigin at Rho=" << outerState->position().perp()
173  << " Z=" << outerState->position().z() << "\n";
174  trackOrigin = SteppingHelixStateInfo(*outerState);
175  } else if (innerState) {
176  LogTrace("TrackAssociator") << "Using innerState as trackOrigin at Rho=" << innerState->position().perp()
177  << " Z=" << innerState->position().z() << "\n";
178  trackOrigin = SteppingHelixStateInfo(*innerState);
179  }
180  }
181 
182  if (trackOrigin.momentum().mag() == 0)
183  return info;
184  if (edm::isNotFinite(trackOrigin.momentum().x()) or edm::isNotFinite(trackOrigin.momentum().y()) or
185  edm::isNotFinite(trackOrigin.momentum().z()))
186  return info;
187  if (!cachedTrajectory_.propagateAll(trackOrigin))
188  return info;
189 
190  // get trajectory in calorimeters
191  cachedTrajectory_.findEcalTrajectory(ecalDetIdAssociator_->volume());
192  cachedTrajectory_.findHcalTrajectory(hcalDetIdAssociator_->volume());
193  cachedTrajectory_.findHOTrajectory(hoDetIdAssociator_->volume());
194  cachedTrajectory_.findPreshowerTrajectory(preshowerDetIdAssociator_->volume());
195 
196  info.trkGlobPosAtEcal = getPoint(cachedTrajectory_.getStateAtEcal().position());
197  info.trkGlobPosAtHcal = getPoint(cachedTrajectory_.getStateAtHcal().position());
198  info.trkGlobPosAtHO = getPoint(cachedTrajectory_.getStateAtHO().position());
199 
200  info.trkMomAtEcal = cachedTrajectory_.getStateAtEcal().momentum();
201  info.trkMomAtHcal = cachedTrajectory_.getStateAtHcal().momentum();
202  info.trkMomAtHO = cachedTrajectory_.getStateAtHO().momentum();
203 
204  if (parameters.useEcal)
205  fillEcal(iEvent, info, parameters);
206  if (parameters.useCalo)
207  fillCaloTowers(iEvent, info, parameters);
208  if (parameters.useHcal)
209  fillHcal(iEvent, info, parameters);
210  if (parameters.useHO)
211  fillHO(iEvent, info, parameters);
212  if (parameters.usePreshower)
213  fillPreshower(iEvent, info, parameters);
214  if (parameters.useMuon)
215  fillMuon(iEvent, info, parameters);
216  if (parameters.truthMatch)
217  fillCaloTruth(iEvent, info, parameters);
218 
219  return info;
220 }
221 
225  const std::vector<SteppingHelixStateInfo>& trajectoryStates = cachedTrajectory_.getEcalTrajectory();
226 
227  for (std::vector<SteppingHelixStateInfo>::const_iterator itr = trajectoryStates.begin();
228  itr != trajectoryStates.end();
229  itr++)
230  LogTrace("TrackAssociator") << "ECAL trajectory point (rho, z, phi): " << itr->position().perp() << ", "
231  << itr->position().z() << ", " << itr->position().phi();
232 
233  std::vector<GlobalPoint> coreTrajectory;
234  for (std::vector<SteppingHelixStateInfo>::const_iterator itr = trajectoryStates.begin();
235  itr != trajectoryStates.end();
236  itr++)
237  coreTrajectory.push_back(itr->position());
238 
239  if (coreTrajectory.empty()) {
240  LogTrace("TrackAssociator") << "ECAL track trajectory is empty; moving on\n";
241  info.isGoodEcal = false;
242  return;
243  }
244  info.isGoodEcal = true;
245 
246  // Find ECAL crystals
248  iEvent.getByToken(parameters.EBRecHitsToken, EBRecHits);
249  if (!EBRecHits.isValid())
250  throw cms::Exception("FatalError") << "Unable to find EBRecHitCollection in the event!\n";
251 
253  iEvent.getByToken(parameters.EERecHitsToken, EERecHits);
254  if (!EERecHits.isValid())
255  throw cms::Exception("FatalError") << "Unable to find EERecHitCollection in event!\n";
256 
257  std::set<DetId> ecalIdsInRegion;
258  if (parameters.accountForTrajectoryChangeCalo) {
259  // get trajectory change with respect to initial state
260  DetIdAssociator::MapRange mapRange =
261  getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::IpToEcal), parameters.dREcalPreselection);
262  ecalIdsInRegion = ecalDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], mapRange);
263  } else
264  ecalIdsInRegion = ecalDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], parameters.dREcalPreselection);
265  LogTrace("TrackAssociator") << "ECAL hits in the region: " << ecalIdsInRegion.size();
266  if (parameters.dREcalPreselection > parameters.dREcal)
267  ecalIdsInRegion = ecalDetIdAssociator_->getDetIdsInACone(ecalIdsInRegion, coreTrajectory, parameters.dREcal);
268  LogTrace("TrackAssociator") << "ECAL hits in the cone: " << ecalIdsInRegion.size();
269  info.crossedEcalIds = ecalDetIdAssociator_->getCrossedDetIds(ecalIdsInRegion, coreTrajectory);
270  const std::vector<DetId>& crossedEcalIds = info.crossedEcalIds;
271  LogTrace("TrackAssociator") << "ECAL crossed hits " << crossedEcalIds.size();
272 
273  // add EcalRecHits
274  for (std::vector<DetId>::const_iterator itr = crossedEcalIds.begin(); itr != crossedEcalIds.end(); itr++) {
275  std::vector<EcalRecHit>::const_iterator ebHit = (*EBRecHits).find(*itr);
276  std::vector<EcalRecHit>::const_iterator eeHit = (*EERecHits).find(*itr);
277  if (ebHit != (*EBRecHits).end())
278  info.crossedEcalRecHits.push_back(&*ebHit);
279  else if (eeHit != (*EERecHits).end())
280  info.crossedEcalRecHits.push_back(&*eeHit);
281  else
282  LogTrace("TrackAssociator") << "Crossed EcalRecHit is not found for DetId: " << itr->rawId();
283  }
284  for (std::set<DetId>::const_iterator itr = ecalIdsInRegion.begin(); itr != ecalIdsInRegion.end(); itr++) {
285  std::vector<EcalRecHit>::const_iterator ebHit = (*EBRecHits).find(*itr);
286  std::vector<EcalRecHit>::const_iterator eeHit = (*EERecHits).find(*itr);
287  if (ebHit != (*EBRecHits).end())
288  info.ecalRecHits.push_back(&*ebHit);
289  else if (eeHit != (*EERecHits).end())
290  info.ecalRecHits.push_back(&*eeHit);
291  else
292  LogTrace("TrackAssociator") << "EcalRecHit from the cone is not found for DetId: " << itr->rawId();
293  }
294 }
295 
299  // use ECAL and HCAL trajectories to match a tower. (HO isn't used for matching).
300  std::vector<GlobalPoint> trajectory;
301  const std::vector<SteppingHelixStateInfo>& ecalTrajectoryStates = cachedTrajectory_.getEcalTrajectory();
302  const std::vector<SteppingHelixStateInfo>& hcalTrajectoryStates = cachedTrajectory_.getHcalTrajectory();
303  for (std::vector<SteppingHelixStateInfo>::const_iterator itr = ecalTrajectoryStates.begin();
304  itr != ecalTrajectoryStates.end();
305  itr++)
306  trajectory.push_back(itr->position());
307  for (std::vector<SteppingHelixStateInfo>::const_iterator itr = hcalTrajectoryStates.begin();
308  itr != hcalTrajectoryStates.end();
309  itr++)
310  trajectory.push_back(itr->position());
311 
312  if (trajectory.empty()) {
313  LogTrace("TrackAssociator") << "HCAL trajectory is empty; moving on\n";
314  info.isGoodCalo = false;
315  return;
316  }
317  info.isGoodCalo = true;
318 
319  // find crossed CaloTowers
321  iEvent.getByToken(parameters.caloTowersToken, caloTowers);
322  if (!caloTowers.isValid())
323  throw cms::Exception("FatalError") << "Unable to find CaloTowers in event!\n";
324 
325  std::set<DetId> caloTowerIdsInRegion;
326  if (parameters.accountForTrajectoryChangeCalo) {
327  // get trajectory change with respect to initial state
328  DetIdAssociator::MapRange mapRange =
329  getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::IpToHcal), parameters.dRHcalPreselection);
330  caloTowerIdsInRegion = caloDetIdAssociator_->getDetIdsCloseToAPoint(trajectory[0], mapRange);
331  } else
332  caloTowerIdsInRegion = caloDetIdAssociator_->getDetIdsCloseToAPoint(trajectory[0], parameters.dRHcalPreselection);
333 
334  LogTrace("TrackAssociator") << "Towers in the region: " << caloTowerIdsInRegion.size();
335 
336  auto caloTowerIdsInAConeBegin = caloTowerIdsInRegion.begin();
337  auto caloTowerIdsInAConeEnd = caloTowerIdsInRegion.end();
338  std::set<DetId> caloTowerIdsInAConeTmp;
339  if (!caloDetIdAssociator_->selectAllInACone(parameters.dRHcal)) {
340  caloTowerIdsInAConeTmp =
341  caloDetIdAssociator_->getDetIdsInACone(caloTowerIdsInRegion, trajectory, parameters.dRHcal);
342  caloTowerIdsInAConeBegin = caloTowerIdsInAConeTmp.begin();
343  caloTowerIdsInAConeEnd = caloTowerIdsInAConeTmp.end();
344  }
345  LogTrace("TrackAssociator") << "Towers in the cone: "
346  << std::distance(caloTowerIdsInAConeBegin, caloTowerIdsInAConeEnd);
347 
348  info.crossedTowerIds = caloDetIdAssociator_->getCrossedDetIds(caloTowerIdsInRegion, trajectory);
349  const std::vector<DetId>& crossedCaloTowerIds = info.crossedTowerIds;
350  LogTrace("TrackAssociator") << "Towers crossed: " << crossedCaloTowerIds.size();
351 
352  // add CaloTowers
353  for (std::vector<DetId>::const_iterator itr = crossedCaloTowerIds.begin(); itr != crossedCaloTowerIds.end(); itr++) {
354  CaloTowerCollection::const_iterator tower = (*caloTowers).find(*itr);
355  if (tower != (*caloTowers).end())
356  info.crossedTowers.push_back(&*tower);
357  else
358  LogTrace("TrackAssociator") << "Crossed CaloTower is not found for DetId: " << (*itr).rawId();
359  }
360 
361  for (std::set<DetId>::const_iterator itr = caloTowerIdsInAConeBegin; itr != caloTowerIdsInAConeEnd; itr++) {
362  CaloTowerCollection::const_iterator tower = (*caloTowers).find(*itr);
363  if (tower != (*caloTowers).end())
364  info.towers.push_back(&*tower);
365  else
366  LogTrace("TrackAssociator") << "CaloTower from the cone is not found for DetId: " << (*itr).rawId();
367  }
368 }
369 
373  std::vector<GlobalPoint> trajectory;
374  const std::vector<SteppingHelixStateInfo>& trajectoryStates = cachedTrajectory_.getPreshowerTrajectory();
375  for (std::vector<SteppingHelixStateInfo>::const_iterator itr = trajectoryStates.begin();
376  itr != trajectoryStates.end();
377  itr++)
378  trajectory.push_back(itr->position());
379 
380  if (trajectory.empty()) {
381  LogTrace("TrackAssociator") << "Preshower trajectory is empty; moving on\n";
382  return;
383  }
384 
385  std::set<DetId> idsInRegion =
386  preshowerDetIdAssociator_->getDetIdsCloseToAPoint(trajectory[0], parameters.dRPreshowerPreselection);
387 
388  LogTrace("TrackAssociator") << "Number of Preshower Ids in the region: " << idsInRegion.size();
389  info.crossedPreshowerIds = preshowerDetIdAssociator_->getCrossedDetIds(idsInRegion, trajectory);
390  LogTrace("TrackAssociator") << "Number of Preshower Ids in crossed: " << info.crossedPreshowerIds.size();
391 }
392 
396  const std::vector<SteppingHelixStateInfo>& trajectoryStates = cachedTrajectory_.getHcalTrajectory();
397 
398  std::vector<GlobalPoint> coreTrajectory;
399  for (std::vector<SteppingHelixStateInfo>::const_iterator itr = trajectoryStates.begin();
400  itr != trajectoryStates.end();
401  itr++)
402  coreTrajectory.push_back(itr->position());
403 
404  if (coreTrajectory.empty()) {
405  LogTrace("TrackAssociator") << "HCAL trajectory is empty; moving on\n";
406  info.isGoodHcal = false;
407  return;
408  }
409  info.isGoodHcal = true;
410 
411  // find crossed Hcals
413  iEvent.getByToken(parameters.HBHEcollToken, collection);
414  if (!collection.isValid())
415  throw cms::Exception("FatalError") << "Unable to find HBHERecHits in event!\n";
416 
417  std::set<DetId> idsInRegion;
418  if (parameters.accountForTrajectoryChangeCalo) {
419  // get trajectory change with respect to initial state
420  DetIdAssociator::MapRange mapRange =
421  getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::IpToHcal), parameters.dRHcalPreselection);
422  idsInRegion = hcalDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], mapRange);
423  } else
424  idsInRegion = hcalDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], parameters.dRHcalPreselection);
425 
426  LogTrace("TrackAssociator") << "HCAL hits in the region: " << idsInRegion.size() << "\n"
427  << DetIdInfo::info(idsInRegion, nullptr);
428 
429  auto idsInAConeBegin = idsInRegion.begin();
430  auto idsInAConeEnd = idsInRegion.end();
431  std::set<DetId> idsInAConeTmp;
432  if (!hcalDetIdAssociator_->selectAllInACone(parameters.dRHcal)) {
433  idsInAConeTmp = hcalDetIdAssociator_->getDetIdsInACone(idsInRegion, coreTrajectory, parameters.dRHcal);
434  idsInAConeBegin = idsInAConeTmp.begin();
435  idsInAConeEnd = idsInAConeTmp.end();
436  }
437  LogTrace("TrackAssociator") << "HCAL hits in the cone: " << std::distance(idsInAConeBegin, idsInAConeEnd) << "\n"
438  << DetIdInfo::info(std::set<DetId>(idsInAConeBegin, idsInAConeEnd), nullptr);
439  info.crossedHcalIds = hcalDetIdAssociator_->getCrossedDetIds(idsInRegion, coreTrajectory);
440  const std::vector<DetId>& crossedIds = info.crossedHcalIds;
441  LogTrace("TrackAssociator") << "HCAL hits crossed: " << crossedIds.size() << "\n"
442  << DetIdInfo::info(crossedIds, nullptr);
443 
444  // add Hcal
445  for (std::vector<DetId>::const_iterator itr = crossedIds.begin(); itr != crossedIds.end(); itr++) {
446  HBHERecHitCollection::const_iterator hit = (*collection).find(*itr);
447  if (hit != (*collection).end())
448  info.crossedHcalRecHits.push_back(&*hit);
449  else
450  LogTrace("TrackAssociator") << "Crossed HBHERecHit is not found for DetId: " << itr->rawId();
451  }
452  for (std::set<DetId>::const_iterator itr = idsInAConeBegin; itr != idsInAConeEnd; itr++) {
453  HBHERecHitCollection::const_iterator hit = (*collection).find(*itr);
454  if (hit != (*collection).end())
455  info.hcalRecHits.push_back(&*hit);
456  else
457  LogTrace("TrackAssociator") << "HBHERecHit from the cone is not found for DetId: " << itr->rawId();
458  }
459 }
460 
464  const std::vector<SteppingHelixStateInfo>& trajectoryStates = cachedTrajectory_.getHOTrajectory();
465 
466  std::vector<GlobalPoint> coreTrajectory;
467  for (std::vector<SteppingHelixStateInfo>::const_iterator itr = trajectoryStates.begin();
468  itr != trajectoryStates.end();
469  itr++)
470  coreTrajectory.push_back(itr->position());
471 
472  if (coreTrajectory.empty()) {
473  LogTrace("TrackAssociator") << "HO trajectory is empty; moving on\n";
474  info.isGoodHO = false;
475  return;
476  }
477  info.isGoodHO = true;
478 
479  // find crossed HOs
481  iEvent.getByToken(parameters.HOcollToken, collection);
482  if (!collection.isValid())
483  throw cms::Exception("FatalError") << "Unable to find HORecHits in event!\n";
484 
485  std::set<DetId> idsInRegion;
486  if (parameters.accountForTrajectoryChangeCalo) {
487  // get trajectory change with respect to initial state
488  DetIdAssociator::MapRange mapRange =
489  getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::IpToHO), parameters.dRHcalPreselection);
490  idsInRegion = hoDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], mapRange);
491  } else
492  idsInRegion = hoDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], parameters.dRHcalPreselection);
493 
494  LogTrace("TrackAssociator") << "idsInRegion.size(): " << idsInRegion.size();
495 
496  auto idsInAConeBegin = idsInRegion.begin();
497  auto idsInAConeEnd = idsInRegion.end();
498  std::set<DetId> idsInAConeTmp;
499  if (!hoDetIdAssociator_->selectAllInACone(parameters.dRHcal)) {
500  idsInAConeTmp = hoDetIdAssociator_->getDetIdsInACone(idsInRegion, coreTrajectory, parameters.dRHcal);
501  idsInAConeBegin = idsInAConeTmp.begin();
502  idsInAConeEnd = idsInAConeTmp.end();
503  }
504  LogTrace("TrackAssociator") << "idsInACone.size(): " << std::distance(idsInAConeBegin, idsInAConeEnd);
505  info.crossedHOIds = hoDetIdAssociator_->getCrossedDetIds(idsInRegion, coreTrajectory);
506  const std::vector<DetId>& crossedIds = info.crossedHOIds;
507 
508  // add HO
509  for (std::vector<DetId>::const_iterator itr = crossedIds.begin(); itr != crossedIds.end(); itr++) {
510  HORecHitCollection::const_iterator hit = (*collection).find(*itr);
511  if (hit != (*collection).end())
512  info.crossedHORecHits.push_back(&*hit);
513  else
514  LogTrace("TrackAssociator") << "Crossed HORecHit is not found for DetId: " << itr->rawId();
515  }
516 
517  for (std::set<DetId>::const_iterator itr = idsInAConeBegin; itr != idsInAConeEnd; itr++) {
518  HORecHitCollection::const_iterator hit = (*collection).find(*itr);
519  if (hit != (*collection).end())
520  info.hoRecHits.push_back(&*hit);
521  else
522  LogTrace("TrackAssociator") << "HORecHit from the cone is not found for DetId: " << itr->rawId();
523  }
524 }
525 
527  const SimTrack& track,
528  const SimVertex& vertex) {
529  GlobalVector vector(track.momentum().x(), track.momentum().y(), track.momentum().z());
530  GlobalPoint point(vertex.position().x(), vertex.position().y(), vertex.position().z());
531 
532  int charge = track.type() > 0 ? -1 : 1; // lepton convention
533  if (abs(track.type()) == 211 || // pion
534  abs(track.type()) == 321 || // kaon
535  abs(track.type()) == 2212)
536  charge = track.type() < 0 ? -1 : 1;
537  return getFreeTrajectoryState(bField, vector, point, charge);
538 }
539 
541  const GlobalVector& momentum,
542  const GlobalPoint& vertex,
543  const int charge) {
544  GlobalTrajectoryParameters tPars(vertex, momentum, charge, bField);
545 
546  ROOT::Math::SMatrixIdentity id;
547  AlgebraicSymMatrix66 covT(id);
548  covT *= 1e-6; // initialize to sigma=1e-3
549  CartesianTrajectoryError tCov(covT);
550 
551  return FreeTrajectoryState(tPars, tCov);
552 }
553 
555  const reco::Track& track) {
556  GlobalVector vector(track.momentum().x(), track.momentum().y(), track.momentum().z());
557 
558  GlobalPoint point(track.vertex().x(), track.vertex().y(), track.vertex().z());
559 
561 
562  // FIX THIS !!!
563  // need to convert from perigee to global or helix (curvilinear) frame
564  // for now just an arbitrary matrix.
565  ROOT::Math::SMatrixIdentity id;
566  AlgebraicSymMatrix66 covT(id);
567  covT *= 1e-6; // initialize to sigma=1e-3
568  CartesianTrajectoryError tCov(covT);
569 
570  return FreeTrajectoryState(tPars, tCov);
571 }
572 
573 DetIdAssociator::MapRange TrackDetectorAssociator::getMapRange(const std::pair<float, float>& delta, const float dR) {
574  DetIdAssociator::MapRange mapRange;
575  mapRange.dThetaPlus = dR;
576  mapRange.dThetaMinus = dR;
577  mapRange.dPhiPlus = dR;
578  mapRange.dPhiMinus = dR;
579  if (delta.first > 0)
580  mapRange.dThetaPlus += delta.first;
581  else
582  mapRange.dThetaMinus += std::abs(delta.first);
583  if (delta.second > 0)
584  mapRange.dPhiPlus += delta.second;
585  else
586  mapRange.dPhiMinus += std::abs(delta.second);
587  LogTrace("TrackAssociator") << "Selection range: (dThetaPlus, dThetaMinus, dPhiPlus, dPhiMinus, dRPreselection): "
588  << mapRange.dThetaPlus << ", " << mapRange.dThetaMinus << ", " << mapRange.dPhiPlus
589  << ", " << mapRange.dPhiMinus << ", " << dR;
590  return mapRange;
591 }
592 
593 void TrackDetectorAssociator::getTAMuonChamberMatches(std::vector<TAMuonChamberMatch>& matches,
595  // Strategy:
596  // Propagate through the whole detector, estimate change in eta and phi
597  // along the trajectory, add this to dRMuon and find DetIds around this
598  // direction using the map. Then propagate fast to each surface and apply
599  // final matching criteria.
600 
601  // get the direction first
602  SteppingHelixStateInfo trajectoryPoint = cachedTrajectory_.getStateAtHcal();
603  // If trajectory point at HCAL is not valid, try to use the outer most state of the
604  // trajectory instead.
605  if (!trajectoryPoint.isValid())
606  trajectoryPoint = cachedTrajectory_.getOuterState();
607  if (!trajectoryPoint.isValid()) {
608  LogTrace("TrackAssociator")
609  << "trajectory position at HCAL is not valid. Assume the track cannot reach muon detectors and skip it";
610  return;
611  }
612 
613  GlobalVector direction = trajectoryPoint.momentum().unit();
614  LogTrace("TrackAssociator") << "muon direction: " << direction
615  << "\n\t and corresponding point: " << trajectoryPoint.position() << "\n";
616 
617  DetIdAssociator::MapRange mapRange =
618  getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::FullTrajectory), parameters.dRMuonPreselection);
619 
620  // and find chamber DetIds
621 
622  std::set<DetId> muonIdsInRegion = muonDetIdAssociator_->getDetIdsCloseToAPoint(trajectoryPoint.position(), mapRange);
623  LogTrace("TrackAssociator") << "Number of chambers to check: " << muonIdsInRegion.size();
624  for (std::set<DetId>::const_iterator detId = muonIdsInRegion.begin(); detId != muonIdsInRegion.end(); detId++) {
625  const GeomDet* geomDet = muonDetIdAssociator_->getGeomDet(*detId);
626  TrajectoryStateOnSurface stateOnSurface = cachedTrajectory_.propagate(&geomDet->surface());
627  if (!stateOnSurface.isValid()) {
628  LogTrace("TrackAssociator") << "Failed to propagate the track; moving on\n\t"
629  << "Element is not crosssed: " << DetIdInfo::info(*detId, nullptr) << "\n";
630  continue;
631  }
632  LocalPoint localPoint = geomDet->surface().toLocal(stateOnSurface.freeState()->position());
633  LocalError localError = stateOnSurface.localError().positionError();
634  float distanceX = 0.f;
635  float distanceY = 0.f;
636  if (const CSCChamber* cscChamber = dynamic_cast<const CSCChamber*>(geomDet)) {
637  const CSCChamberSpecs* chamberSpecs = cscChamber->specs();
638  if (!chamberSpecs) {
639  LogTrace("TrackAssociator") << "Failed to get CSCChamberSpecs from CSCChamber; moving on\n";
640  continue;
641  }
642  const CSCLayerGeometry* layerGeometry = chamberSpecs->oddLayerGeometry(1);
643  if (!layerGeometry) {
644  LogTrace("TrackAssociator") << "Failed to get CSCLayerGeometry from CSCChamberSpecs; moving on\n";
645  continue;
646  }
647  const CSCWireTopology* wireTopology = layerGeometry->wireTopology();
648  if (!wireTopology) {
649  LogTrace("TrackAssociator") << "Failed to get CSCWireTopology from CSCLayerGeometry; moving on\n";
650  continue;
651  }
652 
653  float wideWidth = wireTopology->wideWidthOfPlane();
654  float narrowWidth = wireTopology->narrowWidthOfPlane();
655  float length = wireTopology->lengthOfPlane();
656  // If slanted, there is no y offset between local origin and symmetry center of wire plane
657  float yOfFirstWire = std::abs(wireTopology->wireAngle()) > 1.E-06f ? -0.5 * length : wireTopology->yOfWire(1);
658  // y offset between local origin and symmetry center of wire plane
659  float yCOWPOffset = yOfFirstWire + 0.5f * length;
660 
661  // tangent of the incline angle from inside the trapezoid
662  float tangent = (wideWidth - narrowWidth) / (2.f * length);
663  // y position wrt bottom of trapezoid
664  float yPrime = localPoint.y() + std::abs(yOfFirstWire);
665  // half trapezoid width at y' is 0.5 * narrowWidth + x side of triangle with the above tangent and side y'
666  float halfWidthAtYPrime = 0.5f * narrowWidth + yPrime * tangent;
667  distanceX = std::abs(localPoint.x()) - halfWidthAtYPrime;
668  distanceY = std::abs(localPoint.y() - yCOWPOffset) - 0.5f * length;
669  } else if (dynamic_cast<const GEMChamber*>(geomDet) || dynamic_cast<const GEMSuperChamber*>(geomDet)) {
670  const TrapezoidalPlaneBounds* bounds = dynamic_cast<const TrapezoidalPlaneBounds*>(&geomDet->surface().bounds());
671 
672  float wideWidth = bounds->width();
673  float narrowWidth = 2.f * bounds->widthAtHalfLength() - wideWidth;
674  float length = bounds->length();
675  float tangent = (wideWidth - narrowWidth) / (2.f * length);
676  float halfWidthAtY = tangent * localPoint.y() + 0.25f * (narrowWidth + wideWidth);
677 
678  distanceX = std::abs(localPoint.x()) - halfWidthAtY;
679  distanceY = std::abs(localPoint.y()) - 0.5f * length;
680  } else {
681  distanceX = std::abs(localPoint.x()) - 0.5f * geomDet->surface().bounds().width();
682  distanceY = std::abs(localPoint.y()) - 0.5f * geomDet->surface().bounds().length();
683  }
684  if ((distanceX < parameters.muonMaxDistanceX && distanceY < parameters.muonMaxDistanceY) ||
685  (distanceX * distanceX <
686  localError.xx() * parameters.muonMaxDistanceSigmaX * parameters.muonMaxDistanceSigmaX &&
687  distanceY * distanceY <
688  localError.yy() * parameters.muonMaxDistanceSigmaY * parameters.muonMaxDistanceSigmaY)) {
689  LogTrace("TrackAssociator") << "found a match: " << DetIdInfo::info(*detId, nullptr) << "\n";
690  TAMuonChamberMatch match;
691  match.tState = stateOnSurface;
692  match.localDistanceX = distanceX;
693  match.localDistanceY = distanceY;
694  match.id = *detId;
695  matches.push_back(match);
696  } else {
697  LogTrace("TrackAssociator") << "chamber is too far: " << DetIdInfo::info(*detId, nullptr)
698  << "\n\tdistanceX: " << distanceX << "\t distanceY: " << distanceY
699  << "\t sigmaX: " << distanceX / sqrt(localError.xx())
700  << "\t sigmaY: " << distanceY / sqrt(localError.yy()) << "\n";
701  }
702  }
703 }
704 
708  // Get the segments from the event
710  iEvent.getByToken(parameters.dtSegmentsToken, dtSegments);
711  if (!dtSegments.isValid())
712  throw cms::Exception("FatalError") << "Unable to find DTRecSegment4DCollection in event!\n";
713 
715  iEvent.getByToken(parameters.cscSegmentsToken, cscSegments);
716  if (!cscSegments.isValid())
717  throw cms::Exception("FatalError") << "Unable to find CSCSegmentCollection in event!\n";
718 
720  if (parameters.useGEM)
721  iEvent.getByToken(parameters.gemSegmentsToken, gemSegments);
723  if (parameters.useME0)
724  iEvent.getByToken(parameters.me0SegmentsToken, me0Segments);
725 
727 
728  // check the map of available segments
729  // if there is no segments in a given direction at all,
730  // then there is no point to fly there.
731  //
732  // MISSING
733  // Possible solution: quick search for presence of segments
734  // for the set of DetIds
735 
736  // get a set of matches corresponding to muon chambers
737  std::vector<TAMuonChamberMatch> matchedChambers;
738  LogTrace("TrackAssociator") << "Trying to Get ChamberMatches" << std::endl;
739  getTAMuonChamberMatches(matchedChambers, parameters);
740  LogTrace("TrackAssociator") << "Chambers matched: " << matchedChambers.size() << "\n";
741 
742  // Iterate over all chamber matches and fill segment matching
743  // info if it's available
744  for (std::vector<TAMuonChamberMatch>::iterator matchedChamber = matchedChambers.begin();
745  matchedChamber != matchedChambers.end();
746  matchedChamber++) {
747  const GeomDet* geomDet = muonDetIdAssociator_->getGeomDet((*matchedChamber).id);
748  // DT chamber
749  if (const DTChamber* chamber = dynamic_cast<const DTChamber*>(geomDet)) {
750  // Get the range for the corresponding segments
751  DTRecSegment4DCollection::range range = dtSegments->get(chamber->id());
752  // Loop over the segments of this chamber
753  for (DTRecSegment4DCollection::const_iterator segment = range.first; segment != range.second; segment++) {
754  if (addTAMuonSegmentMatch(*matchedChamber, &(*segment), parameters)) {
755  matchedChamber->segments.back().dtSegmentRef = DTRecSegment4DRef(dtSegments, segment - dtSegments->begin());
756  }
757  }
758  }
759  // CSC Chamber
760  else if (const CSCChamber* chamber = dynamic_cast<const CSCChamber*>(geomDet)) {
761  // Get the range for the corresponding segments
763  // Loop over the segments
764  for (CSCSegmentCollection::const_iterator segment = range.first; segment != range.second; segment++) {
765  if (addTAMuonSegmentMatch(*matchedChamber, &(*segment), parameters)) {
766  matchedChamber->segments.back().cscSegmentRef = CSCSegmentRef(cscSegments, segment - cscSegments->begin());
767  }
768  }
769  } else {
770  // GEM Chamber
771  if (parameters.useGEM) {
772  if (const GEMSuperChamber* chamber = dynamic_cast<const GEMSuperChamber*>(geomDet)) {
773  // Get the range for the corresponding segments
775  // Loop over the segments
776  for (GEMSegmentCollection::const_iterator segment = range.first; segment != range.second; segment++) {
777  if (addTAMuonSegmentMatch(*matchedChamber, &(*segment), parameters)) {
778  matchedChamber->segments.back().gemSegmentRef =
779  GEMSegmentRef(gemSegments, segment - gemSegments->begin());
780  }
781  }
782  }
783  }
784  // ME0 Chamber
785  if (parameters.useME0) {
786  if (const ME0Chamber* chamber = dynamic_cast<const ME0Chamber*>(geomDet)) {
787  // Get the range for the corresponding segments
789  // Loop over the segments
790  for (ME0SegmentCollection::const_iterator segment = range.first; segment != range.second; segment++) {
791  if (addTAMuonSegmentMatch(*matchedChamber, &(*segment), parameters)) {
792  matchedChamber->segments.back().me0SegmentRef =
793  ME0SegmentRef(me0Segments, segment - me0Segments->begin());
794  }
795  }
796  }
797  }
798  }
799  info.chambers.push_back(*matchedChamber);
800  }
801 }
802 
804  const RecSegment* segment,
806  LogTrace("TrackAssociator") << "Segment local position: " << segment->localPosition() << "\n"
807  << std::hex << segment->geographicalId().rawId() << "\n";
808 
809  const GeomDet* chamber = muonDetIdAssociator_->getGeomDet(matchedChamber.id);
810  TrajectoryStateOnSurface trajectoryStateOnSurface = matchedChamber.tState;
811  GlobalPoint segmentGlobalPosition = chamber->toGlobal(segment->localPosition());
812 
813  LogTrace("TrackAssociator") << "Segment global position: " << segmentGlobalPosition
814  << " \t (R_xy,eta,phi): " << segmentGlobalPosition.perp() << ","
815  << segmentGlobalPosition.eta() << "," << segmentGlobalPosition.phi() << "\n";
816 
817  LogTrace("TrackAssociator") << "\teta hit: " << segmentGlobalPosition.eta()
818  << " \tpropagator: " << trajectoryStateOnSurface.freeState()->position().eta() << "\n"
819  << "\tphi hit: " << segmentGlobalPosition.phi()
820  << " \tpropagator: " << trajectoryStateOnSurface.freeState()->position().phi()
821  << std::endl;
822 
823  bool isGood = false;
824  bool isDTWithoutY = false;
825  const DTRecSegment4D* dtseg = dynamic_cast<const DTRecSegment4D*>(segment);
826  if (dtseg && (!dtseg->hasZed()))
827  isDTWithoutY = true;
828 
829  float deltaPhi(std::abs(segmentGlobalPosition.phi() - trajectoryStateOnSurface.freeState()->position().phi()));
830  if (deltaPhi > float(M_PI))
831  deltaPhi = std::abs(deltaPhi - float(M_PI) * 2.f);
832  float deltaEta = std::abs(segmentGlobalPosition.eta() - trajectoryStateOnSurface.freeState()->position().eta());
833 
834  if (isDTWithoutY) {
835  isGood = deltaPhi < parameters.dRMuon;
836  // Be in chamber
837  isGood &= deltaEta < .3f;
838  } else
839  isGood = deltaEta * deltaEta + deltaPhi * deltaPhi < parameters.dRMuon * parameters.dRMuon;
840 
841  if (isGood) {
842  TAMuonSegmentMatch muonSegment;
843  muonSegment.segmentGlobalPosition = getPoint(segmentGlobalPosition);
844  muonSegment.segmentLocalPosition = getPoint(segment->localPosition());
845  muonSegment.segmentLocalDirection = getVector(segment->localDirection());
846  muonSegment.segmentLocalErrorXX = segment->localPositionError().xx();
847  muonSegment.segmentLocalErrorYY = segment->localPositionError().yy();
848  muonSegment.segmentLocalErrorXY = segment->localPositionError().xy();
849  muonSegment.segmentLocalErrorDxDz = segment->localDirectionError().xx();
850  muonSegment.segmentLocalErrorDyDz = segment->localDirectionError().yy();
851 
852  // DANGEROUS - compiler cannot guaranty parameters ordering
853  // AlgebraicSymMatrix segmentCovMatrix = segment->parametersError();
854  // muonSegment.segmentLocalErrorXDxDz = segmentCovMatrix[2][0];
855  // muonSegment.segmentLocalErrorYDyDz = segmentCovMatrix[3][1];
856  muonSegment.segmentLocalErrorXDxDz = -999.f;
857  muonSegment.segmentLocalErrorYDyDz = -999.f;
858  muonSegment.hasZed = true;
859  muonSegment.hasPhi = true;
860 
861  // timing information
862  muonSegment.t0 = 0.f;
863  if (dtseg) {
864  if ((dtseg->hasPhi()) && (!isDTWithoutY)) {
865  int phiHits = dtseg->phiSegment()->specificRecHits().size();
866  // int zHits = dtseg->zSegment()->specificRecHits().size();
867  int hits = 0;
868  double t0 = 0.;
869  // TODO: cuts on hit numbers not optimized in any way yet...
870  if (phiHits > 5 && dtseg->phiSegment()->ist0Valid()) {
871  t0 += dtseg->phiSegment()->t0() * phiHits;
872  hits += phiHits;
873  LogTrace("TrackAssociator") << " Phi t0: " << dtseg->phiSegment()->t0() << " hits: " << phiHits;
874  }
875  // the z segments seem to contain little useful information...
876  // if (zHits>3) {
877  // t0+=s->zSegment()->t0()*zHits;
878  // hits+=zHits;
879  // LogTrace("TrackAssociator") << " Z t0: " << s->zSegment()->t0() << " hits: " << zHits << std::endl;
880  // }
881  if (hits)
882  muonSegment.t0 = t0 / hits;
883  // LogTrace("TrackAssociator") << " --- t0: " << muonSegment.t0 << std::endl;
884  } else {
885  // check and set dimensionality
886  if (isDTWithoutY)
887  muonSegment.hasZed = false;
888  if (!dtseg->hasPhi())
889  muonSegment.hasPhi = false;
890  }
891  }
892  matchedChamber.segments.push_back(muonSegment);
893  }
894 
895  return isGood;
896 }
897 
898 //********************** NON-CORE CODE ******************************//
899 
903  // get list of simulated tracks and their vertices
904  using namespace edm;
906  iEvent.getByToken(parameters.simTracksToken, simTracks);
907  if (!simTracks.isValid())
908  throw cms::Exception("FatalError") << "No simulated tracks found\n";
909 
911  iEvent.getByToken(parameters.simVerticesToken, simVertices);
912  if (!simVertices.isValid())
913  throw cms::Exception("FatalError") << "No simulated vertices found\n";
914 
915  // get sim calo hits
916  Handle<PCaloHitContainer> simEcalHitsEB;
917  iEvent.getByToken(parameters.simEcalHitsEBToken, simEcalHitsEB);
918  if (!simEcalHitsEB.isValid())
919  throw cms::Exception("FatalError") << "No simulated ECAL EB hits found\n";
920 
921  Handle<PCaloHitContainer> simEcalHitsEE;
922  iEvent.getByToken(parameters.simEcalHitsEEToken, simEcalHitsEE);
923  if (!simEcalHitsEE.isValid())
924  throw cms::Exception("FatalError") << "No simulated ECAL EE hits found\n";
925 
926  Handle<PCaloHitContainer> simHcalHits;
927  iEvent.getByToken(parameters.simHcalHitsToken, simHcalHits);
928  if (!simHcalHits.isValid())
929  throw cms::Exception("FatalError") << "No simulated HCAL hits found\n";
930 
931  // find truth partner
932  SimTrackContainer::const_iterator simTrack = simTracks->begin();
933  for (; simTrack != simTracks->end(); ++simTrack) {
934  math::XYZVector simP3(simTrack->momentum().x(), simTrack->momentum().y(), simTrack->momentum().z());
935  math::XYZVector recoP3(info.stateAtIP.momentum().x(), info.stateAtIP.momentum().y(), info.stateAtIP.momentum().z());
936  if (ROOT::Math::VectorUtil::DeltaR(recoP3, simP3) < 0.1)
937  break;
938  }
939  if (simTrack != simTracks->end()) {
940  info.simTrack = &(*simTrack);
941  float ecalTrueEnergy(0);
942  float hcalTrueEnergy(0);
943 
944  // loop over calo hits
945  for (PCaloHitContainer::const_iterator hit = simEcalHitsEB->begin(); hit != simEcalHitsEB->end(); ++hit)
946  if (hit->geantTrackId() == info.simTrack->genpartIndex())
947  ecalTrueEnergy += hit->energy();
948 
949  for (PCaloHitContainer::const_iterator hit = simEcalHitsEE->begin(); hit != simEcalHitsEE->end(); ++hit)
950  if (hit->geantTrackId() == info.simTrack->genpartIndex())
951  ecalTrueEnergy += hit->energy();
952 
953  for (PCaloHitContainer::const_iterator hit = simHcalHits->begin(); hit != simHcalHits->end(); ++hit)
954  if (hit->geantTrackId() == info.simTrack->genpartIndex())
955  hcalTrueEnergy += hit->energy();
956 
957  info.ecalTrueEnergy = ecalTrueEnergy;
958  info.hcalTrueEnergy = hcalTrueEnergy;
959  info.hcalTrueEnergyCorrected = hcalTrueEnergy;
960  if (std::abs(info.trkGlobPosAtHcal.eta()) < 1.3f)
961  info.hcalTrueEnergyCorrected = hcalTrueEnergy * 113.2f;
962  else if (std::abs(info.trkGlobPosAtHcal.eta()) < 3.0f)
963  info.hcalTrueEnergyCorrected = hcalTrueEnergy * 167.2f;
964  }
965 }
966 
968  const edm::EventSetup& iSetup,
969  const reco::Track& track,
971  Direction direction /*= Any*/) {
972  double currentStepSize = cachedTrajectory_.getPropagationStep();
973 
974  const MagneticField* bField = &iSetup.getData(parameters.bFieldToken);
975 
976  if (track.extra().isNull()) {
977  if (direction != InsideOut)
978  throw cms::Exception("FatalError") << "No TrackExtra information is available and association is done with "
979  "something else than InsideOut track.\n"
980  << "Either change the parameter or provide needed data!\n";
981  LogTrace("TrackAssociator") << "Track Extras not found\n";
983  return associate(iEvent, iSetup, parameters, &initialState); // 5th argument is null pointer
984  }
985 
986  LogTrace("TrackAssociator") << "Track Extras found\n";
990 
991  LogTrace("TrackAssociator") << "inner track state (rho, z, phi):" << track.innerPosition().Rho() << ", "
992  << track.innerPosition().z() << ", " << track.innerPosition().phi() << "\n";
993  LogTrace("TrackAssociator") << "innerFreeState (rho, z, phi):" << innerState.position().perp() << ", "
994  << innerState.position().z() << ", " << innerState.position().phi() << "\n";
995 
996  LogTrace("TrackAssociator") << "outer track state (rho, z, phi):" << track.outerPosition().Rho() << ", "
997  << track.outerPosition().z() << ", " << track.outerPosition().phi() << "\n";
998  LogTrace("TrackAssociator") << "outerFreeState (rho, z, phi):" << outerState.position().perp() << ", "
999  << outerState.position().z() << ", " << outerState.position().phi() << "\n";
1000 
1001  // InsideOut first
1002  if (crossedIP(track)) {
1003  switch (direction) {
1004  case InsideOut:
1005  case Any:
1006  return associate(iEvent, iSetup, parameters, &referenceState, &outerState);
1007  break;
1008  case OutsideIn: {
1009  cachedTrajectory_.setPropagationStep(-std::abs(currentStepSize));
1010  TrackDetMatchInfo result = associate(iEvent, iSetup, parameters, &innerState, &referenceState);
1011  cachedTrajectory_.setPropagationStep(currentStepSize);
1012  return result;
1013  break;
1014  }
1015  }
1016  } else {
1017  switch (direction) {
1018  case InsideOut:
1019  return associate(iEvent, iSetup, parameters, &innerState, &outerState);
1020  break;
1021  case OutsideIn: {
1022  cachedTrajectory_.setPropagationStep(-std::abs(currentStepSize));
1023  TrackDetMatchInfo result = associate(iEvent, iSetup, parameters, &outerState, &innerState);
1024  cachedTrajectory_.setPropagationStep(currentStepSize);
1025  return result;
1026  break;
1027  }
1028  case Any: {
1029  // check if we deal with clear outside-in case
1030  if (track.innerPosition().Dot(track.innerMomentum()) < 0 &&
1031  track.outerPosition().Dot(track.outerMomentum()) < 0) {
1032  cachedTrajectory_.setPropagationStep(-std::abs(currentStepSize));
1034  if (track.innerPosition().Mag2() < track.outerPosition().Mag2())
1035  result = associate(iEvent, iSetup, parameters, &innerState, &outerState);
1036  else
1037  result = associate(iEvent, iSetup, parameters, &outerState, &innerState);
1038  cachedTrajectory_.setPropagationStep(currentStepSize);
1039  return result;
1040  }
1041  }
1042  }
1043  }
1044 
1045  // all other cases
1046  return associate(iEvent, iSetup, parameters, &innerState, &outerState);
1047 }
1048 
1050  const edm::EventSetup& iSetup,
1051  const SimTrack& track,
1052  const SimVertex& vertex,
1054  auto const* bField = &iSetup.getData(parameters.bFieldToken);
1055  return associate(iEvent, iSetup, getFreeTrajectoryState(bField, track, vertex), parameters);
1056 }
1057 
1059  const edm::EventSetup& iSetup,
1060  const GlobalVector& momentum,
1061  const GlobalPoint& vertex,
1062  const int charge,
1064  auto const* bField = &iSetup.getData(parameters.bFieldToken);
1065  return associate(iEvent, iSetup, getFreeTrajectoryState(bField, momentum, vertex, charge), parameters);
1066 }
1067 
1069  bool crossed = true;
1070  crossed &= (track.innerPosition().rho() > 3); // something close to active volume
1071  crossed &= (track.outerPosition().rho() > 3); // something close to active volume
1072  crossed &=
1073  ((track.innerPosition().x() * track.innerMomentum().x() + track.innerPosition().y() * track.innerMomentum().y() <
1074  0) !=
1075  (track.outerPosition().x() * track.outerMomentum().x() + track.outerPosition().y() * track.outerMomentum().y() <
1076  0));
1077  return crossed;
1078 }
constexpr double deltaPhi(double phi1, double phi2)
Definition: deltaPhi.h:26
bool ist0Valid() const
static const TGPicture * info(bool iBackgroundIsBlack)
T perp() const
Definition: PV3DBase.h:69
std::pair< const_iterator, const_iterator > range
iterator range
Definition: RangeMap.h:50
virtual float length() const =0
static bool crossedIP(const reco::Track &track)
double wideWidthOfPlane() const
T const & getData(const ESGetToken< T, R > &iToken) const noexcept(false)
Definition: EventSetup.h:119
const LocalTrajectoryError & localError() const
virtual LocalError localDirectionError() const =0
Error on the local direction.
TrajectoryStateOnSurface tState
bool hasPhi() const
Does it have the Phi projection?
float width() const override
const CSCWireTopology * wireTopology() const
double t0() const
Get the segment t0 (if recomputed, 0 is returned otherwise)
T z() const
Definition: PV3DBase.h:61
math::XYZPoint segmentGlobalPosition
constexpr bool isNotFinite(T x)
Definition: isFinite.h:9
edm::Ref< CSCSegmentCollection > CSCSegmentRef
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
edm::Ref< ME0SegmentCollection > ME0SegmentRef
math::XYZPoint segmentLocalPosition
virtual LocalVector localDirection() const =0
Local direction.
T eta() const
Definition: PV3DBase.h:73
int init
Definition: HydjetWrapper.h:64
void useDefaultPropagator()
use the default propagator
std::vector< CaloTower >::const_iterator const_iterator
void fillHcal(const edm::Event &, TrackDetMatchInfo &, const AssociatorParameters &)
bool addTAMuonSegmentMatch(TAMuonChamberMatch &, const RecSegment *, const AssociatorParameters &)
LocalError positionError() const
LocalPoint toLocal(const GlobalPoint &gp) const
math::XYZVector segmentLocalDirection
static const double deltaEta
Definition: CaloConstants.h:8
#define LogTrace(id)
GlobalPoint position() const
void fillCaloTowers(const edm::Event &, TrackDetMatchInfo &, const AssociatorParameters &)
GlobalPoint position() const
void fillCaloTruth(const edm::Event &, TrackDetMatchInfo &, const AssociatorParameters &)
C::const_iterator const_iterator
constant access iterator type
Definition: RangeMap.h:43
float yy() const
Definition: LocalError.h:24
T x() const
Definition: PV3DBase.h:59
T y() const
Definition: PV3DBase.h:60
void init(const edm::EventSetup &, const AssociatorParameters &)
float wireAngle() const override
int iEvent
Definition: GenABIO.cc:224
void setPropagator(const Propagator *)
use a user configured propagator
void getTAMuonChamberMatches(std::vector< TAMuonChamberMatch > &matches, const AssociatorParameters &parameters)
T sqrt(T t)
Definition: SSEVec.h:19
edm::Ref< GEMSegmentCollection > GEMSegmentRef
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventID const &, edm::Timestamp const & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
Definition: Activities.doc:12
std::vector< TAMuonSegmentMatch > segments
distance sign convention: negative - crossed chamber, positive - missed chamber
T mag() const
Definition: PV3DBase.h:64
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double f[11][100]
void fillPreshower(const edm::Event &iEvent, TrackDetMatchInfo &info, const AssociatorParameters &)
#define M_PI
edm::Ref< DTRecSegment4DCollection > DTRecSegment4DRef
ROOT::Math::SMatrix< double, 6, 6, ROOT::Math::MatRepSym< double, 6 > > AlgebraicSymMatrix66
float xy() const
Definition: LocalError.h:23
float yOfWire(float wire, float x=0.) const
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:37
DetId geographicalId() const
const DTChamberRecSegment2D * phiSegment() const
The superPhi segment: 0 if no phi projection available.
const CSCLayerGeometry * oddLayerGeometry(int iendcap) const
Accessors for LayerGeometry&#39;s.
std::vector< DTRecHit1D > specificRecHits() const
Access to specific components.
static std::string info(const DetId &, const TrackerTopology *tTopo)
Definition: DetIdInfo.cc:25
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:31
virtual LocalError localPositionError() const =0
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
void fillMuon(const edm::Event &, TrackDetMatchInfo &, const AssociatorParameters &)
double lengthOfPlane() const
static FreeTrajectoryState getFreeTrajectoryState(const MagneticField *, const reco::Track &)
get FreeTrajectoryState from different track representations
bool isValid() const
Definition: HandleBase.h:70
DetIdAssociator::MapRange getMapRange(const std::pair< float, float > &delta, const float dR)
double narrowWidthOfPlane() const
GlobalVector momentum() const
fixed size matrix
HLT enums.
void fillEcal(const edm::Event &, TrackDetMatchInfo &, const AssociatorParameters &)
TrackDetMatchInfo associate(const edm::Event &, const edm::EventSetup &, const FreeTrajectoryState &, const AssociatorParameters &)
FreeTrajectoryState const * freeState(bool withErrors=true) const
Vector3DBase unit() const
Definition: Vector3DBase.h:54
FreeTrajectoryState innerFreeState(const reco::Track &tk, const MagneticField *field, bool withErr=true)
bool hasZed() const
Does it have the Z projection?
virtual LocalPoint localPosition() const =0
virtual float width() const =0
FreeTrajectoryState outerFreeState(const reco::Track &tk, const MagneticField *field, bool withErr=true)
FreeTrajectoryState initialFreeState(const reco::Track &tk, const MagneticField *field, bool withErr=true)
void fillHO(const edm::Event &, TrackDetMatchInfo &, const AssociatorParameters &)
float xx() const
Definition: LocalError.h:22
def move(src, dest)
Definition: eostools.py:511
*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
const Bounds & bounds() const
Definition: Surface.h:87