CMS 3D CMS Logo

Public Member Functions | Private Member Functions | Private Attributes

TrackCandidateProducer Class Reference

#include <TrackCandidateProducer.h>

Inheritance diagram for TrackCandidateProducer:
edm::EDProducer edm::ProducerBase edm::EDConsumerBase edm::ProductRegistryHelper

List of all members.

Public Member Functions

virtual void beginRun (edm::Run const &run, const edm::EventSetup &es) override
virtual void produce (edm::Event &e, const edm::EventSetup &es) override
 TrackCandidateProducer (const edm::ParameterSet &conf)
virtual ~TrackCandidateProducer ()

Private Member Functions

void addSplitHits (const TrackerRecHit &, std::vector< TrackerRecHit > &)
int findId (const reco::Track &aTrack) const
bool isDuplicateCandidate (const TrackCandidateCollection &candidates, const TrackCandidate &newCand) const
bool sameLocalParameters (const TrackingRecHit *aH, const TrackingRecHit *bH) const

Private Attributes

double estimatorCut_
edm::InputTag hitProducer
bool keepFittedTracks
unsigned int maxNumberOfCrossedLayers
unsigned int minNumberOfCrossedLayers
bool rejectOverlaps
bool seedCleaning
edm::InputTag seedProducer
edm::InputTag simTracks_
bool splitHits
const TrackerGeometrytheGeometry
const MagneticFieldtheMagField
PropagatorWithMaterialthePropagator
std::vector< edm::InputTagtrackProducers

Detailed Description

Definition at line 28 of file TrackCandidateProducer.h.


Constructor & Destructor Documentation

TrackCandidateProducer::TrackCandidateProducer ( const edm::ParameterSet conf) [explicit]

Definition at line 54 of file TrackCandidateProducer.cc.

References gather_cfg::cout, estimatorCut_, edm::ParameterSet::getParameter(), hitProducer, keepFittedTracks, maxNumberOfCrossedLayers, minNumberOfCrossedLayers, rejectOverlaps, seedCleaning, seedProducer, simTracks_, splitHits, and trackProducers.

                                                                         :thePropagator(0) 
{  
#ifdef FAMOS_DEBUG
  std::cout << "TrackCandidateProducer created" << std::endl;
#endif

  // The main product is a track candidate collection.
  produces<TrackCandidateCollection>();

  // These products contain tracks already reconstructed at this level
  // (No need to reconstruct them twice!)
  produces<reco::TrackCollection>();
  produces<TrackingRecHitCollection>();
  produces<reco::TrackExtraCollection>();
  produces<std::vector<Trajectory> >();
  produces<TrajTrackAssociationCollection>();
  
  // The name of the seed producer
  seedProducer = conf.getParameter<edm::InputTag>("SeedProducer");

  // The name of the recHit producer
  hitProducer = conf.getParameter<edm::InputTag>("HitProducer");

  // The name of the track producer (tracks already produced need not be produced again!)
  // trackProducer = conf.getParameter<edm::InputTag>("TrackProducer");
  trackProducers = conf.getParameter<std::vector<edm::InputTag> >("TrackProducers");

  // Copy (or not) the tracks already produced in a new collection
  keepFittedTracks = conf.getParameter<bool>("KeepFittedTracks");

  // The minimum number of crossed layers
  minNumberOfCrossedLayers = conf.getParameter<unsigned int>("MinNumberOfCrossedLayers");

  // The maximum number of crossed layers
  maxNumberOfCrossedLayers = conf.getParameter<unsigned int>("MaxNumberOfCrossedLayers");

  // Reject overlapping hits?
  rejectOverlaps = conf.getParameter<bool>("OverlapCleaning");

  // Split hits ?
  splitHits = conf.getParameter<bool>("SplitHits");

  // Reject tracks with several seeds ?
  // Typically don't do that at HLT for electrons, but do it otherwise
  seedCleaning = conf.getParameter<bool>("SeedCleaning");

  simTracks_ = conf.getParameter<edm::InputTag>("SimTracks");
  estimatorCut_= conf.getParameter<double>("EstimatorCut");
}
TrackCandidateProducer::~TrackCandidateProducer ( ) [virtual]

Definition at line 106 of file TrackCandidateProducer.cc.

References gather_cfg::cout, and thePropagator.

                                                {

  if(thePropagator) delete thePropagator;

  // do nothing
#ifdef FAMOS_DEBUG
  std::cout << "TrackCandidateProducer destructed" << std::endl;
#endif


} 

Member Function Documentation

void TrackCandidateProducer::addSplitHits ( const TrackerRecHit theCurrentRecHit,
std::vector< TrackerRecHit > &  theTrackerRecHits 
) [private]

Definition at line 731 of file TrackCandidateProducer.cc.

References TrackerRecHit::matchedHit(), SiTrackerGSMatchedRecHit2D::monoHit(), SiTrackerGSRecHit2D::simhitId(), and SiTrackerGSMatchedRecHit2D::stereoHit().

