#include <TrackCandidateProducer.h>

Inheritance diagram for TrackCandidateProducer:

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 ()

Public Member Functions inherited from edm::EDProducer
	EDProducer ()

ModuleDescription const &	moduleDescription () const

virtual	~EDProducer ()

Public Member Functions inherited from edm::ProducerBase
	ProducerBase ()

void	registerProducts (ProducerBase , ProductRegistry , ModuleDescription const &)

std::function< void(BranchDescription const &)>	registrationCallback () const
	used by the fwk to register list of products More...

virtual	~ProducerBase ()

Public Member Functions inherited from edm::EDConsumerBase
	EDConsumerBase ()

ProductHolderIndexAndSkipBit	indexFrom (EDGetToken, BranchType, TypeID const &) const

void	itemsMayGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const

void	itemsToGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const

std::vector < ProductHolderIndexAndSkipBit > const &	itemsToGetFromEvent () const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesDependentUpon (const std::string &iProcessName, std::vector< const char * > &oModuleLabels) const

bool	registeredToConsume (ProductHolderIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

void	updateLookup (BranchType iBranchType, ProductHolderIndexHelper const &)

virtual	~EDConsumerBase ()

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
std::vector< edm::EDGetTokenT < TrajTrackAssociationCollection > >	assoMapTokens

double	estimatorCut_

edm::InputTag	hitProducer

bool	keepFittedTracks

unsigned int	maxNumberOfCrossedLayers

unsigned int	minNumberOfCrossedLayers

edm::EDGetTokenT < SiTrackerGSMatchedRecHit2DCollection >	recHitToken

bool	rejectOverlaps

bool	seedCleaning

edm::InputTag	seedProducer

edm::EDGetTokenT< edm::View < TrajectorySeed > >	seedToken

edm::InputTag	simTracks_

edm::EDGetTokenT < edm::SimTrackContainer >	simTrackToken

edm::EDGetTokenT < edm::SimVertexContainer >	simVertexToken

bool	splitHits

const TrackerGeometry *	theGeometry

const MagneticField *	theMagField

PropagatorWithMaterial *	thePropagator

std::vector< edm::InputTag >	trackProducers

std::vector< edm::EDGetTokenT < reco::TrackCollection > >	trackTokens

std::vector< edm::EDGetTokenT < std::vector< Trajectory > > >	trajectoryTokens

Additional Inherited Members
Public Types inherited from edm::EDProducer
typedef EDProducer	ModuleType

Public Types inherited from edm::ProducerBase
typedef ProductRegistryHelper::TypeLabelList	TypeLabelList

Static Public Member Functions inherited from edm::EDProducer
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

Protected Member Functions inherited from edm::ProducerBase
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)

Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	consumes (edm::InputTag const &tag)

EDGetToken	consumes (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	consumes (TypeToGet const &id, edm::InputTag const &tag)

ConsumesCollector	consumesCollector ()
	Use a ConsumesCollector to gather consumes information from helper functions. More...

template<typename ProductType , BranchType B = InEvent>
void	consumesMany ()

void	consumesMany (const TypeToGet &id)

template<BranchType B>
void	consumesMany (const TypeToGet &id)

template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	mayConsume (edm::InputTag const &tag)

EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

Detailed Description

Definition at line 34 of file TrackCandidateProducer.h.

Constructor & Destructor Documentation

TrackCandidateProducer::TrackCandidateProducer ( const edm::ParameterSet & conf )

explicit

Definition at line 48 of file TrackCandidateProducer.cc.

References assoMapTokens, edm::EDConsumerBase::consumes(), gather_cfg::cout, estimatorCut_, edm::ParameterSet::getParameter(), hitProducer, keepFittedTracks, maxNumberOfCrossedLayers, minNumberOfCrossedLayers, recHitToken, rejectOverlaps, seedCleaning, seedProducer, seedToken, simTracks_, simTrackToken, simVertexToken, splitHits, trackProducers, trackTokens, and trajectoryTokens.

                                                                          :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");
 
   // consumes
   seedToken = consumes<edm::View<TrajectorySeed> >(seedProducer);
   recHitToken = consumes<SiTrackerGSMatchedRecHit2DCollection>(hitProducer);
   edm::InputTag _label("famosSimHits");
   simVertexToken = consumes<edm::SimVertexContainer>(_label);
   simTrackToken = consumes<edm::SimTrackContainer>(_label);
   for(unsigned tprod=0; tprod < trackProducers.size(); ++tprod){
     trackTokens.push_back(consumes<reco::TrackCollection>(trackProducers[tprod]));
     trajectoryTokens.push_back(consumes<std::vector<Trajectory> >(trackProducers[tprod]));
     assoMapTokens.push_back(consumes<TrajTrackAssociationCollection>(trackProducers[tprod]));
   }
 }

TrackCandidateProducer::~TrackCandidateProducer ( )

virtual

Definition at line 112 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 737 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
	)

overridevirtual

