RecoTracker/TrackProducer/interface/TrackProducerAlgorithm.icc - master - githubmirror/cmssw

RecoTracker/TrackProducer/interface/TrackProducerAlgorithm.icc
473 строки22 KB
Shahzad Malik Muzaffar
[RECONSTRUCTION] Drop Geometry/CommonDetUnit package
28 май 2026, 12:57
28 май 2026, 12:57685f036
// #include "RecoTracker/TrackProducer/interface/TrackProducerAlgorithm.h"

#include "DataFormats/Common/interface/OrphanHandle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "MagneticField/Engine/interface/MagneticField.h"
#include "Geometry/CommonTopologies/interface/TrackingGeometry.h"

#include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
#include "DataFormats/TrackCandidate/interface/TrackCandidate.h"
#include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"

#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
#include "TrackingTools/TransientTrack/interface/TransientTrack.h"

#include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
#include "FWCore/Utilities/interface/Exception.h"

#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
#include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
#include "RecoTracker/TransientTrackingRecHit/interface/TRecHit5DParamConstraint.h"
#include "RecoTracker/TransientTrackingRecHit/interface/TRecHit2DPosConstraint.h"
#include "RecoTracker/TransientTrackingRecHit/interface/TRecHit1DMomConstraint.h"
// #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
#include "DataFormats/GeometryCommonDetAlgo/interface/ErrorFrameTransformer.h"
#include "TrackingTools/TrackFitters/interface/RecHitSorter.h"

#include "DataFormats/SiStripDetId/interface/SiStripDetId.h"

template <class T>
void TrackProducerAlgorithm<T>::runWithCandidate(const TrackingGeometry* theG,
                                                 const MagneticField* theMF,
                                                 const TrackCandidateCollection& theTCCollection,
                                                 const TrajectoryFitter* theFitter,
                                                 const Propagator* thePropagator,
                                                 const TransientTrackingRecHitBuilder* builder,
                                                 const reco::BeamSpot& bs,
                                                 AlgoProductCollection& algoResults) {
  LogDebug("TrackProducer") << "Number of TrackCandidates: " << theTCCollection.size() << "\n";

  int cont = 0;
  int ntc = 0;
  for (auto const& theTC : theTCCollection) {
    PTrajectoryStateOnDet const& state = theTC.trajectoryStateOnDet();
    const TrajectorySeed& seed = theTC.seed();

    //convert PTrajectoryStateOnDet to TrajectoryStateOnSurface

    DetId detId(state.detId());
    TrajectoryStateOnSurface theTSOS =
        trajectoryStateTransform::transientState(state, &(theG->idToDet(detId)->surface()), theMF);

    LogDebug("TrackProducer") << "Initial TSOS\n" << theTSOS << "\n";

    //convert the TrackingRecHit vector to a TransientTrackingRecHit vector
    TrackingRecHit::RecHitContainer hits;

    float ndof = 0;
    //  we assume are transient...
    for (auto const& recHit : theTC.recHits()) {
      hits.push_back(recHit.cloneSH());  // major waste, will be soon fitted and recloned...
    }

    stopReason_ = theTC.stopReason();
    //build Track
    LogDebug("TrackProducer") << "going to buildTrack"
                              << "\n";
    FitterCloner fc(theFitter, builder);
    bool ok = buildTrack(
        fc.fitter.get(), thePropagator, algoResults, hits, theTSOS, seed, ndof, bs, theTC.seedRef(), 0, theTC.nLoops());
    LogDebug("TrackProducer") << "buildTrack result: " << ok << "\n";
    if (ok) {
      algoResults.back().indexInput = ntc;
      ++cont;
    }
    ++ntc;
  }
  LogDebug("TrackProducer") << "Number of Tracks found: " << cont << "\n";
  // std::cout << "VICkfProducer " << "Number of Tracks found: " << cont << std::endl;
}

// the following are called by the Refitter(s)

