CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CkfTrackCandidateMakerBase.cc
Go to the documentation of this file.
1 #include <memory>
2 #include <string>
3 
8 
12 
14 
19 
24 
25 
30 
33 
34 #include<algorithm>
35 #include<functional>
36 
38 
39 using namespace edm;
40 using namespace std;
41 
42 namespace cms{
43  CkfTrackCandidateMakerBase::CkfTrackCandidateMakerBase(edm::ParameterSet const& conf) :
44 
45  conf_(conf),
46  theTrackCandidateOutput(true),
47  theTrajectoryOutput(false),
48  useSplitting(conf.getParameter<bool>("useHitsSplitting")),
49  doSeedingRegionRebuilding(conf.getParameter<bool>("doSeedingRegionRebuilding")),
50  cleanTrajectoryAfterInOut(conf.getParameter<bool>("cleanTrajectoryAfterInOut")),
51  reverseTrajectories(conf.existsAs<bool>("reverseTrajectories") && conf.getParameter<bool>("reverseTrajectories")),
52  theMaxNSeeds(conf.getParameter<unsigned int>("maxNSeeds")),
53  theTrajectoryBuilderName(conf.getParameter<std::string>("TrajectoryBuilder")),
54  theTrajectoryBuilder(0),
55  theTrajectoryCleanerName(conf.getParameter<std::string>("TrajectoryCleaner")),
56  theTrajectoryCleaner(0),
57  theInitialState(0),
58  theNavigationSchoolName(conf.getParameter<std::string>("NavigationSchool")),
59  theNavigationSchool(0),
60  theSeedCleaner(0),
61  maxSeedsBeforeCleaning_(0)
62  {
63  //produces<TrackCandidateCollection>();
64  // old configuration totally descoped.
65  // if (!conf.exists("src"))
66  // theSeedLabel = InputTag(conf_.getParameter<std::string>("SeedProducer"),conf_.getParameter<std::string>("SeedLabel"));
67  // else
69  if ( conf.exists("maxSeedsBeforeCleaning") )
70  maxSeedsBeforeCleaning_=conf.getParameter<unsigned int>("maxSeedsBeforeCleaning");
71 
72  }
73 
74 
75  // Virtual destructor needed.
77  delete theInitialState;
78  if (theSeedCleaner) delete theSeedCleaner;
79  }
80 
82  {
83  std::string cleaner = conf_.getParameter<std::string>("RedundantSeedCleaner");
84  if (cleaner == "SeedCleanerByHitPosition") {
86  } else if (cleaner == "SeedCleanerBySharedInput") {
88  } else if (cleaner == "CachingSeedCleanerByHitPosition") {
90  } else if (cleaner == "CachingSeedCleanerBySharedInput") {
91  int numHitsForSeedCleaner = conf_.existsAs<int>("numHitsForSeedCleaner") ?
92  conf_.getParameter<int>("numHitsForSeedCleaner") : 4;
93  int onlyPixelHits = conf_.existsAs<bool>("onlyPixelHitsForSeedCleaner") ?
94  conf_.getParameter<bool>("onlyPixelHitsForSeedCleaner") : false;
95  theSeedCleaner = new CachingSeedCleanerBySharedInput(numHitsForSeedCleaner,onlyPixelHits);
96  } else if (cleaner == "none") {
97  theSeedCleaner = 0;
98  } else {
99  throw cms::Exception("RedundantSeedCleaner not found", cleaner);
100  }
101  }
102 
104 
105  //services
108 
109  if (!theInitialState){
110  // constructor uses the EventSetup, it must be in the setEventSetup were it has a proper value.
111  // get nested parameter set for the TransientInitialStateEstimator
112  ParameterSet tise_params = conf_.getParameter<ParameterSet>("TransientInitialStateEstimatorParameters") ;
113  theInitialState = new TransientInitialStateEstimator( es,tise_params);
114  }
115 
117 
118  edm::ESHandle<TrajectoryCleaner> trajectoryCleanerH;
119  es.get<TrajectoryCleaner::Record>().get(theTrajectoryCleanerName, trajectoryCleanerH);
120  theTrajectoryCleaner= trajectoryCleanerH.product();
121 
122  edm::ESHandle<NavigationSchool> navigationSchoolH;
123  es.get<NavigationSchoolRecord>().get(theNavigationSchoolName, navigationSchoolH);
124  theNavigationSchool = navigationSchoolH.product();
125 
126  // set the TrajectoryBuilder
127  edm::ESHandle<TrajectoryBuilder> theTrajectoryBuilderHandle;
128  es.get<CkfComponentsRecord>().get(theTrajectoryBuilderName,theTrajectoryBuilderHandle);
129  theTrajectoryBuilder = theTrajectoryBuilderHandle.product();
130 
131  }
132 
133  // Functions that gets called by framework every event
135  {
136  // getting objects from the EventSetup
137  setEventSetup( es );
138 
139  // set the correct navigation
141 
142  // propagator
143  edm::ESHandle<Propagator> thePropagator;
144  es.get<TrackingComponentsRecord>().get("AnyDirectionAnalyticalPropagator",
145  thePropagator);
146 
147  // method for Debugging
149 
150  // Step A: set Event for the TrajectoryBuilder
152 
153  // Step B: Retrieve seeds
154 
156  e.getByLabel(theSeedLabel, collseed);
157 
158  // Step C: Create empty output collection
159  std::auto_ptr<TrackCandidateCollection> output(new TrackCandidateCollection);
160  std::auto_ptr<std::vector<Trajectory> > outputT (new std::vector<Trajectory>());
161 
162  if ( (*collseed).size()>theMaxNSeeds ) {
163  LogError("TooManySeeds")<<"Exceeded maximum numeber of seeds! theMaxNSeeds="<<theMaxNSeeds<<" nSeed="<<(*collseed).size();
164  if (theTrackCandidateOutput){ e.put(output);}
165  if (theTrajectoryOutput){e.put(outputT);}
167  return;
168  }
169 
170  // Step D: Invoke the building algorithm
171  if ((*collseed).size()>0){
172 
173  unsigned int lastCleanResult=0;
174  vector<Trajectory> rawResult;
175  rawResult.reserve(collseed->size() * 4);
176 
177  if (theSeedCleaner) theSeedCleaner->init( &rawResult );
178 
179  // method for debugging
181 
182  vector<Trajectory> theTmpTrajectories;
183 
184  // Loop over seeds
185  size_t collseed_size = collseed->size();
186  for (size_t j = 0; j < collseed_size; j++){
187 
188  // Check if seed hits already used by another track
189  if (theSeedCleaner && !theSeedCleaner->good( &((*collseed)[j])) ) {
190  LogDebug("CkfTrackCandidateMakerBase")<<" Seed cleaning kills seed "<<j;
191  continue;
192  }
193 
194  // Build trajectory from seed outwards
195  theTmpTrajectories.clear();
196  theTrajectoryBuilder->trajectories( (*collseed)[j], theTmpTrajectories );
197 
198 
199  LogDebug("CkfPattern") << "======== In-out trajectory building found " << theTmpTrajectories.size()
200  << " trajectories from seed " << j << " ========"<<endl
201  <<PrintoutHelper::dumpCandidates(theTmpTrajectories);
202 
204 
205  // Select the best trajectory from this seed (declare others invalid)
206  theTrajectoryCleaner->clean(theTmpTrajectories);
207 
208  LogDebug("CkfPattern") << "======== In-out trajectory cleaning gave the following valid trajectories from seed "
209  << j << " ========"<<endl
210  << PrintoutHelper::dumpCandidates(theTmpTrajectories);
211  }
212 
213  // Optionally continue building trajectory back through
214  // seed and if possible further inwards.
215 
217  theTrajectoryBuilder->rebuildSeedingRegion((*collseed)[j],theTmpTrajectories);
218 
219  LogDebug("CkfPattern") << "======== Out-in trajectory building found " << theTmpTrajectories.size()
220  << " valid/invalid trajectories from seed " << j << " ========"<<endl
221  <<PrintoutHelper::dumpCandidates(theTmpTrajectories);
222  }
223 
224 
225  // Select the best trajectory from this seed (after seed region rebuilding, can be more than one)
226  theTrajectoryCleaner->clean(theTmpTrajectories);
227 
228  LogDebug("CkfPattern") << "======== Trajectory cleaning gave the following valid trajectories from seed "
229  << j << " ========"<<endl
230  <<PrintoutHelper::dumpCandidates(theTmpTrajectories);
231 
232  for(vector<Trajectory>::iterator it=theTmpTrajectories.begin();
233  it!=theTmpTrajectories.end(); it++){
234  if( it->isValid() ) {
235  it->setSeedRef(collseed->refAt(j));
236  // Store trajectory
237  rawResult.push_back(*it);
238  // Tell seed cleaner which hits this trajectory used.
239  //TO BE FIXED: this cut should be configurable via cfi file
240  if (theSeedCleaner && it->foundHits()>3) theSeedCleaner->add( & (*it) );
241  //if (theSeedCleaner ) theSeedCleaner->add( & (*it) );
242  }
243  }
244 
245  theTmpTrajectories.clear();
246 
247  LogDebug("CkfPattern") << "rawResult trajectories found so far = " << rawResult.size();
248 
249  if ( maxSeedsBeforeCleaning_ >0 && rawResult.size() > maxSeedsBeforeCleaning_+lastCleanResult) {
250  theTrajectoryCleaner->clean(rawResult);
251  rawResult.erase(std::remove_if(rawResult.begin(),rawResult.end(),
252  std::not1(std::mem_fun_ref(&Trajectory::isValid))),
253  rawResult.end());
254  lastCleanResult=rawResult.size();
255  }
256 
257  }
258  // end of loop over seeds
259 
261 
262  // Step E: Clean the results to avoid duplicate tracks
263  // Rejected ones just flagged as invalid.
264  theTrajectoryCleaner->clean(rawResult);
265 
266  LogDebug("CkfPattern") << "======== Final cleaning of entire event found " << rawResult.size()
267  << " valid/invalid trajectories ======="<<endl
268  <<PrintoutHelper::dumpCandidates(rawResult);
269 
270  LogDebug("CkfPattern") << "removing invalid trajectories.";
271 
272  vector<Trajectory> & unsmoothedResult(rawResult);
273  unsmoothedResult.erase(std::remove_if(unsmoothedResult.begin(),unsmoothedResult.end(),
274  std::not1(std::mem_fun_ref(&Trajectory::isValid))),
275  unsmoothedResult.end());
276 
277  // If requested, reverse the trajectories creating a new 1-hit seed on the last measurement of the track
278  if (reverseTrajectories) {
279  vector<Trajectory> reversed;
280  reversed.reserve(unsmoothedResult.size());
281  for (vector<Trajectory>::const_iterator it = unsmoothedResult.begin(), ed = unsmoothedResult.end(); it != ed; ++it) {
282  // reverse the trajectory only if it has valid hit on the last measurement (should happen)
283  if (it->lastMeasurement().updatedState().isValid() &&
284  it->lastMeasurement().recHit().get() != 0 &&
285  it->lastMeasurement().recHit()->isValid()) {
286  // I can't use reverse in place, because I want to change the seed
287  // 1) reverse propagation direction
288  PropagationDirection direction = it->direction();
289  if (direction == alongMomentum) direction = oppositeToMomentum;
290  else if (direction == oppositeToMomentum) direction = alongMomentum;
291  // 2) make a seed
292  TrajectoryStateOnSurface initState = it->lastMeasurement().updatedState();
293  DetId initDetId = it->lastMeasurement().recHit()->geographicalId();
296  hits.push_back(*it->lastMeasurement().recHit()->hit());
297  boost::shared_ptr<const TrajectorySeed> seed(new TrajectorySeed(state, hits, direction));
298  // 3) make a trajectory
299  Trajectory trajectory(seed, direction);
300  trajectory.setNLoops(it->nLoops());
301  trajectory.setSeedRef(it->seedRef());
302  // 4) push states in reversed order
303  const Trajectory::DataContainer &meas = it->measurements();
304  for (Trajectory::DataContainer::const_reverse_iterator itmeas = meas.rbegin(), endmeas = meas.rend(); itmeas != endmeas; ++itmeas) {
305  trajectory.push(*itmeas);
306  }
307  reversed.push_back(trajectory);
308  } else {
309  edm::LogWarning("CkfPattern_InvalidLastMeasurement") << "Last measurement of the trajectory is invalid, cannot reverse it";
310  reversed.push_back(*it);
311  }
312  }
313  unsmoothedResult.swap(reversed);
314  }
315 
316  // for (vector<Trajectory>::const_iterator itraw = rawResult.begin();
317  // itraw != rawResult.end(); itraw++) {
318  //if((*itraw).isValid()) unsmoothedResult.push_back( *itraw);
319  //}
320 
321  //analyseCleanedTrajectories(unsmoothedResult);
322 
324  // Step F: Convert to TrackCandidates
325  output->reserve(unsmoothedResult.size());
326  for (vector<Trajectory>::const_iterator it = unsmoothedResult.begin();
327  it != unsmoothedResult.end(); ++it) {
328 
330  //it->recHitsV(thits);
331  LogDebug("CkfPattern") << "retrieving "<<(useSplitting?"splitted":"un-splitted")<<" hits from trajectory";
332  it->recHitsV(thits,useSplitting);
334  recHits.reserve(thits.size());
335  LogDebug("CkfPattern") << "cloning hits into new collection.";
336  for (Trajectory::RecHitContainer::const_iterator hitIt = thits.begin();
337  hitIt != thits.end(); ++hitIt) {
338  recHits.push_back( (**hitIt).hit()->clone());
339  }
340 
341  LogDebug("CkfPattern") << "getting initial state.";
342  const bool doBackFit = !doSeedingRegionRebuilding && !reverseTrajectories;
343  std::pair<TrajectoryStateOnSurface, const GeomDet*> initState =
344  theInitialState->innerState( *it , doBackFit);
345 
346  // temporary protection againt invalid initial states
347  if (! initState.first.isValid() || initState.second == 0 || std::isnan(initState.first.globalPosition().x())) {
348  //cout << "invalid innerState, will not make TrackCandidate" << endl;
349  continue;
350  }
351 
353  if(useSplitting && (initState.second != thits.front()->det()) && thits.front()->det() ){
354  LogDebug("CkfPattern") << "propagating to hit front in case of splitting.";
355  TrajectoryStateOnSurface propagated = thePropagator->propagate(initState.first,thits.front()->det()->surface());
356  if (!propagated.isValid()) continue;
358  thits.front()->det()->geographicalId().rawId());
359  }
360  else state = trajectoryStateTransform::persistentState( initState.first,
361  initState.second->geographicalId().rawId());
362  LogDebug("CkfPattern") << "pushing a TrackCandidate.";
363  output->push_back(TrackCandidate(recHits,it->seed(),state,it->seedRef(),it->nLoops() ) );
364  }
365  }//output trackcandidates
366 
368  es.get<TrackerDigiGeometryRecord>().get(tracker);
369  LogTrace("CkfPattern|TrackingRegressionTest") << "========== CkfTrackCandidateMaker Info =========="
370  << "number of Seed: " << collseed->size()<<endl
371  <<PrintoutHelper::regressionTest(*tracker,unsmoothedResult);
372 
373 
374 
375  if (theTrajectoryOutput){ outputT->swap(unsmoothedResult);}
376 
377  }// end of ((*collseed).size()>0)
378 
379  // method for debugging
381 
382  // Step G: write output to file
383  if (theTrackCandidateOutput){ e.put(output);}
384  if (theTrajectoryOutput){e.put(outputT);}
385 
386  //reset the MT.
388  }
389 
390 }
391 
#define LogDebug(id)
T getParameter(std::string const &) const
static std::string dumpCandidates(collection &candidates)
void setEventSetup(const edm::EventSetup &es)
Call this at each event until this object will come from the EventSetup as it should.
bool existsAs(std::string const &parameterName, bool trackiness=true) const
checks if a parameter exists as a given type
Definition: ParameterSet.h:187
virtual void rebuildSeedingRegion(const TrajectorySeed &, TrajectoryContainer &result) const
const TrajectoryBuilder * theTrajectoryBuilder
std::vector< TrackCandidate > TrackCandidateCollection
void setNLoops(signed char value)
Definition: Trajectory.h:321
PTrajectoryStateOnDet persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid)
bool exists(std::string const &parameterName) const
checks if a parameter exists
virtual void printHitsDebugger(edm::Event &e)
edm::ESHandle< GeometricSearchTracker > theGeomSearchTracker
PropagationDirection
virtual void done()=0
Tells the cleaner that the seeds are finished, and so it can clear any cache it has.
virtual void add(const Trajectory *traj)=0
Informs the cleaner that a new trajectory has been made, in case the cleaner keeps a local collection...
virtual bool good(const TrajectorySeed *seed)=0
Returns true if the seed is not overlapping with another trajectory.
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
static std::string regressionTest(const TrackerGeometry &tracker, std::vector< Trajectory > &unsmoothedResult)
void push_back(D *&d)
Definition: OwnVector.h:273
std::vector< TrajectoryMeasurement > DataContainer
Definition: Trajectory.h:42
const TrajectoryCleaner * theTrajectoryCleaner
const NavigationSchool * theNavigationSchool
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:85
bool isnan(float x)
Definition: math.h:13
virtual void produceBase(edm::Event &e, const edm::EventSetup &es)
virtual void clean(TrajectoryContainer &) const
int j
Definition: DBlmapReader.cc:9
virtual void unset() const
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:356
#define LogTrace(id)
tuple conf
Definition: dbtoconf.py:185
virtual TrajectoryContainer trajectories(const TrajectorySeed &) const =0
virtual void beginRunBase(edm::Run &, edm::EventSetup const &es)
Definition: DetId.h:20
bool isValid() const
Definition: Trajectory.h:259
ConstRecHitContainer RecHitContainer
Definition: Trajectory.h:44
const T & get() const
Definition: EventSetup.h:55
T const * product() const
Definition: ESHandle.h:62
TransientInitialStateEstimator * theInitialState
void setSeedRef(const edm::RefToBase< TrajectorySeed > &seedRef)
Definition: Trajectory.h:298
char state
Definition: procUtils.cc:75
void setEventSetup(const edm::EventSetup &es)
Initialize EventSetup objects at each event.
virtual void setEvent(const edm::Event &event) const =0
void push(const TrajectoryMeasurement &tm)
Definition: Trajectory.cc:35
std::pair< TrajectoryStateOnSurface, const GeomDet * > innerState(const Trajectory &traj, bool doBackFit=true) const
virtual void init(const std::vector< Trajectory > *vect)=0
Provides the cleaner a pointer to the vector where trajectories are stored, in case it does not want ...
void reserve(size_t)
Definition: OwnVector.h:267
Definition: Run.h:33
edm::ESHandle< MagneticField > theMagField