Project CMSSW displayed by LXR CMSSW_15_0_X_2024-12-10-2300/​RecoTracker/​TrackProducer/​interface/​TrackProducerAlgorithm.icc	Source navigation Diff markup Identifier search General search
	Version: CMSSW_15_0_X_2024-12-10-2300 CMSSW_15_0_X_2024-12-09-2300 CMSSW_15_0_X_2024-12-08-2300 CMSSW_15_0_X_2024-12-06-2300 CMSSW_15_0_X_2024-12-05-2300 CMSSW_15_0_X_2024-12-04-2300 Revert   Change

File indexing completed on 2024-09-07 04:38:02

 // #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/CommonDetUnit/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;
   */
 }

This page was automatically generated by the 2.2.1 LXR engine. The LXR team