File indexing completed on 2021-02-14 14:31:38
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
0079 if (avtm.back().recHit()->isValid()) {
0080 {
0081
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
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
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
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
0143
0144
0145
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
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
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 }