CkfPattern/src/TransientInitialStateEstimator.cc

0001 #include "RecoTracker/CkfPattern/interface/TransientInitialStateEstimator.h"
0002
0003 #include "DataFormats/Common/interface/Handle.h"
0004 #include "FWCore/Framework/interface/EventSetup.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006
0007 #include "MagneticField/Engine/interface/MagneticField.h"
0008
0009 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
0010 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0011 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0012
0013 #include "TrackingTools/TrackFitters/interface/KFTrajectoryFitter.h"
0014 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0015
0016 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0017
0018 using namespace std;
0019
0020 TransientInitialStateEstimator::TransientInitialStateEstimator(const edm::ParameterSet& conf, edm::ConsumesCollector iC)
0021     : thePropagatorAlongToken(
0022           iC.esConsumes(edm::ESInputTag("", conf.getParameter<std::string>("propagatorAlongTISE")))),
0023       thePropagatorOppositeToken(
0024           iC.esConsumes(edm::ESInputTag("", conf.getParameter<std::string>("propagatorOppositeTISE")))),
0025       thePropagatorAlong(nullptr),
0026       thePropagatorOpposite(nullptr),
0027       theNumberMeasurementsForFit(conf.getParameter<int>("numberMeasurementsForFit")) {}
0028
0029 void TransientInitialStateEstimator::setEventSetup(const edm::EventSetup& es, const TkClonerImpl& hc) {
0030   theHitCloner = hc;
0031   thePropagatorAlong = &es.getData(thePropagatorAlongToken);
0032   thePropagatorOpposite = &es.getData(thePropagatorOppositeToken);
0033 }
0034
0035 std::pair<TrajectoryStateOnSurface, const GeomDet*> TransientInitialStateEstimator::innerState(const Trajectory& traj,
0036                                                                                                bool doBackFit) const {
0037   if (!doBackFit && traj.firstMeasurement().forwardPredictedState().isValid()) {
0038     LogDebug("TransientInitialStateEstimator")
0039         << "a backward fit will not be done. assuming that the state on first measurement is OK";
0040     TSOS firstStateFromForward = traj.firstMeasurement().forwardPredictedState();
0041     firstStateFromForward.rescaleError(100.);
0042     return std::pair<TrajectoryStateOnSurface, const GeomDet*>(std::move(firstStateFromForward),
0043                                                                traj.firstMeasurement().recHit()->det());
0044   }
0045   if (!doBackFit) {
0046     LogDebug("TransientInitialStateEstimator")
0047         << "told to not do a back fit, but the forward state of the first measurement is not valid. doing a back fit.";
0048   }
0049
0050   int nMeas = traj.measurements().size();
0051   int lastFitted = theNumberMeasurementsForFit >= 0 ? theNumberMeasurementsForFit : nMeas;
0052   if (nMeas - 1 < lastFitted)
0053     lastFitted = nMeas - 1;
0054
0055   auto const& measvec = traj.measurements();
0056   TransientTrackingRecHit::ConstRecHitContainer firstHits;
0057
0058   bool foundLast = false;
0059   int actualLast = -99;
0060
0061   for (int i = lastFitted; i >= 0; i--) {
0062     if (measvec[i].recHit()->det()) {
0063       if (!foundLast) {
0064         actualLast = i;
0065         foundLast = true;
0066       }
0067       firstHits.push_back(measvec[i].recHit());
0068     }
0069   }
0070   TSOS startingState = measvec[actualLast].updatedState();
0071   startingState.rescaleError(100.);
0072
0073   // avoid cloning...
0074   KFUpdator const aKFUpdator;
0075   Chi2MeasurementEstimator const aChi2MeasurementEstimator(100., 3);
0076   KFTrajectoryFitter backFitter(
0077       thePropagatorAlong, &aKFUpdator, &aChi2MeasurementEstimator, firstHits.size(), nullptr, &theHitCloner);
0078
0079   PropagationDirection backFitDirection = traj.direction() == alongMomentum ? oppositeToMomentum : alongMomentum;
0080
0081   // only direction matters in this contest
0082   TrajectorySeed fakeSeed(PTrajectoryStateOnDet(), edm::OwnVector<TrackingRecHit>(), backFitDirection);
0083
0084   Trajectory&& fitres = backFitter.fitOne(
0085       fakeSeed, firstHits, startingState, traj.nLoops() > 0 ? TrajectoryFitter::looper : TrajectoryFitter::standard);
0086
0087   LogDebug("TransientInitialStateEstimator")
0088       << "using a backward fit of :" << firstHits.size() << " hits, starting from:\n"
0089       << startingState << " to get the estimate of the initial state of the track.";
0090
0091   if (!fitres.isValid()) {
0092     LogDebug("TransientInitialStateEstimator") << "FitTester: first hits fit failed!";
0093     return std::pair<TrajectoryStateOnSurface, const GeomDet*>();
0094   }
0095
0096   TrajectoryMeasurement const& firstMeas = fitres.lastMeasurement();
0097
0098   // magnetic field can be different!
0099   TSOS firstState(firstMeas.updatedState().localParameters(),
0100                   firstMeas.updatedState().localError(),
0101                   firstMeas.updatedState().surface(),
0102                   thePropagatorAlong->magneticField());
0103
0104   // TSOS & firstState = const_cast<TSOS&>(firstMeas.updatedState());
0105
0106   // this fails in case of zero field? (for sure for beamhalo reconstruction)
0107   // assert(thePropagatorAlong->magneticField()==firstState.magneticField());
0108
0109   firstState.rescaleError(100.);
0110
0111   LogDebug("TransientInitialStateEstimator")
0112       << "the initial state is found to be:\n:" << firstState
0113       << "\n it's field pointer is: " << firstState.magneticField()
0114       << "\n the pointer from the state of the back fit was: " << firstMeas.updatedState().magneticField();
0115
0116   return std::pair<TrajectoryStateOnSurface, const GeomDet*>(std::move(firstState), firstMeas.recHit()->det());
0117 }