template <class T>
void TrackProducerAlgorithm<T>::runWithTrack(const TrackingGeometry* theG,
                                             const MagneticField* theMF,
                                             const TrackView& theTCollection,
                                             const TrajectoryFitter* theFitter,
                                             const Propagator* thePropagator,
                                             const TransientTrackingRecHitBuilder* gbuilder,
                                             const reco::BeamSpot& bs,
                                             AlgoProductCollection& algoResults) {
  LogDebug("TrackProducer") << "Number of input Tracks: " << theTCollection.size() << "\n";
  const TkTransientTrackingRecHitBuilder* builder = dynamic_cast<TkTransientTrackingRecHitBuilder const*>(gbuilder);
  assert(builder);
  int cont = 0;
  for (auto i = theTCollection.begin(); i != theTCollection.end(); i++) {
    try {
      const T* theT = &(*i);
      float ndof = 0;
      PropagationDirection seedDir = theT->seedDirection();

      // WARNING: here we assume that the hits are correcly sorted according to seedDir
      TransientTrackingRecHit::RecHitContainer hits;
      for (trackingRecHit_iterator i = theT->recHitsBegin(); i != theT->recHitsEnd(); i++) {
        if (reMatchSplitHits_) {
          //re-match hits that belong together
          trackingRecHit_iterator next = i;
          next++;
          if (next != theT->recHitsEnd() && (*i)->isValid()) {
            //check whether hit and nexthit are on glued module
            DetId hitId = (**i).geographicalId();

            if (hitId.det() == DetId::Tracker) {
              if (GeomDetEnumerators::isTrackerStrip((**i).det()->subDetector())) {
                SiStripDetId stripId(hitId);
                if (stripId.partnerDetId() == (*next)->geographicalId().rawId()) {
                  //yes they are parterns in a glued geometry.
                  DetId gluedId = stripId.glued();
                  const SiStripRecHit2D* mono = nullptr;
                  const SiStripRecHit2D* stereo = nullptr;
                  if (stripId.stereo() == 0) {
                    mono = dynamic_cast<const SiStripRecHit2D*>(&**i);
                    stereo = dynamic_cast<const SiStripRecHit2D*>(&**next);
                  } else {
                    mono = dynamic_cast<const SiStripRecHit2D*>(&**next);
                    stereo = dynamic_cast<const SiStripRecHit2D*>(&**i);
                  }
                  if (!mono || !stereo) {
                    edm::LogError("TrackProducerAlgorithm")
                        << "cannot get a SiStripRecHit2D from the rechit." << hitId.rawId() << " " << gluedId.rawId()
                        << " " << stripId.partnerDetId() << " " << (*next)->geographicalId().rawId();
                  }
                  LocalPoint pos;  //null
                  LocalError err;  //null
                  hits.push_back(std::make_shared<SiStripMatchedRecHit2D>(
                      pos, err, *builder->geometry()->idToDet(gluedId), mono, stereo));
                  //the local position and error is dummy but the fitter does not need that anyways
                  i++;
                  continue;  //go to next.
                }  //consecutive hits are on parterns module.
              }
            }  //is a strip module
          }  //next is not the end of hits
        }  //if rematching option is on.

        if ((**i).geographicalId() != 0U)
          hits.push_back((**i).cloneForFit(*builder->geometry()->idToDet((**i).geographicalId())));
      }

      TrajectoryStateOnSurface theInitialStateForRefitting = getInitialState(theT, hits, theG, theMF);

      // the seed has dummy state and hits.What matters for the fitting is the seedDirection;
      //if anyDirection the seed direction is not stored in the root file: keep same order
      if (seedDir == anyDirection) {
        throw cms::Exception("TrackProducer")
            << "ERROR: trying to refit a track which doesn't have a properly filled 'seed direction' data member"
            << std::endl;
      }

      const TrajectorySeed seed({}, {}, seedDir);
      // =========================
      //LogDebug("TrackProducer") << "seed.direction()=" << seed.direction();

      //set the algo_ member in order to propagate the old alog name
      algo_ = theT->algo();
      originalAlgo_ = theT->originalAlgo();
      algoMask_ = theT->algoMask();

      stopReason_ = theT->stopReason();

      //=====  the hits are in the same order as they were in the track::extra.
      FitterCloner fc(theFitter, builder);
      bool ok = buildTrack(fc.fitter.get(),
                           thePropagator,
                           algoResults,
                           hits,
                           theInitialStateForRefitting,
                           seed,
                           ndof,
                           bs,
                           theT->seedRef(),
                           theT->qualityMask(),
                           theT->nLoops());
      if (ok)
        cont++;
    } catch (cms::Exception& e) {
      edm::LogError("TrackProducer") << "Genexception1: " << e.explainSelf() << "\n";
      throw;
    }
  }
  LogDebug("TrackProducer") << "Number of Tracks found: " << cont << "\n";
}

