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 (!predTsos.isValid()) {
0133       return Trajectory();
0134     }
0135
0136     if ((*itm).recHit()->isValid()) {
0137       //update
0138       currTsos = updator()->update(predTsos, *(*itm).recHit());
0139
0140       if (currTsos.isValid() && theConvolutor && !theMatBeforeUpdate)
0141         currTsos = (*theConvolutor)(currTsos, usePropagator->propagationDirection());
0142       if (!currTsos.isValid()) {
0143         edm::LogInfo("GsfTrackFitters") << "GsfTrajectorySmoother: tsos not valid after update / material effects!";
0144         return Trajectory();
0145       }
0146       //3 different possibilities to calculate smoothed state:
0147       //1: update combined predictions with hit
0148       //2: combine fwd-prediction with bwd-filter
0149       //3: combine bwd-prediction with fwd-filter
0150       TSOS combTsos = combiner(predTsos, (*itm).forwardPredictedState());
0151       if (!combTsos.isValid()) {
0152         LogDebug("GsfTrackFitters") << "KFTrajectorySmoother: combined tsos not valid!\n"
0153                                     << "pred Tsos pos: " << predTsos.globalPosition() << "\n"
0154                                     << "pred Tsos mom: " << predTsos.globalMomentum() << "\n"
0155                                     << "TrackingRecHit: "
0156                                     << (*itm).recHit()->surface()->toGlobal((*itm).recHit()->localPosition()) << "\n";
0157         return Trajectory();
0158       }
0159
0160       TSOS smooTsos = combiner((*itm).updatedState(), predTsos);
0161
0162       if (!smooTsos.isValid()) {
0163         LogDebug("GsfTrackFitters") << "KFTrajectorySmoother: smoothed tsos not valid!";
0164         return Trajectory();
0165       }
0166
0167       if (!(*itm).forwardPredictedState().isValid() || !predTsos.isValid() || !smooTsos.isValid()) {
0168         edm::LogError("InvalidState") << "inside hits with combination.";
0169         return Trajectory();
0170       }
0171
0172       auto chi2 = estimator()->estimate(combTsos, *(*itm).recHit()).second;
0173       myTraj.push(TM((*itm).forwardPredictedState(),
0174                      predTsos,
0175                      smooTsos,
0176                      (*itm).recHit(),
0177                      chi2,
0178                      theGeometry->idToLayer((*itm).recHit()->geographicalId())),
0179                   (*itm).estimate());
0180       LogDebug("GsfTrackFitters") << "added measurement #" << hitcounter-- << " with chi2 " << chi2;
0181       dump(predTsos, "predTsos", "GsfTrackFitters");
0182       dump(smooTsos, "smooTsos", "GsfTrackFitters");
0183
0184     } else {
0185       currTsos = predTsos;
0186       TSOS combTsos = combiner(predTsos, (*itm).forwardPredictedState());
0187
0188       if (!combTsos.isValid()) {
0189         LogDebug("GsfTrackFitters") << "KFTrajectorySmoother: combined tsos not valid!";
0190         return Trajectory();
0191       }
0192
0193       if (!(*itm).forwardPredictedState().isValid() || !predTsos.isValid() || !combTsos.isValid()) {
0194         edm::LogError("InvalidState") << "inside hits with invalid rechit.";
0195         return Trajectory();
0196       }
0197
0198       myTraj.push(TM((*itm).forwardPredictedState(),
0199                      predTsos,
0200                      combTsos,
0201                      (*itm).recHit(),
0202                      0.,
0203                      theGeometry->idToLayer((*itm).recHit()->geographicalId())));
0204       LogDebug("GsfTrackFitters") << "added invalid measurement #" << hitcounter--;
0205       dump(predTsos, "predTsos", "GsfTrackFitters");
0206       dump(combTsos, "smooTsos", "GsfTrackFitters");
0207     }
0208     if (theMerger)
0209       currTsos = theMerger->merge(currTsos);
0210
0211     if (!currTsos.isValid()) {
0212       return Trajectory();
0213     }
0214     dump(currTsos, "currTsos", "GsfTrackFitters");
0215   }
0216
0217   //last smoothed tm is last filtered
0218   predTsos = usePropagator->propagate(currTsos, *avtm.front().recHit()->surface());
0219   if (predTsos.isValid() && theConvolutor && theMatBeforeUpdate)
0220     predTsos = (*theConvolutor)(predTsos, usePropagator->propagationDirection());
0221   if (!predTsos.isValid()) {
0222     edm::LogInfo("GsfTrackFitters") << "GsfTrajectorySmoother: predicted tsos not valid!";
0223     return Trajectory();
0224   }
0225   if (theMerger)
0226     predTsos = theMerger->merge(predTsos);
0227
0228   if (avtm.front().recHit()->isValid()) {
0229     //update
0230     currTsos = updator()->update(predTsos, *avtm.front().recHit());
0231     if (currTsos.isValid() && theConvolutor && !theMatBeforeUpdate)
0232       currTsos = (*theConvolutor)(currTsos, usePropagator->propagationDirection());
0233     if (!currTsos.isValid()) {
0234       edm::LogInfo("GsfTrackFitters") << "GsfTrajectorySmoother: tsos not valid after update / material effects!";
0235       return Trajectory();
0236     }
0237
0238     if (!avtm.front().forwardPredictedState().isValid() || !predTsos.isValid() || !currTsos.isValid()) {
0239       edm::LogError("InvalidState") << "last hit";
0240       return Trajectory();
0241     }
0242
0243     auto chi2 = estimator()->estimate(predTsos, *avtm.front().recHit()).second;
0244     myTraj.push(TM(avtm.front().forwardPredictedState(),
0245                    predTsos,
0246                    currTsos,
0247                    avtm.front().recHit(),
0248                    chi2,
0249                    theGeometry->idToLayer(avtm.front().recHit()->geographicalId())),
0250                 avtm.front().estimate());
0251     LogDebug("GsfTrackFitters") << "added measurement #" << hitcounter-- << " with chi2 " << chi2;
0252     dump(predTsos, "predTsos", "GsfTrackFitters");
0253     dump(currTsos, "smooTsos", "GsfTrackFitters");
0254
0255   } else {
0256     if (!avtm.front().forwardPredictedState().isValid()) {
0257       edm::LogError("InvalidState") << "last invalid hit";
0258       return Trajectory();
0259     }
0260     myTraj.push(TM(avtm.front().forwardPredictedState(),
0261                    avtm.front().recHit(),
0262                    0.,
0263                    theGeometry->idToLayer(avtm.front().recHit()->geographicalId())));
0264     LogDebug("GsfTrackFitters") << "added invalid measurement #" << hitcounter--;
0265   }
0266
0267   return myTraj;
0268 }