TrackProducer/interface/TrackProducerAlgorithm.icc

0001 // #include "RecoTracker/TrackProducer/interface/TrackProducerAlgorithm.h"
0002
0003 #include "DataFormats/Common/interface/OrphanHandle.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 #include "MagneticField/Engine/interface/MagneticField.h"
0007 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
0008
0009 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0010 #include "DataFormats/TrackCandidate/interface/TrackCandidate.h"
0011 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
0012
0013 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0014 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0015 #include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
0016 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0017
0018 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
0019 #include "FWCore/Utilities/interface/Exception.h"
0020
0021 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0022 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0023 #include "RecoTracker/TransientTrackingRecHit/interface/TRecHit5DParamConstraint.h"
0024 #include "RecoTracker/TransientTrackingRecHit/interface/TRecHit2DPosConstraint.h"
0025 #include "RecoTracker/TransientTrackingRecHit/interface/TRecHit1DMomConstraint.h"
0026 // #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
0027 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0028 #include "DataFormats/GeometryCommonDetAlgo/interface/ErrorFrameTransformer.h"
0029 #include "TrackingTools/TrackFitters/interface/RecHitSorter.h"
0030
0031 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0032
0033 template <class T>
0034 void TrackProducerAlgorithm<T>::runWithCandidate(const TrackingGeometry* theG,
0035                                                  const MagneticField* theMF,
0036                                                  const TrackCandidateCollection& theTCCollection,
0037                                                  const TrajectoryFitter* theFitter,
0038                                                  const Propagator* thePropagator,
0039                                                  const TransientTrackingRecHitBuilder* builder,
0040                                                  const reco::BeamSpot& bs,
0041                                                  AlgoProductCollection& algoResults) {
0042   LogDebug("TrackProducer") << "Number of TrackCandidates: " << theTCCollection.size() << "\n";
0043
0044   int cont = 0;
0045   int ntc = 0;
0046   for (auto const& theTC : theTCCollection) {
0047     PTrajectoryStateOnDet const& state = theTC.trajectoryStateOnDet();
0048     const TrajectorySeed& seed = theTC.seed();
0049
0050     //convert PTrajectoryStateOnDet to TrajectoryStateOnSurface
0051
0052     DetId detId(state.detId());
0053     TrajectoryStateOnSurface theTSOS =
0054         trajectoryStateTransform::transientState(state, &(theG->idToDet(detId)->surface()), theMF);
0055
0056     LogDebug("TrackProducer") << "Initial TSOS\n" << theTSOS << "\n";
0057
0058     //convert the TrackingRecHit vector to a TransientTrackingRecHit vector
0059     TrackingRecHit::RecHitContainer hits;
0060
0061     float ndof = 0;
0062     //  we assume are transient...
0063     for (auto const& recHit : theTC.recHits()) {
0064       hits.push_back(recHit.cloneSH());  // major waste, will be soon fitted and recloned...
0065     }
0066
0067     stopReason_ = theTC.stopReason();
0068     //build Track
0069     LogDebug("TrackProducer") << "going to buildTrack"
0070                               << "\n";
0071     FitterCloner fc(theFitter, builder);
0072     bool ok = buildTrack(
0073         fc.fitter.get(), thePropagator, algoResults, hits, theTSOS, seed, ndof, bs, theTC.seedRef(), 0, theTC.nLoops());
0074     LogDebug("TrackProducer") << "buildTrack result: " << ok << "\n";
0075     if (ok) {
0076       algoResults.back().indexInput = ntc;
0077       ++cont;
0078     }
0079     ++ntc;
0080   }
0081   LogDebug("TrackProducer") << "Number of Tracks found: " << cont << "\n";
0082   // std::cout << "VICkfProducer " << "Number of Tracks found: " << cont << std::endl;
0083 }
0084
0085 // the following are called by the Refitter(s)
0086
0087 template <class T>
0088 void TrackProducerAlgorithm<T>::runWithTrack(const TrackingGeometry* theG,
0089                                              const MagneticField* theMF,
0090                                              const TrackView& theTCollection,
0091                                              const TrajectoryFitter* theFitter,
0092                                              const Propagator* thePropagator,
0093                                              const TransientTrackingRecHitBuilder* gbuilder,
0094                                              const reco::BeamSpot& bs,
0095                                              AlgoProductCollection& algoResults) {
0096   LogDebug("TrackProducer") << "Number of input Tracks: " << theTCollection.size() << "\n";
0097   const TkTransientTrackingRecHitBuilder* builder = dynamic_cast<TkTransientTrackingRecHitBuilder const*>(gbuilder);
0098   assert(builder);
0099   int cont = 0;
0100   for (auto i = theTCollection.begin(); i != theTCollection.end(); i++) {
0101     try {
0102       const T* theT = &(*i);
0103       float ndof = 0;
0104       PropagationDirection seedDir = theT->seedDirection();
0105
0106       // WARNING: here we assume that the hits are correcly sorted according to seedDir
0107       TransientTrackingRecHit::RecHitContainer hits;
0108       for (trackingRecHit_iterator i = theT->recHitsBegin(); i != theT->recHitsEnd(); i++) {
0109         if (reMatchSplitHits_) {
0110           //re-match hits that belong together
0111           trackingRecHit_iterator next = i;
0112           next++;
0113           if (next != theT->recHitsEnd() && (*i)->isValid()) {
0114             //check whether hit and nexthit are on glued module
0115             DetId hitId = (**i).geographicalId();
0116
0117             if (hitId.det() == DetId::Tracker) {
0118               if (GeomDetEnumerators::isTrackerStrip((**i).det()->subDetector())) {
0119                 SiStripDetId stripId(hitId);
0120                 if (stripId.partnerDetId() == (*next)->geographicalId().rawId()) {
0121                   //yes they are parterns in a glued geometry.
0122                   DetId gluedId = stripId.glued();
0123                   const SiStripRecHit2D* mono = nullptr;
0124                   const SiStripRecHit2D* stereo = nullptr;
0125                   if (stripId.stereo() == 0) {
0126                     mono = dynamic_cast<const SiStripRecHit2D*>(&**i);
0127                     stereo = dynamic_cast<const SiStripRecHit2D*>(&**next);
0128                   } else {
0129                     mono = dynamic_cast<const SiStripRecHit2D*>(&**next);
0130                     stereo = dynamic_cast<const SiStripRecHit2D*>(&**i);
0131                   }
0132                   if (!mono || !stereo) {
0133                     edm::LogError("TrackProducerAlgorithm")
0134                         << "cannot get a SiStripRecHit2D from the rechit." << hitId.rawId() << " " << gluedId.rawId()
0135                         << " " << stripId.partnerDetId() << " " << (*next)->geographicalId().rawId();
0136                   }
0137                   LocalPoint pos;  //null
0138                   LocalError err;  //null
0139                   hits.push_back(std::make_shared<SiStripMatchedRecHit2D>(
0140                       pos, err, *builder->geometry()->idToDet(gluedId), mono, stereo));
0141                   //the local position and error is dummy but the fitter does not need that anyways
0142                   i++;
0143                   continue;  //go to next.
0144                 }            //consecutive hits are on parterns module.
0145               }
0146             }  //is a strip module
0147           }    //next is not the end of hits
0148         }      //if rematching option is on.
0149
0150         if ((**i).geographicalId() != 0U)
0151           hits.push_back((**i).cloneForFit(*builder->geometry()->idToDet((**i).geographicalId())));
0152       }
0153
0154       TrajectoryStateOnSurface theInitialStateForRefitting = getInitialState(theT, hits, theG, theMF);
0155
0156       // the seed has dummy state and hits.What matters for the fitting is the seedDirection;
0157       //if anyDirection the seed direction is not stored in the root file: keep same order
0158       if (seedDir == anyDirection) {
0159         throw cms::Exception("TrackProducer")
0160             << "ERROR: trying to refit a track which doesn't have a properly filled 'seed direction' data member"
0161             << std::endl;
0162       }
0163
0164       const TrajectorySeed seed({}, {}, seedDir);
0165       // =========================
0166       //LogDebug("TrackProducer") << "seed.direction()=" << seed.direction();
0167
0168       //set the algo_ member in order to propagate the old alog name
0169       algo_ = theT->algo();
0170       originalAlgo_ = theT->originalAlgo();
0171       algoMask_ = theT->algoMask();
0172
0173       stopReason_ = theT->stopReason();
0174
0175       //=====  the hits are in the same order as they were in the track::extra.
0176       FitterCloner fc(theFitter, builder);
0177       bool ok = buildTrack(fc.fitter.get(),
0178                            thePropagator,
0179                            algoResults,
0180                            hits,
0181                            theInitialStateForRefitting,
0182                            seed,
0183                            ndof,
0184                            bs,
0185                            theT->seedRef(),
0186                            theT->qualityMask(),
0187                            theT->nLoops());
0188       if (ok)
0189         cont++;
0190     } catch (cms::Exception& e) {
0191       edm::LogError("TrackProducer") << "Genexception1: " << e.explainSelf() << "\n";
0192       throw;
0193     }
0194   }
0195   LogDebug("TrackProducer") << "Number of Tracks found: " << cont << "\n";
0196 }
0197
0198 template <class T>
0199 void TrackProducerAlgorithm<T>::runWithMomentum(
0200     const TrackingGeometry* theG,
0201     const MagneticField* theMF,
0202     const TrackMomConstraintAssociationCollection& theTCollectionWithConstraint,
0203     const TrajectoryFitter* theFitter,
0204     const Propagator* thePropagator,
0205     const TransientTrackingRecHitBuilder* gbuilder,
0206     const reco::BeamSpot& bs,
0207     AlgoProductCollection& algoResults) {
0208   LogDebug("TrackProducer") << "Number of input Tracks: " << theTCollectionWithConstraint.size() << "\n";
0209   const TkTransientTrackingRecHitBuilder* builder = dynamic_cast<TkTransientTrackingRecHitBuilder const*>(gbuilder);
0210
0211   int cont = 0;
0212   for (TrackMomConstraintAssociationCollection::const_iterator i = theTCollectionWithConstraint.begin();
0213        i != theTCollectionWithConstraint.end();
0214        i++) {
0215     try {
0216       const T* theT = i->key.get();
0217
0218       LogDebug("TrackProducer") << "Running Refitter with Momentum Constraint. p=" << i->val->momentum
0219                                 << " err=" << i->val->error;
0220
0221       float ndof = 0;
0222       PropagationDirection seedDir = theT->seedDirection();
0223
0224       // WARNING: here we assume that the hits are correcly sorted according to seedDir
0225       TransientTrackingRecHit::RecHitContainer hits;
0226       for (trackingRecHit_iterator j = theT->recHitsBegin(); j != theT->recHitsEnd(); ++j) {
0227         if ((**j).geographicalId() != 0U)
0228           hits.push_back((**j).cloneForFit(*builder->geometry()->idToDet((**j).geographicalId())));
0229       }
0230
0231       TrajectoryStateOnSurface theInitialStateForRefitting = getInitialState(theT, hits, theG, theMF);
0232
0233       double mom = i->val->momentum;  //10;
0234       double err = i->val->error;     //0.01;
0235       TransientTrackingRecHit::RecHitPointer testhit = TRecHit1DMomConstraint::build(
0236           ((int)(theInitialStateForRefitting.charge())), mom, err, &theInitialStateForRefitting.surface());
0237
0238       //no insert in OwnVector...
0239       TransientTrackingRecHit::RecHitContainer tmpHits;
0240       tmpHits.push_back(testhit);
0241       for (TransientTrackingRecHit::RecHitContainer::const_iterator i = hits.begin(); i != hits.end(); i++) {
0242         tmpHits.push_back(*i);
0243       }
0244       hits.swap(tmpHits);
0245
0246       // the seed has dummy state and hits.What matters for the fitting is the seedDirection;
0247       const TrajectorySeed seed({}, {}, seedDir);
0248       // =========================
0249       //LogDebug("TrackProducer") << "seed.direction()=" << seed.direction();
0250
0251       //=====  the hits are in the same order as they were in the track::extra.
0252       FitterCloner fc(theFitter, builder);
0253       bool ok = buildTrack(fc.fitter.get(),
0254                            thePropagator,
0255                            algoResults,
0256                            hits,
0257                            theInitialStateForRefitting,
0258                            seed,
0259                            ndof,
0260                            bs,
0261                            theT->seedRef(),
0262                            theT->qualityMask(),
0263                            theT->nLoops());
0264       if (ok)
0265         cont++;
0266     } catch (cms::Exception& e) {
0267       edm::LogError("TrackProducer") << "cms::Exception: " << e.explainSelf() << "\n";
0268       throw;
0269     }
0270   }
0271   LogDebug("TrackProducer") << "Number of Tracks found: " << cont << "\n";
0272 }
0273
0274 template <class T>
0275 void TrackProducerAlgorithm<T>::runWithTrackParameters(
0276     const TrackingGeometry* theG,
0277     const MagneticField* theMF,
0278     const TrackParamConstraintAssociationCollection& theTCollectionWithConstraint,
0279     const TrajectoryFitter* theFitter,
0280     const Propagator* thePropagator,
0281     const TransientTrackingRecHitBuilder* gbuilder,
0282     const reco::BeamSpot& bs,
0283     AlgoProductCollection& algoResults) {
0284   LogDebug("TrackProducer") << "Number of input Tracks: " << theTCollectionWithConstraint.size() << "\n";
0285   const TkTransientTrackingRecHitBuilder* builder = dynamic_cast<TkTransientTrackingRecHitBuilder const*>(gbuilder);
0286
0287   int cont = 0;
0288
0289   for (typename TrackParamConstraintAssociationCollection::const_iterator i = theTCollectionWithConstraint.begin();
0290        i != theTCollectionWithConstraint.end();
0291        i++) {
0292     try {
0293       const T* theT = i->key.get();
0294
0295       LogDebug("TrackProducer") << "Running Refitter with Track Parameter Constraint. TSOS = " << *(i->val);
0296
0297       TransientTrackingRecHit::RecHitPointer constraintHit = TRecHit5DParamConstraint::build(*(i->val));
0298
0299       // WARNING: here we assume that the hits are correcly sorted according to seedDir
0300       TransientTrackingRecHit::RecHitContainer hits;
0301       hits.push_back(constraintHit);
0302       for (trackingRecHit_iterator j = theT->recHitsBegin(); j != theT->recHitsEnd(); ++j) {
0303         if ((**j).geographicalId() != 0U)
0304           hits.push_back((**j).cloneForFit(*builder->geometry()->idToDet((**j).geographicalId())));
0305       }
0306
0307       float ndof = 0;
0308       PropagationDirection seedDir = theT->seedDirection();
0309
0310       // the best guess for the initial state for refitting is obviously the constraint tsos itself
0311       TrajectoryStateOnSurface theInitialStateForRefitting = *(i->val);
0312
0313       // the seed has dummy state and hits.What matters for the fitting is the seedDirection;
0314       const TrajectorySeed seed({}, {}, seedDir);
0315
0316       //=====  the hits are in the same order as they were in the track::extra.
0317       FitterCloner fc(theFitter, builder);
0318       bool ok = buildTrack(fc.fitter.get(),
0319                            thePropagator,
0320                            algoResults,
0321                            hits,
0322                            theInitialStateForRefitting,
0323                            seed,
0324                            ndof,
0325                            bs,
0326                            theT->seedRef(),
0327                            theT->qualityMask(),
0328                            theT->nLoops());
0329       if (ok)
0330         cont++;
0331     } catch (cms::Exception& e) {
0332       edm::LogError("TrackProducer") << "cms::Exception: " << e.explainSelf() << "\n";
0333       throw;
0334     }
0335   }
0336   LogDebug("TrackProducer") << "Number of Tracks found: " << cont << "\n";
0337 }
0338
0339 template <class T>
0340 void TrackProducerAlgorithm<T>::runWithVertex(const TrackingGeometry* theG,
0341                                               const MagneticField* theMF,
0342                                               const VtxConstraintAssociationCollection& theTCollectionWithConstraint,
0343                                               const TrajectoryFitter* theFitter,
0344                                               const Propagator* thePropagator,
0345                                               const TransientTrackingRecHitBuilder* gbuilder,
0346                                               const reco::BeamSpot& bs,
0347                                               AlgoProductCollection& algoResults) {
0348   const TkTransientTrackingRecHitBuilder* builder = dynamic_cast<TkTransientTrackingRecHitBuilder const*>(gbuilder);
0349
0350   LogDebug("TrackProducer") << "Number of input Tracks: " << theTCollectionWithConstraint.size() << "\n";
0351
0352   //   PropagatorWithMaterial myPropagator(anyDirection,0.105,theMF);
0353   AnalyticalPropagator myPropagator(theMF, anyDirection);  // no material surface inside 1st hit
0354   TransverseImpactPointExtrapolator extrapolator(myPropagator);
0355
0356   int cont = 0;
0357   for (typename VtxConstraintAssociationCollection::const_iterator i = theTCollectionWithConstraint.begin();
0358        i != theTCollectionWithConstraint.end();
0359        i++) {
0360     try {
0361       const T* theT = i->key.get();
0362
0363       LogDebug("TrackProducer") << "Running Refitter with Vertex Constraint. pos=" << i->val->first
0364                                 << " err=" << i->val->second.matrix();
0365
0366       float ndof = 0;
0367       PropagationDirection seedDir = theT->seedDirection();
0368
0369       // FreeTrajectoryState from track position and momentum
0370       GlobalPoint tkpos(theT->vertex().x(), theT->vertex().y(), theT->vertex().z());
0371       GlobalVector tkmom(theT->momentum().x(), theT->momentum().y(), theT->momentum().z());
0372       FreeTrajectoryState fts(tkpos, tkmom, theT->charge(), theMF);
0373       // Extrapolation to transverse IP plane at vertex constraint position
0374       GlobalPoint pos = i->val->first;
0375       GlobalError err = i->val->second;
0376       TrajectoryStateOnSurface tsosAtVtx = extrapolator.extrapolate(fts, pos);
0377       const Surface& surfAtVtx = tsosAtVtx.surface();
0378       // Conversion of vertex constraint parameters and errors to local system
0379       //   (assumes projection normal to plane - sufficient if transverse components are small)
0380       LocalPoint vtxPosL = surfAtVtx.toLocal(pos);
0381       LocalError vtxErrL = ErrorFrameTransformer().transform(err, surfAtVtx);
0382       // Construction of the virtual hit
0383       TransientTrackingRecHit::RecHitPointer vtxhit = TRecHit2DPosConstraint::build(vtxPosL, vtxErrL, &surfAtVtx);
0384
0385       // WARNING: here we assume that the hits are correcly sorted according to seedDir
0386       TransientTrackingRecHit::RecHitContainer hits;
0387       for (trackingRecHit_iterator j = theT->recHitsBegin(); j != theT->recHitsEnd(); ++j) {
0388         if ((**j).geographicalId() != 0U)
0389           hits.push_back((**j).cloneForFit(*builder->geometry()->idToDet((**j).geographicalId())));
0390       }
0391
0392       TrajectoryStateOnSurface theInitialStateForRefitting = getInitialState(theT, hits, theG, theMF);
0393
0394       //    const LocalPoint testpoint(0,0,0);
0395       //    GlobalPoint pos = i->val->first;
0396       //    GlobalError err = i->val->second;
0397
0398       //    Propagator* myPropagator = new PropagatorWithMaterial(anyDirection,0.105,theMF);
0399       //    TransverseImpactPointExtrapolator extrapolator(*myPropagator);
0400       //    TrajectoryStateOnSurface tsosAtVtx = extrapolator.extrapolate(theInitialStateForRefitting,pos);
0401
0402       //    const Surface * surfAtVtx = &tsosAtVtx.surface();
0403
0404       //    LocalError testerror = ErrorFrameTransformer().transform(err, *surfAtVtx);
0405
0406       //    TransientTrackingRecHit::RecHitPointer testhit = TRecHit2DPosConstraint::build(testpoint,testerror,surfAtVtx);
0407
0408       //push constraining hit and sort along seed direction
0409       //    hits.push_back(testhit);
0410       //    RecHitSorter sorter = RecHitSorter();
0411       //    hits = sorter.sortHits(hits,seedDir);   // warning: I doubt it works for cosmic or beam halo tracks
0412       //insertion of the constraint in the lists of hits
0413       //  (assume hits ordered from inside to outside)
0414       hits.insert(hits.begin(), vtxhit);
0415
0416       //use the state on the surface of the first hit (could be the constraint or not)
0417       theInitialStateForRefitting = myPropagator.propagate(theInitialStateForRefitting, *(hits[0]->surface()));
0418
0419       // the seed has dummy state and hits.What matters for the fitting is the seedDirection;
0420       const TrajectorySeed seed({}, {}, seedDir);
0421       // =========================
0422       //LogDebug("TrackProducer") << "seed.direction()=" << seed.direction();
0423
0424       //=====  the hits are in the same order as they were in the track::extra.
0425       FitterCloner fc(theFitter, builder);
0426       bool ok = buildTrack(fc.fitter.get(),
0427                            thePropagator,
0428                            algoResults,
0429                            hits,
0430                            theInitialStateForRefitting,
0431                            seed,
0432                            ndof,
0433                            bs,
0434                            theT->seedRef(),
0435                            theT->qualityMask(),
0436                            theT->nLoops());
0437       if (ok)
0438         cont++;
0439     } catch (cms::Exception& e) {
0440       edm::LogError("TrackProducer") << "cms::Exception: " << e.explainSelf() << "\n";
0441       throw;
0442     }
0443   }
0444   LogDebug("TrackProducer") << "Number of Tracks found: " << cont << "\n";
0445 }
0446
0447 template <class T>
0448 TrajectoryStateOnSurface TrackProducerAlgorithm<T>::getInitialState(const T* theT,
0449                                                                     TransientTrackingRecHit::RecHitContainer& hits,
0450                                                                     const TrackingGeometry* theG,
0451                                                                     const MagneticField* theMF) {
0452   TrajectoryStateOnSurface theInitialStateForRefitting;
0453   //the starting state is the state closest to the first hit along seedDirection.
0454   //avoiding to use transientTrack, it should be faster;
0455   TrajectoryStateOnSurface innerStateFromTrack = trajectoryStateTransform::innerStateOnSurface(*theT, *theG, theMF);
0456   TrajectoryStateOnSurface outerStateFromTrack = trajectoryStateTransform::outerStateOnSurface(*theT, *theG, theMF);
0457   TrajectoryStateOnSurface initialStateFromTrack =
0458       ((innerStateFromTrack.globalPosition() - hits.front()->det()->position()).mag2() <
0459        (outerStateFromTrack.globalPosition() - hits.front()->det()->position()).mag2())
0460           ? innerStateFromTrack
0461           : outerStateFromTrack;
0462
0463   // error is rescaled, but correlation are kept.
0464   initialStateFromTrack.rescaleError(100);
0465   return initialStateFromTrack;
0466   /*
0467   theInitialStateForRefitting = TrajectoryStateOnSurface(initialStateFromTrack.localParameters(),
0468                              initialStateFromTrack.localError(),
0469                              initialStateFromTrack.surface(),
0470                              theMF);
0471   return theInitialStateForRefitting;
0472   */
0473 }