template <class T>
void TrackProducerAlgorithm<T>::runWithMomentum(
    const TrackingGeometry* theG,
    const MagneticField* theMF,
    const TrackMomConstraintAssociationCollection& theTCollectionWithConstraint,
    const TrajectoryFitter* theFitter,
    const Propagator* thePropagator,
    const TransientTrackingRecHitBuilder* gbuilder,
    const reco::BeamSpot& bs,
    AlgoProductCollection& algoResults) {
  LogDebug("TrackProducer") << "Number of input Tracks: " << theTCollectionWithConstraint.size() << "\n";
  const TkTransientTrackingRecHitBuilder* builder = dynamic_cast<TkTransientTrackingRecHitBuilder const*>(gbuilder);

  int cont = 0;
  for (TrackMomConstraintAssociationCollection::const_iterator i = theTCollectionWithConstraint.begin();
       i != theTCollectionWithConstraint.end();
       i++) {
    try {
      const T* theT = i->key.get();

      LogDebug("TrackProducer") << "Running Refitter with Momentum Constraint. p=" << i->val->momentum
                                << " err=" << i->val->error;

      float ndof = 0;
      PropagationDirection seedDir = theT->seedDirection();

      // WARNING: here we assume that the hits are correcly sorted according to seedDir
      TransientTrackingRecHit::RecHitContainer hits;
      for (trackingRecHit_iterator j = theT->recHitsBegin(); j != theT->recHitsEnd(); ++j) {
        if ((**j).geographicalId() != 0U)
          hits.push_back((**j).cloneForFit(*builder->geometry()->idToDet((**j).geographicalId())));
      }

      TrajectoryStateOnSurface theInitialStateForRefitting = getInitialState(theT, hits, theG, theMF);

      double mom = i->val->momentum;  //10;
      double err = i->val->error;     //0.01;
      TransientTrackingRecHit::RecHitPointer testhit = TRecHit1DMomConstraint::build(
          ((int)(theInitialStateForRefitting.charge())), mom, err, &theInitialStateForRefitting.surface());

      //no insert in OwnVector...
      TransientTrackingRecHit::RecHitContainer tmpHits;
      tmpHits.push_back(testhit);
      for (TransientTrackingRecHit::RecHitContainer::const_iterator i = hits.begin(); i != hits.end(); i++) {
        tmpHits.push_back(*i);
      }
      hits.swap(tmpHits);

      // the seed has dummy state and hits.What matters for the fitting is the seedDirection;
      const TrajectorySeed seed({}, {}, seedDir);
      // =========================
      //LogDebug("TrackProducer") << "seed.direction()=" << seed.direction();

      //=====  the hits are in the same order as they were in the track::extra.
      FitterCloner fc(theFitter, builder);
      bool ok = buildTrack(fc.fitter.get(),
                           thePropagator,
                           algoResults,
                           hits,
                           theInitialStateForRefitting,
                           seed,
                           ndof,
                           bs,
                           theT->seedRef(),
                           theT->qualityMask(),
                           theT->nLoops());
      if (ok)
        cont++;
    } catch (cms::Exception& e) {
      edm::LogError("TrackProducer") << "cms::Exception: " << e.explainSelf() << "\n";
      throw;
    }
  }
  LogDebug("TrackProducer") << "Number of Tracks found: " << cont << "\n";
}

