File indexing completed on 2023-03-17 11:22:07
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
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
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
0099 TSOS firstState(firstMeas.updatedState().localParameters(),
0100 firstMeas.updatedState().localError(),
0101 firstMeas.updatedState().surface(),
0102 thePropagatorAlong->magneticField());
0103
0104
0105
0106
0107
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 }