00001 #include "RecoTracker/CkfPattern/interface/BaseCkfTrajectoryBuilder.h"
00002
00003 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
00004 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
00005 #include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h"
00006 #include "TrackingTools/TrajectoryFiltering/interface/TrajectoryFilter.h"
00007
00008
00009 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
00010 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00011 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
00012 #include "TrackingTools/PatternTools/interface/TrajectoryStateUpdator.h"
00013 #include "TrackingTools/TrajectoryState/interface/BasicSingleTrajectoryState.h"
00014 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimatorBase.h"
00015
00016 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
00017 #include "TrackingTools/PatternTools/interface/TempTrajectory.h"
00018 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00019
00020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00021
00022
00023 BaseCkfTrajectoryBuilder::
00024 BaseCkfTrajectoryBuilder(const edm::ParameterSet& conf,
00025 const TrajectoryStateUpdator* updator,
00026 const Propagator* propagatorAlong,
00027 const Propagator* propagatorOpposite,
00028 const Chi2MeasurementEstimatorBase* estimator,
00029 const TransientTrackingRecHitBuilder* recHitBuilder,
00030 const MeasurementTracker* measurementTracker,
00031 const TrajectoryFilter* filter,
00032 const TrajectoryFilter* inOutFilter):
00033 theUpdator(updator),
00034 thePropagatorAlong(propagatorAlong),thePropagatorOpposite(propagatorOpposite),
00035 theEstimator(estimator),theTTRHBuilder(recHitBuilder),
00036 theMeasurementTracker(measurementTracker),
00037 theLayerMeasurements(new LayerMeasurements(theMeasurementTracker)),
00038 theForwardPropagator(0),theBackwardPropagator(0),
00039 theFilter(filter),
00040 theInOutFilter(inOutFilter)
00041 {}
00042
00043 BaseCkfTrajectoryBuilder::~BaseCkfTrajectoryBuilder(){
00044 delete theLayerMeasurements;
00045 }
00046
00047
00048 void
00049 BaseCkfTrajectoryBuilder::seedMeasurements(const TrajectorySeed& seed, std::vector<TrajectoryMeasurement> & result) const
00050 {
00051 TrajectoryStateTransform tsTransform;
00052
00053 TrajectorySeed::range hitRange = seed.recHits();
00054 for (TrajectorySeed::const_iterator ihit = hitRange.first;
00055 ihit != hitRange.second; ihit++) {
00056 TransientTrackingRecHit::RecHitPointer recHit = theTTRHBuilder->build(&(*ihit));
00057 const GeomDet* hitGeomDet =
00058 theMeasurementTracker->geomTracker()->idToDet( ihit->geographicalId());
00059
00060 const DetLayer* hitLayer =
00061 theMeasurementTracker->geometricSearchTracker()->detLayer(ihit->geographicalId());
00062
00063 TSOS invalidState( new BasicSingleTrajectoryState( hitGeomDet->surface()));
00064 if (ihit == hitRange.second - 1) {
00065
00066 PTrajectoryStateOnDet pState( seed.startingState());
00067 const GeomDet* gdet = theMeasurementTracker->geomTracker()->idToDet( DetId(pState.detId()));
00068 if (&gdet->surface() != &hitGeomDet->surface()) {
00069 edm::LogError("CkfPattern") << "CkfTrajectoryBuilder error: the seed state is not on the surface of the detector of the last seed hit";
00070 return;
00071 }
00072
00073 TSOS updatedState = tsTransform.transientState( pState, &(gdet->surface()),
00074 theForwardPropagator->magneticField());
00075 result.push_back(TM( invalidState, updatedState, recHit, 0, hitLayer));
00076 }
00077 else {
00078 PTrajectoryStateOnDet pState( seed.startingState());
00079
00080 TSOS outerState = tsTransform.transientState(pState,
00081 &((theMeasurementTracker->geomTracker()->idToDet(
00082 (hitRange.second - 1)->geographicalId()))->surface()),
00083 theForwardPropagator->magneticField());
00084 TSOS innerState = theBackwardPropagator->propagate(outerState,hitGeomDet->surface());
00085 if(innerState.isValid()) {
00086 TSOS innerUpdated = theUpdator->update(innerState,*recHit);
00087 result.push_back(TM( invalidState, innerUpdated, recHit, 0, hitLayer));
00088 }
00089 }
00090 }
00091
00092
00093 fillSeedHistoDebugger(result.begin(),result.end());
00094
00095 }
00096
00097
00098 TempTrajectory BaseCkfTrajectoryBuilder::
00099 createStartingTrajectory( const TrajectorySeed& seed) const
00100 {
00101 TempTrajectory result( seed, seed.direction());
00102 if ( seed.direction() == alongMomentum) {
00103 theForwardPropagator = &(*thePropagatorAlong);
00104 theBackwardPropagator = &(*thePropagatorOpposite);
00105 }
00106 else {
00107 theForwardPropagator = &(*thePropagatorOpposite);
00108 theBackwardPropagator = &(*thePropagatorAlong);
00109 }
00110
00111 std::vector<TM> seedMeas;
00112 seedMeasurements(seed, seedMeas);
00113 for (std::vector<TM>::const_iterator i=seedMeas.begin(); i!=seedMeas.end(); i++)
00114 result.push(*i);
00115
00116 LogDebug("CkfPattern")
00117 <<" initial trajectory from the seed: "<<PrintoutHelper::dumpCandidate(result,true);
00118
00119 return result;
00120 }
00121
00122
00123 bool BaseCkfTrajectoryBuilder::toBeContinued (TempTrajectory& traj, bool inOut) const
00124 {
00125 if (traj.measurements().size() > 400) {
00126 edm::LogError("BaseCkfTrajectoryBuilder_InfiniteLoop");
00127 LogTrace("BaseCkfTrajectoryBuilder_InfiniteLoop") <<
00128 "Cropping Track After 400 Measurements:\n" <<
00129 " Last predicted state: " << traj.lastMeasurement().predictedState() << "\n" <<
00130 " Last layer subdetector: " << (traj.lastLayer() ? traj.lastLayer()->subDetector() : -1) << "\n" <<
00131 " Found hits: " << traj.foundHits() << ", lost hits: " << traj.lostHits() << "\n\n";
00132 return false;
00133 }
00134
00135
00136 if (inOut) {
00137 if (theInOutFilter == 0) edm::LogError("CkfPattern") << "CkfTrajectoryBuilder error: trying to use dedicated filter for in-out tracking phase, when none specified";
00138 return theInOutFilter->toBeContinued(traj);
00139 } else {
00140 return theFilter->toBeContinued(traj);
00141 }
00142 }
00143
00144
00145 bool BaseCkfTrajectoryBuilder::qualityFilter( const TempTrajectory& traj, bool inOut) const
00146 {
00147
00148
00149 if (inOut) {
00150 if (theInOutFilter == 0) edm::LogError("CkfPattern") << "CkfTrajectoryBuilder error: trying to use dedicated filter for in-out tracking phase, when none specified";
00151 return theInOutFilter->qualityFilter(traj);
00152 } else {
00153 return theFilter->qualityFilter(traj);
00154 }
00155 }
00156
00157
00158 void
00159 BaseCkfTrajectoryBuilder::addToResult (TempTrajectory& tmptraj,
00160 TrajectoryContainer& result,
00161 bool inOut) const
00162 {
00163
00164 if ( !qualityFilter(tmptraj, inOut) ) return;
00165 Trajectory traj = tmptraj.toTrajectory();
00166
00167 while (!traj.empty() && !traj.lastMeasurement().recHit()->isValid()) traj.pop();
00168 LogDebug("CkfPattern")<<inOut<<"=inOut option. pushing a Trajectory with: "<<traj.foundHits()<<" found hits. "<<traj.lostHits()
00169 <<" lost hits. Popped :"<<(tmptraj.measurements().size())-(traj.measurements().size())<<" hits.";
00170 result.push_back( traj);
00171 }
00172 void
00173 BaseCkfTrajectoryBuilder::addToResult (TempTrajectory& tmptraj,
00174 TempTrajectoryContainer& result,
00175 bool inOut) const
00176 {
00177
00178 if ( !qualityFilter(tmptraj, inOut) ) return;
00179
00180 TempTrajectory traj = tmptraj;
00181 while (!traj.empty() && !traj.lastMeasurement().recHit()->isValid()) traj.pop();
00182 LogDebug("CkfPattern")<<inOut<<"=inOut option. pushing a TempTrajectory with: "<<traj.foundHits()<<" found hits. "<<traj.lostHits()
00183 <<" lost hits. Popped :"<<(tmptraj.measurements().size())-(traj.measurements().size())<<" hits.";
00184 result.push_back( traj );
00185 }
00186
00187
00188
00189 BaseCkfTrajectoryBuilder::StateAndLayers
00190 BaseCkfTrajectoryBuilder::findStateAndLayers(const TempTrajectory& traj) const
00191 {
00192 if (traj.empty())
00193 {
00194
00195 PTrajectoryStateOnDet ptod = traj.seed().startingState();
00196 DetId id(ptod.detId());
00197 const GeomDet * g = theMeasurementTracker->geomTracker()->idToDet(id);
00198 const Surface * surface=&g->surface();
00199 TrajectoryStateTransform tsTransform;
00200
00201 TSOS currentState = TrajectoryStateOnSurface(tsTransform.transientState(ptod,surface,theForwardPropagator->magneticField()));
00202 const DetLayer* lastLayer = theMeasurementTracker->geometricSearchTracker()->detLayer(id);
00203 return StateAndLayers(currentState,lastLayer->nextLayers( *currentState.freeState(), traj.direction()) );
00204 }
00205 else
00206 {
00207 TSOS currentState = traj.lastMeasurement().updatedState();
00208 return StateAndLayers(currentState,traj.lastLayer()->nextLayers( *currentState.freeState(), traj.direction()) );
00209 }
00210 }