CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros 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 // $Id: TrackDetectorAssociator.cc,v 1.42 2010/03/12 13:51:28 gpetrucc Exp $
17 //
18 //
19 
23 
24 // user include files
26 
33 
43 
44 // calorimeter info
51 
55 
58 
61 
63 
64 
67 
70 #include <stack>
71 #include <set>
72 
74 #include "Math/VectorUtil.h"
75 #include <algorithm>
76 
80 // #include "Utilities/Timing/interface/TimerStack.h"
81 
86 
95 
97 
98 using namespace reco;
99 
101 {
102  ivProp_ = 0;
103  defProp_ = 0;
104  useDefaultPropagator_ = false;
105 }
106 
108 {
109  if (defProp_) delete defProp_;
110 }
111 
113 {
114  ivProp_ = ptr;
115  cachedTrajectory_.setPropagator(ivProp_);
116 }
117 
119 {
120  useDefaultPropagator_ = true;
121 }
122 
123 
125 {
126  // access the calorimeter geometry
127  iSetup.get<CaloGeometryRecord>().get(theCaloGeometry_);
128  if (!theCaloGeometry_.isValid())
129  throw cms::Exception("FatalError") << "Unable to find CaloGeometryRecord in event!\n";
130 
131  // get the tracking Geometry
132  iSetup.get<GlobalTrackingGeometryRecord>().get(theTrackingGeometry_);
133  if (!theTrackingGeometry_.isValid())
134  throw cms::Exception("FatalError") << "Unable to find GlobalTrackingGeometryRecord in event!\n";
135 
136  if (useDefaultPropagator_ && (! defProp_ || theMagneticFeildWatcher_.check(iSetup) ) ) {
137  // setup propagator
139  iSetup.get<IdealMagneticFieldRecord>().get(bField);
140 
142  prop->setMaterialMode(false);
143  prop->applyRadX0Correction(true);
144  // prop->setDebug(true); // tmp
145  defProp_ = prop;
146  setPropagator(defProp_);
147  }
148 
149  iSetup.get<DetIdAssociatorRecord>().get("EcalDetIdAssociator", ecalDetIdAssociator_);
150  iSetup.get<DetIdAssociatorRecord>().get("HcalDetIdAssociator", hcalDetIdAssociator_);
151  iSetup.get<DetIdAssociatorRecord>().get("HODetIdAssociator", hoDetIdAssociator_);
152  iSetup.get<DetIdAssociatorRecord>().get("CaloDetIdAssociator", caloDetIdAssociator_);
153  iSetup.get<DetIdAssociatorRecord>().get("MuonDetIdAssociator", muonDetIdAssociator_);
154  iSetup.get<DetIdAssociatorRecord>().get("PreshowerDetIdAssociator", preshowerDetIdAssociator_);
155 }
156 
158  const edm::EventSetup& iSetup,
159  const FreeTrajectoryState& fts,
161 {
162  return associate(iEvent,iSetup,parameters,&fts);
163 }
164 
166  const edm::EventSetup& iSetup,
168  const FreeTrajectoryState* innerState,
169  const FreeTrajectoryState* outerState)
170 {
172  if (! parameters.useEcal && ! parameters.useCalo && ! parameters.useHcal &&
173  ! parameters.useHO && ! parameters.useMuon && !parameters.usePreshower)
174  throw cms::Exception("ConfigurationError") <<
175  "Configuration error! No subdetector was selected for the track association.";
176  // TimerStack timers;
177  // timers.push("TrackDetectorAssociator::associate",TimerStack::DetailedMonitoring);
178 
179  SteppingHelixStateInfo trackOrigin(*innerState);
180  info.stateAtIP = *innerState;
181  cachedTrajectory_.setStateAtIP(trackOrigin);
182 
183  init( iSetup );
184  // get track trajectory
185  // timers.push("TrackDetectorAssociator::fillEcal::propagation");
186  // ECAL points (EB+EE)
187  // If the phi angle between a track entrance and exit points is more
188  // than 2 crystals, it is possible that the track will cross 3 crystals
189  // and therefore one has to check at least 3 points along the track
190  // trajectory inside ECAL. In order to have a chance to cross 4 crystalls
191  // in the barrel, a track should have P_t as low as 3 GeV or smaller
192  // If it's necessary, number of points along trajectory can be increased
193 
194  info.setCaloGeometry(theCaloGeometry_);
195 
196  // timers.push("TrackDetectorAssociator::associate::getTrajectories");
197  cachedTrajectory_.reset_trajectory();
198  // estimate propagation outer boundaries based on
199  // requested sub-detector information. For now limit
200  // propagation region only if muon matching is not
201  // requested.
202  double HOmaxR = hoDetIdAssociator_->volume().maxR();
203  double HOmaxZ = hoDetIdAssociator_->volume().maxZ();
204  double minR = ecalDetIdAssociator_->volume().minR();
205  double minZ = preshowerDetIdAssociator_->volume().minZ();
206  cachedTrajectory_.setMaxHORadius(HOmaxR);
207  cachedTrajectory_.setMaxHOLength(HOmaxZ*2.);
208  cachedTrajectory_.setMinDetectorRadius(minR);
209  cachedTrajectory_.setMinDetectorLength(minZ*2.);
210 
211  double maxR(0);
212  double maxZ(0);
213 
214  if (parameters.useMuon) {
215  maxR = muonDetIdAssociator_->volume().maxR();
216  maxZ = muonDetIdAssociator_->volume().maxZ();
217  cachedTrajectory_.setMaxDetectorRadius(maxR);
218  cachedTrajectory_.setMaxDetectorLength(maxZ*2.);
219  }
220  else {
221  maxR = HOmaxR;
222  maxZ = HOmaxZ;
223  cachedTrajectory_.setMaxDetectorRadius(HOmaxR);
224  cachedTrajectory_.setMaxDetectorLength(HOmaxZ*2.);
225  }
226 
227  // If track extras exist and outerState is before HO maximum, then use outerState
228  if (outerState) {
229  if (outerState->position().perp()<HOmaxR && fabs(outerState->position().z())<HOmaxZ) {
230  LogTrace("TrackAssociator") << "Using outerState as trackOrigin at Rho=" << outerState->position().perp()
231  << " Z=" << outerState->position().z() << "\n";
232  trackOrigin = SteppingHelixStateInfo(*outerState);
233  }
234  else if(innerState) {
235  LogTrace("TrackAssociator") << "Using innerState as trackOrigin at Rho=" << innerState->position().perp()
236  << " Z=" << innerState->position().z() << "\n";
237  trackOrigin = SteppingHelixStateInfo(*innerState);
238  }
239  }
240 
241  if ( ! cachedTrajectory_.propagateAll(trackOrigin) ) return info;
242 
243  // get trajectory in calorimeters
244  cachedTrajectory_.findEcalTrajectory( ecalDetIdAssociator_->volume() );
245  cachedTrajectory_.findHcalTrajectory( hcalDetIdAssociator_->volume() );
246  cachedTrajectory_.findHOTrajectory( hoDetIdAssociator_->volume() );
247  cachedTrajectory_.findPreshowerTrajectory( preshowerDetIdAssociator_->volume() );
248 
249  info.trkGlobPosAtEcal = getPoint( cachedTrajectory_.getStateAtEcal().position() );
250  info.trkGlobPosAtHcal = getPoint( cachedTrajectory_.getStateAtHcal().position() );
251  info.trkGlobPosAtHO = getPoint( cachedTrajectory_.getStateAtHO().position() );
252 
253  info.trkMomAtEcal = cachedTrajectory_.getStateAtEcal().momentum();
254  info.trkMomAtHcal = cachedTrajectory_.getStateAtHcal().momentum();
255  info.trkMomAtHO = cachedTrajectory_.getStateAtHO().momentum();
256 
257  // timers.pop_and_push("TrackDetectorAssociator::associate::fillInfo");
258 
259  if (parameters.useEcal) fillEcal( iEvent, info, parameters);
260  if (parameters.useCalo) fillCaloTowers( iEvent, info, parameters);
261  if (parameters.useHcal) fillHcal( iEvent, info, parameters);
262  if (parameters.useHO) fillHO( iEvent, info, parameters);
263  if (parameters.usePreshower) fillPreshower( iEvent, info, parameters);
264  if (parameters.useMuon) fillMuon( iEvent, info, parameters);
265  if (parameters.truthMatch) fillCaloTruth( iEvent, info, parameters);
266 
267  return info;
268 }
269 
273 {
274  // TimerStack timers;
275  // timers.push("TrackDetectorAssociator::fillEcal");
276 
277  const std::vector<SteppingHelixStateInfo>& trajectoryStates = cachedTrajectory_.getEcalTrajectory();
278 
279  for(std::vector<SteppingHelixStateInfo>::const_iterator itr = trajectoryStates.begin();
280  itr != trajectoryStates.end(); itr++)
281  LogTrace("TrackAssociator") << "ECAL trajectory point (rho, z, phi): " << itr->position().perp() <<
282  ", " << itr->position().z() << ", " << itr->position().phi();
283 
284  std::vector<GlobalPoint> coreTrajectory;
285  for(std::vector<SteppingHelixStateInfo>::const_iterator itr = trajectoryStates.begin();
286  itr != trajectoryStates.end(); itr++) coreTrajectory.push_back(itr->position());
287 
288  if(coreTrajectory.empty()) {
289  LogTrace("TrackAssociator") << "ECAL track trajectory is empty; moving on\n";
290  info.isGoodEcal = 0;
291  return;
292  }
293  info.isGoodEcal = 1;
294 
295  // Find ECAL crystals
296  // timers.push("TrackDetectorAssociator::fillEcal::access::EcalBarrel");
298  iEvent.getByLabel( parameters.theEBRecHitCollectionLabel, EBRecHits );
299  if (!EBRecHits.isValid()) throw cms::Exception("FatalError") << "Unable to find EBRecHitCollection in the event!\n";
300 
301  // timers.pop_and_push("TrackDetectorAssociator::fillEcal::access::EcalEndcaps");
303  iEvent.getByLabel( parameters.theEERecHitCollectionLabel, EERecHits );
304  if (!EERecHits.isValid()) throw cms::Exception("FatalError") << "Unable to find EERecHitCollection in event!\n";
305 
306  // timers.pop_and_push("TrackDetectorAssociator::fillEcal::matching");
307  // timers.push("TrackDetectorAssociator::fillEcal::matching::region");
308  std::set<DetId> ecalIdsInRegion;
309  if (parameters.accountForTrajectoryChangeCalo){
310  // get trajectory change with respect to initial state
311  DetIdAssociator::MapRange mapRange = getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::IpToEcal),
312  parameters.dREcalPreselection);
313  ecalIdsInRegion = ecalDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0],mapRange);
314  } else ecalIdsInRegion = ecalDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], parameters.dREcalPreselection);
315  // timers.pop_and_push("TrackDetectorAssociator::fillEcal::matching::cone");
316  LogTrace("TrackAssociator") << "ECAL hits in the region: " << ecalIdsInRegion.size();
317  if (parameters.dREcalPreselection > parameters.dREcal)
318  ecalIdsInRegion = ecalDetIdAssociator_->getDetIdsInACone(ecalIdsInRegion, coreTrajectory, parameters.dREcal);
319  LogTrace("TrackAssociator") << "ECAL hits in the cone: " << ecalIdsInRegion.size();
320  std::vector<DetId> crossedEcalIds =
321  ecalDetIdAssociator_->getCrossedDetIds(ecalIdsInRegion, coreTrajectory);
322  LogTrace("TrackAssociator") << "ECAL crossed hits " << crossedEcalIds.size();
323 
324  info.crossedEcalIds = crossedEcalIds;
325 
326  // add EcalRecHits
327  // timers.pop_and_push("TrackDetectorAssociator::fillEcal::addCrossedHits");
328  for(std::vector<DetId>::const_iterator itr=crossedEcalIds.begin(); itr!=crossedEcalIds.end();itr++)
329  {
330  std::vector<EcalRecHit>::const_iterator ebHit = (*EBRecHits).find(*itr);
331  std::vector<EcalRecHit>::const_iterator eeHit = (*EERecHits).find(*itr);
332  if(ebHit != (*EBRecHits).end())
333  info.crossedEcalRecHits.push_back(&*ebHit);
334  else if(eeHit != (*EERecHits).end())
335  info.crossedEcalRecHits.push_back(&*eeHit);
336  else
337  LogTrace("TrackAssociator") << "Crossed EcalRecHit is not found for DetId: " << itr->rawId();
338  }
339  // timers.pop_and_push("TrackDetectorAssociator::fillEcal::addHitsInTheRegion");
340  for(std::set<DetId>::const_iterator itr=ecalIdsInRegion.begin(); itr!=ecalIdsInRegion.end();itr++)
341  {
342  std::vector<EcalRecHit>::const_iterator ebHit = (*EBRecHits).find(*itr);
343  std::vector<EcalRecHit>::const_iterator eeHit = (*EERecHits).find(*itr);
344  if(ebHit != (*EBRecHits).end())
345  info.ecalRecHits.push_back(&*ebHit);
346  else if(eeHit != (*EERecHits).end())
347  info.ecalRecHits.push_back(&*eeHit);
348  else
349  LogTrace("TrackAssociator") << "EcalRecHit from the cone is not found for DetId: " << itr->rawId();
350  }
351 }
352 
356 {
357  // TimerStack timers;
358  // timers.push("TrackDetectorAssociator::fillCaloTowers");
359 
360  // use ECAL and HCAL trajectories to match a tower. (HO isn't used for matching).
361  std::vector<GlobalPoint> trajectory;
362  const std::vector<SteppingHelixStateInfo>& ecalTrajectoryStates = cachedTrajectory_.getEcalTrajectory();
363  const std::vector<SteppingHelixStateInfo>& hcalTrajectoryStates = cachedTrajectory_.getHcalTrajectory();
364  for(std::vector<SteppingHelixStateInfo>::const_iterator itr = ecalTrajectoryStates.begin();
365  itr != ecalTrajectoryStates.end(); itr++) trajectory.push_back(itr->position());
366  for(std::vector<SteppingHelixStateInfo>::const_iterator itr = hcalTrajectoryStates.begin();
367  itr != hcalTrajectoryStates.end(); itr++) trajectory.push_back(itr->position());
368 
369  if(trajectory.empty()) {
370  LogTrace("TrackAssociator") << "HCAL trajectory is empty; moving on\n";
371  info.isGoodCalo = 0;
372  return;
373  }
374  info.isGoodCalo = 1;
375 
376  // find crossed CaloTowers
377  // timers.push("TrackDetectorAssociator::fillCaloTowers::access::CaloTowers");
379 
380  iEvent.getByLabel( parameters.theCaloTowerCollectionLabel, caloTowers );
381  if (!caloTowers.isValid()) throw cms::Exception("FatalError") << "Unable to find CaloTowers in event!\n";
382 
383  // timers.pop_and_push("TrackDetectorAssociator::fillCaloTowers::matching");
384  std::set<DetId> caloTowerIdsInRegion;
385  if (parameters.accountForTrajectoryChangeCalo){
386  // get trajectory change with respect to initial state
387  DetIdAssociator::MapRange mapRange = getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::IpToHcal),
388  parameters.dRHcalPreselection);
389  caloTowerIdsInRegion = caloDetIdAssociator_->getDetIdsCloseToAPoint(trajectory[0],mapRange);
390  } else caloTowerIdsInRegion = caloDetIdAssociator_->getDetIdsCloseToAPoint(trajectory[0], parameters.dRHcalPreselection);
391 
392  LogTrace("TrackAssociator") << "Towers in the region: " << caloTowerIdsInRegion.size();
393  std::set<DetId> caloTowerIdsInACone = caloDetIdAssociator_->getDetIdsInACone(caloTowerIdsInRegion, trajectory, parameters.dRHcal);
394  LogTrace("TrackAssociator") << "Towers in the cone: " << caloTowerIdsInACone.size();
395  std::vector<DetId> crossedCaloTowerIds = caloDetIdAssociator_->getCrossedDetIds(caloTowerIdsInRegion, trajectory);
396  LogTrace("TrackAssociator") << "Towers crossed: " << crossedCaloTowerIds.size();
397 
398  info.crossedTowerIds = crossedCaloTowerIds;
399 
400  // add CaloTowers
401  // timers.pop_and_push("TrackDetectorAssociator::fillCaloTowers::addCrossedTowers");
402  for(std::vector<DetId>::const_iterator itr=crossedCaloTowerIds.begin(); itr!=crossedCaloTowerIds.end();itr++)
403  {
404  CaloTowerCollection::const_iterator tower = (*caloTowers).find(*itr);
405  if(tower != (*caloTowers).end())
406  info.crossedTowers.push_back(&*tower);
407  else
408  LogTrace("TrackAssociator") << "Crossed CaloTower is not found for DetId: " << (*itr).rawId();
409  }
410 
411  // timers.pop_and_push("TrackDetectorAssociator::fillCaloTowers::addTowersInTheRegion");
412  for(std::set<DetId>::const_iterator itr=caloTowerIdsInACone.begin(); itr!=caloTowerIdsInACone.end();itr++)
413  {
414  CaloTowerCollection::const_iterator tower = (*caloTowers).find(*itr);
415  if(tower != (*caloTowers).end())
416  info.towers.push_back(&*tower);
417  else
418  LogTrace("TrackAssociator") << "CaloTower from the cone is not found for DetId: " << (*itr).rawId();
419  }
420 
421 }
422 
426 {
427  std::vector<GlobalPoint> trajectory;
428  const std::vector<SteppingHelixStateInfo>& trajectoryStates = cachedTrajectory_.getPreshowerTrajectory();
429  for(std::vector<SteppingHelixStateInfo>::const_iterator itr = trajectoryStates.begin();
430  itr != trajectoryStates.end(); itr++) trajectory.push_back(itr->position());
431 
432  if(trajectory.empty()) {
433  LogTrace("TrackAssociator") << "Preshower trajectory is empty; moving on\n";
434  return;
435  }
436 
437  std::set<DetId> idsInRegion =
438  preshowerDetIdAssociator_->getDetIdsCloseToAPoint(trajectory[0],
439  parameters.dRPreshowerPreselection);
440 
441  LogTrace("TrackAssociator") << "Number of Preshower Ids in the region: " << idsInRegion.size();
442  std::vector<DetId> crossedIds = preshowerDetIdAssociator_->getCrossedDetIds(idsInRegion, trajectory);
443  LogTrace("TrackAssociator") << "Number of Preshower Ids in crossed: " << crossedIds.size();
444  info.crossedPreshowerIds = crossedIds;
445 }
446 
447 
451 {
452  // TimerStack timers;
453  // timers.push("TrackDetectorAssociator::fillHcal");
454 
455  const std::vector<SteppingHelixStateInfo>& trajectoryStates = cachedTrajectory_.getHcalTrajectory();
456 
457  std::vector<GlobalPoint> coreTrajectory;
458  for(std::vector<SteppingHelixStateInfo>::const_iterator itr = trajectoryStates.begin();
459  itr != trajectoryStates.end(); itr++) coreTrajectory.push_back(itr->position());
460 
461  if(coreTrajectory.empty()) {
462  LogTrace("TrackAssociator") << "HCAL trajectory is empty; moving on\n";
463  info.isGoodHcal = 0;
464  return;
465  }
466  info.isGoodHcal = 1;
467 
468  // find crossed Hcals
469  // timers.push("TrackDetectorAssociator::fillHcal::access::Hcal");
471 
472  iEvent.getByLabel( parameters.theHBHERecHitCollectionLabel, collection );
473  if ( ! collection.isValid() ) throw cms::Exception("FatalError") << "Unable to find HBHERecHits in event!\n";
474 
475  // timers.pop_and_push("TrackDetectorAssociator::fillHcal::matching");
476  std::set<DetId> idsInRegion;
477  if (parameters.accountForTrajectoryChangeCalo){
478  // get trajectory change with respect to initial state
479  DetIdAssociator::MapRange mapRange = getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::IpToHcal),
480  parameters.dRHcalPreselection);
481  idsInRegion = hcalDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], mapRange);
482  } else idsInRegion = hcalDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], parameters.dRHcalPreselection);
483 
484  LogTrace("TrackAssociator") << "HCAL hits in the region: " << idsInRegion.size() << "\n" << DetIdInfo::info(idsInRegion);
485  std::set<DetId> idsInACone = hcalDetIdAssociator_->getDetIdsInACone(idsInRegion, coreTrajectory, parameters.dRHcal);
486  LogTrace("TrackAssociator") << "HCAL hits in the cone: " << idsInACone.size() << "\n" << DetIdInfo::info(idsInACone);
487  std::vector<DetId> crossedIds =
488  hcalDetIdAssociator_->getCrossedDetIds(idsInRegion, coreTrajectory);
489  LogTrace("TrackAssociator") << "HCAL hits crossed: " << crossedIds.size() << "\n" << DetIdInfo::info(crossedIds);
490 
491  info.crossedHcalIds = crossedIds;
492  // timers.pop_and_push("TrackDetectorAssociator::fillHcal::addCrossedHits");
493  // add Hcal
494  // timers.push("TrackDetectorAssociator::fillHcal::addHcal");
495  for(std::vector<DetId>::const_iterator itr=crossedIds.begin(); itr!=crossedIds.end();itr++)
496  {
497  HBHERecHitCollection::const_iterator hit = (*collection).find(*itr);
498  if( hit != (*collection).end() )
499  info.crossedHcalRecHits.push_back(&*hit);
500  else
501  LogTrace("TrackAssociator") << "Crossed HBHERecHit is not found for DetId: " << itr->rawId();
502  }
503  // timers.pop_and_push("TrackDetectorAssociator::fillHcal::addHitsInTheRegion");
504  for(std::set<DetId>::const_iterator itr=idsInACone.begin(); itr!=idsInACone.end();itr++)
505  {
506  HBHERecHitCollection::const_iterator hit = (*collection).find(*itr);
507  if( hit != (*collection).end() )
508  info.hcalRecHits.push_back(&*hit);
509  else
510  LogTrace("TrackAssociator") << "HBHERecHit from the cone is not found for DetId: " << itr->rawId();
511  }
512 }
513 
517 {
518  // TimerStack timers;
519  // timers.push("TrackDetectorAssociator::fillHO");
520 
521  const std::vector<SteppingHelixStateInfo>& trajectoryStates = cachedTrajectory_.getHOTrajectory();
522 
523  std::vector<GlobalPoint> coreTrajectory;
524  for(std::vector<SteppingHelixStateInfo>::const_iterator itr = trajectoryStates.begin();
525  itr != trajectoryStates.end(); itr++) coreTrajectory.push_back(itr->position());
526 
527  if(coreTrajectory.empty()) {
528  LogTrace("TrackAssociator") << "HO trajectory is empty; moving on\n";
529  info.isGoodHO = 0;
530  return;
531  }
532  info.isGoodHO = 1;
533 
534  // find crossed HOs
535  // timers.pop_and_push("TrackDetectorAssociator::fillHO::access::HO");
537 
538  iEvent.getByLabel( parameters.theHORecHitCollectionLabel, collection );
539  if ( ! collection.isValid() ) throw cms::Exception("FatalError") << "Unable to find HORecHits in event!\n";
540 
541  // timers.pop_and_push("TrackDetectorAssociator::fillHO::matching");
542  std::set<DetId> idsInRegion;
543  if (parameters.accountForTrajectoryChangeCalo){
544  // get trajectory change with respect to initial state
545  DetIdAssociator::MapRange mapRange = getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::IpToHO),
546  parameters.dRHcalPreselection);
547  idsInRegion = hoDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], mapRange);
548  } else idsInRegion = hoDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], parameters.dRHcalPreselection);
549 
550  LogTrace("TrackAssociator") << "idsInRegion.size(): " << idsInRegion.size();
551  std::set<DetId> idsInACone = hoDetIdAssociator_->getDetIdsInACone(idsInRegion, coreTrajectory, parameters.dRHcal);
552  LogTrace("TrackAssociator") << "idsInACone.size(): " << idsInACone.size();
553  std::vector<DetId> crossedIds =
554  hoDetIdAssociator_->getCrossedDetIds(idsInRegion, coreTrajectory);
555 
556  info.crossedHOIds = crossedIds;
557 
558  // add HO
559  // timers.pop_and_push("TrackDetectorAssociator::fillHO::addCrossedHits");
560  for(std::vector<DetId>::const_iterator itr=crossedIds.begin(); itr!=crossedIds.end();itr++)
561  {
562  HORecHitCollection::const_iterator hit = (*collection).find(*itr);
563  if( hit != (*collection).end() )
564  info.crossedHORecHits.push_back(&*hit);
565  else
566  LogTrace("TrackAssociator") << "Crossed HORecHit is not found for DetId: " << itr->rawId();
567  }
568 
569  // timers.pop_and_push("TrackDetectorAssociator::fillHO::addHitsInTheRegion");
570  for(std::set<DetId>::const_iterator itr=idsInACone.begin(); itr!=idsInACone.end();itr++)
571  {
572  HORecHitCollection::const_iterator hit = (*collection).find(*itr);
573  if( hit != (*collection).end() )
574  info.hoRecHits.push_back(&*hit);
575  else
576  LogTrace("TrackAssociator") << "HORecHit from the cone is not found for DetId: " << itr->rawId();
577  }
578 }
579 
581  const SimTrack& track,
582  const SimVertex& vertex )
583 {
584  GlobalVector vector( track.momentum().x(), track.momentum().y(), track.momentum().z() );
585  GlobalPoint point( vertex.position().x(), vertex.position().y(), vertex.position().z() );
586 
587  int charge = track.type( )> 0 ? -1 : 1; // lepton convention
588  if ( abs(track.type( )) == 211 || // pion
589  abs(track.type( )) == 321 || // kaon
590  abs(track.type( )) == 2212 )
591  charge = track.type( )< 0 ? -1 : 1;
592  return getFreeTrajectoryState(iSetup, vector, point, charge);
593 }
594 
596  const GlobalVector& momentum,
597  const GlobalPoint& vertex,
598  const int charge)
599 {
601  iSetup.get<IdealMagneticFieldRecord>().get(bField);
602 
603  GlobalTrajectoryParameters tPars(vertex, momentum, charge, &*bField);
604 
605  CLHEP::HepSymMatrix covT(6,1); covT *= 1e-6; // initialize to sigma=1e-3
606  CartesianTrajectoryError tCov(covT);
607 
608  return FreeTrajectoryState(tPars, tCov);
609 }
610 
611 
613  const reco::Track& track )
614 {
616  iSetup.get<IdealMagneticFieldRecord>().get(bField);
617 
618  GlobalVector vector( track.momentum().x(), track.momentum().y(), track.momentum().z() );
619 
620  GlobalPoint point( track.vertex().x(), track.vertex().y(), track.vertex().z() );
621 
622  GlobalTrajectoryParameters tPars(point, vector, track.charge(), &*bField);
623 
624  // FIX THIS !!!
625  // need to convert from perigee to global or helix (curvilinear) frame
626  // for now just an arbitrary matrix.
627  CLHEP::HepSymMatrix covT(6,1); covT *= 1e-6; // initialize to sigma=1e-3
628  CartesianTrajectoryError tCov(covT);
629 
630  return FreeTrajectoryState(tPars, tCov);
631 }
632 
634  const float dR )
635 {
636  DetIdAssociator::MapRange mapRange;
637  mapRange.dThetaPlus = dR;
638  mapRange.dThetaMinus = dR;
639  mapRange.dPhiPlus = dR;
640  mapRange.dPhiMinus = dR;
641  if ( delta.first > 0 )
642  mapRange.dThetaPlus += delta.first;
643  else
644  mapRange.dThetaMinus += fabs(delta.first);
645  if ( delta.second > 0 )
646  mapRange.dPhiPlus += delta.second;
647  else
648  mapRange.dPhiMinus += fabs(delta.second);
649  LogTrace("TrackAssociator") << "Selection range: (dThetaPlus, dThetaMinus, dPhiPlus, dPhiMinus, dRPreselection): " <<
650  mapRange.dThetaPlus << ", " << mapRange.dThetaMinus << ", " <<
651  mapRange.dPhiPlus << ", " << mapRange.dPhiMinus << ", " << dR;
652  return mapRange;
653 }
654 
655 void TrackDetectorAssociator::getTAMuonChamberMatches(std::vector<TAMuonChamberMatch>& matches,
657 {
658  // Strategy:
659  // Propagate through the whole detector, estimate change in eta and phi
660  // along the trajectory, add this to dRMuon and find DetIds around this
661  // direction using the map. Then propagate fast to each surface and apply
662  // final matching criteria.
663 
664  // TimerStack timers(TimerStack::Disableable);
665  // timers.push("MuonDetIdAssociator::getTrajectoryInMuonDetector");
666  // timers.push("MuonDetIdAssociator::getTrajectoryInMuonDetector::propagation",TimerStack::FastMonitoring);
667  // timers.pop();
668  // get the direction first
669  SteppingHelixStateInfo trajectoryPoint = cachedTrajectory_.getStateAtHcal();
670  // If trajectory point at HCAL is not valid, try to use the outer most state of the
671  // trajectory instead.
672  if(! trajectoryPoint.isValid() ) trajectoryPoint = cachedTrajectory_.getOuterState();
673  if(! trajectoryPoint.isValid() ) {
674  LogTrace("TrackAssociator") <<
675  "trajectory position at HCAL is not valid. Assume the track cannot reach muon detectors and skip it";
676  return;
677  }
678 
679  GlobalVector direction = trajectoryPoint.momentum().unit();
680  LogTrace("TrackAssociator") << "muon direction: " << direction << "\n\t and corresponding point: " <<
681  trajectoryPoint.position() <<"\n";
682 
683  DetIdAssociator::MapRange mapRange = getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::FullTrajectory),
684  parameters.dRMuonPreselection);
685 
686  // and find chamber DetIds
687 
688  // timers.push("MuonDetIdAssociator::getTrajectoryInMuonDetector::getDetIdsCloseToAPoint",TimerStack::FastMonitoring);
689  std::set<DetId> muonIdsInRegion =
690  muonDetIdAssociator_->getDetIdsCloseToAPoint(trajectoryPoint.position(), mapRange);
691  // timers.pop_and_push("MuonDetIdAssociator::getTrajectoryInMuonDetector::matching",TimerStack::FastMonitoring);
692  LogTrace("TrackAssociator") << "Number of chambers to check: " << muonIdsInRegion.size();
693  for(std::set<DetId>::const_iterator detId = muonIdsInRegion.begin(); detId != muonIdsInRegion.end(); detId++)
694  {
695  const GeomDet* geomDet = muonDetIdAssociator_->getGeomDet(*detId);
696  // timers.push("MuonDetIdAssociator::getTrajectoryInMuonDetector::matching::localPropagation",TimerStack::FastMonitoring);
697  TrajectoryStateOnSurface stateOnSurface = cachedTrajectory_.propagate( &geomDet->surface() );
698  if (! stateOnSurface.isValid()) {
699  LogTrace("TrackAssociator") << "Failed to propagate the track; moving on\n\t"<<
700  detId->rawId() << " not crossed\n";
701  continue;
702  }
703  // timers.pop_and_push("MuonDetIdAssociator::getTrajectoryInMuonDetector::matching::geometryAccess",TimerStack::FastMonitoring);
704  LocalPoint localPoint = geomDet->surface().toLocal(stateOnSurface.freeState()->position());
705  LocalError localError = stateOnSurface.localError().positionError();
706  float distanceX = 0;
707  float distanceY = 0;
708  float sigmaX = 0.0;
709  float sigmaY = 0.0;
710  if(const CSCChamber* cscChamber = dynamic_cast<const CSCChamber*>(geomDet) ) {
711  const CSCChamberSpecs* chamberSpecs = cscChamber->specs();
712  if(! chamberSpecs) {
713  LogTrace("TrackAssociator") << "Failed to get CSCChamberSpecs from CSCChamber; moving on\n";
714  continue;
715  }
716  const CSCLayerGeometry* layerGeometry = chamberSpecs->oddLayerGeometry(1);
717  if(! layerGeometry) {
718  LogTrace("TrackAssociator") << "Failed to get CSCLayerGeometry from CSCChamberSpecs; moving on\n";
719  continue;
720  }
721  const CSCWireTopology* wireTopology = layerGeometry->wireTopology();
722  if(! wireTopology) {
723  LogTrace("TrackAssociator") << "Failed to get CSCWireTopology from CSCLayerGeometry; moving on\n";
724  continue;
725  }
726 
727  float wideWidth = wireTopology->wideWidthOfPlane();
728  float narrowWidth = wireTopology->narrowWidthOfPlane();
729  float length = wireTopology->lengthOfPlane();
730  // If slanted, there is no y offset between local origin and symmetry center of wire plane
731  float yOfFirstWire = fabs(wireTopology->wireAngle())>1.E-06 ? -0.5*length : wireTopology->yOfWire(1);
732  // y offset between local origin and symmetry center of wire plane
733  float yCOWPOffset = yOfFirstWire+0.5*length;
734 
735  // tangent of the incline angle from inside the trapezoid
736  float tangent = (wideWidth-narrowWidth)/(2.*length);
737  // y position wrt bottom of trapezoid
738  float yPrime = localPoint.y()+fabs(yOfFirstWire);
739  // half trapezoid width at y' is 0.5 * narrowWidth + x side of triangle with the above tangent and side y'
740  float halfWidthAtYPrime = 0.5*narrowWidth+yPrime*tangent;
741  distanceX = fabs(localPoint.x()) - halfWidthAtYPrime;
742  distanceY = fabs(localPoint.y()-yCOWPOffset) - 0.5*length;
743  sigmaX = distanceX/sqrt(localError.xx());
744  sigmaY = distanceY/sqrt(localError.yy());
745  } else {
746  distanceX = fabs(localPoint.x()) - geomDet->surface().bounds().width()/2.;
747  distanceY = fabs(localPoint.y()) - geomDet->surface().bounds().length()/2.;
748  sigmaX = distanceX/sqrt(localError.xx());
749  sigmaY = distanceY/sqrt(localError.yy());
750  }
751  // timers.pop_and_push("MuonDetIdAssociator::getTrajectoryInMuonDetector::matching::checking",TimerStack::FastMonitoring);
752  if ( (distanceX < parameters.muonMaxDistanceX && distanceY < parameters.muonMaxDistanceY) ||
753  (sigmaX < parameters.muonMaxDistanceSigmaX && sigmaY < parameters.muonMaxDistanceSigmaY) ) {
754  LogTrace("TrackAssociator") << "found a match, DetId: " << detId->rawId();
756  match.tState = stateOnSurface;
757  match.localDistanceX = distanceX;
758  match.localDistanceY = distanceY;
759  match.id = *detId;
760  matches.push_back(match);
761  }
762  //timers.pop();
763  }
764  //timers.pop();
765 
766 }
767 
771 {
772  // TimerStack timers;
773  // timers.push("TrackDetectorAssociator::fillMuon");
774 
775  // Get the segments from the event
776  // timers.push("TrackDetectorAssociator::fillMuon::access");
778  iEvent.getByLabel( parameters.theDTRecSegment4DCollectionLabel, dtSegments );
779  if (! dtSegments.isValid())
780  throw cms::Exception("FatalError") << "Unable to find DTRecSegment4DCollection in event!\n";
781 
783  iEvent.getByLabel( parameters.theCSCSegmentCollectionLabel, cscSegments );
784  if (! cscSegments.isValid())
785  throw cms::Exception("FatalError") << "Unable to find CSCSegmentCollection in event!\n";
786 
788 
789  // check the map of available segments
790  // if there is no segments in a given direction at all,
791  // then there is no point to fly there.
792  //
793  // MISSING
794  // Possible solution: quick search for presence of segments
795  // for the set of DetIds
796 
797  // timers.pop_and_push("TrackDetectorAssociator::fillMuon::matchChembers");
798 
799  // get a set of matches corresponding to muon chambers
800  std::vector<TAMuonChamberMatch> matchedChambers;
801  LogTrace("TrackAssociator") << "Trying to Get ChamberMatches" << std::endl;
802  getTAMuonChamberMatches(matchedChambers, parameters);
803  LogTrace("TrackAssociator") << "Chambers matched: " << matchedChambers.size() << "\n";
804 
805  // Iterate over all chamber matches and fill segment matching
806  // info if it's available
807  // timers.pop_and_push("TrackDetectorAssociator::fillMuon::findSemgents");
808  for(std::vector<TAMuonChamberMatch>::iterator matchedChamber = matchedChambers.begin();
809  matchedChamber != matchedChambers.end(); matchedChamber++)
810  {
811  const GeomDet* geomDet = muonDetIdAssociator_->getGeomDet((*matchedChamber).id);
812  // DT chamber
813  if(const DTChamber* chamber = dynamic_cast<const DTChamber*>(geomDet) ) {
814  // Get the range for the corresponding segments
815  DTRecSegment4DCollection::range range = dtSegments->get(chamber->id());
816  // Loop over the segments of this chamber
817  for (DTRecSegment4DCollection::const_iterator segment = range.first; segment!=range.second; segment++) {
818  if (addTAMuonSegmentMatch(*matchedChamber, &(*segment), parameters)) {
819  matchedChamber->segments.back().dtSegmentRef = DTRecSegment4DRef(dtSegments, segment - dtSegments->begin());
820  }
821  }
822  }else{
823  // CSC Chamber
824  if(const CSCChamber* chamber = dynamic_cast<const CSCChamber*>(geomDet) ) {
825  // Get the range for the corresponding segments
826  CSCSegmentCollection::range range = cscSegments->get(chamber->id());
827  // Loop over the segments
828  for (CSCSegmentCollection::const_iterator segment = range.first; segment!=range.second; segment++) {
829  if (addTAMuonSegmentMatch(*matchedChamber, &(*segment), parameters)) {
830  matchedChamber->segments.back().cscSegmentRef = CSCSegmentRef(cscSegments, segment - cscSegments->begin());
831  }
832  }
833  }else{
834  throw cms::Exception("FatalError") << "Failed to cast GeomDet object to either DTChamber or CSCChamber. Who is this guy anyway?\n";
835  }
836  }
837  info.chambers.push_back(*matchedChamber);
838  }
839 }
840 
841 
843  const RecSegment* segment,
845 {
846  LogTrace("TrackAssociator")
847  << "Segment local position: " << segment->localPosition() << "\n"
848  << std::hex << segment->geographicalId().rawId() << "\n";
849 
850  const GeomDet* chamber = muonDetIdAssociator_->getGeomDet(matchedChamber.id);
851  TrajectoryStateOnSurface trajectoryStateOnSurface = matchedChamber.tState;
852  GlobalPoint segmentGlobalPosition = chamber->toGlobal(segment->localPosition());
853 
854  LogTrace("TrackAssociator")
855  << "Segment global position: " << segmentGlobalPosition << " \t (R_xy,eta,phi): "
856  << segmentGlobalPosition.perp() << "," << segmentGlobalPosition.eta() << "," << segmentGlobalPosition.phi() << "\n";
857 
858  LogTrace("TrackAssociator")
859  << "\teta hit: " << segmentGlobalPosition.eta() << " \tpropagator: " << trajectoryStateOnSurface.freeState()->position().eta() << "\n"
860  << "\tphi hit: " << segmentGlobalPosition.phi() << " \tpropagator: " << trajectoryStateOnSurface.freeState()->position().phi() << std::endl;
861 
862  bool isGood = false;
863  bool isDTWithoutY = false;
864  const DTRecSegment4D* dtseg = dynamic_cast<const DTRecSegment4D*>(segment);
865  if ( dtseg && (! dtseg->hasZed()) )
866  isDTWithoutY = true;
867 
868  double deltaPhi(fabs(segmentGlobalPosition.phi()-trajectoryStateOnSurface.freeState()->position().phi()));
869  if(deltaPhi>M_PI) deltaPhi = fabs(deltaPhi-M_PI*2.);
870 
871  if( isDTWithoutY )
872  {
873  isGood = deltaPhi < parameters.dRMuon;
874  // Be in chamber
875  isGood &= fabs(segmentGlobalPosition.eta()-trajectoryStateOnSurface.freeState()->position().eta()) < .3;
876  } else isGood = sqrt( pow(segmentGlobalPosition.eta()-trajectoryStateOnSurface.freeState()->position().eta(),2) +
877  deltaPhi*deltaPhi) < parameters.dRMuon;
878 
879  if(isGood) {
880  TAMuonSegmentMatch muonSegment;
881  muonSegment.segmentGlobalPosition = getPoint(segmentGlobalPosition);
882  muonSegment.segmentLocalPosition = getPoint( segment->localPosition() );
883  muonSegment.segmentLocalDirection = getVector( segment->localDirection() );
884  muonSegment.segmentLocalErrorXX = segment->localPositionError().xx();
885  muonSegment.segmentLocalErrorYY = segment->localPositionError().yy();
886  muonSegment.segmentLocalErrorXY = segment->localPositionError().xy();
887  muonSegment.segmentLocalErrorDxDz = segment->localDirectionError().xx();
888  muonSegment.segmentLocalErrorDyDz = segment->localDirectionError().yy();
889 
890  // DANGEROUS - compiler cannot guaranty parameters ordering
891  // AlgebraicSymMatrix segmentCovMatrix = segment->parametersError();
892  // muonSegment.segmentLocalErrorXDxDz = segmentCovMatrix[2][0];
893  // muonSegment.segmentLocalErrorYDyDz = segmentCovMatrix[3][1];
894  muonSegment.segmentLocalErrorXDxDz = -999;
895  muonSegment.segmentLocalErrorYDyDz = -999;
896  muonSegment.hasZed = true;
897  muonSegment.hasPhi = true;
898 
899  // timing information
900  muonSegment.t0 = 0;
901  if ( dtseg ) {
902  if ( (dtseg->hasPhi()) && (! isDTWithoutY) ) {
903  int phiHits = dtseg->phiSegment()->specificRecHits().size();
904  // int zHits = dtseg->zSegment()->specificRecHits().size();
905  int hits=0;
906  double t0=0.;
907  // TODO: cuts on hit numbers not optimized in any way yet...
908  if (phiHits>5 && dtseg->phiSegment()->ist0Valid()) {
909  t0+=dtseg->phiSegment()->t0()*phiHits;
910  hits+=phiHits;
911  LogTrace("TrackAssociator") << " Phi t0: " << dtseg->phiSegment()->t0() << " hits: " << phiHits;
912  }
913  // the z segments seem to contain little useful information...
914 // if (zHits>3) {
915 // t0+=s->zSegment()->t0()*zHits;
916 // hits+=zHits;
917 // std::cout << " Z t0: " << s->zSegment()->t0() << " hits: " << zHits << std::endl;
918 // }
919  if (hits) muonSegment.t0 = t0/hits;
920 // std::cout << " --- t0: " << muonSegment.t0 << std::endl;
921  } else {
922  // check and set dimensionality
923  if (isDTWithoutY) muonSegment.hasZed = false;
924  if (! dtseg->hasPhi()) muonSegment.hasPhi = false;
925  }
926  }
927  matchedChamber.segments.push_back(muonSegment);
928  }
929 
930  return isGood;
931 }
932 
933 //********************** NON-CORE CODE ******************************//
934 
938 {
939  // get list of simulated tracks and their vertices
940  using namespace edm;
941  Handle<SimTrackContainer> simTracks;
942  iEvent.getByType<SimTrackContainer>(simTracks);
943  if (! simTracks.isValid() ) throw cms::Exception("FatalError") << "No simulated tracks found\n";
944 
945  Handle<SimVertexContainer> simVertices;
946  iEvent.getByType<SimVertexContainer>(simVertices);
947  if (! simVertices.isValid() ) throw cms::Exception("FatalError") << "No simulated vertices found\n";
948 
949  // get sim calo hits
950  Handle<PCaloHitContainer> simEcalHitsEB;
951  iEvent.getByLabel("g4SimHits","EcalHitsEB",simEcalHitsEB);
952  if (! simEcalHitsEB.isValid() ) throw cms::Exception("FatalError") << "No simulated ECAL EB hits found\n";
953 
954  Handle<PCaloHitContainer> simEcalHitsEE;
955  iEvent.getByLabel("g4SimHits","EcalHitsEE",simEcalHitsEE);
956  if (! simEcalHitsEE.isValid() ) throw cms::Exception("FatalError") << "No simulated ECAL EE hits found\n";
957 
958  Handle<PCaloHitContainer> simHcalHits;
959  iEvent.getByLabel("g4SimHits","HcalHits",simHcalHits);
960  if (! simHcalHits.isValid() ) throw cms::Exception("FatalError") << "No simulated HCAL hits found\n";
961 
962  // find truth partner
963  SimTrackContainer::const_iterator simTrack = simTracks->begin();
964  for( ; simTrack != simTracks->end(); ++simTrack){
965  math::XYZVector simP3( simTrack->momentum().x(), simTrack->momentum().y(), simTrack->momentum().z() );
966  math::XYZVector recoP3( info.stateAtIP.momentum().x(), info.stateAtIP.momentum().y(), info.stateAtIP.momentum().z() );
967  if ( ROOT::Math::VectorUtil::DeltaR(recoP3, simP3) < 0.1 ) break;
968  }
969  if ( simTrack != simTracks->end() ) {
970  info.simTrack = &(*simTrack);
971  double ecalTrueEnergy(0);
972  double hcalTrueEnergy(0);
973 
974  // loop over calo hits
975  for( PCaloHitContainer::const_iterator hit = simEcalHitsEB->begin(); hit != simEcalHitsEB->end(); ++hit )
976  if ( hit->geantTrackId() == info.simTrack->genpartIndex() ) ecalTrueEnergy += hit->energy();
977 
978  for( PCaloHitContainer::const_iterator hit = simEcalHitsEE->begin(); hit != simEcalHitsEE->end(); ++hit )
979  if ( hit->geantTrackId() == info.simTrack->genpartIndex() ) ecalTrueEnergy += hit->energy();
980 
981  for( PCaloHitContainer::const_iterator hit = simHcalHits->begin(); hit != simHcalHits->end(); ++hit )
982  if ( hit->geantTrackId() == info.simTrack->genpartIndex() ) hcalTrueEnergy += hit->energy();
983 
984  info.ecalTrueEnergy = ecalTrueEnergy;
985  info.hcalTrueEnergy = hcalTrueEnergy;
986  info.hcalTrueEnergyCorrected = hcalTrueEnergy;
987  if ( fabs(info.trkGlobPosAtHcal.eta()) < 1.3 )
988  info.hcalTrueEnergyCorrected = hcalTrueEnergy*113.2;
989  else
990  if ( fabs(info.trkGlobPosAtHcal.eta()) < 3.0 )
991  info.hcalTrueEnergyCorrected = hcalTrueEnergy*167.2;
992  }
993 }
994 
996  const edm::EventSetup& iSetup,
997  const reco::Track& track,
999  Direction direction /*= Any*/ )
1000 {
1001  double currentStepSize = cachedTrajectory_.getPropagationStep();
1002  TrajectoryStateTransform tsTransform;
1004  iSetup.get<IdealMagneticFieldRecord>().get(bField);
1005 
1006  if(track.extra().isNull()) {
1007  if ( direction != InsideOut )
1008  throw cms::Exception("FatalError") <<
1009  "No TrackExtra information is available and association is done with something else than InsideOut track.\n" <<
1010  "Either change the parameter or provide needed data!\n";
1011  LogTrace("TrackAssociator") << "Track Extras not found\n";
1012  FreeTrajectoryState initialState = tsTransform.initialFreeState(track,&*bField);
1013  return associate(iEvent, iSetup, parameters, &initialState); // 5th argument is null pointer
1014  }
1015 
1016  LogTrace("TrackAssociator") << "Track Extras found\n";
1017  FreeTrajectoryState innerState = tsTransform.innerFreeState(track,&*bField);
1018  FreeTrajectoryState outerState = tsTransform.outerFreeState(track,&*bField);
1019  FreeTrajectoryState referenceState = tsTransform.initialFreeState(track,&*bField);
1020 
1021  LogTrace("TrackAssociator") << "inner track state (rho, z, phi):" <<
1022  track.innerPosition().Rho() << ", " << track.innerPosition().z() <<
1023  ", " << track.innerPosition().phi() << "\n";
1024  LogTrace("TrackAssociator") << "innerFreeState (rho, z, phi):" <<
1025  innerState.position().perp() << ", " << innerState.position().z() <<
1026  ", " << innerState.position().phi() << "\n";
1027 
1028  LogTrace("TrackAssociator") << "outer track state (rho, z, phi):" <<
1029  track.outerPosition().Rho() << ", " << track.outerPosition().z() <<
1030  ", " << track.outerPosition().phi() << "\n";
1031  LogTrace("TrackAssociator") << "outerFreeState (rho, z, phi):" <<
1032  outerState.position().perp() << ", " << outerState.position().z() <<
1033  ", " << outerState.position().phi() << "\n";
1034 
1035  // InsideOut first
1036  if ( crossedIP( track ) ) {
1037  switch ( direction ) {
1038  case InsideOut:
1039  case Any:
1040  return associate(iEvent, iSetup, parameters, &referenceState, &outerState);
1041  break;
1042  case OutsideIn:
1043  {
1044  cachedTrajectory_.setPropagationStep( -fabs(currentStepSize) );
1045  TrackDetMatchInfo result = associate(iEvent, iSetup, parameters, &innerState, &referenceState);
1046  cachedTrajectory_.setPropagationStep( currentStepSize );
1047  return result;
1048  break;
1049  }
1050  }
1051  } else {
1052  switch ( direction ) {
1053  case InsideOut:
1054  return associate(iEvent, iSetup, parameters, &innerState, &outerState);
1055  break;
1056  case OutsideIn:
1057  {
1058  cachedTrajectory_.setPropagationStep( -fabs(currentStepSize) );
1059  TrackDetMatchInfo result = associate(iEvent, iSetup, parameters, &outerState, &innerState);
1060  cachedTrajectory_.setPropagationStep( currentStepSize );
1061  return result;
1062  break;
1063  }
1064  case Any:
1065  {
1066  // check if we deal with clear outside-in case
1067  if ( track.innerPosition().Dot( track.innerMomentum() ) < 0 &&
1068  track.outerPosition().Dot( track.outerMomentum() ) < 0 )
1069  {
1070  cachedTrajectory_.setPropagationStep( -fabs(currentStepSize) );
1072  if ( track.innerPosition().R() < track.outerPosition().R() )
1073  result = associate(iEvent, iSetup, parameters, &innerState, &outerState);
1074  else
1075  result = associate(iEvent, iSetup, parameters, &outerState, &innerState);
1076  cachedTrajectory_.setPropagationStep( currentStepSize );
1077  return result;
1078  }
1079  }
1080  }
1081  }
1082 
1083  // all other cases
1084  return associate(iEvent, iSetup, parameters, &innerState, &outerState);
1085 }
1086 
1088  const edm::EventSetup& iSetup,
1089  const SimTrack& track,
1090  const SimVertex& vertex,
1092 {
1093  return associate(iEvent, iSetup, getFreeTrajectoryState(iSetup, track, vertex), parameters);
1094 }
1095 
1097  const edm::EventSetup& iSetup,
1098  const GlobalVector& momentum,
1099  const GlobalPoint& vertex,
1100  const int charge,
1102 {
1103  return associate(iEvent, iSetup, getFreeTrajectoryState(iSetup, momentum, vertex, charge), parameters);
1104 }
1105 
1107 {
1108  bool crossed = true;
1109  crossed &= (track.innerPosition().rho() > 3 ); // something close to active volume
1110  crossed &= (track.outerPosition().rho() > 3 ); // something close to active volume
1111  crossed &= ( ( track.innerPosition().x()*track.innerMomentum().x() +
1112  track.innerPosition().y()*track.innerMomentum().y() < 0 ) !=
1113  ( track.outerPosition().x()*track.outerMomentum().x() +
1114  track.outerPosition().y()*track.outerMomentum().y() < 0 ) );
1115  return crossed;
1116 }
dbl * delta
Definition: mlp_gen.cc:36
std::vector< DetId > crossedPreshowerIds
math::XYZPoint trkGlobPosAtHO
float xx() const
Definition: LocalError.h:19
std::pair< const_iterator, const_iterator > range
iterator range
Definition: RangeMap.h:52
const CSCWireTopology * wireTopology() const
virtual float length() const =0
T perp() const
Definition: PV3DBase.h:66
static bool crossedIP(const reco::Track &track)
const Vector & momentum() const
track momentum vector
Definition: TrackBase.h:149
std::vector< const CaloTower * > crossedTowers
virtual LocalError localDirectionError() const =0
Error on the local direction.
std::vector< const CaloTower * > towers
TrajectoryStateOnSurface tState
const TrackExtraRef & extra() const
reference to &quot;extra&quot; object
Definition: Track.h:97
std::vector< const HBHERecHit * > crossedHcalRecHits
std::vector< DetId > crossedTowerIds
math::XYZPoint segmentGlobalPosition
static FreeTrajectoryState getFreeTrajectoryState(const edm::EventSetup &, const reco::Track &)
get FreeTrajectoryState from different track representations
edm::Ref< CSCSegmentCollection > CSCSegmentRef
const DTChamberRecSegment2D * phiSegment() const
The superPhi segment: 0 if no phi projection available.
std::vector< const EcalRecHit * > ecalRecHits
hits in the cone
std::vector< DetId > crossedEcalIds
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:49
math::XYZPoint segmentLocalPosition
void setCaloGeometry(edm::ESHandle< CaloGeometry > geometry)
virtual LocalVector localDirection() const =0
Local direction.
int init
Definition: HydjetWrapper.h:63
void useDefaultPropagator()
use the default propagator
Geom::Phi< T > phi() const
Definition: PV3DBase.h:63
std::vector< CaloTower >::const_iterator const_iterator
T y() const
Definition: PV3DBase.h:57
#define abs(x)
Definition: mlp_lapack.h:159
void fillHcal(const edm::Event &, TrackDetMatchInfo &, const AssociatorParameters &)
bool getByType(Handle< PROD > &result) const
Definition: Event.h:403
bool addTAMuonSegmentMatch(TAMuonChamberMatch &, const RecSegment *, const AssociatorParameters &)
const math::XYZPoint & outerPosition() const
position of the outermost hit
Definition: Track.h:47
std::vector< DetId > crossedHcalIds
GlobalVector momentum() const
std::vector< const EcalRecHit * > crossedEcalRecHits
hits in detector elements crossed by a track
double charge(const std::vector< uint8_t > &Ampls)
math::XYZVector segmentLocalDirection
LocalError positionError() const
math::XYZPoint trkGlobPosAtHcal
FreeTrajectoryState innerFreeState(const reco::Track &tk, const MagneticField *field) const
FreeTrajectoryState stateAtIP
track info
void fillCaloTowers(const edm::Event &, TrackDetMatchInfo &, const AssociatorParameters &)
const math::XYZPoint & innerPosition() const
position of the innermost hit
Definition: Track.h:42
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
void fillCaloTruth(const edm::Event &, TrackDetMatchInfo &, const AssociatorParameters &)
C::const_iterator const_iterator
constant access iterator type
Definition: RangeMap.h:45
float xy() const
Definition: LocalError.h:20
float wireAngle() const
int iEvent
Definition: GenABIO.cc:243
bool isNull() const
Checks for null.
Definition: Ref.h:244
GlobalVector trkMomAtHO
void setPropagator(const Propagator *)
use a user configured propagator
void getTAMuonChamberMatches(std::vector< TAMuonChamberMatch > &matches, const AssociatorParameters &parameters)
FreeTrajectoryState * freeState(bool withErrors=true) const
float yy() const
Definition: LocalError.h:21
void applyRadX0Correction(bool applyRadX0Correction)
GlobalPoint position() const
T sqrt(T t)
Definition: SSEVec.h:28
LocalPoint toLocal(const GlobalPoint &gp) const
T z() const
Definition: PV3DBase.h:58
std::vector< TAMuonSegmentMatch > segments
distance sign convention: negative - crossed chamber, positive - missed chamber
tuple result
Definition: query.py:137
FreeTrajectoryState outerFreeState(const reco::Track &tk, const MagneticField *field) const
int genpartIndex() const
index of the corresponding Generator particle in the Event container (-1 if no Genpart) ...
Definition: SimTrack.h:33
std::vector< TAMuonChamberMatch > chambers
std::vector< const HBHERecHit * > hcalRecHits
std::vector< DetId > crossedHOIds
const math::XYZTLorentzVectorD & position() const
Definition: CoreSimVertex.h:26
static std::string info(const DetId &)
Definition: DetIdInfo.cc:28
bool hasPhi() const
Does it have the Phi projection?
double lengthOfPlane() const
const LocalTrajectoryError & localError() const
float yOfWire(float wire, float x=0.) const
bool isValid() const
Definition: HandleBase.h:76
void fillPreshower(const edm::Event &iEvent, TrackDetMatchInfo &info, const AssociatorParameters &)
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:359
GlobalVector momentum() const
#define LogTrace(id)
std::vector< DTRecHit1D > specificRecHits() const
Access to specific components.
bool ist0Valid() const
double narrowWidthOfPlane() const
edm::Ref< DTRecSegment4DCollection > DTRecSegment4DRef
const Point & vertex() const
reference point on the track. This method is DEPRECATED, please use referencePoint() instead ...
Definition: TrackBase.h:155
double deltaPhi(double phi1, double phi2)
Definition: deltaPhi.h:12
Vector3DBase unit() const
Definition: Vector3DBase.h:57
bool hasZed() const
Does it have the Z projection?
GlobalPoint position() const
const Bounds & bounds() const
Definition: BoundSurface.h:89
#define M_PI
Definition: BFit3D.cc:3
const math::XYZVector & outerMomentum() const
momentum vector at the outermost hit position
Definition: Track.h:49
double wideWidthOfPlane() const
GlobalVector trkMomAtEcal
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:31
virtual LocalError localPositionError() const =0
const T & get() const
Definition: EventSetup.h:55
std::vector< SimVertex > SimVertexContainer
void fillMuon(const edm::Event &, TrackDetMatchInfo &, const AssociatorParameters &)
void init(const edm::EventSetup &)
const CSCLayerGeometry *const oddLayerGeometry(int iendcap) const
Accessors for LayerGeometry&#39;s.
DetIdAssociator::MapRange getMapRange(const std::pair< float, float > &delta, const float dR)
T eta() const
Definition: PV3DBase.h:70
GlobalVector trkMomAtHcal
int type() const
particle type (HEP PDT convension)
Definition: CoreSimTrack.h:40
void fillEcal(const edm::Event &, TrackDetMatchInfo &, const AssociatorParameters &)
const math::XYZTLorentzVectorD & momentum() const
particle info...
Definition: CoreSimTrack.h:36
const math::XYZVector & innerMomentum() const
momentum vector at the innermost hit position
Definition: Track.h:45
std::vector< const HORecHit * > crossedHORecHits
TrackDetMatchInfo associate(const edm::Event &, const edm::EventSetup &, const FreeTrajectoryState &, const AssociatorParameters &)
math::XYZPoint trkGlobPosAtEcal
Track position at different parts of the calorimeter.
int charge() const
track electric charge
Definition: TrackBase.h:112
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
Definition: Utils.h:6
DetId geographicalId() const
edm::InputTag theEBRecHitCollectionLabel
Labels of the detector EDProducts.
double t0() const
Get the segment t0 (if recomputed, 0 is returned otherwise)
std::vector< const HORecHit * > hoRecHits
T x() const
Definition: PV3DBase.h:56
virtual LocalPoint localPosition() const =0
virtual float width() const =0
const SimTrack * simTrack
MC truth info.
void fillHO(const edm::Event &, TrackDetMatchInfo &, const AssociatorParameters &)
std::vector< SimTrack > SimTrackContainer
virtual const BoundPlane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:37
void setMaterialMode(bool noMaterial)
Switch for material effects mode: no material effects if true.
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
*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
FreeTrajectoryState initialFreeState(const reco::Track &tk, const MagneticField *field) const