CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TrackCandidateProducer.cc
Go to the documentation of this file.
1 #include <memory>
2 
7 
17 
20 
23 
25 //#include "FastSimulation/Tracking/interface/TrackerRecHitSplit.h"
26 
28 
29 #include <vector>
30 #include <map>
31 
36 
42 
43 //Propagator withMaterial
45 
48 
49 //
50 
51 //for debug only
52 //#define FAMOS_DEBUG
53 
55 {
56 #ifdef FAMOS_DEBUG
57  std::cout << "TrackCandidateProducer created" << std::endl;
58 #endif
59 
60  // The main product is a track candidate collection.
61  produces<TrackCandidateCollection>();
62 
63  // These products contain tracks already reconstructed at this level
64  // (No need to reconstruct them twice!)
65  produces<reco::TrackCollection>();
66  produces<TrackingRecHitCollection>();
67  produces<reco::TrackExtraCollection>();
68  produces<std::vector<Trajectory> >();
69  produces<TrajTrackAssociationCollection>();
70 
71  // The name of the seed producer
72  seedProducer = conf.getParameter<edm::InputTag>("SeedProducer");
73 
74  // The name of the recHit producer
75  hitProducer = conf.getParameter<edm::InputTag>("HitProducer");
76 
77  // The name of the track producer (tracks already produced need not be produced again!)
78  // trackProducer = conf.getParameter<edm::InputTag>("TrackProducer");
79  trackProducers = conf.getParameter<std::vector<edm::InputTag> >("TrackProducers");
80 
81  // Copy (or not) the tracks already produced in a new collection
82  keepFittedTracks = conf.getParameter<bool>("KeepFittedTracks");
83 
84  // The minimum number of crossed layers
85  minNumberOfCrossedLayers = conf.getParameter<unsigned int>("MinNumberOfCrossedLayers");
86 
87  // The maximum number of crossed layers
88  maxNumberOfCrossedLayers = conf.getParameter<unsigned int>("MaxNumberOfCrossedLayers");
89 
90  // Reject overlapping hits?
91  rejectOverlaps = conf.getParameter<bool>("OverlapCleaning");
92 
93  // Split hits ?
94  splitHits = conf.getParameter<bool>("SplitHits");
95 
96  // Reject tracks with several seeds ?
97  // Typically don't do that at HLT for electrons, but do it otherwise
98  seedCleaning = conf.getParameter<bool>("SeedCleaning");
99 
100  simTracks_ = conf.getParameter<edm::InputTag>("SimTracks");
101  estimatorCut_= conf.getParameter<double>("EstimatorCut");
102 }
103 
104 
105 // Virtual destructor needed.
107 
108  if(thePropagator) delete thePropagator;
109 
110  // do nothing
111 #ifdef FAMOS_DEBUG
112  std::cout << "TrackCandidateProducer destructed" << std::endl;
113 #endif
114 
115 
116 }
117 
118 void
120 
121  //services
122 
125 
126  es.get<IdealMagneticFieldRecord>().get(magField);
127  es.get<TrackerDigiGeometryRecord>().get(geometry);
128 
129  theMagField = &(*magField);
130  theGeometry = &(*geometry);
131 
133 }
134 
135  // Functions that gets called by framework every event
136 void
138 
139 #ifdef FAMOS_DEBUG
140  std::cout << "################################################################" << std::endl;
141  std::cout << " TrackCandidateProducer produce init " << std::endl;
142 #endif
143 
144  // Useful typedef's to avoid retyping
145  typedef std::pair<reco::TrackRef,edm::Ref<std::vector<Trajectory> > > TrackPair;
146  typedef std::map<unsigned,TrackPair> TrackMap;
147 
148  // The produced objects
149  std::auto_ptr<TrackCandidateCollection> output(new TrackCandidateCollection);
150  std::auto_ptr<reco::TrackCollection> recoTracks(new reco::TrackCollection);
151  std::auto_ptr<TrackingRecHitCollection> recoHits(new TrackingRecHitCollection);
152  std::auto_ptr<reco::TrackExtraCollection> recoTrackExtras(new reco::TrackExtraCollection);
153  std::auto_ptr<std::vector<Trajectory> > recoTrajectories(new std::vector<Trajectory> );
154  std::auto_ptr<TrajTrackAssociationCollection> recoTrajTrackMap( new TrajTrackAssociationCollection() );
155 
156  // Get the seeds
157  // edm::Handle<TrajectorySeedCollection> theSeeds;
159  e.getByLabel(seedProducer,theSeeds);
160 
161  // No seed -> output an empty track collection
162  if(theSeeds->size() == 0) {
163  e.put(output);
164  e.put(recoTracks);
165  e.put(recoHits);
166  e.put(recoTrackExtras);
167  e.put(recoTrajectories);
168  e.put(recoTrajTrackMap);
169  return;
170  }
171 
172  // Get the GS RecHits
173  // edm::Handle<SiTrackerGSRecHit2DCollection> theGSRecHits;
175  e.getByLabel(hitProducer, theGSRecHits);
176 
177  //get other general things
178  const std::vector<unsigned> theSimTrackIds = theGSRecHits->ids();
179  // SimTracks and SimVertices
180 
182  e.getByLabel("famosSimHits",theSimVtx);
184  e.getByLabel("famosSimHits",theSTC);
185 
186  const edm::SimTrackContainer* theSimTracks = &(*theSTC);
187  LogDebug("FastTracking")<<"looking at: "<< theSimTrackIds.size()<<" simtracks.";
188 
189 
190 
191 
192  // The input track collection + extra's
193  /*
194  edm::Handle<reco::TrackCollection> theTrackCollection;
195  edm:: Handle<std::vector<Trajectory> > theTrajectoryCollection;
196  edm::Handle<TrajTrackAssociationCollection> theAssoMap;
197  bool isTrackCollection = e.getByLabel(trackProducer,theTrackCollection);
198  */
199  std::vector<edm::Handle<reco::TrackCollection> > theTrackCollections;
200  std::vector<edm:: Handle<std::vector<Trajectory> > > theTrajectoryCollections;
201  std::vector<edm::Handle<TrajTrackAssociationCollection> > theAssoMaps;
202  std::vector<bool> isTrackCollections;
205  TrackMap theTrackMap;
206  unsigned nCollections = trackProducers.size();
207  unsigned nRecoHits = 0;
208 
209  if ( nCollections ) {
210  theTrackCollections.resize(nCollections);
211  theTrajectoryCollections.resize(nCollections);
212  theAssoMaps.resize(nCollections);
213  isTrackCollections.resize(nCollections);
214  for ( unsigned tprod=0; tprod < nCollections; ++tprod ) {
215  isTrackCollections[tprod] = e.getByLabel(trackProducers[tprod],theTrackCollections[tprod]);
216 
217  if ( isTrackCollections[tprod] ) {
218  // The track collection
219  reco::TrackCollection::const_iterator aTrack = theTrackCollections[tprod]->begin();
220  reco::TrackCollection::const_iterator lastTrack = theTrackCollections[tprod]->end();
221  // The numbers of hits
222  for ( ; aTrack!=lastTrack; ++aTrack ) nRecoHits+= aTrack->recHitsSize();
223  e.getByLabel(trackProducers[tprod],theTrajectoryCollections[tprod]);
224  e.getByLabel(trackProducers[tprod],theAssoMaps[tprod]);
225  // The association between trajectories and tracks
226  anAssociation = theAssoMaps[tprod]->begin();
227  lastAssociation = theAssoMaps[tprod]->end();
228 #ifdef FAMOS_DEBUG
229  std::cout << "Input Track Producer " << tprod << " : " << trackProducers[tprod] << std::endl;
230  std::cout << "List of tracks already reconstructed " << std::endl;
231 #endif
232  // Build the map of correspondance between reco tracks and sim tracks
233  for ( ; anAssociation != lastAssociation; ++anAssociation ) {
234  edm::Ref<std::vector<Trajectory> > aTrajectoryRef = anAssociation->key;
235  reco::TrackRef aTrackRef = anAssociation->val;
236  // Find the simtrack id of the reconstructed track
237  int recoTrackId = findId(*aTrackRef);
238  if ( recoTrackId < 0 ) continue;
239 #ifdef FAMOS_DEBUG
240  std::cout << recoTrackId << " ";
241 #endif
242  // And store it.
243  theTrackMap[recoTrackId] = TrackPair(aTrackRef,aTrajectoryRef);
244  }
245 #ifdef FAMOS_DEBUG
246  std::cout << std::endl;
247 #endif
248  }
249  }
250  // This is to save some time at push_back.
251  recoHits->reserve(nRecoHits);
252  }
253 
254 
255  // Loop over the seeds
256  int currentTrackId = -1;
257  /*
258  TrajectorySeedCollection::const_iterator aSeed = theSeeds->begin();
259  TrajectorySeedCollection::const_iterator lastSeed = theSeeds->end();
260  for ( ; aSeed!=lastSeed; ++aSeed ) {
261  // The first hit of the seed and its simtrack id
262  */
263  /* */
264 #ifdef FAMOS_DEBUG
265  std::cout << "Input seed Producer : " << seedProducer << std::endl;
266  std::cout << "Number of seeds : " << theSeeds->size() << std::endl;
267 #endif
268  unsigned seed_size = theSeeds->size();
269  for (unsigned seednr = 0; seednr < seed_size; ++seednr){
270 
271  LogDebug("FastTracking")<<"looking at seed #:"<<seednr;
272 
273  // The seed
274  const BasicTrajectorySeed* aSeed = &((*theSeeds)[seednr]);
275 
276  std::vector<int> simTrackIds;
277  std::map<int,TrajectoryStateOnSurface> seedStates;
278  std::map<int,TrajectoryStateOnSurface> simtkStates;
279 
280  TrackerRecHit theFirstSeedingTrackerRecHit;
281  if (theSeeds->at(seednr).nHits()==0){
282  //new stuff for no hits on seed
283 
284  LogDebug("FastTracking")<<" seed with no hits to be considered.";
285 
286  //moved out of the loop
287  //edm::ESHandle<MagneticField> field;
288  //es.get<IdealMagneticFieldRecord>().get(field);
289 
290  PTrajectoryStateOnDet ptod =theSeeds->at(seednr).startingState();
291  DetId id(ptod.detId());
292  const GeomDet * g = theGeometry->idToDet(id);
293  const Surface * surface=&g->surface();
294  TrajectoryStateTransform tsTransform;
295  TrajectoryStateOnSurface seedState(tsTransform.transientState(ptod,surface,theMagField));
296 
298  es.get<TrackingComponentsRecord>().get("AnyDirectionAnalyticalPropagator",propagator);
299 
300  //moved out of the loop
301  // const std::vector<unsigned> theSimTrackIds = theGSRecHits->ids();
302  // edm::Handle<edm::SimTrackContainer> theSTC;
303  // e.getByLabel(simTracks_,theSTC);
304  // const edm::SimTrackContainer* theSimTracks = &(*theSTC);
305 
306  double minimunEst=1000000;
307  LogDebug("FastTracking")<<"looking at: "<< theSimTrackIds.size()<<" simtracks.";
308  for ( unsigned tkId=0; tkId != theSimTrackIds.size(); ++tkId ) {
309 
310  const SimTrack & simtrack = (*theSimTracks)[theSimTrackIds[tkId]];
311 
313  simtrack.trackerSurfacePosition().y(),
314  simtrack.trackerSurfacePosition().z());
315 
316  GlobalVector momentum(simtrack.trackerSurfaceMomentum().x(),
317  simtrack.trackerSurfaceMomentum().y(),
318  simtrack.trackerSurfaceMomentum().z());
319 
320  if (position.basicVector().dot( momentum.basicVector() ) * seedState.globalPosition().basicVector().dot(seedState.globalMomentum().basicVector()) <0. ){
321  LogDebug("FastTracking")<<"not on the same direction.";
322  continue;
323  }
324 
325  //no charge mis-identification ... FIXME
326  int charge = (int) simtrack.charge();
327  GlobalTrajectoryParameters glb_parameters(position,
328  momentum,
329  charge,
330  theMagField);
331  FreeTrajectoryState simtrack_trackerstate(glb_parameters);
332 
333  TrajectoryStateOnSurface simtrack_comparestate = propagator->propagate(simtrack_trackerstate,*surface);
334 
335 
336  if (!simtrack_comparestate.isValid()){
337  LogDebug("FastTracking")<<" ok this is a state-based seed. simtrack state does not propagate to the seed surface. skipping.";
338  continue;}
339 
340  if (simtrack_comparestate.globalPosition().basicVector().dot(simtrack_comparestate.globalMomentum().basicVector()) * seedState.globalPosition().basicVector().dot(seedState.globalMomentum().basicVector()) <0. ){
341  LogDebug("FastTracking")<<"not on the same direction.";
342  continue;
343  }
344 
345 
346  AlgebraicVector5 v(seedState.localParameters().vector() - simtrack_comparestate.localParameters().vector());
347  AlgebraicSymMatrix55 m(seedState.localError().matrix());
348  bool ierr = !m.Invert();
349  if ( ierr ){
350  edm::LogWarning("FastTracking") <<" Candidate Producer cannot invert the error matrix! - Skipping...";
351  continue;
352  }
353  double est = ROOT::Math::Similarity(v,m);
354  LogDebug("FastTracking")<<"comparing two state on the seed surface:\n"
355  <<"seed: "<<seedState
356  <<"sim: "<<simtrack_comparestate
357  <<"\n estimator is: "<<est;
358 
359  if (est<minimunEst) minimunEst=est;
360  if (est<estimatorCut_){
361  simTrackIds.push_back(theSimTrackIds[tkId]);
362  //making a state with exactly the sim track parameters
363  //the initial errors are set to unity just for kicks
364  // AlgebraicSymMatrix C(5,1);// C*=50;
365  //new attempt!!!!
366  AlgebraicSymMatrix55 C = seedState.curvilinearError().matrix();
367  C *= 0.0000001;
368 
369  seedStates[theSimTrackIds[tkId]] = TrajectoryStateOnSurface(simtrack_comparestate.globalParameters(),
371  seedState.surface());
372  LogDebug("FastTracking")<<"the compatibility estimator is: "<<est<<" for track id: "<<simTrackIds.back();
373  }
374  }//SimTrack loop
375  if (simTrackIds.size()==0) LogDebug("FastTracking")<<"could not find any simtrack within errors, closest was at: "<<minimunEst;
376  }//seed has 0 hit.
377  else{
378  //same old stuff
379  // Find the first hit of the Seed
380  TrajectorySeed::range theSeedingRecHitRange = aSeed->recHits();
381  const SiTrackerGSMatchedRecHit2D * theFirstSeedingRecHit = (const SiTrackerGSMatchedRecHit2D*) (&(*(theSeedingRecHitRange.first)));
382  theFirstSeedingTrackerRecHit = TrackerRecHit(theFirstSeedingRecHit,theGeometry);
383  // The SimTrack id associated to that recHit
384  simTrackIds.push_back( theFirstSeedingRecHit->simtrackId() );
385  }
386 
387  //from then on, only the simtrack IDs are usefull.
388  //now loop over all possible trackid for this seed.
389  //an actual seed can be shared by two tracks in dense envirronement, and also for hit-less seeds.
390  for (unsigned int iToMake=0;iToMake!=simTrackIds.size();++iToMake){
391  int simTrackId = simTrackIds[iToMake];
392 
393  // Don't consider seeds belonging to a track already considered
394  // (Equivalent to seed cleaning)
395  if ( seedCleaning && currentTrackId == simTrackId ) continue;
396  currentTrackId = simTrackId;
397 
398  // A vector of TrackerRecHits belonging to the track and the number of crossed layers
399  std::vector<TrackerRecHit> theTrackerRecHits;
400  unsigned theNumberOfCrossedLayers = 0;
401 
402  // The track has indeed been reconstructed already -> Save the pertaining info
403  TrackMap::const_iterator theTrackIt = theTrackMap.find(simTrackId);
404  if ( nCollections && theTrackIt != theTrackMap.end() ) {
405 
406  if ( keepFittedTracks ) {
407  LogDebug("FastTracking") << "Track " << simTrackId << " already reconstructed -> copy it";
408  // The track and trajectroy references
409  reco::TrackRef aTrackRef = theTrackIt->second.first;
410  edm::Ref<std::vector<Trajectory> > aTrajectoryRef = theTrackIt->second.second;
411 
412  // A copy of the track
413  reco::Track aRecoTrack(*aTrackRef);
414  recoTracks->push_back(aRecoTrack);
415 
416  // A copy of the hits
417  unsigned nh = aRecoTrack.recHitsSize();
418  for ( unsigned ih=0; ih<nh; ++ih ) {
419  TrackingRecHit *hit = aRecoTrack.recHit(ih)->clone();
420  recoHits->push_back(hit);
421  }
422 
423  // A copy of the trajectories
424  recoTrajectories->push_back(*aTrajectoryRef);
425 
426  }// keepFitterTracks
427  else {
428  LogDebug("FastTracking") << "Track " << simTrackId << " already reconstructed -> ignore it";
429  }
430 
431  // The track was not saved -> create a track candidate.
432 
433  } //already existing collection of tracks
434  else{//no collection of tracks already exists
435 
436  LogDebug("FastTracking")<<"Track " << simTrackId << " is considered to return a track candidate" ;
437 
438  // Get all the rechits associated to this track
439  SiTrackerGSMatchedRecHit2DCollection::range theRecHitRange = theGSRecHits->get(simTrackId);
440  SiTrackerGSMatchedRecHit2DCollection::const_iterator theRecHitRangeIteratorBegin = theRecHitRange.first;
441  SiTrackerGSMatchedRecHit2DCollection::const_iterator theRecHitRangeIteratorEnd = theRecHitRange.second;
443 
444  LogDebug("FastTracking")<<"counting: "<<theRecHitRangeIteratorEnd-theRecHitRangeIteratorBegin<<" hits to be considered.";
445 
446  bool firstRecHit = true;
447  TrackerRecHit theCurrentRecHit, thePreviousRecHit;
448  TrackerRecHit theFirstHitComp, theSecondHitComp;
449 
450  for ( iterRecHit = theRecHitRangeIteratorBegin;
451  iterRecHit != theRecHitRangeIteratorEnd;
452  ++iterRecHit) {
453 
454  // Check the number of crossed layers
455  if ( theNumberOfCrossedLayers >= maxNumberOfCrossedLayers ) continue;
456 
457  // Get current and previous rechits
458  if(!firstRecHit) thePreviousRecHit = theCurrentRecHit;
459  theCurrentRecHit = TrackerRecHit(&(*iterRecHit),theGeometry);
460 
461  //>>>>>>>>>BACKBUILDING CHANGE: DO NOT STAT FROM THE FIRST HIT OF THE SEED
462 
463  // NOTE: checking the direction --> specific for OIHit only
464  // if( aSeed->direction()!=oppositeToMomentum ) {
465  // // Check that the first rechit is indeed the first seeding hit
466  // if ( firstRecHit && theCurrentRecHit != theFirstSeedingTrackerRecHit && theSeeds->at(seednr).nHits()!=0 ) continue;
467  // }
468  //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
469 
470  // Count the number of crossed layers
471  if ( !theCurrentRecHit.isOnTheSameLayer(thePreviousRecHit) )
472  ++theNumberOfCrossedLayers;
473 
474  // Add all rechits (Grouped Trajectory Builder) from this hit onwards
475  // Always add the first seeding rechit anyway
476  if ( !rejectOverlaps || firstRecHit ) {
477  // Split matched hits (if requested / possible )
478  if ( splitHits && theCurrentRecHit.matchedHit()->isMatched() ) addSplitHits(theCurrentRecHit,theTrackerRecHits);
479  else theTrackerRecHits.push_back(theCurrentRecHit); // No splitting
480  firstRecHit = false;
481 
482  // And now treat the following RecHits if hits in the same layer
483  // have to be rejected - The split option is not
484  } else {
485 
486  // Not the same layer : Add the current hit
487  if ( theCurrentRecHit.subDetId() != thePreviousRecHit.subDetId() ||
488  theCurrentRecHit.layerNumber() != thePreviousRecHit.layerNumber() ) {
489 
490  // Split matched hits (if requested / possible )
491  if ( splitHits && theCurrentRecHit.matchedHit()->isMatched() ) addSplitHits(theCurrentRecHit,theTrackerRecHits);
492  else theTrackerRecHits.push_back(theCurrentRecHit); // No splitting
493 
494  // Same layer : keep the current hit if better, and drop the other - otherwise do nothing
495  } else if ( theCurrentRecHit.localError() < thePreviousRecHit.localError() ) {
496 
497  // Split matched hits (if requested / possible )
498  if( splitHits && theCurrentRecHit.matchedHit()->isMatched() ){
499  // Remove the previous hit(s)
500  theTrackerRecHits.pop_back();
501  if ( thePreviousRecHit.matchedHit()->isMatched() ) theTrackerRecHits.pop_back();
502 
503  // Replace by the new hits
504  addSplitHits(theCurrentRecHit,theTrackerRecHits);
505  }
506  // No splitting
507  else {
508  theTrackerRecHits.back() = theCurrentRecHit; // Replace the previous hit by the current hit
509  }
510 
511  } else {
512 
513  //keep the older rechit as a reference if the error of the new one is worse
514  theCurrentRecHit = thePreviousRecHit;
515  }
516  }
517  }
518 
519  // End of loop over the track rechits
520  }//no collection of track already existed. adding the hits by hand.
521 
522  LogDebug("FastTracking")<<" number of hits: " << theTrackerRecHits.size()<<" after counting overlaps and splitting.";
523 
524  // 1) Create the OwnWector of TrackingRecHits
526  unsigned nTrackerHits = theTrackerRecHits.size();
527  recHits.reserve(nTrackerHits); // To save some time at push_back
528 
529  if (aSeed->direction()==oppositeToMomentum){
530  LogDebug("FastTracking")<<"reversing the order of the hits";
531  std::reverse(theTrackerRecHits.begin(),theTrackerRecHits.end());
532  }
533 
534  for ( unsigned ih=0; ih<nTrackerHits; ++ih ) {
535  TrackingRecHit* aTrackingRecHit = theTrackerRecHits[ih].hit()->clone();
536  recHits.push_back(aTrackingRecHit);
537 
538  const DetId& detId = theTrackerRecHits[ih].hit()->geographicalId();
539  LogDebug("FastTracking")
540  << "Added RecHit from detid " << detId.rawId()
541  << " subdet = " << theTrackerRecHits[ih].subDetId()
542  << " layer = " << theTrackerRecHits[ih].layerNumber()
543  << " ring = " << theTrackerRecHits[ih].ringNumber()
544  << " error = " << theTrackerRecHits[ih].localError()
545  << std::endl
546 
547  << "Track/z/r : "
548  << simTrackId << " "
549  << theTrackerRecHits[ih].globalPosition().z() << " "
550  << theTrackerRecHits[ih].globalPosition().perp() << std::endl;
551  if ( theTrackerRecHits[ih].matchedHit() && theTrackerRecHits[ih].matchedHit()->isMatched() )
552  LogTrace("FastTracking") << "Matched : " << theTrackerRecHits[ih].matchedHit()->isMatched()
553  << "Rphi Hit = " << theTrackerRecHits[ih].matchedHit()->monoHit()->simhitId()
554  << "Stereo Hit = " << theTrackerRecHits[ih].matchedHit()->stereoHit()->simhitId()
555  <<std::endl;
556  }//loop over the rechits
557 
558  // Check the number of crossed layers
559  if ( theNumberOfCrossedLayers < minNumberOfCrossedLayers ) {
560  LogDebug("FastTracking")<<"not enough layer crossed ("<<theNumberOfCrossedLayers<<")";
561  continue;
562  }
563 
564 
565  //>>>>>>>>>BACKBUILDING CHANGE: REPLACE THE STARTING STATE
566 
567  // Create a track Candidate (now with the reference to the seed!) .
568  //PTrajectoryStateOnDet PTSOD = aSeed->startingState();
569  PTrajectoryStateOnDet PTSOD;
570 
571  if (aSeed->nHits()==0){
572  //stabilize the fit with the true simtrack state
573  //in case of zero hits
574  TrajectoryStateTransform tsTransform;
575  //PTSOD = *tsTransform.persistentState(seedStates[simTrackId],aSeed->startingState().detId());
576  PTrajectoryStateOnDet * aPointer = tsTransform.persistentState(seedStates[simTrackId],aSeed->startingState().detId());
577  PTSOD = *aPointer;
578  delete aPointer;
579 
580  } else {
581  //create the initial state from the SimTrack
582  int vertexIndex = (*theSimTracks)[currentTrackId].vertIndex();
583  // a) origin vertex
584  GlobalPoint position((*theSimVtx)[vertexIndex].position().x(),
585  (*theSimVtx)[vertexIndex].position().y(),
586  (*theSimVtx)[vertexIndex].position().z());
587 
588  // b) initial momentum
589  GlobalVector momentum( (*theSimTracks)[currentTrackId].momentum().x() ,
590  (*theSimTracks)[currentTrackId].momentum().y() ,
591  (*theSimTracks)[currentTrackId].momentum().z() );
592  // c) electric charge
593  float charge = (*theSimTracks)[simTrackId].charge();
594  // -> inital parameters
595  GlobalTrajectoryParameters initialParams(position,momentum,(int)charge,theMagField);
596  // -> large initial errors
597  AlgebraicSymMatrix55 errorMatrix= AlgebraicMatrixID();
598  CurvilinearTrajectoryError initialError(errorMatrix);
599  // -> initial state
600  FreeTrajectoryState initialFTS(initialParams, initialError);
601 #ifdef FAMOS_DEBUG
602  std::cout << "TrajectorySeedProducer: FTS momentum " << initialFTS.momentum() << std::endl;
603 #endif
604  const GeomDet* initialLayer = theGeometry->idToDet(recHits.front().geographicalId());
605  //this is wrong because the FTS is defined at vertex, and it need to be properly propagated to the first rechit
606  // const TrajectoryStateOnSurface initialTSOS(initialFTS, initialLayer->surface());
607  const TrajectoryStateOnSurface initialTSOS = thePropagator->propagate(initialFTS,initialLayer->surface()) ;
608  if (!initialTSOS.isValid()) continue;
609  TrajectoryStateTransform tsTransform;
610 
611  PTrajectoryStateOnDet * aPointer = tsTransform.persistentState(initialTSOS,recHits.front().geographicalId().rawId());
612  PTSOD = *aPointer;
613  delete aPointer;
614  }
615 
617  newTrackCandidate(recHits,
618  *aSeed,
619  PTSOD,
620  edm::RefToBase<TrajectorySeed>(theSeeds,seednr));
621 
622  LogDebug("FastTracking")<< "\tSeed Information " << std::endl
623  << "\tSeed Direction = " << aSeed->direction() << std::endl
624  << "\tSeed StartingDet = " << aSeed->startingState().detId() << std::endl
625  << "\tTrajectory Parameters " << std::endl
626  << "\t\t detId = " << newTrackCandidate.trajectoryStateOnDet().detId() << std::endl
627  << "\t\t loc.px = "
628  << newTrackCandidate.trajectoryStateOnDet().parameters().momentum().x()
629  << std::endl
630  << "\t\t loc.py = "
631  << newTrackCandidate.trajectoryStateOnDet().parameters().momentum().y()
632  << std::endl
633  << "\t\t loc.pz = "
634  << newTrackCandidate.trajectoryStateOnDet().parameters().momentum().z()
635  << std::endl
636  << "\t\t error = ";
637  //<< newTrackCandidate.trajectoryStateOnDet().errorMatrix()<<std::endl;
638  // for(std::vector< float >::const_iterator iElement = newTrackCandidate.trajectoryStateOnDet().errorMatrix().begin();
639  // iElement < newTrackCandidate.trajectoryStateOnDet().errorMatrix().end();
640  // ++iElement) {
641  // std::cout << "\t" << *iElement;
642  // }
643 
644  output->push_back(newTrackCandidate);
645  LogDebug("FastTracking")<<"filling a track candidate into the collection, now having: "<<output->size();
646 
647  }//loop over possible simtrack associated.
648  }//loop over all possible seeds.
649 
650  // Save the track candidates in the event
651  LogDebug("FastTracking") << "Saving "
652  << output->size() << " track candidates and "
653  << recoTracks->size() << " reco::Tracks ";
654  // Save the track candidates
655  e.put(output);
656 
657 
658 
659  // Save the tracking recHits
660 
661  edm::OrphanHandle <TrackingRecHitCollection> theRecoHits = e.put(recoHits );
662 
663  // Create the track extras and add the references to the rechits
664  unsigned hits=0;
665  unsigned nTracks = recoTracks->size();
666  recoTrackExtras->reserve(nTracks); // To save some time at push_back
667  for ( unsigned index = 0; index < nTracks; ++index ) {
668  //reco::TrackExtra aTrackExtra;
669  reco::Track& aTrack = recoTracks->at(index);
670  reco::TrackExtra aTrackExtra(aTrack.outerPosition(),
671  aTrack.outerMomentum(),
672  aTrack.outerOk(),
673  aTrack.innerPosition(),
674  aTrack.innerMomentum(),
675  aTrack.innerOk(),
676  aTrack.outerStateCovariance(),
677  aTrack.outerDetId(),
678  aTrack.innerStateCovariance(),
679  aTrack.innerDetId(),
680  aTrack.seedDirection(),
681  aTrack.seedRef());
682 
683  unsigned nHits = aTrack.recHitsSize();
684  for ( unsigned int ih=0; ih<nHits; ++ih) {
685  aTrackExtra.add(TrackingRecHitRef(theRecoHits,hits++));
686  }
687  recoTrackExtras->push_back(aTrackExtra);
688  }
689 
690 
691  // Save the track extras
692  edm::OrphanHandle<reco::TrackExtraCollection> theRecoTrackExtras = e.put(recoTrackExtras);
693 
694  // Add the reference to the track extra in the tracks
695  for ( unsigned index = 0; index<nTracks; ++index ) {
696  const reco::TrackExtraRef theTrackExtraRef(theRecoTrackExtras,index);
697  (recoTracks->at(index)).setExtra(theTrackExtraRef);
698  }
699 
700  // Save the tracks
701  edm::OrphanHandle<reco::TrackCollection> theRecoTracks = e.put(recoTracks);
702 
703  // Save the trajectories
704  edm::OrphanHandle<std::vector<Trajectory> > theRecoTrajectories = e.put(recoTrajectories);
705 
706  // Create and set the trajectory/track association map
707  for ( unsigned index = 0; index<nTracks; ++index ) {
708  edm::Ref<std::vector<Trajectory> > trajRef( theRecoTrajectories, index );
709  edm::Ref<reco::TrackCollection> tkRef( theRecoTracks, index );
710  recoTrajTrackMap->insert(trajRef,tkRef);
711  }
712 
713 
714  // Save the association map.
715  e.put(recoTrajTrackMap);
716 
717 }
718 
719 int
721  int trackId = -1;
722  trackingRecHit_iterator aHit = aTrack.recHitsBegin();
723  trackingRecHit_iterator lastHit = aTrack.recHitsEnd();
724  for ( ; aHit!=lastHit; ++aHit ) {
725  if ( !aHit->get()->isValid() ) continue;
726  // const SiTrackerGSRecHit2D * rechit = (const SiTrackerGSRecHit2D*) (aHit->get());
727  const SiTrackerGSMatchedRecHit2D * rechit = (const SiTrackerGSMatchedRecHit2D*) (aHit->get());
728  trackId = rechit->simtrackId();
729  break;
730  }
731  return trackId;
732 }
733 
734 void
736  std::vector<TrackerRecHit>& theTrackerRecHits) {
737 
738  const SiTrackerGSRecHit2D* mHit = theCurrentRecHit.matchedHit()->monoHit();
739  const SiTrackerGSRecHit2D* sHit = theCurrentRecHit.matchedHit()->stereoHit();
740 
741  // Add the new hits
742  if( mHit->simhitId() < sHit->simhitId() ) {
743 
744  theTrackerRecHits.push_back(TrackerRecHit(mHit,theCurrentRecHit));
745  theTrackerRecHits.push_back(TrackerRecHit(sHit,theCurrentRecHit));
746 
747  } else {
748 
749  theTrackerRecHits.push_back(TrackerRecHit(sHit,theCurrentRecHit));
750  theTrackerRecHits.push_back(TrackerRecHit(mHit,theCurrentRecHit));
751 
752  }
753 
754 }
#define LogDebug(id)
PropagationDirection direction() const
const math::XYZVectorD & trackerSurfacePosition() const
Definition: SimTrack.h:36
T getParameter(std::string const &) const
std::pair< const_iterator, const_iterator > range
iterator range
Definition: RangeMap.h:53
const int & simhitId() const
const LocalTrajectoryParameters & localParameters() const
size_t recHitsSize() const
Get number of RecHits. (Warning, this includes invalid hits, which are not physical hits)...
Definition: Track.h:69
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
TrackCandidateProducer(const edm::ParameterSet &conf)
int findId(const reco::Track &aTrack) const
std::vector< TrackCandidate > TrackCandidateCollection
size_type size() const
Definition: OwnVector.h:262
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:10
T y() const
Definition: PV3DBase.h:57
GlobalPoint globalPosition() const
bool innerOk() const
return true if the innermost hit is valid
Definition: Track.h:40
const SiTrackerGSRecHit2D * stereoHit() const
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
const math::XYZPoint & outerPosition() const
position of the outermost hit
Definition: Track.h:47
float charge() const
charge
Definition: CoreSimTrack.cc:3
const MagneticField * theMagField
double charge(const std::vector< uint8_t > &Ampls)
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
double double double z
static int position[TOTALCHAMBERS][3]
Definition: ReadPGInfo.cc:509
bool isOnTheSameLayer(const TrackerRecHit &other) const
Check if two hits are on the same layer of the same subdetector.
const math::XYZPoint & innerPosition() const
position of the innermost hit
Definition: Track.h:42
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
PTrajectoryStateOnDet const & trajectoryStateOnDet() const
C::const_iterator const_iterator
constant access iterator type
Definition: RangeMap.h:46
AlgebraicVector5 vector() const
void push_back(D *&d)
Definition: OwnVector.h:290
virtual void beginRun(edm::Run &run, const edm::EventSetup &es)
const SiTrackerGSMatchedRecHit2D * matchedHit() const
The Hit itself.
Definition: TrackerRecHit.h:70
const SiTrackerGSRecHit2D * monoHit() const
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:84
PropagatorWithMaterial * thePropagator
edm::Ref< TrackingRecHitCollection > TrackingRecHitRef
persistent reference to a TrackingRecHit
CovarianceMatrix outerStateCovariance() const
outermost trajectory state curvilinear errors
Definition: Track.h:51
T z() const
Definition: PV3DBase.h:58
std::vector< edm::InputTag > trackProducers
unsigned int outerDetId() const
DetId of the detector on which surface the outermost state is located.
Definition: Track.h:59
std::pair< const_iterator, const_iterator > range
TrajectoryStateOnSurface transientState(const PTrajectoryStateOnDet &ts, const Surface *surface, const MagneticField *field) const
PTrajectoryStateOnDet * persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid) const
LocalVector momentum() const
Momentum vector in the local frame.
trackingRecHit_iterator recHitsBegin() const
Iterator to first hit on the track.
Definition: Track.h:63
edm::AssociationMap< edm::OneToOne< std::vector< Trajectory >, reco::TrackCollection, unsigned short > > TrajTrackAssociationCollection
virtual const GeomDet * idToDet(DetId) const
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:355
GlobalVector momentum() const
#define LogTrace(id)
tuple conf
Definition: dbtoconf.py:185
virtual TrackingRecHit * clone() const =0
virtual TrajectoryStateOnSurface propagate(const TrajectoryStateOnSurface &tsos, const Plane &plane) const
Definition: DetId.h:20
std::vector< TrackExtra > TrackExtraCollection
collection of TrackExtra objects
Definition: TrackExtraFwd.h:9
edm::RefToBase< TrajectorySeed > seedRef() const
Definition: Track.h:112
PTrajectoryStateOnDet const & startingState() const
const math::XYZVector & outerMomentum() const
momentum vector at the outermost hit position
Definition: Track.h:49
ROOT::Math::SVector< double, 5 > AlgebraicVector5
bool outerOk() const
return true if the outermost hit is valid
Definition: Track.h:38
const GlobalTrajectoryParameters & globalParameters() const
const math::XYZTLorentzVectorD & trackerSurfaceMomentum() const
Definition: SimTrack.h:38
const T & get() const
Definition: EventSetup.h:55
const unsigned int detId() const
CovarianceMatrix innerStateCovariance() const
innermost trajectory state curvilinear errors
Definition: Track.h:53
range recHits() const
unsigned int layerNumber() const
The Layer Number.
Definition: TrackerRecHit.h:81
ESHandle< TrackerGeometry > geometry
unsigned int nHits() const
GlobalVector globalMomentum() const
virtual void produce(edm::Event &e, const edm::EventSetup &es)
const math::XYZVector & innerMomentum() const
momentum vector at the innermost hit position
Definition: Track.h:45
TrackingRecHitRef recHit(size_t i) const
Get i-th hit on the track.
Definition: Track.h:67
unsigned int subDetId() const
The subdet Id.
Definition: TrackerRecHit.h:78
double localError()
PropagationDirection seedDirection() const
direction of how the hits were sorted in the original seed
Definition: Track.h:105
tuple cout
Definition: gather_cfg.py:41
const TrackerGeometry * theGeometry
DetId geographicalId() const
Definition: DDAxes.h:10
T x() const
Definition: PV3DBase.h:56
const BasicVectorType & basicVector() const
Definition: PV3DBase.h:54
unsigned int innerDetId() const
DetId of the detector on which surface the innermost state is located.
Definition: Track.h:61
std::vector< SimTrack > SimTrackContainer
reference front()
Definition: OwnVector.h:373
virtual const BoundPlane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:37
mathSSE::Vec4< T > v
void reserve(size_t)
Definition: OwnVector.h:284
const LocalTrajectoryParameters & parameters() const
Definition: Run.h:32
void addSplitHits(const TrackerRecHit &, std::vector< TrackerRecHit > &)
T dot(const Basic3DVector &rh) const
Scalar product, or &quot;dot&quot; product, with a vector of same type.
trackingRecHit_iterator recHitsEnd() const
Iterator to last hit on the track.
Definition: Track.h:65