Reimplemented from edm::EDProducer.

Definition at line 125 of file TrackCandidateProducer.cc.

References alongMomentum, geometry, edm::EventSetup::get(), theGeometry, theMagField, and thePropagator.

                                                                         {
 
   //services
 
   edm::ESHandle<MagneticField>          magField;
   edm::ESHandle<TrackerGeometry>        geometry;
 
   es.get<IdealMagneticFieldRecord>().get(magField);
   es.get<TrackerDigiGeometryRecord>().get(geometry);
 
   theMagField = &(*magField);
   theGeometry = &(*geometry);
 
   thePropagator = new PropagatorWithMaterial(alongMomentum,0.105,&(*theMagField)); 
 }

int TrackCandidateProducer::findId ( const reco::Track & aTrack ) const

private

Definition at line 722 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 758 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
	)

overridevirtual

Implements edm::EDProducer.

Definition at line 143 of file TrackCandidateProducer.cc.

References addSplitHits(), assoMapTokens, PV3DBase< T, PVType, FrameType >::basicVector(), funct::C, DeDxDiscriminatorTools::charge(), CoreSimTrack::charge(), TrackingRecHit::clone(), edm::AssociationMap< edm::OneToOne< std::vector< Trajectory >, reco::TrackCollection, unsigned short > >::const_iterator, gather_cfg::cout, TrajectoryStateOnSurface::curvilinearError(), PTrajectoryStateOnDet::detId(), TrajectorySeed::direction(), Basic3DVector< T >::dot(), estimatorCut_, findId(), edm::OwnVector< T, P >::front(), g, TrackingRecHit::geographicalId(), edm::EventSetup::get(), edm::Event::getByToken(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalParameters(), TrajectoryStateOnSurface::globalPosition(), TrackingRecHit::hit(), TrackerGeometry::idToDet(), cmsHarvester::index, reco::Track::innerDetId(), reco::Track::innerMomentum(), reco::Track::innerOk(), reco::Track::innerPosition(), reco::Track::innerStateCovariance(), isDuplicateCandidate(), trackerHitRTTI::isMatched(), SiTrackerGSMatchedRecHit2D::isMatched(), TrackerRecHit::isOnTheSameLayer(), TrajectoryStateOnSurface::isValid(), keepFittedTracks, edm::Ref< C, T, F >::key(), TrackerRecHit::layerNumber(), TrajectoryStateOnSurface::localError(), TrackerRecHit::localError(), TrajectoryStateOnSurface::localParameters(), LogDebug, LogTrace, m, TrackerRecHit::matchedHit(), LocalTrajectoryError::matrix(), CurvilinearTrajectoryError::matrix(), maxNumberOfCrossedLayers, minNumberOfCrossedLayers, FreeTrajectoryState::momentum(), LocalTrajectoryParameters::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< class >::product(), Propagator::propagate(), edm::OwnVector< T, P >::push_back(), edm::Event::put(), DetId::rawId(), reco::Track::recHit(), TrajectorySeed::recHits(), reco::Track::recHitsSize(), recHitToken, rejectOverlaps, edm::OwnVector< T, P >::reserve(), seedCleaning, reco::Track::seedDirection(), seedProducer, reco::Track::seedRef(), seedToken, SiTrackerGSMatchedRecHit2D::simtrackId(), simTrackToken, simVertexToken, edm::OwnVector< T, P >::size(), splitHits, TrajectorySeed::startingState(), TrackerRecHit::subDetId(), GeomDet::surface(), TrajectoryStateOnSurface::surface(), theGeometry, theMagField, thePropagator, SimTrack::trackerSurfaceMomentum(), SimTrack::trackerSurfacePosition(), trackProducers, trackTokens, TrackCandidate::trajectoryStateOnDet(), trajectoryTokens, trajectoryStateTransform::transientState(), findQualityFiles::v, LocalTrajectoryParameters::vector(), x, PV3DBase< T, PVType, FrameType >::x(), detailsBasic3DVector::y, PV3DBase< T, PVType, FrameType >::y(), detailsBasic3DVector::z, and PV3DBase< T, PVType, FrameType >::z().

Referenced by JSONExport.JsonExport::export(), HTMLExport.HTMLExport::export(), and HTMLExport.HTMLExportStatic::export().

                                                                     {        
 
 #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.getByToken(seedToken,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.getByToken(recHitToken, theGSRecHits);
 
   //get other general things
   const std::vector<unsigned> theSimTrackIds = theGSRecHits->ids(); 
   // SimTracks and SimVertices
 
   edm::Handle<edm::SimVertexContainer> theSimVtx;
   e.getByToken(simVertexToken,theSimVtx);
   edm::Handle<edm::SimTrackContainer> theSTC;
   e.getByToken(simTrackToken,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.getByToken(trackTokens[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.getByToken(trajectoryTokens[tprod],theTrajectoryCollections[tprod]);
     e.getByToken(assoMapTokens[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 790 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;
 }