GsfTracking/src/GsfTrajectorySmoother.cc

0001 #include "TrackingTools/GsfTracking/interface/GsfTrajectorySmoother.h"
0002
0003 #include "TrackingTools/GsfTracking/interface/MultiTrajectoryStateMerger.h"
0004 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0005 #include "DataFormats/GeometrySurface/interface/BoundPlane.h"
0006 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008
0009 #include "TrackingTools/TrackFitters/interface/DebugHelpers.h"
0010
0011 GsfTrajectorySmoother::GsfTrajectorySmoother(const GsfPropagatorWithMaterial& aPropagator,
0012                                              const TrajectoryStateUpdator& aUpdator,
0013                                              const MeasurementEstimator& aEstimator,
0014                                              const MultiTrajectoryStateMerger& aMerger,
0015                                              float errorRescaling,
0016                                              const bool materialBeforeUpdate,
0017                                              const DetLayerGeometry* detLayerGeometry)
0018     : theAlongPropagator(nullptr),
0019       theOppositePropagator(nullptr),
0020       theGeomPropagator(nullptr),
0021       theConvolutor(nullptr),
0022       theUpdator(aUpdator.clone()),
0023       theEstimator(aEstimator.clone()),
0024       theMerger(aMerger.clone()),
0025       theMatBeforeUpdate(materialBeforeUpdate),
0026       theErrorRescaling(errorRescaling),
0027       theGeometry(detLayerGeometry) {
0028   auto p = aPropagator.clone();
0029   p->setPropagationDirection(alongMomentum);
0030   theAlongPropagator = p;
0031   p = aPropagator.clone();
0032   p->setPropagationDirection(oppositeToMomentum);
0033   theOppositePropagator = p;
0034   if (!theMatBeforeUpdate) {
0035     theGeomPropagator = new GsfPropagatorAdapter(aPropagator.geometricalPropagator());
0036     theConvolutor = aPropagator.convolutionWithMaterial().clone();
0037   }
0038
0039   if (!theGeometry)
0040     theGeometry = &dummyGeometry;
0041 }
0042
0043 GsfTrajectorySmoother::~GsfTrajectorySmoother() {
0044   delete theAlongPropagator;
0045   delete theOppositePropagator;
0046   delete theGeomPropagator;
0047   delete theConvolutor;
0048   delete theUpdator;
0049   delete theEstimator;
0050   delete theMerger;
0051 }
0052
0053 Trajectory GsfTrajectorySmoother::trajectory(const Trajectory& aTraj) const {
0054   if (aTraj.empty())
0055     return Trajectory();
0056
0057   const Propagator* usePropagator = theAlongPropagator;
0058   if (aTraj.direction() == alongMomentum) {
0059     usePropagator = theOppositePropagator;
0060   }
0061   if (not usePropagator) {
0062     usePropagator = theGeomPropagator;
0063   }
0064
0065   Trajectory myTraj(aTraj.seed(), usePropagator->propagationDirection());
0066
0067   std::vector<TM> const& avtm = aTraj.measurements();
0068
0069   TSOS predTsos = avtm.back().forwardPredictedState();
0070   predTsos.rescaleError(theErrorRescaling);
0071
0072   if (!predTsos.isValid()) {
0073     edm::LogInfo("GsfTrackFitters") << "GsfTrajectorySmoother: predicted tsos of last measurement not valid!";
0074     return Trajectory();
0075   }
0076   TSOS currTsos;
0077
0078   //first smoothed tm is last fitted
0079   if (avtm.back().recHit()->isValid()) {
0080     {
0081       //       TimeMe t(*updateTimer,false);
0082       currTsos = updator()->update(predTsos, *avtm.back().recHit());
0083     }
0084
0085     if (!currTsos.isValid()) {
0086       edm::LogInfo("GsfTrajectorySmoother") << "GsfTrajectorySmoother: tsos not valid after update!";
0087       return Trajectory();
0088     }
0089
0090     //check validity
0091     if (!avtm.back().forwardPredictedState().isValid() || !predTsos.isValid() ||
0092         !avtm.back().updatedState().isValid()) {
0093       edm::LogError("InvalidState") << "first hit";
0094       return Trajectory();
0095     }
0096
0097     myTraj.push(TM(avtm.back().forwardPredictedState(),
0098                    predTsos,
0099                    avtm.back().updatedState(),
0100                    avtm.back().recHit(),
0101                    avtm.back().estimate(),
0102                    theGeometry->idToLayer(avtm.back().recHit()->geographicalId())),
0103                 avtm.back().estimate());
0104   } else {
0105     currTsos = predTsos;
0106     //check validity
0107     if (!avtm.back().forwardPredictedState().isValid()) {
0108       edm::LogError("InvalidState") << "first hit on invalid hit";
0109       return Trajectory();
0110     }
0111
0112     myTraj.push(TM(avtm.back().forwardPredictedState(),
0113                    avtm.back().recHit(),
0114                    0.,
0115                    theGeometry->idToLayer(avtm.back().recHit()->geographicalId())));
0116   }
0117
0118   TrajectoryStateCombiner combiner;
0119
0120   int hitcounter = avtm.size() - 1;
0121   for (std::vector<TM>::const_reverse_iterator itm = avtm.rbegin() + 1; itm < avtm.rend() - 1; ++itm) {
0122     predTsos = usePropagator->propagate(currTsos, *(*itm).recHit()->surface());
0123     if (predTsos.isValid() && theConvolutor && theMatBeforeUpdate)
0124       predTsos = (*theConvolutor)(predTsos, usePropagator->propagationDirection());
0125     if (!predTsos.isValid()) {
0126       edm::LogInfo("GsfTrackFitters") << "GsfTrajectorySmoother: predicted tsos not valid!";
0127       return Trajectory();
0128     }
0129     if (theMerger)
0130       predTsos = theMerger->merge(predTsos);
0131
0132     if ((*itm).recHit()->isValid()) {
0133       //update
0134       currTsos = updator()->update(predTsos, *(*itm).recHit());
0135
0136       if (currTsos.isValid() && theConvolutor && !theMatBeforeUpdate)
0137         currTsos = (*theConvolutor)(currTsos, usePropagator->propagationDirection());
0138       if (!currTsos.isValid()) {
0139         edm::LogInfo("GsfTrackFitters") << "GsfTrajectorySmoother: tsos not valid after update / material effects!";
0140         return Trajectory();
0141       }
0142       //3 different possibilities to calculate smoothed state:
0143       //1: update combined predictions with hit
0144       //2: combine fwd-prediction with bwd-filter
0145       //3: combine bwd-prediction with fwd-filter
0146       TSOS combTsos = combiner(predTsos, (*itm).forwardPredictedState());
0147       if (!combTsos.isValid()) {
0148         LogDebug("GsfTrackFitters") << "KFTrajectorySmoother: combined tsos not valid!\n"
0149                                     << "pred Tsos pos: " << predTsos.globalPosition() << "\n"
0150                                     << "pred Tsos mom: " << predTsos.globalMomentum() << "\n"
0151                                     << "TrackingRecHit: "
0152                                     << (*itm).recHit()->surface()->toGlobal((*itm).recHit()->localPosition()) << "\n";
0153         return Trajectory();
0154       }
0155
0156       TSOS smooTsos = combiner((*itm).updatedState(), predTsos);
0157
0158       if (!smooTsos.isValid()) {
0159         LogDebug("GsfTrackFitters") << "KFTrajectorySmoother: smoothed tsos not valid!";
0160         return Trajectory();
0161       }
0162
0163       if (!(*itm).forwardPredictedState().isValid() || !predTsos.isValid() || !smooTsos.isValid()) {
0164         edm::LogError("InvalidState") << "inside hits with combination.";
0165         return Trajectory();
0166       }
0167
0168       auto chi2 = estimator()->estimate(combTsos, *(*itm).recHit()).second;
0169       myTraj.push(TM((*itm).forwardPredictedState(),
0170                      predTsos,
0171                      smooTsos,
0172                      (*itm).recHit(),
0173                      chi2,
0174                      theGeometry->idToLayer((*itm).recHit()->geographicalId())),
0175                   (*itm).estimate());
0176       LogDebug("GsfTrackFitters") << "added measurement #" << hitcounter-- << " with chi2 " << chi2;
0177       dump(predTsos, "predTsos", "GsfTrackFitters");
0178       dump(smooTsos, "smooTsos", "GsfTrackFitters");
0179
0180     } else {
0181       currTsos = predTsos;
0182       TSOS combTsos = combiner(predTsos, (*itm).forwardPredictedState());
0183
0184       if (!combTsos.isValid()) {
0185         LogDebug("GsfTrackFitters") << "KFTrajectorySmoother: combined tsos not valid!";
0186         return Trajectory();
0187       }
0188
0189       if (!(*itm).forwardPredictedState().isValid() || !predTsos.isValid() || !combTsos.isValid()) {
0190         edm::LogError("InvalidState") << "inside hits with invalid rechit.";
0191         return Trajectory();
0192       }
0193
0194       myTraj.push(TM((*itm).forwardPredictedState(),
0195                      predTsos,
0196                      combTsos,
0197                      (*itm).recHit(),
0198                      0.,
0199                      theGeometry->idToLayer((*itm).recHit()->geographicalId())));
0200       LogDebug("GsfTrackFitters") << "added invalid measurement #" << hitcounter--;
0201       dump(predTsos, "predTsos", "GsfTrackFitters");
0202       dump(combTsos, "smooTsos", "GsfTrackFitters");
0203     }
0204     if (theMerger)
0205       currTsos = theMerger->merge(currTsos);
0206
0207     dump(currTsos, "currTsos", "GsfTrackFitters");
0208   }
0209
0210   //last smoothed tm is last filtered
0211   predTsos = usePropagator->propagate(currTsos, *avtm.front().recHit()->surface());
0212   if (predTsos.isValid() && theConvolutor && theMatBeforeUpdate)
0213     predTsos = (*theConvolutor)(predTsos, usePropagator->propagationDirection());
0214   if (!predTsos.isValid()) {
0215     edm::LogInfo("GsfTrackFitters") << "GsfTrajectorySmoother: predicted tsos not valid!";
0216     return Trajectory();
0217   }
0218   if (theMerger)
0219     predTsos = theMerger->merge(predTsos);
0220
0221   if (avtm.front().recHit()->isValid()) {
0222     //update
0223     currTsos = updator()->update(predTsos, *avtm.front().recHit());
0224     if (currTsos.isValid() && theConvolutor && !theMatBeforeUpdate)
0225       currTsos = (*theConvolutor)(currTsos, usePropagator->propagationDirection());
0226     if (!currTsos.isValid()) {
0227       edm::LogInfo("GsfTrackFitters") << "GsfTrajectorySmoother: tsos not valid after update / material effects!";
0228       return Trajectory();
0229     }
0230
0231     if (!avtm.front().forwardPredictedState().isValid() || !predTsos.isValid() || !currTsos.isValid()) {
0232       edm::LogError("InvalidState") << "last hit";
0233       return Trajectory();
0234     }
0235
0236     auto chi2 = estimator()->estimate(predTsos, *avtm.front().recHit()).second;
0237     myTraj.push(TM(avtm.front().forwardPredictedState(),
0238                    predTsos,
0239                    currTsos,
0240                    avtm.front().recHit(),
0241                    chi2,
0242                    theGeometry->idToLayer(avtm.front().recHit()->geographicalId())),
0243                 avtm.front().estimate());
0244     LogDebug("GsfTrackFitters") << "added measurement #" << hitcounter-- << " with chi2 " << chi2;
0245     dump(predTsos, "predTsos", "GsfTrackFitters");
0246     dump(currTsos, "smooTsos", "GsfTrackFitters");
0247
0248   } else {
0249     if (!avtm.front().forwardPredictedState().isValid()) {
0250       edm::LogError("InvalidState") << "last invalid hit";
0251       return Trajectory();
0252     }
0253     myTraj.push(TM(avtm.front().forwardPredictedState(),
0254                    avtm.front().recHit(),
0255                    0.,
0256                    theGeometry->idToLayer(avtm.front().recHit()->geographicalId())));
0257     LogDebug("GsfTrackFitters") << "added invalid measurement #" << hitcounter--;
0258   }
0259
0260   return myTraj;
0261 }