Referenced by produce().

                                                                                  { 
  
  const SiTrackerGSRecHit2D* mHit = theCurrentRecHit.matchedHit()->monoHit();
  const SiTrackerGSRecHit2D* sHit = theCurrentRecHit.matchedHit()->stereoHit();
  
  // Add the new hits
  if( mHit->simhitId() < sHit->simhitId() ) {
    
    theTrackerRecHits.push_back(TrackerRecHit(mHit,theCurrentRecHit));
    theTrackerRecHits.push_back(TrackerRecHit(sHit,theCurrentRecHit));
    
  } else {
    
    theTrackerRecHits.push_back(TrackerRecHit(sHit,theCurrentRecHit));
    theTrackerRecHits.push_back(TrackerRecHit(mHit,theCurrentRecHit));
    
  }

}
void TrackCandidateProducer::beginRun ( edm::Run const &  run,
const edm::EventSetup es 
) [override, virtual]
int TrackCandidateProducer::findId ( const reco::Track aTrack) const [private]

Definition at line 716 of file TrackCandidateProducer.cc.

References reco::Track::recHitsBegin(), reco::Track::recHitsEnd(), and SiTrackerGSMatchedRecHit2D::simtrackId().

Referenced by produce().

                                                            {
  int trackId = -1;
  trackingRecHit_iterator aHit = aTrack.recHitsBegin();
  trackingRecHit_iterator lastHit = aTrack.recHitsEnd();
  for ( ; aHit!=lastHit; ++aHit ) {
    if ( !aHit->get()->isValid() ) continue;
    //    const SiTrackerGSRecHit2D * rechit = (const SiTrackerGSRecHit2D*) (aHit->get());
    const SiTrackerGSMatchedRecHit2D * rechit = (const SiTrackerGSMatchedRecHit2D*) (aHit->get());
    trackId = rechit->simtrackId();
    break;
  }
  return trackId;
}
bool TrackCandidateProducer::isDuplicateCandidate ( const TrackCandidateCollection candidates,
const TrackCandidate newCand 
) const [private]

Definition at line 752 of file TrackCandidateProducer.cc.

References LocalTrajectoryParameters::dxdz(), LogDebug, PTrajectoryStateOnDet::parameters(), LocalTrajectoryParameters::qbp(), TrackCandidate::recHits(), sameLocalParameters(), and TrackCandidate::trajectoryStateOnDet().

Referenced by produce().

                                                                                                                                {
  typedef TrackCandidateCollection::const_iterator TCCI;
  
  TCCI candsEnd = candidates.end();
  double newQbp = newCand.trajectoryStateOnDet().parameters().qbp();
  double newDxdz = newCand.trajectoryStateOnDet().parameters().dxdz();
  TrackCandidate::range newHits = newCand.recHits();

  for (TCCI iCand = candidates.begin(); iCand!= candsEnd; ++iCand){
    //pick some first traits of duplication: qbp and dxdz
    double iQbp = iCand->trajectoryStateOnDet().parameters().qbp();
    double iDxdz = iCand->trajectoryStateOnDet().parameters().dxdz();
    if (newQbp == iQbp && newDxdz == iDxdz){
      LogDebug("isDuplicateCandidate")<<"Probably a duplicate "<<iQbp <<" : "<<iDxdz;
      TrackCandidate::range iHits = iCand->recHits();
      unsigned int nHits = 0;
      unsigned int nShared = 0;
      unsigned int nnewHits = 0;
      for (TrackCandidate::const_iterator iiHit = iHits.first; iiHit != iHits.second; ++iiHit){
        nHits++;
        for (TrackCandidate::const_iterator inewHit = newHits.first; inewHit != newHits.second; ++inewHit){
          if (iiHit == iHits.first) nnewHits++;
          if (sameLocalParameters(&*iiHit,&*inewHit)) nShared++;
        }
      }
      LogDebug("isDuplicateCandidate") <<"nHits  "<<nHits<<" nnewHits "<< nnewHits<< " nShared "<<nShared;
      if (nHits == nShared && nHits == nnewHits) return true;
    }
  }
  return false;
}
void TrackCandidateProducer::produce ( edm::Event e,
const edm::EventSetup es 
) [override, virtual]

Implements edm::EDProducer.

Definition at line 137 of file TrackCandidateProducer.cc.