template <class T>
void TrackProducerAlgorithm<T>::runWithTrackParameters(
    const TrackingGeometry* theG,
    const MagneticField* theMF,
    const TrackParamConstraintAssociationCollection& theTCollectionWithConstraint,
    const TrajectoryFitter* theFitter,
    const Propagator* thePropagator,
    const TransientTrackingRecHitBuilder* gbuilder,
    const reco::BeamSpot& bs,
    AlgoProductCollection& algoResults) {
  LogDebug("TrackProducer") << "Number of input Tracks: " << theTCollectionWithConstraint.size() << "\n";
  const TkTransientTrackingRecHitBuilder* builder = dynamic_cast<TkTransientTrackingRecHitBuilder const*>(gbuilder);

  int cont = 0;

  for (typename TrackParamConstraintAssociationCollection::const_iterator i = theTCollectionWithConstraint.begin();
       i != theTCollectionWithConstraint.end();
       i++) {
    try {
      const T* theT = i->key.get();

      LogDebug("TrackProducer") << "Running Refitter with Track Parameter Constraint. TSOS = " << *(i->val);

      TransientTrackingRecHit::RecHitPointer constraintHit = TRecHit5DParamConstraint::build(*(i->val));

      // WARNING: here we assume that the hits are correcly sorted according to seedDir
      TransientTrackingRecHit::RecHitContainer hits;
      hits.push_back(constraintHit);
      for (trackingRecHit_iterator j = theT->recHitsBegin(); j != theT->recHitsEnd(); ++j) {
        if ((**j).geographicalId() != 0U)
          hits.push_back((**j).cloneForFit(*builder->geometry()->idToDet((**j).geographicalId())));
      }

      float ndof = 0;
      PropagationDirection seedDir = theT->seedDirection();

      // the best guess for the initial state for refitting is obviously the constraint tsos itself
      TrajectoryStateOnSurface theInitialStateForRefitting = *(i->val);

      // the seed has dummy state and hits.What matters for the fitting is the seedDirection;
      const TrajectorySeed seed({}, {}, seedDir);

      //=====  the hits are in the same order as they were in the track::extra.
      FitterCloner fc(theFitter, builder);
      bool ok = buildTrack(fc.fitter.get(),
                           thePropagator,
                           algoResults,
                           hits,
                           theInitialStateForRefitting,
                           seed,
                           ndof,
                           bs,
                           theT->seedRef(),
                           theT->qualityMask(),
                           theT->nLoops());
      if (ok)
        cont++;
    } catch (cms::Exception& e) {
      edm::LogError("TrackProducer") << "cms::Exception: " << e.explainSelf() << "\n";
      throw;
    }
  }
  LogDebug("TrackProducer") << "Number of Tracks found: " << cont << "\n";
}

