SeedMatcher.cc

Go to the documentation of this file.

#include "FastSimulation/Tracking/interface/SeedMatcher.h"

#include "SimDataFormats/Track/interface/SimTrackContainer.h"
#include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
#include "TrackingTools/GeomPropagators/interface/Propagator.h"
#include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHit.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "iostream"

std::vector<int32_t> SeedMatcher::matchRecHitCombinations(
    const TrajectorySeed& seed,
    const FastTrackerRecHitCombinationCollection& recHitCombinationCollection,
    const edm::SimTrackContainer& simTrackCollection,
    double maxMatchEstimator,
    const Propagator& propagator,
    const MagneticField& magneticField,
    const TrackerGeometry& trackerGeometry) {
  // container for result
  std::vector<int32_t> result;

  // seed state
  PTrajectoryStateOnDet ptod = seed.startingState();
  DetId seedState_detId(ptod.detId());
  const GeomDet* seedState_det = trackerGeometry.idToDet(seedState_detId);
  const Surface* seedState_surface = &seedState_det->surface();
  TrajectoryStateOnSurface seedState(trajectoryStateTransform::transientState(ptod, seedState_surface, &magneticField));

  // find matches
  int nSimTracks = simTrackCollection.size();
  for (unsigned recHitCombinationIndex = 0; recHitCombinationIndex < recHitCombinationCollection.size();
       recHitCombinationIndex++) {
    const auto& recHitCombination = recHitCombinationCollection[recHitCombinationIndex];
    int simTrackIndex = recHitCombination.back()->simTrackId(0);
    if (simTrackIndex < 0 || simTrackIndex >= nSimTracks) {
      throw cms::Exception("SeedMatcher") << "SimTrack index out of range: " << simTrackIndex << std::endl;
    }
    const auto& simTrack = simTrackCollection[recHitCombination.back()->simTrackId(0)];
    double matchEstimator = matchSimTrack(seedState, simTrack, propagator, magneticField);
    if (matchEstimator < maxMatchEstimator) {
      result.push_back(recHitCombinationIndex);
    }
  }
  return result;
}

double SeedMatcher::matchSimTrack(const TrajectoryStateOnSurface& seedState,
                                  const SimTrack& simTrack,
                                  const Propagator& propagator,
                                  const MagneticField& magneticField) {
  // simtrack and simvertex at tracker surface
  GlobalPoint simTrack_atTrackerSurface_position(simTrack.trackerSurfacePosition().x(),
                                                 simTrack.trackerSurfacePosition().y(),
                                                 simTrack.trackerSurfacePosition().z());
  GlobalVector simTrack_atTrackerSurface_momentum(simTrack.trackerSurfaceMomentum().x(),
                                                  simTrack.trackerSurfaceMomentum().y(),
                                                  simTrack.trackerSurfaceMomentum().z());

  // no match if seedstate and simtrack in oposite direction
  if (simTrack_atTrackerSurface_position.basicVector().dot(simTrack_atTrackerSurface_momentum.basicVector()) *
          seedState.globalPosition().basicVector().dot(seedState.globalMomentum().basicVector()) <
      0.) {
    return 9999.;
  }

  // find simtrack state on surface of seed state
  GlobalTrajectoryParameters simTrack_atTrackerSurface_parameters(
      simTrack_atTrackerSurface_position, simTrack_atTrackerSurface_momentum, simTrack.charge(), &magneticField);
  FreeTrajectoryState simtrack_atTrackerSurface_state(simTrack_atTrackerSurface_parameters);
  TrajectoryStateOnSurface simtrack_atSeedStateSurface_state =
      propagator.propagate(simtrack_atTrackerSurface_state, seedState.surface());

  // simtrack does not propagate to surface of seed state
  if (!simtrack_atSeedStateSurface_state.isValid()) {
    return 9999.;
  }

  // simtrack and seed state have opposite direction
  if (simtrack_atSeedStateSurface_state.globalPosition().basicVector().dot(
          simtrack_atSeedStateSurface_state.globalMomentum().basicVector()) *
          seedState.globalPosition().basicVector().dot(seedState.globalMomentum().basicVector()) <
      0.) {
    return 9999.;
  }

  AlgebraicVector5 difference(seedState.localParameters().vector() -
                              simtrack_atSeedStateSurface_state.localParameters().vector());
  AlgebraicSymMatrix55 error(seedState.localError().matrix());
  if (!error.Invert()) {
    edm::LogWarning("FastSim SeedToSimTrackMatcher") << "Cannot invert seed state error matrix => no match...";
    return 9999.;
  }

  return ROOT::Math::Similarity(difference, error);
}