References addSplitHits(), PV3DBase< T, PVType, FrameType >::basicVector(), funct::C, CoreSimTrack::charge(), DeDxDiscriminatorTools::charge(), TrackingRecHit::clone(), gather_cfg::cout, TrajectoryStateOnSurface::curvilinearError(), PTrajectoryStateOnDet::detId(), TrajectorySeed::direction(), estimatorCut_, findId(), edm::OwnVector< T, P >::front(), g, edm::EventSetup::get(), edm::Event::getByLabel(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalParameters(), TrajectoryStateOnSurface::globalPosition(), hitProducer, TrackerGeometry::idToDet(), getHLTprescales::index, reco::Track::innerDetId(), reco::Track::innerMomentum(), reco::Track::innerOk(), reco::Track::innerPosition(), reco::Track::innerStateCovariance(), isDuplicateCandidate(), SiTrackerGSMatchedRecHit2D::isMatched(), trackerHitRTTI::isMatched(), TrackerRecHit::isOnTheSameLayer(), TrajectoryStateOnSurface::isValid(), keepFittedTracks, TrackerRecHit::layerNumber(), TrajectoryStateOnSurface::localError(), TrackerRecHit::localError(), TrajectoryStateOnSurface::localParameters(), LogDebug, LogTrace, m, TrackerRecHit::matchedHit(), CurvilinearTrajectoryError::matrix(), LocalTrajectoryError::matrix(), maxNumberOfCrossedLayers, minNumberOfCrossedLayers, LocalTrajectoryParameters::momentum(), FreeTrajectoryState::momentum(), ecaldqm::nCollections, TrajectorySeed::nHits(), oppositeToMomentum, reco::Track::outerDetId(), reco::Track::outerMomentum(), reco::Track::outerOk(), reco::Track::outerPosition(), reco::Track::outerStateCovariance(), convertSQLitetoXML_cfg::output, PTrajectoryStateOnDet::parameters(), trajectoryStateTransform::persistentState(), position, edm::ESHandle< T >::product(), LargeD0_PixelPairStep_cff::propagator, edm::OwnVector< T, P >::push_back(), edm::Event::put(), DetId::rawId(), reco::Track::recHit(), TrajectorySeed::recHits(), reco::Track::recHitsSize(), rejectOverlaps, edm::OwnVector< T, P >::reserve(), seedCleaning, reco::Track::seedDirection(), seedProducer, reco::Track::seedRef(), SiTrackerGSMatchedRecHit2D::simtrackId(), edm::OwnVector< T, P >::size(), splitHits, TrajectorySeed::startingState(), TrackerRecHit::subDetId(), GeomDet::surface(), TrajectoryStateOnSurface::surface(), theGeometry, theMagField, thePropagator, SimTrack::trackerSurfaceMomentum(), SimTrack::trackerSurfacePosition(), trackProducers, TrackCandidate::trajectoryStateOnDet(), trajectoryStateTransform::transientState(), findQualityFiles::v, LocalTrajectoryParameters::vector(), PV3DBase< T, PVType, FrameType >::x(), x, PV3DBase< T, PVType, FrameType >::y(), detailsBasic3DVector::y, PV3DBase< T, PVType, FrameType >::z(), and z.

                                                                    {        

#ifdef FAMOS_DEBUG
  std::cout << "################################################################" << std::endl;
  std::cout << " TrackCandidateProducer produce init " << std::endl;
#endif

  // Useful typedef's to avoid retyping
  typedef std::pair<reco::TrackRef,edm::Ref<std::vector<Trajectory> > > TrackPair;
  typedef std::map<unsigned,TrackPair> TrackMap;

  // The produced objects
  std::auto_ptr<TrackCandidateCollection> output(new TrackCandidateCollection);    
  std::auto_ptr<reco::TrackCollection> recoTracks(new reco::TrackCollection);    
  std::auto_ptr<TrackingRecHitCollection> recoHits(new TrackingRecHitCollection);
  std::auto_ptr<reco::TrackExtraCollection> recoTrackExtras(new reco::TrackExtraCollection);
  std::auto_ptr<std::vector<Trajectory> > recoTrajectories(new std::vector<Trajectory> );
  std::auto_ptr<TrajTrackAssociationCollection> recoTrajTrackMap( new TrajTrackAssociationCollection() );
  
  // Get the seeds
  // edm::Handle<TrajectorySeedCollection> theSeeds;
  edm::Handle<edm::View<TrajectorySeed> > theSeeds;
  e.getByLabel(seedProducer,theSeeds);

  //Retrieve tracker topology from geometry
  edm::ESHandle<TrackerTopology> tTopoHand;
  es.get<IdealGeometryRecord>().get(tTopoHand);
  const TrackerTopology *tTopo=tTopoHand.product();


  // No seed -> output an empty track collection
  if(theSeeds->size() == 0) {
    e.put(output);
    e.put(recoTracks);
    e.put(recoHits);
    e.put(recoTrackExtras);
    e.put(recoTrajectories);
    e.put(recoTrajTrackMap);
    return;
  }

  // Get the GS RecHits
  //  edm::Handle<SiTrackerGSRecHit2DCollection> theGSRecHits;
  edm::Handle<SiTrackerGSMatchedRecHit2DCollection> theGSRecHits;
  e.getByLabel(hitProducer, theGSRecHits);

  //get other general things
  const std::vector<unsigned> theSimTrackIds = theGSRecHits->ids(); 
  // SimTracks and SimVertices

  edm::Handle<edm::SimVertexContainer> theSimVtx;
  e.getByLabel("famosSimHits",theSimVtx);
  edm::Handle<edm::SimTrackContainer> theSTC;
  e.getByLabel("famosSimHits",theSTC);

  const edm::SimTrackContainer* theSimTracks = &(*theSTC);
  LogDebug("FastTracking")<<"looking at: "<< theSimTrackIds.size()<<" simtracks.";
 

     

  // The input track collection + extra's
  /*
  edm::Handle<reco::TrackCollection> theTrackCollection;
  edm:: Handle<std::vector<Trajectory> > theTrajectoryCollection;
  edm::Handle<TrajTrackAssociationCollection> theAssoMap;  
  bool isTrackCollection = e.getByLabel(trackProducer,theTrackCollection);
  */
  std::vector<edm::Handle<reco::TrackCollection> > theTrackCollections;
  std::vector<edm:: Handle<std::vector<Trajectory> > > theTrajectoryCollections;
  std::vector<edm::Handle<TrajTrackAssociationCollection> > theAssoMaps;
  std::vector<bool> isTrackCollections;
  TrajTrackAssociationCollection::const_iterator anAssociation;  
  TrajTrackAssociationCollection::const_iterator lastAssociation;
  TrackMap theTrackMap;
  unsigned nCollections = trackProducers.size();
  unsigned nRecoHits = 0;

  if ( nCollections ) { 
    theTrackCollections.resize(nCollections);
    theTrajectoryCollections.resize(nCollections);
    theAssoMaps.resize(nCollections);
    isTrackCollections.resize(nCollections);
    for ( unsigned tprod=0; tprod < nCollections; ++tprod ) { 
      isTrackCollections[tprod] = e.getByLabel(trackProducers[tprod],theTrackCollections[tprod]); 

      if ( isTrackCollections[tprod] ) { 
        // The track collection
        reco::TrackCollection::const_iterator aTrack = theTrackCollections[tprod]->begin();
        reco::TrackCollection::const_iterator lastTrack = theTrackCollections[tprod]->end();
        // The numbers of hits
        for ( ; aTrack!=lastTrack; ++aTrack ) nRecoHits+= aTrack->recHitsSize();
        e.getByLabel(trackProducers[tprod],theTrajectoryCollections[tprod]);
        e.getByLabel(trackProducers[tprod],theAssoMaps[tprod]);
        // The association between trajectories and tracks
        anAssociation = theAssoMaps[tprod]->begin();
        lastAssociation = theAssoMaps[tprod]->end(); 
#ifdef FAMOS_DEBUG
        std::cout << "Input Track Producer " << tprod << " : " << trackProducers[tprod] << std::endl;
        std::cout << "List of tracks already reconstructed " << std::endl;
#endif
        // Build the map of correspondance between reco tracks and sim tracks
        for ( ; anAssociation != lastAssociation; ++anAssociation ) { 
          edm::Ref<std::vector<Trajectory> > aTrajectoryRef = anAssociation->key;
          reco::TrackRef aTrackRef = anAssociation->val;
          // Find the simtrack id of the reconstructed track
          int recoTrackId = findId(*aTrackRef);
          if ( recoTrackId < 0 ) continue;
#ifdef FAMOS_DEBUG
          std::cout << recoTrackId << " ";
#endif
          // And store it.
          theTrackMap[recoTrackId] = TrackPair(aTrackRef,aTrajectoryRef);
        }
#ifdef FAMOS_DEBUG
        std::cout << std::endl;
#endif
      }
    }
    // This is to save some time at push_back.
    recoHits->reserve(nRecoHits); 
  }


  // Loop over the seeds
  int currentTrackId = -1;
  /*
  TrajectorySeedCollection::const_iterator aSeed = theSeeds->begin();
  TrajectorySeedCollection::const_iterator lastSeed = theSeeds->end();
  for ( ; aSeed!=lastSeed; ++aSeed ) { 
    // The first hit of the seed  and its simtrack id
  */
  /* */
#ifdef FAMOS_DEBUG
  std::cout << "Input seed Producer : " << seedProducer << std::endl;
  std::cout << "Number of seeds : " << theSeeds->size() << std::endl;
#endif
  unsigned seed_size = theSeeds->size(); 
  for (unsigned seednr = 0; seednr < seed_size; ++seednr){
    
    LogDebug("FastTracking")<<"looking at seed #:"<<seednr;

    // The seed
    const BasicTrajectorySeed* aSeed = &((*theSeeds)[seednr]);

    std::vector<int> simTrackIds;
    std::map<int,TrajectoryStateOnSurface> seedStates;
    std::map<int,TrajectoryStateOnSurface> simtkStates;

    TrackerRecHit theFirstSeedingTrackerRecHit;
    if (theSeeds->at(seednr).nHits()==0){
      //new stuff for no hits on seed

      LogDebug("FastTracking")<<" seed with no hits to be considered.";

      //moved out of the loop
      //edm::ESHandle<MagneticField> field;
      //es.get<IdealMagneticFieldRecord>().get(field);

      PTrajectoryStateOnDet ptod =theSeeds->at(seednr).startingState();
      DetId id(ptod.detId());
      const GeomDet * g = theGeometry->idToDet(id);
      const Surface * surface=&g->surface();
      
      TrajectoryStateOnSurface seedState(trajectoryStateTransform::transientState(ptod,surface,theMagField));
      
      edm::ESHandle<Propagator> propagator;
      es.get<TrackingComponentsRecord>().get("AnyDirectionAnalyticalPropagator",propagator);
      
      //moved out of the loop 
      //      const std::vector<unsigned> theSimTrackIds = theGSRecHits->ids(); 
      //      edm::Handle<edm::SimTrackContainer> theSTC;
      //      e.getByLabel(simTracks_,theSTC);
      //      const edm::SimTrackContainer* theSimTracks = &(*theSTC);

      double minimunEst=1000000;
      LogDebug("FastTracking")<<"looking at: "<< theSimTrackIds.size()<<" simtracks.";
      for ( unsigned tkId=0;  tkId != theSimTrackIds.size(); ++tkId ) {
        
        const SimTrack & simtrack = (*theSimTracks)[theSimTrackIds[tkId]];

        GlobalPoint position(simtrack.trackerSurfacePosition().x(),
                             simtrack.trackerSurfacePosition().y(),
                             simtrack.trackerSurfacePosition().z());
        
        GlobalVector momentum(simtrack.trackerSurfaceMomentum().x(),
                              simtrack.trackerSurfaceMomentum().y(),
                              simtrack.trackerSurfaceMomentum().z());

        if (position.basicVector().dot( momentum.basicVector() ) * seedState.globalPosition().basicVector().dot(seedState.globalMomentum().basicVector()) <0. ){
          LogDebug("FastTracking")<<"not on the same direction.";
          continue;
        }

        //no charge mis-identification ... FIXME
        int charge = (int) simtrack.charge();
        GlobalTrajectoryParameters glb_parameters(position,
                                                  momentum,
                                                  charge,
                                                  theMagField);
        FreeTrajectoryState simtrack_trackerstate(glb_parameters);
        
        TrajectoryStateOnSurface simtrack_comparestate = propagator->propagate(simtrack_trackerstate,*surface);

          
        if (!simtrack_comparestate.isValid()){
          LogDebug("FastTracking")<<" ok this is a state-based seed. simtrack state does not propagate to the seed surface. skipping.";
          continue;}
        
        if (simtrack_comparestate.globalPosition().basicVector().dot(simtrack_comparestate.globalMomentum().basicVector()) * seedState.globalPosition().basicVector().dot(seedState.globalMomentum().basicVector()) <0. ){
          LogDebug("FastTracking")<<"not on the same direction.";
          continue;
        }

        
        AlgebraicVector5 v(seedState.localParameters().vector() - simtrack_comparestate.localParameters().vector());
        AlgebraicSymMatrix55 m(seedState.localError().matrix());
        bool ierr = !m.Invert();
        if ( ierr ){
          edm::LogWarning("FastTracking") <<" Candidate Producer cannot invert the error matrix! - Skipping...";
          continue;
        }
        double est = ROOT::Math::Similarity(v,m);
        LogDebug("FastTracking")<<"comparing two state on the seed surface:\n"
                                          <<"seed: "<<seedState
                                          <<"sim: "<<simtrack_comparestate
                                          <<"\n estimator is: "<<est;

        if (est<minimunEst)       minimunEst=est;
        if (est<estimatorCut_){
          simTrackIds.push_back(theSimTrackIds[tkId]);
          //making a state with exactly the sim track parameters
          //the initial errors are set to unity just for kicks
          //      AlgebraicSymMatrix C(5,1);// C*=50;
          //new attempt!!!!
          AlgebraicSymMatrix55 C = seedState.curvilinearError().matrix();
          C *= 0.0000001;

          seedStates[theSimTrackIds[tkId]] = TrajectoryStateOnSurface(simtrack_comparestate.globalParameters(),
                                                                      CurvilinearTrajectoryError(C),
                                                                      seedState.surface());
          LogDebug("FastTracking")<<"the compatibility estimator is: "<<est<<" for track id: "<<simTrackIds.back();
        }
      }//SimTrack loop
      if (simTrackIds.size()==0) LogDebug("FastTracking")<<"could not find any simtrack within errors, closest was at: "<<minimunEst;
    }//seed has 0 hit.
    else{
      //same old stuff
      // Find the first hit of the Seed
      TrajectorySeed::range theSeedingRecHitRange = aSeed->recHits();
      const SiTrackerGSMatchedRecHit2D * theFirstSeedingRecHit = (const SiTrackerGSMatchedRecHit2D*) (&(*(theSeedingRecHitRange.first)));
      theFirstSeedingTrackerRecHit = TrackerRecHit(theFirstSeedingRecHit,theGeometry,tTopo);
      // The SimTrack id associated to that recHit
      simTrackIds.push_back( theFirstSeedingRecHit->simtrackId() );
    }

    //from then on, only the simtrack IDs are usefull.
    //now loop over all possible trackid for this seed.
    //an actual seed can be shared by two tracks in dense envirronement, and also for hit-less seeds.
    for (unsigned int iToMake=0;iToMake!=simTrackIds.size();++iToMake){
      int simTrackId = simTrackIds[iToMake];
      
      // Don't consider seeds belonging to a track already considered 
      // (Equivalent to seed cleaning)
      if ( seedCleaning && currentTrackId == simTrackId ) continue;
      currentTrackId = simTrackId;
      
      // A vector of TrackerRecHits belonging to the track and the number of crossed layers
      std::vector<TrackerRecHit> theTrackerRecHits;
      unsigned theNumberOfCrossedLayers = 0;
      
      // The track has indeed been reconstructed already -> Save the pertaining info
      TrackMap::const_iterator theTrackIt = theTrackMap.find(simTrackId);
      if ( nCollections && theTrackIt != theTrackMap.end() ) { 
        
        if ( keepFittedTracks ) { 
          LogDebug("FastTracking") << "Track " << simTrackId << " already reconstructed -> copy it";
          // The track and trajectroy references
          reco::TrackRef aTrackRef = theTrackIt->second.first;
          edm::Ref<std::vector<Trajectory> > aTrajectoryRef = theTrackIt->second.second;
          
          // A copy of the track
          reco::Track aRecoTrack(*aTrackRef);
          recoTracks->push_back(aRecoTrack);      
          
          // A copy of the hits
          unsigned nh = aRecoTrack.recHitsSize();
          for ( unsigned ih=0; ih<nh; ++ih ) {
            TrackingRecHit *hit = aRecoTrack.recHit(ih)->clone();
            recoHits->push_back(hit);
          }
          
          // A copy of the trajectories
          recoTrajectories->push_back(*aTrajectoryRef);
          
        }// keepFitterTracks
        else {    
          LogDebug("FastTracking") << "Track " << simTrackId << " already reconstructed -> ignore it";
        }

        // The track was not saved -> create a track candidate.
        
      } //already existing collection of tracks
      else{//no collection of tracks already exists

        LogDebug("FastTracking")<<"Track " << simTrackId << " is considered to return a track candidate" ;

        // Get all the rechits associated to this track
        SiTrackerGSMatchedRecHit2DCollection::range theRecHitRange = theGSRecHits->get(simTrackId);
        SiTrackerGSMatchedRecHit2DCollection::const_iterator theRecHitRangeIteratorBegin = theRecHitRange.first;
        SiTrackerGSMatchedRecHit2DCollection::const_iterator theRecHitRangeIteratorEnd   = theRecHitRange.second;
        SiTrackerGSMatchedRecHit2DCollection::const_iterator iterRecHit;

        LogDebug("FastTracking")<<"counting: "<<theRecHitRangeIteratorEnd-theRecHitRangeIteratorBegin<<" hits to be considered.";

        bool firstRecHit = true;
        TrackerRecHit theCurrentRecHit, thePreviousRecHit;
        TrackerRecHit theFirstHitComp, theSecondHitComp;
        
        for ( iterRecHit = theRecHitRangeIteratorBegin; 
              iterRecHit != theRecHitRangeIteratorEnd; 
              ++iterRecHit) {
          
          // Check the number of crossed layers
          if ( theNumberOfCrossedLayers >= maxNumberOfCrossedLayers ) continue;
          
          // Get current and previous rechits
          if(!firstRecHit) thePreviousRecHit = theCurrentRecHit;
          theCurrentRecHit = TrackerRecHit(&(*iterRecHit),theGeometry,tTopo);
          
          //>>>>>>>>>BACKBUILDING CHANGE: DO NOT STAT FROM THE FIRST HIT OF THE SEED

          // NOTE: checking the direction --> specific for OIHit only
          //      if( aSeed->direction()!=oppositeToMomentum ) { 
          //  // Check that the first rechit is indeed the first seeding hit
          //  if ( firstRecHit && theCurrentRecHit != theFirstSeedingTrackerRecHit && theSeeds->at(seednr).nHits()!=0 ) continue;
          // }
          //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

          // Count the number of crossed layers
          if ( !theCurrentRecHit.isOnTheSameLayer(thePreviousRecHit) ) 
            ++theNumberOfCrossedLayers;
          
          // Add all rechits (Grouped Trajectory Builder) from this hit onwards
          // Always add the first seeding rechit anyway
          if ( !rejectOverlaps || firstRecHit ) {  
            // Split matched hits (if requested / possible )
            if ( splitHits && theCurrentRecHit.matchedHit()->isMatched() ) addSplitHits(theCurrentRecHit,theTrackerRecHits);
            else theTrackerRecHits.push_back(theCurrentRecHit);       // No splitting   
            firstRecHit = false;
            
            // And now treat the following RecHits if hits in the same layer 
            // have to be rejected - The split option is not 
          } else { 
            
            // Not the same layer : Add the current hit
            if ( theCurrentRecHit.subDetId()    != thePreviousRecHit.subDetId() || 
                 theCurrentRecHit.layerNumber() != thePreviousRecHit.layerNumber() ) {
              
              // Split matched hits (if requested / possible )
              if ( splitHits && theCurrentRecHit.matchedHit()->isMatched() ) addSplitHits(theCurrentRecHit,theTrackerRecHits);
              else              theTrackerRecHits.push_back(theCurrentRecHit);          // No splitting               
              
              // Same layer : keep the current hit if better, and drop the other - otherwise do nothing  
            } else if ( theCurrentRecHit.localError() < thePreviousRecHit.localError() ) { 
              
              // Split matched hits (if requested / possible )
              if( splitHits && theCurrentRecHit.matchedHit()->isMatched() ){
                // Remove the previous hit(s)
                theTrackerRecHits.pop_back();
                if ( thePreviousRecHit.matchedHit()->isMatched() ) theTrackerRecHits.pop_back();
                
                // Replace by the new hits
                addSplitHits(theCurrentRecHit,theTrackerRecHits);
              }
              // No splitting   
              else {
                theTrackerRecHits.back() = theCurrentRecHit; // Replace the previous hit by the current hit
              }
              
            } else {
              
              //keep the older rechit as a reference if the error of the new one is worse
              theCurrentRecHit = thePreviousRecHit;
            }  
          }
        }
        
        // End of loop over the track rechits
      }//no collection of track already existed. adding the hits by hand.
    
      LogDebug("FastTracking")<<" number of hits: " << theTrackerRecHits.size()<<" after counting overlaps and splitting.";

      // 1) Create the OwnWector of TrackingRecHits
      edm::OwnVector<TrackingRecHit> recHits;
      unsigned nTrackerHits = theTrackerRecHits.size();
      recHits.reserve(nTrackerHits); // To save some time at push_back

      if (aSeed->direction()==oppositeToMomentum){
        LogDebug("FastTracking")<<"reversing the order of the hits";
        std::reverse(theTrackerRecHits.begin(),theTrackerRecHits.end());
      }

      for ( unsigned ih=0; ih<nTrackerHits; ++ih ) {
        TrackingRecHit* aTrackingRecHit = theTrackerRecHits[ih].hit()->clone();
        recHits.push_back(aTrackingRecHit);
        
        const DetId& detId = theTrackerRecHits[ih].hit()->geographicalId();
        LogDebug("FastTracking")
          << "Added RecHit from detid " << detId.rawId() 
          << " subdet = " << theTrackerRecHits[ih].subDetId() 
          << " layer = " << theTrackerRecHits[ih].layerNumber()
          << " ring = " << theTrackerRecHits[ih].ringNumber()
          << " error = " << theTrackerRecHits[ih].localError()
          << std::endl
          
          << "Track/z/r : "
          << simTrackId << " " 
          << theTrackerRecHits[ih].globalPosition().z() << " " 
          << theTrackerRecHits[ih].globalPosition().perp() << std::endl;
        if ( theTrackerRecHits[ih].matchedHit() && theTrackerRecHits[ih].matchedHit()->isMatched() ) 
          LogTrace("FastTracking") << "Matched : " << theTrackerRecHits[ih].matchedHit()->isMatched() 
                                             << "Rphi Hit = " <<  theTrackerRecHits[ih].matchedHit()->monoHit()->simhitId()              
                                             << "Stereo Hit = " <<  theTrackerRecHits[ih].matchedHit()->stereoHit()->simhitId()
                                             <<std::endl;
      }//loop over the rechits

    // Check the number of crossed layers
    if ( theNumberOfCrossedLayers < minNumberOfCrossedLayers ) {
      LogDebug("FastTracking")<<"not enough layer crossed ("<<theNumberOfCrossedLayers<<")";
      continue;
    }


    //>>>>>>>>>BACKBUILDING CHANGE: REPLACE THE STARTING STATE

    // Create a track Candidate (now with the reference to the seed!) .
    //PTrajectoryStateOnDet PTSOD = aSeed->startingState();
    PTrajectoryStateOnDet PTSOD;

    if (aSeed->nHits()==0){
      //stabilize the fit with the true simtrack state
      //in case of zero hits
      
      PTSOD = trajectoryStateTransform::persistentState(seedStates[simTrackId],aSeed->startingState().detId());
 
    } else {
      //create the initial state from the SimTrack
      int vertexIndex = (*theSimTracks)[currentTrackId].vertIndex();
      //   a) origin vertex
      GlobalPoint  position((*theSimVtx)[vertexIndex].position().x(),
                            (*theSimVtx)[vertexIndex].position().y(),
                            (*theSimVtx)[vertexIndex].position().z());
      
      //   b) initial momentum
      GlobalVector momentum( (*theSimTracks)[currentTrackId].momentum().x() , 
                             (*theSimTracks)[currentTrackId].momentum().y() , 
                             (*theSimTracks)[currentTrackId].momentum().z() );
      //   c) electric charge
      float        charge   = (*theSimTracks)[simTrackId].charge();
      //  -> inital parameters
      GlobalTrajectoryParameters initialParams(position,momentum,(int)charge,theMagField);
 //  -> large initial errors
      AlgebraicSymMatrix55 errorMatrix= AlgebraicMatrixID();    
      CurvilinearTrajectoryError initialError(errorMatrix);
      // -> initial state
      FreeTrajectoryState initialFTS(initialParams, initialError);      
#ifdef FAMOS_DEBUG
      std::cout << "TrajectorySeedProducer: FTS momentum " << initialFTS.momentum() << std::endl;
#endif
      const GeomDet* initialLayer = theGeometry->idToDet(recHits.front().geographicalId());
      //this is wrong because the FTS is defined at vertex, and it need to be properly propagated to the first rechit
      //      const TrajectoryStateOnSurface initialTSOS(initialFTS, initialLayer->surface());      
       const TrajectoryStateOnSurface initialTSOS = thePropagator->propagate(initialFTS,initialLayer->surface()) ;
       if (!initialTSOS.isValid()) continue; 
       

       PTSOD = trajectoryStateTransform::persistentState(initialTSOS,recHits.front().geographicalId().rawId()); 
    }
    
    TrackCandidate  
      newTrackCandidate(recHits, 
                        *aSeed, 
                        PTSOD, 
                        edm::RefToBase<TrajectorySeed>(theSeeds,seednr));

    LogDebug("FastTracking")<< "\tSeed Information " << std::endl
                                      << "\tSeed Direction = " << aSeed->direction() << std::endl
                                      << "\tSeed StartingDet = " << aSeed->startingState().detId() << std::endl
                                      << "\tTrajectory Parameters "           << std::endl
                                      << "\t\t detId  = "             << newTrackCandidate.trajectoryStateOnDet().detId()             << std::endl
                                      << "\t\t loc.px = "
                                      << newTrackCandidate.trajectoryStateOnDet().parameters().momentum().x()    
                                      << std::endl
                                      << "\t\t loc.py = "
                                      << newTrackCandidate.trajectoryStateOnDet().parameters().momentum().y()    
                                      << std::endl
                                      << "\t\t loc.pz = "
                                      << newTrackCandidate.trajectoryStateOnDet().parameters().momentum().z()    
                                      << std::endl
                                      << "\t\t error  = ";

    bool newTrackCandidateIsDuplicate = isDuplicateCandidate(*output,newTrackCandidate);
    if (!newTrackCandidateIsDuplicate) output->push_back(newTrackCandidate);
    LogDebug("FastTracking")<<"filling a track candidate into the collection, now having: "<<output->size();
    
    }//loop over possible simtrack associated.
  }//loop over all possible seeds.
  
  // Save the track candidates in the event
  LogDebug("FastTracking") << "Saving " 
                                     << output->size() << " track candidates and " 
                                     << recoTracks->size() << " reco::Tracks ";
  // Save the track candidates
  e.put(output);



  // Save the tracking recHits

  edm::OrphanHandle <TrackingRecHitCollection> theRecoHits = e.put(recoHits );

  // Create the track extras and add the references to the rechits
  unsigned hits=0;
  unsigned nTracks = recoTracks->size();
  recoTrackExtras->reserve(nTracks); // To save some time at push_back
  for ( unsigned index = 0; index < nTracks; ++index ) { 
    //reco::TrackExtra aTrackExtra;
    reco::Track& aTrack = recoTracks->at(index);
    reco::TrackExtra aTrackExtra(aTrack.outerPosition(),
                                 aTrack.outerMomentum(),
                                 aTrack.outerOk(),
                                 aTrack.innerPosition(),
                                 aTrack.innerMomentum(),
                                 aTrack.innerOk(),
                                 aTrack.outerStateCovariance(),
                                 aTrack.outerDetId(),
                                 aTrack.innerStateCovariance(),
                                 aTrack.innerDetId(),
                                 aTrack.seedDirection(),
                                 aTrack.seedRef());

    unsigned nHits = aTrack.recHitsSize();
    for ( unsigned int ih=0; ih<nHits; ++ih) {
      aTrackExtra.add(TrackingRecHitRef(theRecoHits,hits++));
    }
    recoTrackExtras->push_back(aTrackExtra);
  }
  

  // Save the track extras
  edm::OrphanHandle<reco::TrackExtraCollection> theRecoTrackExtras = e.put(recoTrackExtras);

  // Add the reference to the track extra in the tracks
  for ( unsigned index = 0; index<nTracks; ++index ) { 
    const reco::TrackExtraRef theTrackExtraRef(theRecoTrackExtras,index);
    (recoTracks->at(index)).setExtra(theTrackExtraRef);
  }

  // Save the tracks
  edm::OrphanHandle<reco::TrackCollection> theRecoTracks = e.put(recoTracks);

  // Save the trajectories
  edm::OrphanHandle<std::vector<Trajectory> > theRecoTrajectories = e.put(recoTrajectories);
  
  // Create and set the trajectory/track association map 
  for ( unsigned index = 0; index<nTracks; ++index ) { 
    edm::Ref<std::vector<Trajectory> > trajRef( theRecoTrajectories, index );
    edm::Ref<reco::TrackCollection>    tkRef( theRecoTracks, index );
    recoTrajTrackMap->insert(trajRef,tkRef);
  }


  // Save the association map.
  e.put(recoTrajTrackMap);

}
bool TrackCandidateProducer::sameLocalParameters ( const TrackingRecHit aH,
const TrackingRecHit bH 
) const [private]

Definition at line 784 of file TrackCandidateProducer.cc.

References TrackingRecHit::localPosition(), TrackingRecHit::localPositionError(), LocalError::xx(), LocalError::xy(), and LocalError::yy().

Referenced by isDuplicateCandidate().

                                                                                                         {
  bool localPointSame = aH->localPosition() == bH->localPosition();
  bool localErrorSame = aH->localPositionError().xx() == bH->localPositionError().xx();
  localErrorSame &= aH->localPositionError().yy() == bH->localPositionError().yy();
  localErrorSame &= aH->localPositionError().xy() == bH->localPositionError().xy();
  return localPointSame && localErrorSame;
}

Member Data Documentation

Definition at line 69 of file TrackCandidateProducer.h.

Referenced by produce(), and TrackCandidateProducer().

Definition at line 56 of file TrackCandidateProducer.h.

Referenced by produce(), and TrackCandidateProducer().

Definition at line 66 of file TrackCandidateProducer.h.

Referenced by produce(), and TrackCandidateProducer().

Definition at line 61 of file TrackCandidateProducer.h.

Referenced by produce(), and TrackCandidateProducer().

Definition at line 60 of file TrackCandidateProducer.h.

Referenced by produce(), and TrackCandidateProducer().

Definition at line 63 of file TrackCandidateProducer.h.

Referenced by produce(), and TrackCandidateProducer().

Definition at line 65 of file TrackCandidateProducer.h.

Referenced by produce(), and TrackCandidateProducer().

Definition at line 55 of file TrackCandidateProducer.h.

Referenced by produce(), and TrackCandidateProducer().

Definition at line 68 of file TrackCandidateProducer.h.

Referenced by TrackCandidateProducer().

Definition at line 64 of file TrackCandidateProducer.h.

Referenced by produce(), and TrackCandidateProducer().

Definition at line 50 of file TrackCandidateProducer.h.

Referenced by beginRun(), and produce().

Definition at line 51 of file TrackCandidateProducer.h.

Referenced by beginRun(), and produce().

Definition at line 52 of file TrackCandidateProducer.h.

Referenced by beginRun(), produce(), and ~TrackCandidateProducer().

Definition at line 58 of file TrackCandidateProducer.h.

Referenced by produce(), and TrackCandidateProducer().