template <class T>
void TrackProducerAlgorithm<T>::runWithVertex(const TrackingGeometry* theG,
                                              const MagneticField* theMF,
                                              const VtxConstraintAssociationCollection& theTCollectionWithConstraint,
                                              const TrajectoryFitter* theFitter,
                                              const Propagator* thePropagator,
                                              const TransientTrackingRecHitBuilder* gbuilder,
                                              const reco::BeamSpot& bs,
                                              AlgoProductCollection& algoResults) {
  const TkTransientTrackingRecHitBuilder* builder = dynamic_cast<TkTransientTrackingRecHitBuilder const*>(gbuilder);

  LogDebug("TrackProducer") << "Number of input Tracks: " << theTCollectionWithConstraint.size() << "\n";

  //   PropagatorWithMaterial myPropagator(anyDirection,0.105,theMF);
  AnalyticalPropagator myPropagator(theMF, anyDirection);  // no material surface inside 1st hit
  TransverseImpactPointExtrapolator extrapolator(myPropagator);

  int cont = 0;
  for (typename VtxConstraintAssociationCollection::const_iterator i = theTCollectionWithConstraint.begin();
       i != theTCollectionWithConstraint.end();
       i++) {
    try {
      const T* theT = i->key.get();

      LogDebug("TrackProducer") << "Running Refitter with Vertex Constraint. pos=" << i->val->first
                                << " err=" << i->val->second.matrix();

      float ndof = 0;
      PropagationDirection seedDir = theT->seedDirection();

      // FreeTrajectoryState from track position and momentum
      GlobalPoint tkpos(theT->vertex().x(), theT->vertex().y(), theT->vertex().z());
      GlobalVector tkmom(theT->momentum().x(), theT->momentum().y(), theT->momentum().z());
      FreeTrajectoryState fts(tkpos, tkmom, theT->charge(), theMF);
      // Extrapolation to transverse IP plane at vertex constraint position
      GlobalPoint pos = i->val->first;
      GlobalError err = i->val->second;
      TrajectoryStateOnSurface tsosAtVtx = extrapolator.extrapolate(fts, pos);
      const Surface& surfAtVtx = tsosAtVtx.surface();
      // Conversion of vertex constraint parameters and errors to local system
      //   (assumes projection normal to plane - sufficient if transverse components are small)
      LocalPoint vtxPosL = surfAtVtx.toLocal(pos);
      LocalError vtxErrL = ErrorFrameTransformer().transform(err, surfAtVtx);
      // Construction of the virtual hit
      TransientTrackingRecHit::RecHitPointer vtxhit = TRecHit2DPosConstraint::build(vtxPosL, vtxErrL, &surfAtVtx);

      // WARNING: here we assume that the hits are correcly sorted according to seedDir
      TransientTrackingRecHit::RecHitContainer hits;
      for (trackingRecHit_iterator j = theT->recHitsBegin(); j != theT->recHitsEnd(); ++j) {
        if ((**j).geographicalId() != 0U)
          hits.push_back((**j).cloneForFit(*builder->geometry()->idToDet((**j).geographicalId())));
      }

      TrajectoryStateOnSurface theInitialStateForRefitting = getInitialState(theT, hits, theG, theMF);

      // 	const LocalPoint testpoint(0,0,0);
      // 	GlobalPoint pos = i->val->first;
      // 	GlobalError err = i->val->second;

      // 	Propagator* myPropagator = new PropagatorWithMaterial(anyDirection,0.105,theMF);
      // 	TransverseImpactPointExtrapolator extrapolator(*myPropagator);
      // 	TrajectoryStateOnSurface tsosAtVtx = extrapolator.extrapolate(theInitialStateForRefitting,pos);

      // 	const Surface * surfAtVtx = &tsosAtVtx.surface();

      // 	LocalError testerror = ErrorFrameTransformer().transform(err, *surfAtVtx);

      // 	TransientTrackingRecHit::RecHitPointer testhit = TRecHit2DPosConstraint::build(testpoint,testerror,surfAtVtx);

      //push constraining hit and sort along seed direction
      // 	hits.push_back(testhit);
      // 	RecHitSorter sorter = RecHitSorter();
      // 	hits = sorter.sortHits(hits,seedDir);   // warning: I doubt it works for cosmic or beam halo tracks
      //insertion of the constraint in the lists of hits
      //  (assume hits ordered from inside to outside)
      hits.insert(hits.begin(), vtxhit);

      //use the state on the surface of the first hit (could be the constraint or not)
      theInitialStateForRefitting = myPropagator.propagate(theInitialStateForRefitting, *(hits[0]->surface()));

      // the seed has dummy state and hits.What matters for the fitting is the seedDirection;
      const TrajectorySeed seed({}, {}, seedDir);
      // =========================
      //LogDebug("TrackProducer") << "seed.direction()=" << seed.direction();

      //=====  the hits are in the same order as they were in the track::extra.
      FitterCloner fc(theFitter, builder);
      bool ok = buildTrack(fc.fitter.get(),
                           thePropagator,
                           algoResults,
                           hits,
                           theInitialStateForRefitting,
                           seed,
                           ndof,
                           bs,
                           theT->seedRef(),
                           theT->qualityMask(),
                           theT->nLoops());
      if (ok)
        cont++;
    } catch (cms::Exception& e) {
      edm::LogError("TrackProducer") << "cms::Exception: " << e.explainSelf() << "\n";
      throw;
    }
  }
  LogDebug("TrackProducer") << "Number of Tracks found: " << cont << "\n";
}

template <class T>
TrajectoryStateOnSurface TrackProducerAlgorithm<T>::getInitialState(const T* theT,
                                                                    TransientTrackingRecHit::RecHitContainer& hits,
                                                                    const TrackingGeometry* theG,
                                                                    const MagneticField* theMF) {
  TrajectoryStateOnSurface theInitialStateForRefitting;
  //the starting state is the state closest to the first hit along seedDirection.
  //avoiding to use transientTrack, it should be faster;
  TrajectoryStateOnSurface innerStateFromTrack = trajectoryStateTransform::innerStateOnSurface(*theT, *theG, theMF);
  TrajectoryStateOnSurface outerStateFromTrack = trajectoryStateTransform::outerStateOnSurface(*theT, *theG, theMF);
  TrajectoryStateOnSurface initialStateFromTrack =
      ((innerStateFromTrack.globalPosition() - hits.front()->det()->position()).mag2() <
       (outerStateFromTrack.globalPosition() - hits.front()->det()->position()).mag2())
          ? innerStateFromTrack
          : outerStateFromTrack;

  // error is rescaled, but correlation are kept.
  initialStateFromTrack.rescaleError(100);
  return initialStateFromTrack;
  /*
  theInitialStateForRefitting = TrajectoryStateOnSurface(initialStateFromTrack.localParameters(),
							 initialStateFromTrack.localError(), 		      
							 initialStateFromTrack.surface(),
							 theMF);
  return theInitialStateForRefitting;
  */
}