File indexing completed on 2024-04-06 12:11:23
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
0023 std::vector<int32_t> result;
0024
0025
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
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
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
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
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
0077 if (!simtrack_atSeedStateSurface_state.isValid()) {
0078 return 9999.;
0079 }
0080
0081
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 }