Tracking/src/SeedMatcher.cc

0001 #include "FastSimulation/Tracking/interface/SeedMatcher.h"
0002
0003 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0004 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0005 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0006 #include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h"
0007 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0008 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0009 #include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHit.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011
0012 #include "iostream"
0013
0014 std::vector<int32_t> SeedMatcher::matchRecHitCombinations(
0015     const TrajectorySeed& seed,
0016     const FastTrackerRecHitCombinationCollection& recHitCombinationCollection,
0017     const edm::SimTrackContainer& simTrackCollection,
0018     double maxMatchEstimator,
0019     const Propagator& propagator,
0020     const MagneticField& magneticField,
0021     const TrackerGeometry& trackerGeometry) {
0022   // container for result
0023   std::vector<int32_t> result;
0024
0025   // seed state
0026   PTrajectoryStateOnDet ptod = seed.startingState();
0027   DetId seedState_detId(ptod.detId());
0028   const GeomDet* seedState_det = trackerGeometry.idToDet(seedState_detId);
0029   const Surface* seedState_surface = &seedState_det->surface();
0030   TrajectoryStateOnSurface seedState(trajectoryStateTransform::transientState(ptod, seedState_surface, &magneticField));
0031
0032   // find matches
0033   int nSimTracks = simTrackCollection.size();
0034   for (unsigned recHitCombinationIndex = 0; recHitCombinationIndex < recHitCombinationCollection.size();
0035        recHitCombinationIndex++) {
0036     const auto& recHitCombination = recHitCombinationCollection[recHitCombinationIndex];
0037     int simTrackIndex = recHitCombination.back()->simTrackId(0);
0038     if (simTrackIndex < 0 || simTrackIndex >= nSimTracks) {
0039       throw cms::Exception("SeedMatcher") << "SimTrack index out of range: " << simTrackIndex << std::endl;
0040     }
0041     const auto& simTrack = simTrackCollection[recHitCombination.back()->simTrackId(0)];
0042     double matchEstimator = matchSimTrack(seedState, simTrack, propagator, magneticField);
0043     if (matchEstimator < maxMatchEstimator) {
0044       result.push_back(recHitCombinationIndex);
0045     }
0046   }
0047   return result;
0048 }
0049
0050 double SeedMatcher::matchSimTrack(const TrajectoryStateOnSurface& seedState,
0051                                   const SimTrack& simTrack,
0052                                   const Propagator& propagator,
0053                                   const MagneticField& magneticField) {
0054   // simtrack and simvertex at tracker surface
0055   GlobalPoint simTrack_atTrackerSurface_position(simTrack.trackerSurfacePosition().x(),
0056                                                  simTrack.trackerSurfacePosition().y(),
0057                                                  simTrack.trackerSurfacePosition().z());
0058   GlobalVector simTrack_atTrackerSurface_momentum(simTrack.trackerSurfaceMomentum().x(),
0059                                                   simTrack.trackerSurfaceMomentum().y(),
0060                                                   simTrack.trackerSurfaceMomentum().z());
0061
0062   // no match if seedstate and simtrack in oposite direction
0063   if (simTrack_atTrackerSurface_position.basicVector().dot(simTrack_atTrackerSurface_momentum.basicVector()) *
0064           seedState.globalPosition().basicVector().dot(seedState.globalMomentum().basicVector()) <
0065       0.) {
0066     return 9999.;
0067   }
0068
0069   // find simtrack state on surface of seed state
0070   GlobalTrajectoryParameters simTrack_atTrackerSurface_parameters(
0071       simTrack_atTrackerSurface_position, simTrack_atTrackerSurface_momentum, simTrack.charge(), &magneticField);
0072   FreeTrajectoryState simtrack_atTrackerSurface_state(simTrack_atTrackerSurface_parameters);
0073   TrajectoryStateOnSurface simtrack_atSeedStateSurface_state =
0074       propagator.propagate(simtrack_atTrackerSurface_state, seedState.surface());
0075
0076   // simtrack does not propagate to surface of seed state
0077   if (!simtrack_atSeedStateSurface_state.isValid()) {
0078     return 9999.;
0079   }
0080
0081   // simtrack and seed state have opposite direction
0082   if (simtrack_atSeedStateSurface_state.globalPosition().basicVector().dot(
0083           simtrack_atSeedStateSurface_state.globalMomentum().basicVector()) *
0084           seedState.globalPosition().basicVector().dot(seedState.globalMomentum().basicVector()) <
0085       0.) {
0086     return 9999.;
0087   }
0088
0089   AlgebraicVector5 difference(seedState.localParameters().vector() -
0090                               simtrack_atSeedStateSurface_state.localParameters().vector());
0091   AlgebraicSymMatrix55 error(seedState.localError().matrix());
0092   if (!error.Invert()) {
0093     edm::LogWarning("FastSim SeedToSimTrackMatcher") << "Cannot invert seed state error matrix => no match...";
0094     return 9999.;
0095   }
0096
0097   return ROOT::Math::Similarity(difference, error);
0098 }