CMSSW Reference Manual

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

Definition at line 73 of file RoadSearchTrackCandidateMakerAlgorithm.cc.

References Chi2MeasurementEstimatorESProducer_cfi::Chi2MeasurementEstimator, conf_, CosmicReco_, cosmicSeedPt_, CosmicTrackMerging_, debug_, debugCosmics_, edm::ParameterSet::getParameter(), initialVertexErrorXY_, int, maxPropagationDistance, measurementTrackerName_, MinChunkLength_, nFoundMin_, splitMatchedHits_, theChi2Cut, theEstimator, theNumHitCut, theTrajectoryCleaner, theTransformer, and theUpdator.

00073                                                                                                           : conf_(conf) { 
00074   
00075   theNumHitCut = (unsigned int)conf_.getParameter<int>("NumHitCut");
00076   theChi2Cut   = conf_.getParameter<double>("HitChi2Cut");
00077   
00078   theEstimator = new Chi2MeasurementEstimator(theChi2Cut);
00079   theUpdator = new KFUpdator();
00080   theTransformer = new TrajectoryStateTransform;
00081   theTrajectoryCleaner = new TrajectoryCleanerBySharedHits;
00082   
00083   CosmicReco_  = conf_.getParameter<bool>("StraightLineNoBeamSpotCloud");
00084   CosmicTrackMerging_ = conf_.getParameter<bool>("CosmicTrackMerging");
00085   MinChunkLength_ = conf_.getParameter<int>("MinimumChunkLength");
00086   nFoundMin_      = conf_.getParameter<int>("nFoundMin");
00087   
00088   initialVertexErrorXY_  = conf_.getParameter<double>("InitialVertexErrorXY");
00089   splitMatchedHits_  = conf_.getParameter<bool>("SplitMatchedHits");
00090   cosmicSeedPt_  = conf_.getParameter<double>("CosmicSeedPt");
00091 
00092   measurementTrackerName_ = conf_.getParameter<std::string>("MeasurementTrackerName");
00093   
00094   debug_ = false;
00095   debugCosmics_ = false;
00096 
00097   maxPropagationDistance = 1000.0; // 10m
00098 }

RoadSearchTrackCandidateMakerAlgorithm::~RoadSearchTrackCandidateMakerAlgorithm ( )

Definition at line 100 of file RoadSearchTrackCandidateMakerAlgorithm.cc.

References theEstimator, theTrajectoryCleaner, theTransformer, and theUpdator.

00100                                                                                 {
00101   delete theEstimator;
00102   delete theUpdator;
00103   delete theTransformer;
00104   delete theTrajectoryCleaner;
00105   // delete theMeasurementTracker;
00106 
00107 }

bool RoadSearchTrackCandidateMakerAlgorithm::chooseStartingLayers	(	std::vector< std::pair< const DetLayer *, RoadSearchCloud::RecHitVector > > &	RecHitsByLayer,
		std::vector< std::pair< const DetLayer *, RoadSearchCloud::RecHitVector > >::iterator	ilyr0,
		const std::multimap< int, const DetLayer * > &	layer_map,
		std::set< const DetLayer * > &	good_layers,
		std::vector< const DetLayer * > &	middle_layers,
		RoadSearchCloud::RecHitVector &	recHits_middle
	)

Definition at line 1190 of file RoadSearchTrackCandidateMakerAlgorithm.cc.

References CosmicReco_, lstereo, and nFoundMin_.

Referenced by run().

01197 {
01198       const uint max_middle_layers = 2;
01199 
01200       // consider the best nFoundMin layers + other layers with only one hit      
01201       // This has implications, based on the way we locate the hits.  
01202       // For now, use only the low occupancy layers in the first pass
01203       //const int nfound_min = min_chunk_length-1;
01204       //const int nfound_min = 4;
01205       std::multimap<int, const DetLayer*>::const_iterator ilm = layer_map.begin();
01206       int ngoodlayers = 0;
01207       while (ilm != layer_map.end()) {
01208         if (ngoodlayers >= nFoundMin_ && ilm->first > 1) break;
01209         //if (ilm->first > 1) break;
01210         good_layers.insert(ilm->second);
01211         ++ngoodlayers;
01212         ++ilm;
01213       }
01214 
01215       // choose intermediate layers
01216       for (std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >::iterator ilayer = ilyr0+1;
01217            ilayer != recHitsByLayer.end(); ++ilayer) {
01218         // only use useful layers
01219         if (good_layers.find(ilayer->first) == good_layers.end()) continue;
01220         // only use stereo layers
01221         if (!CosmicReco_ && !lstereo[ilayer-recHitsByLayer.begin()]) continue;
01222         middle_layers.push_back(ilayer->first);
01223         if (middle_layers.size() >= max_middle_layers) break;
01224       }
01225       
01226       for (std::vector<const DetLayer*>::iterator ml = middle_layers.begin();
01227            ml!=middle_layers.end();++ml){
01228         uint middle_layers_found = 0;
01229         for (std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >::iterator ilyr = recHitsByLayer.begin();
01230              ilyr != recHitsByLayer.end(); ++ilyr) {
01231           if (ilyr->first == *ml){
01232             for (RoadSearchCloud::RecHitVector::const_iterator ih = ilyr->second.begin();
01233                  ih != ilyr->second.end(); ++ih) {
01234               recHits_middle.push_back(*ih);
01235             }
01236             ++middle_layers_found;
01237           }
01238           if (middle_layers_found == middle_layers.size()) continue;
01239         }
01240 
01241       }
01242       return (recHits_middle.size()>0);
01243 }

Trajectory RoadSearchTrackCandidateMakerAlgorithm::createSeedTrajectory	(	FreeTrajectoryState &	fts,
		const TrackingRecHit *	InnerHit,
		const DetLayer *	innerHitLayer
	)

Definition at line 1406 of file RoadSearchTrackCandidateMakerAlgorithm.cc.

References alongMomentum, TransientTrackingRecHitBuilder::build(), GenMuonPlsPt100GeV_cfg::cout, debug_, TrajectorySeed::direction(), lat::endl(), MeasurementEstimator::estimate(), TrackerGeometry::idToDet(), TrajectoryStateOnSurface::isValid(), TrajectoryStateOnSurface::localDirection(), SiStripRecHitMatcher::match(), oppositeToMomentum, TrajectoryStateTransform::persistentState(), PropagatorWithMaterial::propagate(), Propagator::propagationDirection(), Trajectory::push(), edm::OwnVector< T, P >::push_back(), GeomDet::surface(), theEstimator, theHitMatcher, thePropagator, theTransformer, theUpdator, trackerGeom, ttrhBuilder, and TrajectoryStateUpdator::update().

Referenced by run().

01410 {
01411   Trajectory theSeedTrajectory;
01412 
01413   // Need to put the first hit on the trajectory
01414   const GeomDet* innerDet = trackerGeom->idToDet((theInnerHit)->geographicalId());
01415   const TrajectoryStateOnSurface innerState = 
01416     thePropagator->propagate(fts,innerDet->surface());
01417   if ( !innerState.isValid()) {
01418     if (debug_) std::cout<<"*******DISASTER ********* seed doesn't make it to first hit!!!!!" << std::endl;
01419     return theSeedTrajectory;
01420   }
01421   TransientTrackingRecHit::RecHitPointer intrhit = ttrhBuilder->build(theInnerHit);
01422   // if this first hit is a matched hit, it should be updated for the trajectory
01423   const SiStripMatchedRecHit2D *origHit = dynamic_cast<const SiStripMatchedRecHit2D *>(theInnerHit);
01424   if (origHit !=0){
01425     const GluedGeomDet *gdet = dynamic_cast<const GluedGeomDet*>(innerDet);
01426     const SiStripMatchedRecHit2D *corrHit = theHitMatcher->match(origHit,gdet,innerState.localDirection());
01427     if (corrHit!=0){
01428       intrhit = ttrhBuilder->build(&(*corrHit));
01429       delete corrHit;
01430     }
01431   }
01432   
01433   MeasurementEstimator::HitReturnType est = theEstimator->estimate(innerState, *intrhit);
01434   if (!est.first) return theSeedTrajectory;         
01435   TrajectoryStateOnSurface innerUpdated= theUpdator->update( innerState,*intrhit);                         
01436   TrajectoryMeasurement tm = TrajectoryMeasurement(innerState, innerUpdated, &(*intrhit),est.second,theInnerHitLayer);
01437   
01438   PTrajectoryStateOnDet* pFirstStateTwo = theTransformer->persistentState(innerUpdated,
01439                                                                           intrhit->geographicalId().rawId());
01440   edm::OwnVector<TrackingRecHit> newHitsTwo;
01441   newHitsTwo.push_back(intrhit->hit()->clone());
01442   
01443   TrajectorySeed tmpseedTwo = TrajectorySeed(*pFirstStateTwo, 
01444                                              newHitsTwo,
01445                                              alongMomentum);
01446   if (thePropagator->propagationDirection()==oppositeToMomentum) {
01447     tmpseedTwo = TrajectorySeed(*pFirstStateTwo, 
01448                                 newHitsTwo,
01449                                 oppositeToMomentum);
01450   }
01451 
01452   delete pFirstStateTwo;
01453   
01454   //Trajectory seedTraj(tmpseedTwo, alongMomentum);
01455   theSeedTrajectory = Trajectory(tmpseedTwo, tmpseedTwo.direction());
01456   
01457   theSeedTrajectory.push(tm,est.second);
01458   
01459   return theSeedTrajectory;
01460 }

std::vector< Trajectory > RoadSearchTrackCandidateMakerAlgorithm::extrapolateTrajectory	(	const Trajectory &	inputTrajectory,
		RoadSearchCloud::RecHitVector &	theLayerHits,
		const DetLayer *	innerHitLayer,
		const TrackingRecHit *	outerHit,
		const DetLayer *	outerHitLayer
	)

Definition at line 1464 of file RoadSearchTrackCandidateMakerAlgorithm.cc.

References TransientTrackingRecHitBuilder::build(), GenMuonPlsPt100GeV_cfg::cout, debug_, GeometricSearchTracker::detLayer(), lat::endl(), MeasurementEstimator::estimate(), MeasurementTracker::geometricSearchTracker(), TrajectoryStateOnSurface::globalPosition(), TrackerGeometry::idToDet(), TrajectoryStateOnSurface::isValid(), TrajectoryStateOnSurface::localDirection(), TrajectoryStateOnSurface::localError(), TrackingRecHit::localPosition(), TrajectoryStateOnSurface::localPosition(), LogDebug, PV3DBase< T, PVType, FrameType >::mag(), SiStripRecHitMatcher::match(), maxPropagationDistance, Trajectory::measurements(), p1, p2, GloballyPositioned< T >::position(), LocalTrajectoryError::positionError(), PropagatorWithMaterial::propagate(), Trajectory::push(), funct::sqrt(), GeomDet::surface(), theEstimator, theHitMatcher, theMeasurementTracker, thePropagator, theUpdator, trackerGeom, ttrhBuilder, TrajectoryStateUpdator::update(), PV3DBase< T, PVType, FrameType >::x(), LocalError::xx(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by run().

01469 {
01470   std::vector<Trajectory> newTrajectories;
01471 
01472               for(RoadSearchCloud::RecHitVector::const_iterator ihit = theLayerHits.begin();
01473                   ihit != theLayerHits.end(); ihit++) {
01474                 Trajectory traj = theTrajectory;
01475                 const DetLayer* thisLayer =
01476                   theMeasurementTracker->geometricSearchTracker()->detLayer((*ihit)->geographicalId());
01477                 if (thisLayer == innerHitLayer){
01478                   // Right now we are assuming that ONLY single hit layers are in this initial collection
01479                   //if (thisLayer == innerHitLayer && !(ihit->recHit() == innerHit->recHit())){
01480                   //  if (debug_) std::cout<<"On inner hit layer, but have wrong hit?!?!?" << std::endl;
01481                   continue;
01482                 }
01483                 if (thisLayer == outerHitLayer){
01484                   LocalPoint p1 = (*ihit)->localPosition();
01485                   LocalPoint p2 = outerHit->localPosition();
01486                   if (p1.x()!=p2.x() || p1.y()!=p2.y()) continue;
01487                 }
01488                 // extrapolate
01489                 
01490                 TransientTrackingRecHit::RecHitPointer rhit = ttrhBuilder->build(*ihit);
01491                 
01492                 if (debug_){
01493                   if (rhit->isValid()) {
01494                     LogDebug("RoadSearch") << "RecHit " << ihit-theLayerHits.begin()
01495                                            << ", det " << rhit->det() << ", r/phi/z = "
01496                                            << rhit->globalPosition().perp() << " "
01497                                            << rhit->globalPosition().phi() << " "
01498                                            << rhit->globalPosition().z();
01499                   } else {
01500                     LogDebug("RoadSearch") << "RecHit " << ihit-theLayerHits.begin()
01501                                            << " (invalid)";
01502                   }
01503                 }
01504                 
01505                 const GeomDet* det = trackerGeom->idToDet(rhit->geographicalId());
01506                 
01507                 TrajectoryStateOnSurface predTsos;
01508                 TrajectoryStateOnSurface currTsos;
01509                 
01510                 if (traj.measurements().empty()) {
01511                   //predTsos = thePropagator->propagate(fts, det->surface());
01512                   if (debug_) std::cout<<"BIG ERROR!!! How did we make it to here with no trajectory measurements?!?!?"<<std::endl;
01513                 } else {
01514                   currTsos = traj.measurements().back().updatedState();
01515                   predTsos = thePropagator->propagate(currTsos, det->surface());
01516                 }
01517                 if (!predTsos.isValid()){
01518                   continue;
01519                 }
01520                 GlobalVector propagationDistance = predTsos.globalPosition() - currTsos.globalPosition();
01521                 if (propagationDistance.mag() > maxPropagationDistance) continue;
01522                 if (debug_) {
01523                   std::cout << "currTsos (x/y/z): " 
01524                             << currTsos.globalPosition().x() << " / "
01525                             << currTsos.globalPosition().y() << " / "
01526                             << currTsos.globalPosition().z() << std::endl;
01527                   std::cout << "predTsos (x/y/z): " 
01528                             << predTsos.globalPosition().x() << " / "
01529                             << predTsos.globalPosition().y() << " / "
01530                             << predTsos.globalPosition().z() << std::endl;
01531                   std::cout << "Propagation distance1 is " << propagationDistance.mag() << std::endl;
01532                 }
01533                 TrajectoryMeasurement tm;
01534                 if (debug_){
01535                   std::cout<< "trajectory at r/z=" <<  det->surface().position().perp() 
01536                            << "  " <<  det->surface().position().z() 
01537                            << ", hit " << ihit-theLayerHits.begin()
01538                            << " local prediction " << predTsos.localPosition().x() 
01539                            << " +- " << sqrt(predTsos.localError().positionError().xx()) 
01540                            << ", hit at " << rhit->localPosition().x() << " +- " << sqrt(rhit->localPositionError().xx())
01541                            << std::endl;
01542                 }
01543                 
01544                 // update
01545                 // first correct for angle
01546 
01547                 const SiStripMatchedRecHit2D *origHit = dynamic_cast<const SiStripMatchedRecHit2D *>(*ihit);
01548                 if (origHit !=0){
01549                   const GluedGeomDet *gdet = dynamic_cast<const GluedGeomDet*>(rhit->det());
01550                   const SiStripMatchedRecHit2D *corrHit = theHitMatcher->match(origHit,gdet,predTsos.localDirection());
01551                   if (corrHit!=0){
01552                     rhit = ttrhBuilder->build(&(*corrHit));
01553                     delete corrHit;
01554                   }
01555                 }
01556 
01557                 MeasurementEstimator::HitReturnType est = theEstimator->estimate(predTsos, *rhit);
01558                 if (debug_) std::cout << "estimation: " << est.first << " " << est.second << std::endl;
01559                 if (!est.first) continue;
01560                 currTsos = theUpdator->update(predTsos, *rhit);
01561                 tm = TrajectoryMeasurement(predTsos, currTsos, &(*rhit),est.second,thisLayer);
01562                 traj.push(tm,est.second);
01563                 newTrajectories.push_back(traj);
01564 
01565 
01566               }
01567 
01568             return newTrajectories;
01569 }

std::vector< TrajectoryMeasurement > RoadSearchTrackCandidateMakerAlgorithm::FindBestHit	(	const TrajectoryStateOnSurface &	tsosBefore,
		const std::set< const GeomDet * > &	theDets,
		edm::OwnVector< TrackingRecHit > &	theHits
	)

Definition at line 896 of file RoadSearchTrackCandidateMakerAlgorithm.cc.

References edm::OwnVector< T, P >::begin(), TransientTrackingRecHitBuilder::build(), GenMuonPlsPt100GeV_cfg::cout, debug_, GeometricSearchTracker::detLayer(), edm::OwnVector< T, P >::end(), lat::endl(), MeasurementEstimator::estimate(), MeasurementTracker::geometricSearchTracker(), TrajectoryStateOnSurface::globalPosition(), TrackerGeometry::idToDet(), TrajectoryStateOnSurface::isValid(), PV3DBase< T, PVType, FrameType >::mag(), maxPropagationDistance, PropagatorWithMaterial::propagate(), TrajectoryStateOnSurface::surface(), theEstimator, theMeasurementTracker, thePropagator, theUpdator, trackerGeom, ttrhBuilder, and TrajectoryStateUpdator::update().

00899 {
00900   
00901   std::vector<TrajectoryMeasurement> theBestHits;
00902   
00903   double bestchi = 10000.0;
00904   // extrapolate to all detectors from the list
00905   std::map<const GeomDet*, TrajectoryStateOnSurface> dmmap;
00906   for (std::set<const GeomDet*>::iterator idet = theDets.begin();
00907        idet != theDets.end(); ++idet) {
00908     TrajectoryStateOnSurface predTsos = thePropagator->propagate(tsosBefore, (**idet).surface());
00909     if (predTsos.isValid()) {
00910       GlobalVector propagationDistance = predTsos.globalPosition() - tsosBefore.globalPosition();
00911       if (propagationDistance.mag() > maxPropagationDistance) continue; 
00912       dmmap.insert(std::make_pair(*idet, predTsos));
00913     }
00914   }
00915   // evaluate hit residuals
00916   std::map<const GeomDet*, TrajectoryMeasurement> dtmmap;
00917   for (edm::OwnVector<TrackingRecHit>::const_iterator ih = theHits.begin();
00918        ih != theHits.end(); ++ih) {
00919     const GeomDet* det = trackerGeom->idToDet(ih->geographicalId());
00920     
00921     std::map<const GeomDet*, TrajectoryStateOnSurface>::iterator idm = dmmap.find(det);
00922     if (idm == dmmap.end()) continue;
00923     TrajectoryStateOnSurface predTsos = idm->second;
00924     TransientTrackingRecHit::RecHitPointer rhit = ttrhBuilder->build(&(*ih));
00925     MeasurementEstimator::HitReturnType est = theEstimator->estimate(predTsos, *rhit);
00926     
00927     // Take the best hit on any Det
00928     if (est.first) {
00929       TrajectoryStateOnSurface currTsos = theUpdator->update(predTsos, *rhit);
00930       if (est.second < bestchi){
00931         if(!theBestHits.empty()){
00932           theBestHits.erase(theBestHits.begin());
00933         }
00934         bestchi = est.second;
00935         TrajectoryMeasurement tm(predTsos, currTsos, &(*rhit),est.second,
00936                                  theMeasurementTracker->geometricSearchTracker()->detLayer(ih->geographicalId()));      
00937         theBestHits.push_back(tm);
00938       }
00939     }
00940   }
00941   
00942   if (theBestHits.empty()){
00943     if (debug_) std::cout<< "no hits to add! " <<std::endl;
00944   }
00945   else{
00946     for (std::vector<TrajectoryMeasurement>::const_iterator im=theBestHits.begin();im!=theBestHits.end();++im)
00947       if (debug_) std::cout<<" Measurement on layer "
00948                            << theMeasurementTracker->geometricSearchTracker()->detLayer(im->recHit()->geographicalId())
00949                            << " with estimate " << im->estimate()<<std::endl ;
00950   }
00951   
00952   return theBestHits;
00953 }

std::vector< TrajectoryMeasurement > RoadSearchTrackCandidateMakerAlgorithm::FindBestHits	(	const TrajectoryStateOnSurface &	tsosBefore,
		const std::set< const GeomDet * > &	theDets,
		const SiStripRecHitMatcher *	theHitMatcher,
		edm::OwnVector< TrackingRecHit > &	theHits
	)

Definition at line 956 of file RoadSearchTrackCandidateMakerAlgorithm.cc.

References edm::OwnVector< T, P >::begin(), TransientTrackingRecHitBuilder::build(), GenMuonPlsPt100GeV_cfg::cout, debug_, GeometricSearchTracker::detLayer(), edm::OwnVector< T, P >::end(), lat::endl(), MeasurementEstimator::estimate(), TrajectoryMeasurement::estimate(), MeasurementTracker::geometricSearchTracker(), TrajectoryStateOnSurface::globalPosition(), TrackerGeometry::idToDet(), TrajectoryStateOnSurface::isValid(), TrajectoryStateOnSurface::localDirection(), TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localPosition(), PV3DBase< T, PVType, FrameType >::mag(), SiStripRecHitMatcher::match(), maxPropagationDistance, LocalTrajectoryError::positionError(), PropagatorWithMaterial::propagate(), TrajectoryMeasurement::recHit(), python::multivaluedict::sort(), funct::sqrt(), GeomDet::surface(), TrajectoryStateOnSurface::surface(), theEstimator, theMeasurementTracker, thePropagator, theUpdator, trackerGeom, ttrhBuilder, TrajectoryStateUpdator::update(), PV3DBase< T, PVType, FrameType >::x(), LocalError::xx(), PV3DBase< T, PVType, FrameType >::y(), LocalError::yy(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by run().

00961 {
00962   
00963 
00964   std::vector<TrajectoryMeasurement> theBestHits;
00965   //TrajectoryMeasurement* theBestTM = 0;
00966   TrajectoryMeasurement theBestTM;
00967   bool firstTM = true;
00968   
00969   // extrapolate to all detectors from the list
00970   std::map<const GeomDet*, TrajectoryStateOnSurface> dmmap;
00971   for (std::set<const GeomDet*>::iterator idet = theDets.begin();
00972        idet != theDets.end(); ++idet) {
00973     TrajectoryStateOnSurface predTsos = thePropagator->propagate(tsosBefore, (**idet).surface());
00974     if (predTsos.isValid()) {
00975       GlobalVector propagationDistance = predTsos.globalPosition() - tsosBefore.globalPosition();
00976       if (propagationDistance.mag() > maxPropagationDistance) continue; 
00977       dmmap.insert(std::make_pair(*idet, predTsos));
00978     }
00979   }
00980 
00981   if(debug_) std::cout << "TRAJECTORY INTERSECTS " << dmmap.size() << " DETECTORS." << std::endl;
00982 
00983   // evaluate hit residuals
00984   std::map<const GeomDet*, TrajectoryMeasurement> dtmmap;
00985   for (edm::OwnVector<TrackingRecHit>::const_iterator ih = theHits.begin(); ih != theHits.end(); ++ih) {
00986     const GeomDet* det = trackerGeom->idToDet(ih->geographicalId());
00987     //if (*isl != theMeasurementTracker->geometricSearchTracker()->detLayer(ih->geographicalId())) 
00988     //  std::cout <<" You don't know what you're doing !!!!" << std::endl;
00989     
00990     std::map<const GeomDet*, TrajectoryStateOnSurface>::iterator idm = dmmap.find(det);
00991     if (idm == dmmap.end()) continue;
00992     TrajectoryStateOnSurface predTsos = idm->second;
00993     TransientTrackingRecHit::RecHitPointer rhit = ttrhBuilder->build(&(*ih));
00994     
00995     const SiStripMatchedRecHit2D *origHit = dynamic_cast<const SiStripMatchedRecHit2D *>(&(*ih));
00996     if (origHit !=0){
00997       const GluedGeomDet *gdet = dynamic_cast<const GluedGeomDet*>(det);
00998       const SiStripMatchedRecHit2D *corrHit = theHitMatcher->match(origHit,gdet,predTsos.localDirection());
00999       if (corrHit!=0){
01000         rhit = ttrhBuilder->build(&(*corrHit));
01001         delete corrHit;
01002       }
01003     }
01004 
01005     if (debug_) {
01006       std::cout << "predTSOS (x/y/z): " << predTsos.globalPosition().x()<< " / " << predTsos.globalPosition().y()<< " / " << predTsos.globalPosition().z() << std::endl;
01007       std::cout << "local position x: " << predTsos.localPosition().x() << "+-" << sqrt(predTsos.localError().positionError().xx()) << std::endl;
01008       std::cout << "local position y: " << predTsos.localPosition().y() << "+-" << sqrt(predTsos.localError().positionError().yy()) << std::endl;
01009       std::cout << "rhit local position x: " << rhit->localPosition().x() << "+-" << sqrt(rhit->localPositionError().xx()) << std::endl;
01010       std::cout << "rhit local position y: " << rhit->localPosition().y() << "+-" << sqrt(rhit->localPositionError().yy()) << std::endl;
01011     }
01012  
01013     MeasurementEstimator::HitReturnType est = theEstimator->estimate(predTsos, *rhit);
01014     if (debug_) std::cout<< "hit " << ih-theHits.begin() 
01015                          << ": est = " << est.first << " " << est.second  <<std::endl;
01016     
01017     
01018     // Take the best hit on a given Det
01019     if (est.first) {
01020       TrajectoryMeasurement tm;
01021       TrajectoryStateOnSurface currTsos = theUpdator->update(predTsos, *rhit);
01022       std::map<const GeomDet*, TrajectoryMeasurement>::iterator idtm = dtmmap.find(det);
01023       if (idtm == dtmmap.end()) {
01024         tm = TrajectoryMeasurement (predTsos, currTsos, &(*rhit),est.second,
01025                                     theMeasurementTracker->geometricSearchTracker()->detLayer(ih->geographicalId()));
01026         dtmmap.insert(std::make_pair(det, tm));
01027       } else if (idtm->second.estimate() > est.second) {
01028         dtmmap.erase(idtm);
01029         tm = TrajectoryMeasurement(predTsos, currTsos, &(*rhit),est.second,
01030                                    theMeasurementTracker->geometricSearchTracker()->detLayer(ih->geographicalId()));
01031         dtmmap.insert(std::make_pair(det, tm));
01032       }
01033       if ((firstTM)){
01034         theBestTM = tm; 
01035         if (debug_) std::cout <<"Initialize best to " << theBestTM.estimate() << std::endl;
01036         firstTM = false;
01037       }
01038       else if (!firstTM) {
01039         if (debug_) std::cout << "Current best is " << theBestTM.estimate() << " while this hit is " << est.second;
01040         if (est.second < theBestTM.estimate()) {
01041           if (debug_) std::cout << " so replace it " ;
01042           theBestTM = tm;
01043         }
01044         if (debug_) std::cout << std::endl;
01045       }
01046     }
01047   }
01048   if (debug_) std::cout<< "Hits(Dets) to add: " << dtmmap.size() <<std::endl;
01049   if (!dtmmap.empty()) {
01050     
01051     std::vector<std::pair<TransientTrackingRecHit::ConstRecHitPointer, TrajectoryMeasurement*> > OverlapHits;
01052     for (std::map<const GeomDet*, TrajectoryMeasurement>::iterator idtm = dtmmap.begin();
01053          idtm != dtmmap.end(); ++idtm) {
01054       OverlapHits.push_back(std::make_pair(idtm->second.recHit(),&idtm->second));
01055       
01056       if (debug_) std::cout<<" Measurement on layer "
01057                            << theMeasurementTracker->geometricSearchTracker()->detLayer(idtm->second.recHit()->geographicalId())
01058                            << " with estimate " << idtm->second.estimate()<<std::endl ;
01059     }
01060     if (debug_)
01061       std::cout<<" Best  Measurement is on layer "
01062                << theMeasurementTracker->geometricSearchTracker()->detLayer(theBestTM.recHit()->geographicalId())
01063                << " with estimate " << theBestTM.estimate()<<std::endl ;
01064     
01065     
01066     if (dtmmap.size()==1){  // only one hit so we can just return that one
01067       for (std::map<const GeomDet*, TrajectoryMeasurement>::iterator idtm = dtmmap.begin();
01068            idtm != dtmmap.end(); ++idtm) {
01069         TrajectoryMeasurement itm = idtm->second;
01070         if (debug_) std::cout<<" Measurement on layer "
01071                              << theMeasurementTracker->geometricSearchTracker()->detLayer(itm.recHit()->geographicalId())
01072                              << " with estimate " << itm.estimate()<<std::endl ;
01073         theBestHits.push_back(itm);
01074       }
01075     }
01076     else if (dtmmap.size()>=2) { // try for the overlaps -- first have to sort inside out
01077       
01078       if (debug_) std::cout<<"Unsorted OverlapHits has size " <<OverlapHits.size() << std::endl;
01079       
01080       for (std::vector<std::pair<TransientTrackingRecHit::ConstRecHitPointer,TrajectoryMeasurement*> >::iterator irh =OverlapHits.begin();
01081            irh!=OverlapHits.end();++irh){
01082         if (debug_) std::cout << "Hit " << irh-OverlapHits.begin()
01083                               << " on det " << irh->first->det() 
01084                               << " detLayer " 
01085                               << theMeasurementTracker->geometricSearchTracker()->detLayer(irh->first->geographicalId())
01086                               << ", r/phi/z = "
01087                               << irh->first->globalPosition().perp() << " "
01088                               << irh->first->globalPosition().phi() << " "
01089                               << irh->first->globalPosition().z()
01090                               << std::endl;
01091       }
01092       
01093       std::sort( OverlapHits.begin(),OverlapHits.end(),SortHitTrajectoryPairsByGlobalPosition());
01094       if (debug_) std::cout<<"Sorted OverlapHits has size " <<OverlapHits.size() << std::endl;
01095       
01096       float workingBestChi2 = 1000000.0;
01097       std::vector<TrajectoryMeasurement> workingBestHits;
01098       
01099       std::vector<std::pair<TransientTrackingRecHit::ConstRecHitPointer,TrajectoryMeasurement*> >::iterator irh1;
01100       std::vector<std::pair<TransientTrackingRecHit::ConstRecHitPointer,TrajectoryMeasurement*> >::iterator irh2;
01101       for (irh1 =OverlapHits.begin(); irh1!=--OverlapHits.end(); ++irh1){
01102         theBestHits.clear();
01103         float running_chi2=0;
01104         if (debug_) std::cout << "Hit " << irh1-OverlapHits.begin()
01105                               << " on det " << irh1->first->det() 
01106                               << " detLayer " 
01107                               << theMeasurementTracker->geometricSearchTracker()->detLayer(irh1->first->geographicalId())
01108                               << ", r/phi/z = "
01109           
01110                               << irh1->first->globalPosition().perp() << " "
01111                               << irh1->first->globalPosition().phi() << " "
01112                               << irh1->first->globalPosition().z()
01113                               << std::endl;
01114         
01115         TrajectoryStateOnSurface currTsos = irh1->second->updatedState();
01116         TransientTrackingRecHit::ConstRecHitPointer rhit = irh1->first;
01117         theBestHits.push_back(*(irh1->second));
01118         if (debug_)  std::cout<<"Added first hit with chi2 = " << irh1->second->estimate() << std::endl;
01119         running_chi2 += irh1->second->estimate();
01120         for (irh2 = irh1; irh2!=OverlapHits.end(); ++irh2){
01121           if (irh2 == irh1) continue;
01122           TransientTrackingRecHit::ConstRecHitPointer rh = irh2->first;
01123           const GeomDet* det = irh2->first->det();
01124           // extrapolate the trajectory to the next hit
01125           TrajectoryStateOnSurface predTsos = thePropagator->propagate(currTsos, det->surface());
01126           // test if matches
01127           if (predTsos.isValid()){
01128             MeasurementEstimator::HitReturnType est = theEstimator->estimate(predTsos, *rh);
01129             if (debug_)  std::cout<<"Added overlap hit with est = " << est.first << "   " << est.second << std::endl;
01130             if (est.first){
01131               TrajectoryMeasurement tm(predTsos, currTsos, &(*rh),est.second,
01132                                        theMeasurementTracker->geometricSearchTracker()->detLayer(rh->geographicalId()));
01133               theBestHits.push_back(tm);
01134               running_chi2 += est.second ;
01135             }
01136             else { // couldn't add 2nd hit so return best single hit
01137             }
01138           }
01139           
01140         }
01141         if (theBestHits.size()==dtmmap.size()){ // added the best hit in every layer
01142           if (debug_) std::cout<<"Added all "<<theBestHits.size()<<" hits out of " << dtmmap.size() << std::endl;
01143           break;
01144         }
01145         // Didn't add hits from every Det
01146         if (theBestHits.size() < dtmmap.size()){
01147           if (debug_) std::cout<<"Added only "<<theBestHits.size()<<" hits out of " << dtmmap.size() << std::endl;
01148           // Take the combination with the most hits
01149           if (theBestHits.size() > workingBestHits.size()){
01150             if (debug_) std::cout<<"Current combo has more hits so replace best" << std::endl;
01151             workingBestHits = theBestHits;
01152           }
01153           // has same number of hits as best, so check chi2
01154           else if (theBestHits.size() == workingBestHits.size()){ 
01155             if (running_chi2< workingBestChi2){
01156               if (debug_) std::cout<<"Current combo has same # of hits but lower chi2 so replace best" << std::endl;
01157               workingBestHits = theBestHits;
01158               workingBestChi2 = running_chi2;
01159             }
01160           }
01161         }
01162       }
01163       if (theBestHits.size()<2){
01164         if (debug_) std::cout<<"Only one good hit in overlap"<<std::endl;
01165         if (debug_) std::cout<<" Added hit on layer on det " 
01166                              << theBestTM.recHit()->det() 
01167                              << " detLayer " 
01168                              << theMeasurementTracker->geometricSearchTracker()->detLayer(theBestTM.recHit()->geographicalId())
01169                              << ", r/phi/z = "
01170                              << theBestTM.recHit()->globalPosition().perp() << " "
01171                              << theBestTM.recHit()->globalPosition().phi() << " "
01172                              << theBestTM.recHit()->globalPosition().z()
01173                              << " with estimate " << theBestTM.estimate()<<std::endl ;
01174         theBestHits.clear();
01175         theBestHits.push_back(theBestTM);
01176       }
01177       
01178     }
01179     else {
01180       if (debug_)std::cout << "ERROR: Unexpected size from DTMMAP = " << dtmmap.size() << std::endl;
01181       theBestHits.push_back(theBestTM);
01182     }
01183   }
01184   
01185   return theBestHits;
01186 }

std::vector< TrajectoryMeasurement > RoadSearchTrackCandidateMakerAlgorithm::FindBestHitsByDet	(	const TrajectoryStateOnSurface &	tsosBefore,
		const std::set< const GeomDet * > &	theDets,
		edm::OwnVector< TrackingRecHit > &	theHits
	)

Definition at line 834 of file RoadSearchTrackCandidateMakerAlgorithm.cc.

References edm::OwnVector< T, P >::begin(), TransientTrackingRecHitBuilder::build(), GeometricSearchTracker::detLayer(), edm::OwnVector< T, P >::end(), MeasurementEstimator::estimate(), MeasurementTracker::geometricSearchTracker(), TrajectoryStateOnSurface::globalPosition(), TrackerGeometry::idToDet(), TrajectoryStateOnSurface::isValid(), PV3DBase< T, PVType, FrameType >::mag(), maxPropagationDistance, PropagatorWithMaterial::propagate(), TrajectoryStateOnSurface::surface(), theEstimator, theMeasurementTracker, thePropagator, theUpdator, trackerGeom, ttrhBuilder, and TrajectoryStateUpdator::update().

00838 {
00839   
00840   //edm::OwnVector<TrackingRecHit> theBestHits;
00841   std::vector<TrajectoryMeasurement> theBestHits;
00842   
00843   // extrapolate to all detectors from the list
00844   std::map<const GeomDet*, TrajectoryStateOnSurface> dmmap;
00845   for (std::set<const GeomDet*>::iterator idet = theDets.begin();
00846        idet != theDets.end(); ++idet) {
00847     TrajectoryStateOnSurface predTsos = thePropagator->propagate(tsosBefore, (**idet).surface());
00848     if (predTsos.isValid()) {
00849       GlobalVector propagationDistance = predTsos.globalPosition() - tsosBefore.globalPosition();
00850       if (propagationDistance.mag() > maxPropagationDistance) continue; 
00851       dmmap.insert(std::make_pair(*idet, predTsos));
00852     }
00853   }
00854   // evaluate hit residuals
00855   std::map<const GeomDet*, TrajectoryMeasurement> dtmmap;
00856   for (edm::OwnVector<TrackingRecHit>::const_iterator ih = theHits.begin();
00857        ih != theHits.end(); ++ih) {
00858     const GeomDet* det = trackerGeom->idToDet(ih->geographicalId());
00859     
00860     std::map<const GeomDet*, TrajectoryStateOnSurface>::iterator idm = dmmap.find(det);
00861     if (idm == dmmap.end()) continue;
00862     TrajectoryStateOnSurface predTsos = idm->second;
00863     TransientTrackingRecHit::RecHitPointer rhit = ttrhBuilder->build(&(*ih));
00864     MeasurementEstimator::HitReturnType est = theEstimator->estimate(predTsos, *rhit);
00865     
00866     // Take the best hit on a given Det
00867     if (est.first) {
00868       TrajectoryStateOnSurface currTsos = theUpdator->update(predTsos, *rhit);
00869       std::map<const GeomDet*, TrajectoryMeasurement>::iterator idtm = dtmmap.find(det);
00870       if (idtm == dtmmap.end()) {
00871         TrajectoryMeasurement tm(predTsos, currTsos, &(*rhit),est.second,
00872                                  theMeasurementTracker->geometricSearchTracker()->detLayer(ih->geographicalId()));
00873         dtmmap.insert(std::make_pair(det, tm));
00874       } else if (idtm->second.estimate() > est.second) {
00875         dtmmap.erase(idtm);
00876         TrajectoryMeasurement tm(predTsos, currTsos, &(*rhit),est.second,
00877                                  theMeasurementTracker->geometricSearchTracker()->detLayer(ih->geographicalId()));
00878         dtmmap.insert(std::make_pair(det, tm));
00879       }
00880     }
00881   }
00882 
00883   if (!dtmmap.empty()) {
00884     for (std::map<const GeomDet*, TrajectoryMeasurement>::iterator idtm = dtmmap.begin();
00885          idtm != dtmmap.end(); ++idtm) {
00886       TrajectoryMeasurement itm = idtm->second;
00887       theBestHits.push_back(itm);
00888     }
00889   }
00890   
00891   return theBestHits;
00892 }

FreeTrajectoryState RoadSearchTrackCandidateMakerAlgorithm::initialTrajectory	(	const edm::EventSetup &	es,
		const TrackingRecHit *	InnerHit,
		const TrackingRecHit *	OuterHit
	)

Definition at line 1245 of file RoadSearchTrackCandidateMakerAlgorithm.cc.

References funct::abs(), funct::C, FastLine::c(), funct::cos(), CosmicReco_, cosmicSeedPt_, GenMuonPlsPt100GeV_cfg::cout, GlobalErrorBase< T, ErrorWeightType >::cxx(), GlobalErrorBase< T, ErrorWeightType >::czz(), debug_cff::d0, d1, debug_, lat::endl(), TrackingRecHit::geographicalId(), TrackerGeometry::idToDet(), initialError, initialVertexErrorXY_, SurfaceOrientation::inner, FastHelix::isValid(), l1, l2, TrackingRecHit::localPosition(), TrackingRecHit::localPositionError(), LogDebug, magField, FastLine::n1(), FastLine::n2(), NoFieldCosmic_, SurfaceOrientation::outer, FreeTrajectoryState::parameters(), PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), funct::sin(), funct::sqrt(), FastHelix::stateAtVertex(), theAloPropagator, thePropagator, theRevPropagator, GeomDet::toGlobal(), trackerGeom, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), LocalError::yy(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by run().

01248 {
01249   FreeTrajectoryState fts;
01250 
01251           GlobalPoint inner = trackerGeom->idToDet(theInnerHit->geographicalId())->surface().toGlobal(theInnerHit->localPosition());
01252           GlobalPoint outer = trackerGeom->idToDet(theOuterHit->geographicalId())->surface().toGlobal(theOuterHit->localPosition());
01253           
01254           LogDebug("RoadSearch") << "inner hit: r/phi/z = "<< inner.perp() << " " << inner.phi() << " " << inner.z() ;
01255           LogDebug("RoadSearch") << "outer hit: r/phi/z = "<< outer.perp() << " " << outer.phi() << " " << outer.z() ;
01256           
01257           // hits should be reasonably separated in r
01258           const double dRmin = 0.1; // cm
01259           if (outer.perp() - inner.perp() < dRmin) return fts;
01260           //GlobalPoint vertexPos(0,0,0);
01261           const double dr2 = initialVertexErrorXY_*initialVertexErrorXY_;
01262           //const double dr2 = 0.0015*0.0015;
01263           //const double dr2 = 0.2*0.2;
01264           const double dz2 = 5.3*5.3;
01265 
01266           // linear z extrapolation of two hits have to be inside tracker ( |z| < 275 cm)
01267           FastLine linearFit(outer, inner);
01268           double z_0 = -linearFit.c()/linearFit.n2();
01269           if ( std::abs(z_0) > 275 && !CosmicReco_ ) return fts;
01270 
01271           GlobalError vertexErr(dr2,
01272                                 0, dr2,
01273                                 0, 0, dz2);
01274           //TrivialVertex vtx( vertexPos, vertexErr);
01275           //FastHelix helix(outerHit.globalPosition(),
01276           //              (*innerHit).globalPosition(),
01277           //              vtx.position());
01278           
01279           double x0=0.0,y0=0.0,z0=0.0;
01280           double phi0 = -999.0;
01281           if (NoFieldCosmic_){
01282             phi0=atan2(outer.y()-inner.y(),outer.x()-inner.x());
01283             thePropagator = theAloPropagator;//GC
01284             if (inner.y()<outer.y()){
01285               if (debug_) std::cout<<"Flipping direction!!!" << std::endl;
01286               thePropagator = theRevPropagator;
01287               phi0=phi0-M_PI;
01288             } 
01289             double alpha=atan2(inner.y(),inner.x());
01290             double d1=sqrt(inner.x()*inner.x()+inner.y()*inner.y());
01291             double d0=-d1*sin(alpha-phi0); x0=d0*sin(phi0); y0=-d0*cos(phi0);
01292             double l1=0.0,l2=0.0;
01293             if (fabs(cos(phi0))>0.1){
01294               l1=(inner.x()-x0)/cos(phi0);l2=(outer.x()-x0)/cos(phi0);
01295             }else{
01296               l1=(inner.y()-y0)/sin(phi0);l2=(outer.y()-y0)/sin(phi0);
01297             }
01298             z0=(l2*inner.z()-l1*outer.z())/(l2-l1);
01299           }
01300           //FastHelix helix(outer, inner, vertexPos, es);
01301           FastHelix helix(outer, inner, GlobalPoint(x0,y0,z0), es);
01302           if (!NoFieldCosmic_ && !helix.isValid()) return fts;
01303           
01304           AlgebraicSymMatrix55 C = AlgebraicMatrixID();
01305           float zErr = vertexErr.czz();
01306           float transverseErr = vertexErr.cxx(); // assume equal cxx cyy
01307           C(3, 3) = transverseErr;
01308           C(4, 4) = zErr;
01309           CurvilinearTrajectoryError initialError(C);
01310 
01311           if (NoFieldCosmic_) {
01312             TrackCharge q = 1;      
01313             FastLine flfit(outer, inner);
01314             double dzdr = -flfit.n1()/flfit.n2();
01315             if (inner.y()<outer.y()) dzdr*=-1;
01316 
01317             GlobalPoint XYZ0(x0,y0,z0);
01318             if (debug_) std::cout<<"Initial Point (x0/y0/z0): " << x0 <<'\t'<< y0 <<'\t'<< z0 << std::endl;
01319             GlobalVector PXYZ(cosmicSeedPt_*cos(phi0),
01320                               cosmicSeedPt_*sin(phi0),
01321                               cosmicSeedPt_*dzdr);
01322             GlobalTrajectoryParameters thePars(XYZ0,PXYZ,q,magField);
01323             AlgebraicSymMatrix66 CErr = AlgebraicMatrixID();
01324             CErr *= 5.0;
01325             // CErr(3,3) = (theInnerHit->localPositionError().yy()*theInnerHit->localPositionError().yy() +
01326             //           theOuterHit->localPositionError().yy()*theOuterHit->localPositionError().yy());
01327             fts = FreeTrajectoryState(thePars,
01328                                       CartesianTrajectoryError(CErr));
01329             if (debug_) std::cout<<"\nInitial CError (dx/dy/dz): " << CErr(1,1) <<'\t'<< CErr(2,2) <<'\t'<< CErr(3,3) << std::endl;
01330             if (debug_) std::cout<<"\n\ninner dy = " << theInnerHit->localPositionError().yy() <<"\t\touter dy = " << theOuterHit->localPositionError().yy() << std::endl;
01331           }
01332           else {
01333             fts = FreeTrajectoryState( helix.stateAtVertex().parameters(), initialError);
01334           }
01335           //                       RoadSearchSeedFinderAlgorithm::initialError( *outerHit, *(*innerHit),
01336           //                                  vertexPos, vertexErr));
01337 
01338           return fts;
01339 }

FreeTrajectoryState RoadSearchTrackCandidateMakerAlgorithm::initialTrajectoryFromTriplet	(	const edm::EventSetup &	es,
		const TrackingRecHit *	InnerHit,
		const TrackingRecHit *	MiddleHit,
		const TrackingRecHit *	OuterHit
	)

Definition at line 1341 of file RoadSearchTrackCandidateMakerAlgorithm.cc.

References funct::abs(), funct::C, FastLine::c(), GlobalTrajectoryParameters::charge(), CosmicReco_, GenMuonPlsPt100GeV_cfg::cout, debug_, lat::endl(), TrackingRecHit::geographicalId(), TrackerGeometry::idToDet(), initialError, initialVertexErrorXY_, SurfaceOrientation::inner, edm::isnan(), FastHelix::isValid(), TrackingRecHit::localPosition(), LogDebug, magField, middle, GlobalTrajectoryParameters::momentum(), FastLine::n2(), SurfaceOrientation::outer, FreeTrajectoryState::parameters(), PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), GlobalTrajectoryParameters::position(), FastHelix::stateAtVertex(), theAloPropagator, thePropagator, theRevPropagator, GeomDet::toGlobal(), trackerGeom, PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by run().

01345 {
01346   FreeTrajectoryState fts;
01347 
01348           GlobalPoint inner = trackerGeom->idToDet(theInnerHit->geographicalId())->surface().toGlobal(theInnerHit->localPosition());
01349           GlobalPoint middle= trackerGeom->idToDet(theMiddleHit->geographicalId())->surface().toGlobal(theMiddleHit->localPosition());
01350           GlobalPoint outer = trackerGeom->idToDet(theOuterHit->geographicalId())->surface().toGlobal(theOuterHit->localPosition());
01351           
01352           LogDebug("RoadSearch") << "inner hit: r/phi/z = "<< inner.perp() << " " << inner.phi() << " " << inner.z() ;
01353           LogDebug("RoadSearch") << "middlehit: r/phi/z = "<< inner.perp() << " " << inner.phi() << " " << inner.z() ;
01354           LogDebug("RoadSearch") << "outer hit: r/phi/z = "<< outer.perp() << " " << outer.phi() << " " << outer.z() ;
01355           
01356           // hits should be reasonably separated in r
01357           const double dRmin = 0.1; // cm
01358           if (outer.perp() - inner.perp() < dRmin) return fts;
01359           const double dr2 = initialVertexErrorXY_*initialVertexErrorXY_;
01360           const double dz2 = 5.3*5.3;
01361 
01362           // linear z extrapolation of two hits have to be inside tracker ( |z| < 275 cm)
01363           FastLine linearFit(outer, inner);
01364           double z_0 = -linearFit.c()/linearFit.n2();
01365           if ( std::abs(z_0) > 275 && !CosmicReco_ ) return fts;
01366 
01367 
01368           FastHelix helix(outer, middle, inner, es);
01369           // check that helix is OK
01370           if (!helix.isValid() || 
01371               isnan(helix.stateAtVertex().parameters().momentum().perp())) return fts;
01372           // simple cuts on pt and dz
01373           if (helix.stateAtVertex().parameters().momentum().perp() < 0.5 ||
01374               std::abs(helix.stateAtVertex().parameters().position().z()) > 550.0)
01375             return fts;
01376 
01377           AlgebraicSymMatrix55 C = AlgebraicMatrixID();
01378           float zErr = dz2;
01379           float transverseErr = dr2; // assume equal cxx cyy
01380           C(3, 3) = transverseErr;
01381           C(4, 4) = zErr;
01382           CurvilinearTrajectoryError initialError(C);
01383 
01384 
01385           thePropagator = theAloPropagator;//GC
01386           GlobalVector gv=helix.stateAtVertex().parameters().momentum();
01387           float charge=helix.stateAtVertex().parameters().charge();
01388 
01389           if (CosmicReco_ && gv.y()>0){
01390             if (debug_) std::cout<<"Flipping direction!!!" << std::endl;
01391             thePropagator = theRevPropagator;
01392             gv=-1.*gv;
01393             charge=-1.*charge;
01394           }
01395 
01396           GlobalTrajectoryParameters Gtp(inner,gv,int(charge),&(*magField));
01397           fts = FreeTrajectoryState(Gtp, initialError);
01398 
01399          //fts = FreeTrajectoryState( helix.stateAtVertex().parameters(), initialError);
01400 
01401           return fts;
01402 }

TrackCandidateCollection RoadSearchTrackCandidateMakerAlgorithm::PrepareTrackCandidates ( std::vector< Trajectory > & theTrajectories )

Definition at line 1572 of file RoadSearchTrackCandidateMakerAlgorithm.cc.

References alongMomentum, TrajectoryMeasurement::backwardPredictedState(), TrajectoryCleanerBySharedHits::clean(), clone(), CosmicTrackMerging_, GenMuonPlsPt100GeV_cfg::cout, debug_, debugCosmics_, GeometricSearchTracker::detLayer(), lat::endl(), MeasurementEstimator::estimate(), TrajectoryMeasurement::forwardPredictedState(), MeasurementTracker::geometricSearchTracker(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), i, TrackerGeometry::idToDet(), TrajectoryStateOnSurface::isValid(), it, j, k, Trajectory::lastMeasurement(), PV3DBase< T, PVType, FrameType >::mag(), maxPropagationDistance, Trajectory::measurements(), TrajectoryStateTransform::persistentState(), PV3DBase< T, PVType, FrameType >::phi(), AnalyticalPropagator::propagate(), PropagatorWithMaterial::propagate(), Trajectory::push(), edm::OwnVector< T, P >::push_back(), DetId::rawId(), TrajectoryMeasurement::recHit(), Trajectory::recHits(), Trajectory::seed(), splitMatchedHits_, state, GeomDet::surface(), t, theAloPropagator, theAnalyticalPropagator, theEstimator, theMeasurementTracker, theSmoother, theTrajectoryCleaner, theTransformer, theUpdator, trackerGeom, KFTrajectorySmoother::trajectories(), TrajectoryStateUpdator::update(), TrajectoryMeasurement::updatedState(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by run().

01573 {
01574 
01575   TrackCandidateCollection theCollection;
01576 
01577   theTrajectoryCleaner->clean(theTrajectories);
01578   
01579   //==========NEW CODE ==========
01580   
01581   if(CosmicTrackMerging_) {
01582     
01583     //generate vector of *valid* trajectories -> traj
01584     std::vector<Trajectory> traj;
01585 
01586     //keep track of trajectories which are used during merging
01587     std::vector<bool> trajUsed;
01588 
01589     for (std::vector<Trajectory>::iterator it = theTrajectories.begin(); it != theTrajectories.end(); ++it) {
01590       if (it->isValid()) {
01591         traj.push_back(*it);
01592         trajUsed.push_back(false);
01593       }
01594     }
01595     
01596     if(debugCosmics_) {
01597       std::cout << "==========ENTERING COSMIC MODE===========" << std::endl;
01598       int t=0;
01599       for (std::vector<Trajectory>::iterator it = traj.begin(); it != traj.end(); it++) {
01600         std::cout << "Trajectory " << it-traj.begin() << " has "<<it->recHits().size()<<" hits and is valid: " << it->isValid() << std::endl;
01601         TransientTrackingRecHit::ConstRecHitContainer itHits = it->recHits();
01602         
01603 
01604         for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator rhit=itHits.begin(); rhit!=itHits.end(); ++rhit)
01605           std::cout << "-->good hit position: " << (*rhit)->globalPosition().x() << ", " 
01606                     << (*rhit)->globalPosition().y() << ", "<< (*rhit)->globalPosition().z() << std::endl;
01607 
01608       }
01609     }
01610 
01611     //double nested looop to find trajectories that match in phi
01612     for ( unsigned int i = 0; i < traj.size(); ++i) {
01613       if (trajUsed[i]) continue;    
01614       for ( unsigned int j = i+1; j != traj.size(); ++j) {
01615         if (trajUsed[j]) continue;    
01616         
01617         if (debugCosmics_) std::cout<< "Trajectory " <<i<< " has "<<traj[i].recHits().size()<<" hits with chi2=" << traj[i].chiSquared() << " and is valid"<<std::endl;
01618         if (debugCosmics_) std::cout<< "Trajectory " <<j<< " has "<<traj[j].recHits().size()<<" hits with chi2=" << traj[j].chiSquared() << " and is valid"<<std::endl;           
01619         
01620         TrajectoryMeasurement firstTraj1 = traj[i].firstMeasurement();
01621         TrajectoryMeasurement firstTraj2 = traj[j].firstMeasurement();
01622         TrajectoryStateOnSurface firstTraj1TSOS = firstTraj1.updatedState();
01623         TrajectoryStateOnSurface firstTraj2TSOS = firstTraj2.updatedState();
01624 
01625         
01626         if(debugCosmics_) 
01627           std::cout << "phi1: " << firstTraj1TSOS.globalMomentum().phi() 
01628                     << " phi2: " << firstTraj2TSOS.globalMomentum().phi() 
01629                     << " --> delta_phi: " << firstTraj1TSOS.globalMomentum().phi()-firstTraj2TSOS.globalMomentum().phi() << std::endl;
01630         
01631         //generate new trajectory if delta_phi<0.3
01632         //use phi of momentum vector associated to *innermost* hit of trajectories
01633         if( fabs(firstTraj1TSOS.globalMomentum().phi()-firstTraj2TSOS.globalMomentum().phi())<0.3 ) {
01634           if(debugCosmics_) std::cout << "-->match successful" << std::endl;
01635         } else {
01636           if(debugCosmics_) std::cout << "-->match not successful" << std::endl;
01637         }
01638         if( fabs(firstTraj1TSOS.globalMomentum().phi()-firstTraj2TSOS.globalMomentum().phi())>0.3 ) continue;
01639         
01640         
01641         //choose starting trajectory: use trajectory in lower hemisphere (with y<0) to start new combined trajectory
01642         //use y position of outermost hit
01643         TrajectoryMeasurement lastTraj1 = traj[i].lastMeasurement();
01644         TrajectoryMeasurement lastTraj2 = traj[j].lastMeasurement();
01645         TrajectoryStateOnSurface lastTraj1TSOS = lastTraj1.updatedState();
01646         TrajectoryStateOnSurface lastTraj2TSOS = lastTraj2.updatedState();
01647         
01648         if(debugCosmics_){
01649           std::cout<<"Traj1 first (x/y/z): " 
01650                    << firstTraj1TSOS.globalPosition().x() <<" / "
01651                    << firstTraj1TSOS.globalPosition().y() <<" / "
01652                    << firstTraj1TSOS.globalPosition().z() 
01653                    << "   phi: " << firstTraj1TSOS.globalMomentum().phi() << std::endl;
01654           std::cout<<"Traj1  last (x/y/z): " 
01655                    << lastTraj1TSOS.globalPosition().x() <<" / "
01656                    << lastTraj1TSOS.globalPosition().y() <<" / "
01657                    << lastTraj1TSOS.globalPosition().z() 
01658                    << "   phi: " << lastTraj1TSOS.globalMomentum().phi() << std::endl;
01659 
01660           std::cout<<"Traj2 first (x/y/z): " 
01661                    << firstTraj2TSOS.globalPosition().x() <<" / "
01662                    << firstTraj2TSOS.globalPosition().y() <<" / "
01663                    << firstTraj2TSOS.globalPosition().z()
01664                    << "   phi: " << firstTraj2TSOS.globalMomentum().phi() << std::endl;
01665           std::cout<<"Traj2  last (x/y/z): " 
01666                    << lastTraj2TSOS.globalPosition().x() <<" / "
01667                    << lastTraj2TSOS.globalPosition().y() <<" / "
01668                    << lastTraj2TSOS.globalPosition().z()
01669                    << "   phi: " << lastTraj2TSOS.globalMomentum().phi() << std::endl;
01670 
01671         }
01672 
01673         Trajectory *upperTrajectory, *lowerTrajectory;
01674         
01675         TrajectoryStateOnSurface lowerTSOS1,upperTSOS1;
01676         if (lastTraj1TSOS.globalPosition().y()<firstTraj1TSOS.globalPosition().y()) {
01677           lowerTSOS1=lastTraj1TSOS; upperTSOS1=firstTraj1TSOS;
01678         }
01679         else {
01680           lowerTSOS1=firstTraj1TSOS; upperTSOS1=lastTraj1TSOS;
01681         }
01682         TrajectoryStateOnSurface lowerTSOS2;
01683         if (lastTraj2TSOS.globalPosition().y()<firstTraj2TSOS.globalPosition().y()) lowerTSOS2=lastTraj2TSOS;
01684         else lowerTSOS2=firstTraj2TSOS;
01685         if(lowerTSOS1.globalPosition().y() > lowerTSOS2.globalPosition().y()) {
01686           if(debugCosmics_) 
01687             std::cout << "-->case A: "<< lowerTSOS1.globalPosition().y() << " > " << lowerTSOS2.globalPosition().y() << std::endl;
01688           
01689           upperTrajectory = &(traj[i]);
01690           lowerTrajectory = &(traj[j]);
01691           
01692         } else {
01693           if(debugCosmics_) 
01694             std::cout << "-->case B: "<< lowerTSOS1.globalPosition().y() << " < " << lowerTSOS2.globalPosition().y() << std::endl;
01695           
01696           upperTrajectory = &(traj[j]);
01697           lowerTrajectory = &(traj[i]);
01698         } 
01699         
01700         std::vector<Trajectory> freshStartUpperTrajectory = theSmoother->trajectories(*upperTrajectory);
01701         std::vector<Trajectory> freshStartLowerTrajectory = theSmoother->trajectories(*lowerTrajectory);
01702         //--JR
01703         if (freshStartUpperTrajectory.empty() || freshStartLowerTrajectory .empty()){
01704           if (debugCosmics_) std::cout << " the smoother has failed."<<std::endl;
01705           continue;
01706         }
01707         //--JR
01708         TrajectoryStateOnSurface NewFirstTsos = freshStartUpperTrajectory.begin()->lastMeasurement().updatedState();
01709         TrajectoryStateOnSurface forwardTsos  = freshStartUpperTrajectory.begin()->firstMeasurement().forwardPredictedState();
01710         TrajectoryStateOnSurface backwardTsos = freshStartUpperTrajectory.begin()->lastMeasurement().backwardPredictedState();
01711 
01712         Trajectory freshStartTrajectory = *freshStartUpperTrajectory.begin();
01713         if(debugCosmics_) std::cout << "seed hit updatedState: " << NewFirstTsos.globalMomentum().x() << ", " << NewFirstTsos.globalMomentum().y() << ", " << NewFirstTsos.globalMomentum().z()  <<  std::endl;
01714         if(debugCosmics_) std::cout << "seed hit updatedState (pos x/y/z): " << NewFirstTsos.globalPosition().x() << ", " << NewFirstTsos.globalPosition().y() << ", " << NewFirstTsos.globalPosition().z()  <<  std::endl;
01715         if(debugCosmics_) std::cout << "seed hit forwardPredictedState: " << forwardTsos.globalMomentum().x() << ", " << forwardTsos.globalMomentum().y() << ", " << forwardTsos.globalMomentum().z()  <<  std::endl;
01716         if(debugCosmics_) std::cout << "seed hit forwardPredictedState (pos x/y/z): " << forwardTsos.globalPosition().x() << ", " << forwardTsos.globalPosition().y() << ", " << forwardTsos.globalPosition().z()  <<  std::endl;
01717         if(debugCosmics_) std::cout << "seed hit backwardPredictedState: " << backwardTsos.globalMomentum().x() << ", " << backwardTsos.globalMomentum().y() << ", " << backwardTsos.globalMomentum().z() <<  std::endl;
01718         if(debugCosmics_) std::cout << "seed hit backwardPredictedState (pos x/y/z): " << backwardTsos.globalPosition().x() << ", " << backwardTsos.globalPosition().y() << ", " << backwardTsos.globalPosition().z() <<  std::endl;
01719         
01720         if(debugCosmics_) std::cout<<"#hits for upper trajectory: " << freshStartTrajectory.measurements().size() << std::endl;
01721         
01722         //loop over hits in upper trajectory and add them to lower trajectory
01723         TransientTrackingRecHit::ConstRecHitContainer ttHits = lowerTrajectory->recHits(splitMatchedHits_);
01724         
01725         if(debugCosmics_) std::cout << "loop over hits in lower trajectory..." << std::endl;
01726         
01727         bool addHitToFreshStartTrajectory = false;
01728         bool propagationFailed = false;
01729         int lostHits = 0;
01730         for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator rhit=ttHits.begin(); rhit!=ttHits.end(); ++rhit){
01731 
01732           if(debugCosmics_ && lostHits>0){
01733             std::cout << " Lost " << lostHits << " of " << ttHits.size() << " on lower trajectory " << std::endl;
01734             std::cout << " Lost " << ((float)lostHits/(float)ttHits.size()) << " of hits of on lower trajectory " << std::endl;
01735           }
01736           if ((float)lostHits/(float)ttHits.size() > 0.5) {
01737               propagationFailed = true;
01738               break;
01739           }
01740           if(debugCosmics_) std::cout << "-->hit position: " << (*rhit)->globalPosition().x() << ", " << (*rhit)->globalPosition().y() << ", "<< (*rhit)->globalPosition().z() << std::endl;
01741           
01742           TrajectoryStateOnSurface predTsos;
01743           TrajectoryStateOnSurface currTsos;
01744           TrajectoryMeasurement tm;
01745           
01746           TransientTrackingRecHit::ConstRecHitPointer rh = (*rhit);
01747 
01748 
01749 
01750           /*
01751           if( rh->isValid() ) { 
01752             if(debugCosmics_) std::cout << "-->hit is valid"<<std::endl;
01753             const GeomDet* det = trackerGeom->idToDet(rh->geographicalId());
01754             
01755             std::vector<TrajectoryMeasurement> measL = freshStartTrajectory.measurements();
01756             bool alreadyUsed = false;         
01757             for (std::vector<TrajectoryMeasurement>::iterator mh=measL.begin();mh!=measL.end();++mh) {
01758               if (rh->geographicalId().rawId()==mh->recHit()->geographicalId().rawId()) {
01759                 if (debugCosmics_) std::cout << "this hit is already in the trajectory, skip it" << std::endl;
01760                 alreadyUsed = true;
01761                 break;
01762               }
01763             }
01764             if (alreadyUsed) continue;
01765             //std::vector<TrajectoryMeasurement> measU = freshStartUpperTrajectory[0].measurements();
01766             if (freshStartTrajectory.direction()==0) {
01767               std::vector<TrajectoryMeasurement>::iterator ml;
01768               for (ml=measL.begin();ml!=measL.end();++ml) {
01769                 if (debugCosmics_)  std::cout << "hit y="<<ml->recHit()->globalPosition().y()
01770                                               << " tsos validity fwd="<<ml->forwardPredictedState().isValid() 
01771                                               << " bwd="<<ml->backwardPredictedState().isValid() 
01772                                               << " upd="<<ml->updatedState().isValid() 
01773                                               <<std::endl;
01774                 if (ml->recHit()->globalPosition().y()>(*rhit)->globalPosition().y() && ml!=measL.begin()) {
01775                   break;
01776                 }
01777               }
01778               if ((ml-1)->forwardPredictedState().isValid()) currTsos = (ml-1)->forwardPredictedState();
01779               else currTsos = (ml-1)->backwardPredictedState();
01780               
01781               if (debugCosmics_) std::cout << "currTsos y=" << currTsos.globalPosition().y() << std::endl;
01782             }
01783             else {
01784               std::vector<TrajectoryMeasurement>::reverse_iterator ml;
01785               for (ml=measL.rbegin();ml!=measL.rend();++ml) {
01786                 if (debugCosmics_) std::cout << "hit y="<<ml->recHit()->globalPosition().y()
01787                                              << " tsos validity fwd="<<ml->forwardPredictedState().isValid() 
01788                                              << " bwd="<<ml->backwardPredictedState().isValid() 
01789                                              << " upd="<<ml->updatedState().isValid() 
01790                                              <<std::endl;
01791                 if (ml->recHit()->globalPosition().y()>(*rhit)->globalPosition().y() && ml!=measL.rbegin()) {
01792                   break;
01793                 }
01794               }
01795               if ((ml-1)->forwardPredictedState().isValid()) {
01796                 currTsos = (ml-1)->forwardPredictedState();
01797               }
01798               else {
01799                 currTsos = (ml-1)->backwardPredictedState();
01800               }
01801               
01802               if (debugCosmics_) std::cout << "reverse. currTsos y=" << currTsos.globalPosition().y() << std::endl;
01803             }
01804             
01805             
01806           }
01807           */
01808 
01809 
01810           if( rh->isValid() ) { 
01811             if(debugCosmics_) std::cout << "-->hit is valid"<<std::endl;
01812             const GeomDet* det = trackerGeom->idToDet(rh->geographicalId());
01813             if( addHitToFreshStartTrajectory==false ) {
01814               //first hit from upper trajectory that is being added to lower trajectory requires usage of backwardPredictedState (of lower trajectory)
01815               currTsos = freshStartTrajectory.lastMeasurement().backwardPredictedState();
01816             } else {
01817               //remaining hits from upper trajectory that are being added to lower trajectory require usage of forwardPredictedState
01818               currTsos = freshStartTrajectory.lastMeasurement().forwardPredictedState();
01819             }
01820             
01821             if(debugCosmics_) std::cout << "current TSOS: " << currTsos.globalPosition().x() << ", " << currTsos.globalPosition().y() << ", " << currTsos.globalPosition().z() << ", " << std::endl;
01822             
01823             predTsos = theAloPropagator->propagate(currTsos, det->surface());
01824             
01825             if (!predTsos.isValid()) {
01826               if(debugCosmics_) std::cout<<"predTsos is not valid!" <<std::endl;
01827               //propagationFailed = true;
01828               ++lostHits;
01829               //break;
01830               continue;
01831             }
01832             GlobalVector propagationDistance = predTsos.globalPosition() - currTsos.globalPosition();
01833             if (propagationDistance.mag() > maxPropagationDistance) continue;
01834             
01835             if(debugCosmics_) std::cout << "predicted TSOS: " << predTsos.globalPosition().x() << ", " << predTsos.globalPosition().y() << ", " << predTsos.globalPosition().z() << ", " << std::endl;
01836             MeasurementEstimator::HitReturnType est = theEstimator->estimate(predTsos, *rh);
01837             if (!est.first) {
01838               if(debugCosmics_) std::cout<<"Failed to add one of the original hits on a low occupancy layer!!!!" << std::endl;
01839               //propagationFailed = true;
01840               ++lostHits;
01841               //break;
01842               continue;
01843             }
01844             
01845             currTsos = theUpdator->update(predTsos, *rh);
01846             if(debugCosmics_) std::cout << "current updated TSOS: " << currTsos.globalPosition().x() << ", " << currTsos.globalPosition().y() << ", " << currTsos.globalPosition().z() << ", " << std::endl;
01847             tm = TrajectoryMeasurement(predTsos, currTsos,&(*rh),est.second,theMeasurementTracker->geometricSearchTracker()->detLayer(rh->geographicalId()));
01848             freshStartTrajectory.push(tm,est.second);
01849             addHitToFreshStartTrajectory=true;
01850             
01851           }
01852 
01853           if(debugCosmics_) std::cout<<"#hits for new trajectory (his from upper trajectory added): " << freshStartTrajectory.measurements().size() << std::endl;
01854         }
01855 
01856         if (propagationFailed) {
01857           if (debugCosmics_) std::cout<<"Propagation failed so go to next trajectory" << std::endl;
01858           continue;
01859         }
01860 
01861         //put final trajectory together
01862         if(debugCosmics_) std::cout << "put final trajectory together..." << std::endl;
01863         edm::OwnVector<TrackingRecHit> goodHits;
01864         TransientTrackingRecHit::ConstRecHitContainer tttempHits = freshStartTrajectory.recHits(splitMatchedHits_);
01865         
01866         for (int k=tttempHits.size()-1; k>=0; k--) {
01867           if(debugCosmics_) std::cout << "-->good hit position: " << tttempHits[k]->globalPosition().x() << ", " << tttempHits[k]->globalPosition().y() << ", "<< tttempHits[k]->globalPosition().z() << std::endl;
01868           goodHits.push_back(tttempHits[k]->hit()->clone());
01869         }
01870         TrajectoryStateOnSurface firstState;
01871         
01872         // check if Trajectory from seed is on first hit of the cloud, if not, propagate
01873         // exclude if first state on first hit is not valid
01874         DetId FirstHitId = (*(freshStartTrajectory.recHits().end()-1))->geographicalId(); //was begin
01875         
01876         TrajectoryMeasurement maxYMeasurement = freshStartTrajectory.lastMeasurement();
01877         const GeomDet* det = trackerGeom->idToDet(FirstHitId);
01878         firstState = theAnalyticalPropagator->propagate(maxYMeasurement.updatedState(),det->surface());
01879         if (firstState.isValid() == false) continue;    
01880         PTrajectoryStateOnDet *state = theTransformer->persistentState(firstState,FirstHitId.rawId());
01881         
01882         //generate new trajectory seed
01883         TrajectoryStateOnSurface firstTSOS = freshStartTrajectory.lastMeasurement().updatedState();
01884         if(debugCosmics_) std::cout << "generate new trajectory seed with hit (x/y/z): " << firstTSOS.globalPosition().x() << ", " << firstTSOS.globalPosition().y() << ", " << firstTSOS.globalPosition().z() << ", " << std::endl;
01885         TransientTrackingRecHit::ConstRecHitPointer rhit = freshStartTrajectory.lastMeasurement().recHit();
01886         PTrajectoryStateOnDet* pFirstState = TrajectoryStateTransform().persistentState(NewFirstTsos,rhit->geographicalId().rawId());
01887         edm::OwnVector<TrackingRecHit> newHits;
01888         newHits.push_back(rhit->hit()->clone());
01889         TrajectorySeed tmpseed = TrajectorySeed(*pFirstState,newHits,alongMomentum);
01890         
01891         theCollection.push_back(TrackCandidate(goodHits,freshStartTrajectory.seed(),*state));
01892         delete state;
01893         delete pFirstState;
01894         
01895         //trajectory usage
01896         trajUsed[i]=true;
01897         trajUsed[j]=true;
01898         
01899       } //for loop trajectory2
01900       
01901     } //for loop trajectory1
01902 
01903     //add all trajectories to the resulting vector if they have *not* been used by the trajectory merging algorithm
01904     for ( unsigned int i = 0; i < traj.size(); ++i) {
01905       
01906       if (trajUsed[i]==true) continue;
01907 
01908       if (debugCosmics_) std::cout<< "Trajectory (not merged) has "<<traj[i].recHits().size()<<" hits with chi2=" << traj[i].chiSquared() << " and is valid? "<< traj[i].isValid()<<std::endl;
01909       edm::OwnVector<TrackingRecHit> goodHits;
01910       TransientTrackingRecHit::ConstRecHitContainer ttHits = traj[i].recHits(splitMatchedHits_);
01911       for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator rhit=ttHits.begin(); rhit!=ttHits.end(); ++rhit){
01912         goodHits.push_back((*rhit)->hit()->clone());
01913       }
01914       TrajectoryStateOnSurface firstState;
01915       
01916       // check if Trajectory from seed is on first hit of the cloud, if not, propagate
01917       // exclude if first state on first hit is not valid
01918       DetId FirstHitId = (*(traj[i].recHits().begin()))->geographicalId();
01919       
01920       // propagate back to first hit    
01921       TrajectoryMeasurement firstMeasurement = traj[i].measurements()[0];
01922       const GeomDet* det = trackerGeom->idToDet(FirstHitId);
01923       firstState = theAnalyticalPropagator->propagate(firstMeasurement.updatedState(), det->surface());   
01924       if (firstState.isValid()) {
01925         PTrajectoryStateOnDet *state = theTransformer->persistentState(firstState,FirstHitId.rawId());
01926         theCollection.push_back(TrackCandidate(goodHits,traj[i].seed(),*state));
01927         delete state;
01928       }
01929     }
01930     if (debugCosmics_) std::cout << "Original collection had " << theTrajectories.size() 
01931                                  << " candidates.  Merged collection has " << theCollection.size() << std::endl;
01932   } //if(CosmicTrackMerging_)
01933   
01934   
01935   if(!CosmicTrackMerging_) {
01936      for (std::vector<Trajectory>::const_iterator it = theTrajectories.begin(); it != theTrajectories.end(); it++) {
01937         if (debug_) std::cout<< " Trajectory has "<<it->recHits().size()<<" hits with chi2=" << it->chiSquared() << " and is valid? "<<it->isValid()<<std::endl;
01938         if (it->isValid()){
01939 
01940            edm::OwnVector<TrackingRecHit> goodHits;
01941            TransientTrackingRecHit::ConstRecHitContainer ttHits = it->recHits(splitMatchedHits_);
01942            for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator rhit=ttHits.begin(); rhit!=ttHits.end(); ++rhit){
01943               goodHits.push_back((*rhit)->hit()->clone());
01944            }
01945            TrajectoryStateOnSurface firstState;
01946 
01947            // check if Trajectory from seed is on first hit of the cloud, if not, propagate
01948            // exclude if first state on first hit is not valid
01949            DetId FirstHitId = (*(it->recHits().begin()))->geographicalId();
01950 
01951            // propagate back to first hit    
01952            TrajectoryMeasurement firstMeasurement = it->measurements()[0];
01953            const GeomDet* det = trackerGeom->idToDet(FirstHitId);
01954            firstState = theAnalyticalPropagator->propagate(firstMeasurement.updatedState(), det->surface());      
01955            if (firstState.isValid() == false) continue;    
01956            PTrajectoryStateOnDet *state = theTransformer->persistentState(firstState,FirstHitId.rawId());
01957            theCollection.push_back(TrackCandidate(goodHits,it->seed(),*state));
01958            delete state;
01959 
01960         }
01961      }
01962   }
01963 
01964   return theCollection;
01965 }

void RoadSearchTrackCandidateMakerAlgorithm::run	(	const RoadSearchCloudCollection *	input,
		const edm::Event &	e,
		const edm::EventSetup &	es,
		TrackCandidateCollection &	output
	)

Runs the algorithm.

Definition at line 109 of file RoadSearchTrackCandidateMakerAlgorithm.cc.

References alongMomentum, AnalyticalPropagator_cfi::AnalyticalPropagator, anyDirection, chooseStartingLayers(), TrajectoryCleanerBySharedHits::clean(), conf_, CosmicReco_, GenMuonPlsPt100GeV_cfg::cout, createSeedTrajectory(), debug_, debugCosmics_, GeometricSearchTracker::detLayer(), TrajectorySeed::direction(), lat::endl(), MeasurementEstimator::estimate(), extrapolateTrajectory(), FindBestHits(), MeasurementTracker::geometricSearchTracker(), edm::EventSetup::get(), edm::ParameterSet::getParameter(), TrajectoryStateOnSurface::globalPosition(), FreeTrajectoryState::hasError(), i, TrackerGeometry::idToDet(), initialTrajectory(), initialTrajectoryFromTriplet(), MagneticField::inTesla(), Trajectory::isValid(), TrajectoryStateOnSurface::isValid(), it, KFTrajectorySmootherESProducer_cfi::KFTrajectorySmoother, Trajectory::lastMeasurement(), TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localPosition(), LogDebug, lstereo, PV3DBase< T, PVType, FrameType >::mag(), magField, maxPropagationDistance, Trajectory::measurements(), measurementTrackerName_, MinChunkLength_, nFoundMin_, NoFieldCosmic_, oppositeToMomentum, PV3DBase< T, PVType, FrameType >::perp(), LocalTrajectoryError::positionError(), PrepareTrackCandidates(), edm::ESHandle< T >::product(), edm::OwnVector< T, P >::push_back(), TrajectoryMeasurement::recHit(), edm::second(), edm::OwnVector< T, P >::size(), python::multivaluedict::sort(), funct::sqrt(), GeomDet::surface(), theAloPropagator, theAnalyticalPropagator, theEstimator, theHitMatcher, theMeasurementTracker, theNumHitCut, thePropagator, theRevPropagator, theSmoother, theTrajectoryCleaner, theUpdator, GeomDet::toGlobal(), trackerGeom, KFTrajectorySmoother::trajectories(), ttrhBuilder, MeasurementTracker::update(), TrajectoryStateUpdator::update(), TrajectoryMeasurement::updatedState(), PV3DBase< T, PVType, FrameType >::x(), LocalError::xx(), PV3DBase< T, PVType, FrameType >::y(), LocalError::yy(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by RoadSearchTrackCandidateMaker::produce().

00113 {
00114   
00115   //
00116   // right now, track candidates are just filled from cleaned
00117   // clouds. The trajectory of the seed is taken as the initial
00118   // trajectory for the final fit
00119   //
00120   
00121   //
00122   // get the transient builder
00123   //
00124   edm::ESHandle<TransientTrackingRecHitBuilder> theBuilder;
00125   std::string builderName = conf_.getParameter<std::string>("TTRHBuilder");   
00126   es.get<TransientRecHitRecord>().get(builderName,theBuilder);
00127   ttrhBuilder = theBuilder.product();
00128   
00129   edm::ESHandle<MeasurementTracker>    measurementTrackerHandle;
00130   es.get<CkfComponentsRecord>().get(measurementTrackerName_, measurementTrackerHandle);
00131   theMeasurementTracker = measurementTrackerHandle.product();
00132   
00133   std::vector<Trajectory> FinalTrajectories;
00134   
00135   
00136   // need this to sort recHits, sorting done after getting seed because propagationDirection is needed
00137   // get tracker geometry
00138   edm::ESHandle<TrackerGeometry> tracker;
00139   es.get<TrackerDigiGeometryRecord>().get(tracker);
00140   trackerGeom = tracker.product();
00141   
00142   edm::ESHandle<MagneticField> magField_;
00143   es.get<IdealMagneticFieldRecord>().get(magField_);
00144   magField = magField_.product();
00145   
00146   NoFieldCosmic_ = (CosmicReco_ && (magField->inTesla(GlobalPoint(0,0,0)).mag() < 0.01));
00147 
00148   theMeasurementTracker->update(e);
00149   //const MeasurementTracker*  theMeasurementTracker = new MeasurementTracker(es,mt_params); // will need this later
00150   
00151   theAloPropagator = new PropagatorWithMaterial(alongMomentum,.1057,&(*magField)); 
00152   theRevPropagator = new PropagatorWithMaterial(oppositeToMomentum,.1057,&(*magField)); 
00153   theAnalyticalPropagator = new AnalyticalPropagator(magField,anyDirection);
00154 
00155   thePropagator = theAloPropagator;
00156 
00157   //KFTrajectorySmoother theSmoother(*theRevPropagator, *theUpdator, *theEstimator);
00158   theSmoother = new KFTrajectorySmoother(*theRevPropagator, *theUpdator, *theEstimator);
00159   
00160   // get hit matcher
00161   theHitMatcher = new SiStripRecHitMatcher(3.0);
00162 
00163   //debug_ = true;
00164   //if (input->size()>0) debug_ = true;
00165 
00166   LogDebug("RoadSearch") << "Clean Clouds input size: " << input->size();
00167   if (debug_) std::cout << std::endl << std::endl
00168                         << "*** NEW EVENT: Clean Clouds input size: " << input->size() << std::endl;
00169   
00170   int i_c = 0;
00171   int nchit = 0;
00172   for ( RoadSearchCloudCollection::const_iterator cloud = input->begin(); cloud != input->end(); ++cloud ) {
00173     
00174     // fill rechits from cloud into new
00175     RoadSearchCloud::RecHitVector recHits = cloud->recHits();
00176     nchit = recHits.size();
00177     
00178     std::vector<Trajectory> CloudTrajectories;
00179     
00180     if (!CosmicReco_){
00181       std::sort(recHits.begin(),recHits.end(),SortHitPointersByGlobalPosition(tracker.product(),alongMomentum));
00182     }
00183     else {
00184       std::sort(recHits.begin(),recHits.end(),SortHitPointersByY(*tracker));
00185     }
00186 
00187     const uint nlost_max = 2;
00188             
00189     // make a list of layers in cloud and mark stereo layers
00190     
00191     const uint max_layers = 128;
00192 
00193     // collect hits in cloud by layer
00194     std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > > RecHitsByLayer;
00195     std::map<const DetLayer*, int> cloud_layer_reference; // for debugging
00196     std::map<const DetLayer*, int>::iterator hiter;
00197     for(RoadSearchCloud::RecHitVector::const_iterator ihit = recHits.begin();
00198         ihit != recHits.end(); ihit++) {
00199       // only use useful layers
00200       const DetLayer* thisLayer =
00201         theMeasurementTracker->geometricSearchTracker()->detLayer((*ihit)->geographicalId());
00202 
00203       std::map<const DetLayer*, int>::const_iterator ilyr = cloud_layer_reference.find(thisLayer);
00204       if (ilyr==cloud_layer_reference.end())
00205         cloud_layer_reference.insert(std::make_pair( thisLayer, RecHitsByLayer.size()));
00206 
00207       if (!RecHitsByLayer.empty() && RecHitsByLayer.back().first == thisLayer) { // Same as previous layer
00208         RecHitsByLayer.back().second.push_back(*ihit);
00209       }
00210       else { // check if this is a new layer
00211         if (ilyr != cloud_layer_reference.end()){// Not a New Layer
00212           int ilayer = ilyr->second;
00213           (RecHitsByLayer.begin()+ilayer)->second.push_back(*ihit);
00214         }
00215         else{// New Layer
00216           if (RecHitsByLayer.size() >= max_layers) break; // should never happen
00217           lstereo[RecHitsByLayer.size()] = false;
00218           if ((*ihit)->localPositionError().yy()<1.) lstereo[RecHitsByLayer.size()] = true;
00219           RoadSearchCloud::RecHitVector rhc;
00220           rhc.push_back(*ihit);
00221           RecHitsByLayer.push_back(std::make_pair(thisLayer, rhc));
00222         }       
00223       }
00224     }
00225 
00226     LogDebug("RoadSearch")<<"Cloud #"<<i_c<<" has "<<recHits.size()<<" hits in "<<RecHitsByLayer.size()<<" layers ";
00227     if (debug_) std::cout <<"Cloud "<<i_c<<" has "<<recHits.size()<<" hits in " <<RecHitsByLayer.size() << " layers " <<std::endl;;
00228     ++i_c;
00229 
00230     if (debug_){
00231       int ntothit = 0;
00232 
00233       for (std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >::iterator ilhv = RecHitsByLayer.begin();
00234            ilhv != RecHitsByLayer.end(); ++ilhv) {
00235         std::cout<<"   Layer " << ilhv-RecHitsByLayer.begin() << " has " << ilhv->second.size() << " hits " << std::endl;
00236       }
00237       std::cout<<std::endl;
00238       for (std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >::iterator ilhv = RecHitsByLayer.begin();
00239            ilhv != RecHitsByLayer.end(); ++ilhv) {
00240         RoadSearchCloud::RecHitVector theLayerHits = ilhv->second;
00241         for (RoadSearchCloud::RecHitVector::const_iterator ihit = theLayerHits.begin();
00242            ihit != theLayerHits.end(); ++ihit) {
00243         
00244           GlobalPoint gp = trackerGeom->idToDet((*ihit)->geographicalId())->surface().toGlobal((*ihit)->localPosition());
00245           if (CosmicReco_){
00246             std::cout << "   Hit "<< ntothit
00247                       << " x/y/z = "
00248                       << gp.x() << " " << gp.y() << " " << gp.z()
00249                       <<" in layer " << ilhv-RecHitsByLayer.begin()
00250                       << " is hit " << (ihit-theLayerHits.begin())+1 
00251                       << " of " << theLayerHits.size() << std::endl;
00252           }
00253           else {
00254             std::cout << "   Hit "<< ntothit
00255                       << " r/z = "
00256                       << gp.perp() << " " << gp.z()
00257                       <<" in layer " << ilhv-RecHitsByLayer.begin()
00258                       << " is hit " << (ihit-theLayerHits.begin())+1 
00259                       << " of " << theLayerHits.size() << std::endl;
00260           }
00261           ntothit++;
00262         }
00263       }
00264       std::cout<<std::endl;
00265     }
00266 
00267     // try to start from all layers until the chunk is too short
00268     //
00269     
00270     for (std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >::iterator ilyr0 = RecHitsByLayer.begin();
00271          ilyr0 != RecHitsByLayer.end(); ++ilyr0) {
00272 
00273       uint ilayer0 = (uint)(ilyr0-RecHitsByLayer.begin());
00274       if (ilayer0 > RecHitsByLayer.size()-MinChunkLength_) continue;      
00275 
00276       std::vector<Trajectory> ChunkTrajectories;
00277       std::vector<Trajectory> CleanChunks;
00278       bool all_chunk_layers_used = false;
00279       
00280       if (debug_) std::cout  << "*** START NEW CHUNK --> layer range (" << ilyr0-RecHitsByLayer.begin() 
00281                              << "-" << RecHitsByLayer.size()-1 << ")";
00282 
00283       // collect hits from the starting layer
00284       RoadSearchCloud::RecHitVector recHits_start = ilyr0->second;
00285 
00286       //
00287       // Step 1: find small tracks (chunks) made of hits
00288       // in layers with low occupancy
00289       //
00290       
00291       // find layers with small number of hits
00292       // TODO: try to keep earliest layers + at least one stereo layer
00293       std::multimap<int, const DetLayer*> layer_map;
00294       std::map<const DetLayer*, int> layer_reference; // for debugging
00295                                                       // skip starting layer, as it is always included
00296       for (std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >::iterator ilayer = ilyr0+1;
00297            ilayer != RecHitsByLayer.end(); ++ilayer) {
00298         layer_map.insert(std::make_pair(ilayer->second.size(), ilayer->first));
00299         layer_reference.insert(std::make_pair(ilayer->first, ilayer-RecHitsByLayer.begin()));
00300       }
00301       
00302       if (debug_) {
00303         std::cout<<std::endl<<"   Available layers are: " << std::endl;
00304         for (std::multimap<int, const DetLayer*>::iterator ilm1 = layer_map.begin();
00305              ilm1 != layer_map.end(); ++ilm1) {
00306           std::map<const DetLayer*, int>::iterator ilr = layer_reference.find(ilm1->second);
00307           if (ilr != layer_reference.end() && debug_) 
00308             std::cout << "Layer " << ilr->second << " with " << ilm1->first <<" hits" <<std::endl;;
00309         }
00310       }
00311 
00312       //const int max_middle_layers = 2;
00313       std::set<const DetLayer*> the_good_layers;
00314       std::vector<const DetLayer*> the_middle_layers;
00315       RoadSearchCloud::RecHitVector the_recHits_middle;
00316 
00317       bool StartLayers = chooseStartingLayers(RecHitsByLayer,ilyr0,layer_map,the_good_layers,the_middle_layers,the_recHits_middle);
00318       if (debug_) {
00319         std::cout << " From new code... With " << the_good_layers.size() << " useful layers: ";
00320         for (std::set<const DetLayer*>::iterator igl = the_good_layers.begin();
00321              igl!= the_good_layers.end(); ++igl){
00322           std::map<const DetLayer*, int>::iterator ilr = layer_reference.find(*igl);
00323           if (ilr != layer_reference.end()) std::cout << " " << ilr->second;
00324         }
00325         std::cout << std::endl;
00326         std::cout << " From new code... and middle layers: ";
00327         for (std::vector<const DetLayer*>::iterator iml = the_middle_layers.begin();
00328              iml!= the_middle_layers.end(); ++iml){
00329           std::map<const DetLayer*, int>::iterator ilr = layer_reference.find(*iml);
00330           if (ilr != layer_reference.end()) std::cout << " " << ilr->second;
00331         }
00332         std::cout << std::endl;
00333       }
00334       RoadSearchCloud::RecHitVector recHits_inner = recHits_start;
00335       RoadSearchCloud::RecHitVector recHits_outer = the_recHits_middle;
00336       std::set<const DetLayer*> good_layers = the_good_layers;
00337       uint ngoodlayers = good_layers.size();
00338 
00339       if (debug_)
00340         std::cout<<"Found " << recHits_inner.size() << " inner hits and " << recHits_outer.size() << " outer hits" << std::endl;
00341 
00342       // collect hits in useful layers
00343       std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > > goodHits;
00344       // mark layers that will be skipped in first pass
00345       std::set<const DetLayer*> skipped_layers;
00346       std::map<int, const DetLayer*> skipped_layer_detmap;
00347 
00348       
00349       goodHits.push_back(*ilyr0); // save hits from starting layer
00350       // save hits from other good layers
00351       for (std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >::iterator ilyr = ilyr0+1;
00352          ilyr != RecHitsByLayer.end(); ++ilyr) {
00353         if (good_layers.find(ilyr->first) != good_layers.end()){
00354           goodHits.push_back(*ilyr);
00355         }
00356         else {
00357           skipped_layers.insert(ilyr->first);
00358           std::map<const DetLayer*, int>::iterator ilr = layer_reference.find(ilyr->first);
00359           if (ilr != layer_reference.end())
00360             skipped_layer_detmap.insert(std::make_pair(ilr->second,ilyr->first));
00361           else
00362             if (debug_) std::cout<<"Couldn't find thisLayer to insert into map..."<<std::endl;
00363         }
00364       }
00365 
00366       // try various hit combinations
00367       for (RoadSearchCloud::RecHitVector::const_iterator innerHit = recHits_inner.begin();
00368            innerHit != recHits_inner.end(); ++innerHit) {
00369 
00370         const DetLayer* innerHitLayer =
00371           theMeasurementTracker->geometricSearchTracker()->detLayer((*innerHit)->geographicalId());
00372 
00373 
00374         RoadSearchCloud::RecHitVector::iterator middleHit, outerHit;
00375         RoadSearchCloud::RecHitVector::iterator firstHit, lastHit;
00376 
00377         bool triplets = (CosmicReco_ && (magField->inTesla(GlobalPoint(0,0,0)).mag() > 0.01));
00378         
00379         if (!triplets){
00380           firstHit = recHits_outer.begin();
00381           lastHit  = recHits_outer.end();
00382         }
00383         else if (triplets){
00384           firstHit = recHits_outer.begin()+1;
00385           lastHit = recHits_outer.end();
00386         }
00387 
00388         for (RoadSearchCloud::RecHitVector::iterator outerHit = firstHit; outerHit != lastHit; ++outerHit) {
00389           
00390           const DetLayer* middleHitLayer = 0;
00391           if (triplets){
00392             middleHit = outerHit-1;
00393             middleHitLayer = theMeasurementTracker->geometricSearchTracker()->detLayer((*middleHit)->geographicalId());
00394           }
00395           const DetLayer* outerHitLayer =
00396             theMeasurementTracker->geometricSearchTracker()->detLayer((*outerHit)->geographicalId());
00397           if (middleHitLayer == outerHitLayer) continue;
00398 
00399           FreeTrajectoryState fts;
00400           if (!triplets){
00401             if (debug_){
00402               std::map<const DetLayer*, int>::iterator ilro = layer_reference.find(outerHitLayer);
00403               if (ilro != layer_reference.end()) {
00404                 std::cout << "Try trajectory with Inner Hit on Layer " << ilayer0 << " and  " ;
00405                 std::cout << "Outer Hit on Layer " << ilro->second << std::endl;
00406               }
00407             }
00408             fts = initialTrajectory(es,*innerHit,*outerHit);
00409           }
00410           else if (triplets){
00411             if (debug_){
00412               std::map<const DetLayer*, int>::iterator ilrm = layer_reference.find(middleHitLayer);
00413               std::map<const DetLayer*, int>::iterator ilro = layer_reference.find(outerHitLayer);
00414               if (ilro != layer_reference.end() && ilrm != layer_reference.end()) {
00415                 std::cout << "Try trajectory with Hits on Layers " << ilayer0 << " , "
00416                           << ilrm->second <<  " and  " << ilro->second << std::endl;
00417               }
00418             }
00419             fts = initialTrajectoryFromTriplet(es,*innerHit,*middleHit,*outerHit);
00420           }
00421 
00422           if (!fts.hasError()) continue;
00423           if (debug_) std::cout<<"FTS: " << fts << std::endl;
00424 
00425           Trajectory seedTraj = createSeedTrajectory(fts,*innerHit,innerHitLayer);
00426 
00427           std::vector<Trajectory> rawTrajectories;          
00428           if (seedTraj.isValid() && !seedTraj.measurements().empty() ) rawTrajectories.push_back(seedTraj);//GC
00429           //rawTrajectories.push_back(seedTraj);
00430 
00431           int ntested = 0;
00432           // now loop on hits
00433           std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >::iterator ilyr_start = (goodHits.begin()+1);
00434           for (std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >::iterator ilhv = ilyr_start;
00435                ilhv != goodHits.end(); ++ilhv) {
00436             RoadSearchCloud::RecHitVector& hits = ilhv->second;
00437             //std::vector<Trajectory> newTrajectories;
00438             ++ntested;
00439             if (debug_){
00440               std::map<const DetLayer*, int>::iterator ilr = cloud_layer_reference.find(ilhv->first);
00441               if (ilr != cloud_layer_reference.end())
00442                 std::cout << "extrapolating " << rawTrajectories.size() 
00443                           << " trajectories to layer " << ilr->second 
00444                           << " which has  " << hits.size() << " hits " << std::endl;
00445             }
00446             
00447             std::vector<Trajectory>newTrajectories;
00448             for (std::vector<Trajectory>::const_iterator it = rawTrajectories.begin();
00449                      it != rawTrajectories.end(); it++) {
00450               if (debug_) std::cout << "extrapolating Trajectory #" << it-rawTrajectories.begin() << std::endl;
00451               if (it->direction()==alongMomentum) thePropagator = theAloPropagator;//GC
00452               else thePropagator = theRevPropagator;
00453 
00454               std::vector<Trajectory> theTrajectories = extrapolateTrajectory(*it,hits,
00455                                                                               innerHitLayer, *outerHit, outerHitLayer);
00456               if (theTrajectories.empty()) {
00457                 if (debug_) std::cout<<" Could not add the hit in this layer " << std::endl;
00458                 if (debug_){
00459                   std::cout << " --> trajectory " << it-rawTrajectories.begin() 
00460                             << " has "<<it->recHits().size()<<" hits after "
00461                             << (ilhv-ilyr_start+1) << " tested (ntested=" <<ntested<<") "
00462                             << " --> misses="<< (ilhv-ilyr_start+1)-(it->recHits().size()-1)
00463                             << " but there are " << (goodHits.end() - ilhv)
00464                             <<" more layers in first pass and "<< skipped_layers.size() <<" skipped layers " <<std::endl;
00465                   
00466                   
00467                 }
00468                 // layer missed
00469                 if ((ilhv-ilyr_start+1)-(it->recHits().size()-1) <= nlost_max){
00470                   newTrajectories.push_back(*it);
00471                 }
00472               }
00473               else{ // added hits in this layers
00474                 for (std::vector<Trajectory>::const_iterator it = theTrajectories.begin();
00475                      it != theTrajectories.end(); it++) {
00476                   newTrajectories.push_back(*it);
00477                 }
00478               }
00479             } // end loop over rawTrajectories
00480             rawTrajectories = newTrajectories;
00481             if (newTrajectories.empty()) break;
00482           }
00483           if (rawTrajectories.size()==0){
00484             continue;
00485             if (debug_) std::cout<<" --> yields ZERO raw trajectories!" << std::endl;
00486           }
00487           if (debug_){
00488             for (std::vector<Trajectory>::const_iterator it = rawTrajectories.begin();
00489                  it != rawTrajectories.end(); it++) {
00490               std::cout << " --> yields trajectory with "<<it->recHits().size()<<" hits with chi2="
00491                         <<it->chiSquared()<<" and is valid? "<<it->isValid() <<std::endl;
00492             }
00493           }
00494           std::vector<Trajectory> rawCleaned;
00495           theTrajectoryCleaner->clean(rawTrajectories);
00496           for (std::vector<Trajectory>::const_iterator itr = rawTrajectories.begin();
00497                itr != rawTrajectories.end(); ++itr) {
00498             // see how many layers have been found
00499             if (!itr->isValid()) continue;
00500             std::set<const DetLayer*> used_layers;
00501             Trajectory::DataContainer tmv = itr->measurements();
00502             for (Trajectory::DataContainer::iterator itm = tmv.begin();
00503                  itm != tmv.end(); ++itm) {
00504               TransientTrackingRecHit::ConstRecHitPointer rh = itm->recHit();
00505               if (!rh->isValid()) continue;
00506               used_layers.insert(theMeasurementTracker->geometricSearchTracker()->detLayer(rh->geographicalId()));
00507             }
00508 
00509             // need to subtract 1 from used_layers since it includes the starting layer
00510             if (debug_) std::cout<<"Used " << (used_layers.size()-1) << " layers out of " << ngoodlayers
00511                                  << " good layers, so " << ngoodlayers - (used_layers.size()-1) << " missed "
00512                                  << std::endl;
00513             if ((int)used_layers.size() < nFoundMin_) continue;
00514             uint nlostlayers = ngoodlayers - (used_layers.size()-1);
00515             if (nlostlayers > nlost_max) continue;
00516             
00517             rawCleaned.push_back(*itr);
00518             
00519           }
00520           if (!rawCleaned.empty()) {
00521             ChunkTrajectories.insert(ChunkTrajectories.end(), rawCleaned.begin(), rawCleaned.end());
00522           }
00523         } // END LOOP OVER OUTER HITS
00524       } // END LOOP OVER INNER HITS
00525       // At this point we have made all the trajectories from the low occupancy layers
00526       // We clean these trajectories first, and then try to add hits from the skipped layers
00527       
00528       //    }
00529       if (debug_) std::cout << "Clean the " << ChunkTrajectories.size()<<" trajectories for this chunk" << std::endl;
00530       // clean the intermediate result
00531       theTrajectoryCleaner->clean(ChunkTrajectories);
00532       for (std::vector<Trajectory>::const_iterator it = ChunkTrajectories.begin();
00533            it != ChunkTrajectories.end(); it++) {
00534         if (it->isValid())  CleanChunks.push_back(*it);
00535       }
00536       if (debug_) std::cout <<"After cleaning there are " << CleanChunks.size() << " trajectories for this chunk" << std::endl;
00537 
00538       
00539       // *********************  BEGIN NEW ADDITION      
00540       
00541       //
00542       // Step 2: recover measurements from busy layers
00543       //
00544       
00545       std::vector<Trajectory> extendedChunks;
00546       
00547 
00548       // see if there are layers that we skipped
00549 
00550       if (debug_){
00551         if (skipped_layers.empty()) {
00552           std::cout << "all layers were used in first pass" << std::endl;
00553         } else {
00554           std::cout << "There are " << skipped_layer_detmap.size() << " skipped layers:";
00555           for (std::map<int, const DetLayer*>::const_iterator imap = skipped_layer_detmap.begin();
00556                imap!=skipped_layer_detmap.end(); imap++){
00557             std::cout<< " " <<imap->first;
00558           }
00559           std::cout << std::endl;
00560         }
00561       }
00562       
00563       for (std::vector<Trajectory>::const_iterator i = CleanChunks.begin();
00564            i != CleanChunks.end(); i++) {
00565         if (!(*i).isValid()) continue;
00566         if (debug_) std::cout<< "Now process CleanChunk trajectory " << i-CleanChunks.begin() << std::endl;
00567         bool all_layers_used = false;
00568         if (skipped_layers.empty() && i->measurements().size() >= theNumHitCut) {
00569           if (debug_) std::cout<<"The trajectory has " << i->measurements().size() 
00570                                << " hits from a cloud of " << RecHitsByLayer.size() 
00571                                << " layers and a chunk of " << (RecHitsByLayer.size() - ilayer0) << " layers " << std::endl;
00572           extendedChunks.insert(extendedChunks.end(), *i);
00573           if (i->measurements().size() >= (RecHitsByLayer.size() - ilayer0)){
00574             all_layers_used = true;
00575             break;
00576           }
00577         } 
00578         else {
00579           
00580           Trajectory temptraj = *i;
00581           Trajectory::DataContainer tmv = (*i).measurements();
00582           if (tmv.size()+skipped_layer_detmap.size() < theNumHitCut) continue;          
00583 
00584           // Debug dump of hits
00585           if (debug_){
00586             for (Trajectory::DataContainer::const_iterator ih=tmv.begin();
00587                  ih!=tmv.end();++ih){
00588               TransientTrackingRecHit::ConstRecHitPointer rh = ih->recHit();
00589               const DetLayer* Layer =
00590                 theMeasurementTracker->geometricSearchTracker()->detLayer(rh->geographicalId());      
00591               std::map<const DetLayer*, int>::iterator ilr = cloud_layer_reference.find(Layer);
00592               if (ilr != cloud_layer_reference.end())
00593                 std::cout << "   Hit #"<<ih-tmv.begin() << " of " << tmv.size()
00594                           <<" is on Layer " << ilr->second << std::endl;
00595               else 
00596                 std::cout << "   Layer for Hit #"<<ih-tmv.begin() <<" can't be found " << std::endl;
00597               std::cout<<" updatedState:\n" << ih->updatedState() << std::endl;
00598               std::cout<<" predictState:\n" << ih->predictedState() << std::endl;
00599             }
00600           }
00601           
00602           // Loop over the layers in the cloud
00603           
00604           std::set<const DetLayer*> final_layers;
00605           Trajectory::DataContainer::const_iterator im = tmv.begin();
00606           Trajectory::DataContainer::const_iterator im2 = tmv.begin();
00607           
00608           TrajectoryMeasurement firstMeasurement = i->firstMeasurement();
00609           const DetLayer* firstDetLayer = 
00610             theMeasurementTracker->geometricSearchTracker()->detLayer(firstMeasurement.recHit()->geographicalId());
00611           
00612           std::vector<Trajectory> freshStartv = theSmoother->trajectories(*i);
00613           
00614           if(debug_) std::cout<<"Smoothing has returned " << freshStartv.size() <<" trajectory " << std::endl;
00615           if (!freshStartv.empty()){
00616              if(debug_)  std::cout<< "Smoothing of trajectory " <<i-CleanChunks.begin() << " has succeeded with " << freshStartv.begin()->measurements().size() << " hits and a chi2 of " << freshStartv.begin()->chiSquared() <<" for " << freshStartv.begin()->ndof() << " DOF.  Now add hits." <<std::endl;
00617           } else {
00618              if (debug_) std::cout<< "Smoothing of trajectory " <<i-CleanChunks.begin() <<" has failed"<<std::endl;
00619              continue;
00620           }
00621           
00622           Trajectory freshStart = *freshStartv.begin();
00623           std::vector<TrajectoryMeasurement> freshStartTM = freshStart.measurements();
00624           
00625           if (debug_) {
00626              for (std::vector<TrajectoryMeasurement>::const_iterator itm = freshStartTM.begin();itm != freshStartTM.end(); ++itm){
00627                 std::cout<<"Trajectory hit " << itm-freshStartTM.begin() <<" updatedState:\n" << itm->updatedState() << std::endl;
00628              }
00629           }
00630 
00631           TrajectoryStateOnSurface NewFirstTsos = freshStart.lastMeasurement().updatedState();
00632           if(debug_) std::cout<<"NewFirstTSOS is valid? " << NewFirstTsos.isValid() << std::endl;
00633           if(debug_) std::cout << " NewFirstTSOS:\n " << NewFirstTsos << std::endl;
00634           TransientTrackingRecHit::ConstRecHitPointer rh = freshStart.lastMeasurement().recHit();
00635 
00636           if(debug_) {
00637              std::cout<< "First hit for fresh start on det " << rh->det() << ", r/phi/z = " << rh->globalPosition().perp() << " " << rh->globalPosition().phi() << " " << rh->globalPosition().z();
00638           }
00639           
00640           PTrajectoryStateOnDet* pFirstState = TrajectoryStateTransform().persistentState(NewFirstTsos,
00641                                                                                           rh->geographicalId().rawId());
00642           edm::OwnVector<TrackingRecHit> newHits;
00643           newHits.push_back(rh->hit()->clone());
00644           
00645           TrajectorySeed tmpseed = TrajectorySeed(*pFirstState, 
00646                                                   newHits,
00647                                                   i->direction());
00648 
00649           thePropagator = theAloPropagator;
00650           if (i->direction()==oppositeToMomentum) thePropagator = theRevPropagator;
00651 
00652           delete pFirstState;
00653           
00654           
00655           Trajectory newTrajectory(tmpseed,tmpseed.direction());
00656           
00657           const GeomDet* det = trackerGeom->idToDet(rh->geographicalId());
00658           TrajectoryStateOnSurface invalidState(new BasicSingleTrajectoryState(det->surface()));
00659           newTrajectory.push(TrajectoryMeasurement(invalidState, NewFirstTsos, rh, 0, firstDetLayer));
00660           final_layers.insert(firstDetLayer);
00661 
00662           if(debug_) std::cout << "TRAJ is valid: " << newTrajectory.isValid() <<std::endl;
00663 
00664           TrajectoryStateOnSurface testTsos = newTrajectory.measurements().back().updatedState();
00665           
00666           if(debug_) {
00667             std::cout << "testTSOS is valid!!" << testTsos.isValid() << std::endl;
00668             std::cout << " testTSOS (x/y/z): " << testTsos.globalPosition().x()<< " / " << testTsos.globalPosition().y()<< " / " << testTsos.globalPosition().z() << std::endl;
00669             std::cout << " local position x: " << testTsos.localPosition().x() << "+-" << sqrt(testTsos.localError().positionError().xx()) << std::endl;
00670           }
00671 
00672           if (firstDetLayer != ilyr0->first){
00673             if (debug_) std::cout<<"!!! ERROR !!! firstDetLayer ("<<firstDetLayer<<") != ilyr0 ( " <<ilyr0->first <<")"<< std::endl;
00674           }          
00675           ++im;
00676 
00677           if (debug_){
00678             std::map<const DetLayer*, int>::iterator ilr = cloud_layer_reference.find(firstDetLayer);
00679             if (ilr != cloud_layer_reference.end() ){
00680               std::cout << "   First hit is on layer  " << ilr->second << std::endl;
00681             }
00682           }
00683           for (std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >::iterator ilyr = ilyr0+1;
00684                ilyr != RecHitsByLayer.end(); ++ilyr) {
00685 
00686             TrajectoryStateOnSurface predTsos;
00687             TrajectoryStateOnSurface currTsos;
00688             TrajectoryMeasurement tm;
00689 
00690             if(debug_)
00691               std::cout<<"Trajectory has " << newTrajectory.measurements().size() << " measurements with " << (RecHitsByLayer.end()-ilyr)
00692                        << " remaining layers " << std::endl;
00693 
00694             if (im != tmv.end()) im2 = im;
00695             TransientTrackingRecHit::ConstRecHitPointer rh = im2->recHit(); // Current 
00696             if (rh->isValid() && 
00697                 theMeasurementTracker->geometricSearchTracker()->detLayer(rh->geographicalId()) == ilyr->first) {
00698               
00699                 if (debug_) std::cout<<"   Layer " << ilyr-RecHitsByLayer.begin() <<" has a good hit " << std::endl;
00700                 ++im;
00701                 
00702                 const GeomDet* det = trackerGeom->idToDet(rh->geographicalId());
00703 
00704                 currTsos = newTrajectory.measurements().back().updatedState();
00705                 predTsos = thePropagator->propagate(currTsos, det->surface());
00706                 if (!predTsos.isValid()) continue;
00707                 GlobalVector propagationDistance = predTsos.globalPosition() - currTsos.globalPosition();
00708                 if (propagationDistance.mag() > maxPropagationDistance) continue;
00709                 MeasurementEstimator::HitReturnType est = theEstimator->estimate(predTsos, *rh);
00710                 if(debug_) {
00711                   std::cout << "Propagation distance2 is " << propagationDistance.mag() << std::endl;
00712                   std::cout << "predTSOS is valid!!" << std::endl;
00713                   std::cout << " predTSOS (x/y/z): " << predTsos.globalPosition().x()<< " / " << predTsos.globalPosition().y()<< " / " << predTsos.globalPosition().z() << std::endl;
00714                   std::cout << " local position x: " << predTsos.localPosition().x() << "+-" << sqrt(predTsos.localError().positionError().xx()) << std::endl;
00715                   std::cout << " local position y: " << predTsos.localPosition().y() << "+-" << sqrt(predTsos.localError().positionError().yy()) << std::endl;
00716                   std::cout << "currTSOS is valid!! " << currTsos.isValid() <<  std::endl;
00717                   std::cout << " currTSOS (x/y/z): " << currTsos.globalPosition().x()<< " / " << currTsos.globalPosition().y()<< " / " << currTsos.globalPosition().z() << std::endl;
00718                   std::cout << " local position x: " << currTsos.localPosition().x() << "+-" << sqrt(currTsos.localError().positionError().xx()) << std::endl;
00719                   std::cout << " local position y: " << currTsos.localPosition().y() << "+-" << sqrt(currTsos.localError().positionError().yy()) << std::endl;
00720                 }
00721 
00722                 if (!est.first) {
00723                   if (debug_) std::cout<<"Failed to add one of the original hits on a low occupancy layer!!!!" << std::endl;
00724                   continue;
00725                 }
00726                 currTsos = theUpdator->update(predTsos, *rh);
00727                 tm = TrajectoryMeasurement(predTsos, currTsos, &(*rh),est.second,ilyr->first);
00728 
00729                 const TrajectoryStateOnSurface theTSOS = newTrajectory.lastMeasurement().updatedState();
00730                 
00731                 //std::cout << "11TsosBefore (x/y/z): " << theTSOS.globalPosition().x()<< " / " << theTSOS.globalPosition().y()<< " / " << theTSOS.globalPosition().z() << std::endl;
00732                 //std::cout << " 11local position x: " << theTSOS.localPosition().x() << "+-" << sqrt(theTSOS.localError().positionError().xx()) << std::endl;
00733                 //std::cout << " 11local position y: " << theTSOS.localPosition().y() << "+-" << sqrt(theTSOS.localError().positionError().yy()) << std::endl;
00734 
00735 
00736 
00737                 newTrajectory.push(tm,est.second);
00738                 final_layers.insert(ilyr->first);                
00739             }       
00740             else{
00741               if (debug_) std::cout<<"   Layer " << ilyr-RecHitsByLayer.begin() <<" is one of the skipped layers " << std::endl;
00742               
00743               //collect hits in the skipped layer
00744               edm::OwnVector<TrackingRecHit> skipped_hits;
00745               std::set<const GeomDet*> dets;
00746               for (RoadSearchCloud::RecHitVector::const_iterator ih = ilyr->second.begin();
00747                    ih != ilyr->second.end(); ++ih) {
00748                 skipped_hits.push_back((*ih)->clone());
00749                 dets.insert(trackerGeom->idToDet((*ih)->geographicalId()));
00750               }
00751               
00752               if(debug_){
00753                 std::cout<<"   ---> probing missing hits (nh="<< skipped_hits.size() << ", nd=" << dets.size() 
00754                          << ")  in layer " <<  ilyr-RecHitsByLayer.begin() <<std::endl;
00755               }
00756 
00757               const TrajectoryStateOnSurface theTSOS = newTrajectory.lastMeasurement().updatedState();
00758               
00759               //std::cout << "TsosBefore (x/y/z): " << theTSOS.globalPosition().x()<< " / " << theTSOS.globalPosition().y()<< " / " << theTSOS.globalPosition().z() << std::endl;
00760               //std::cout << " local position x: " << theTSOS.localPosition().x() << "+-" << sqrt(theTSOS.localError().positionError().xx()) << std::endl;
00761               //std::cout << " local position y: " << theTSOS.localPosition().y() << "+-" << sqrt(theTSOS.localError().positionError().yy()) << std::endl;
00762 
00763               std::vector<TrajectoryMeasurement> theGoodHits = FindBestHits(theTSOS,dets,theHitMatcher,skipped_hits);
00764               if (!theGoodHits.empty()){
00765                 final_layers.insert(ilyr->first);
00766                 if (debug_) std::cout<<"Found " << theGoodHits.size() << " good hits to add" << std::endl;
00767                 for (std::vector<TrajectoryMeasurement>::const_iterator im=theGoodHits.begin();im!=theGoodHits.end();++im){
00768                   newTrajectory.push(*im,im->estimate());
00769                 }
00770               }
00771             }
00772           } // END 2nd PASS LOOP OVER LAYERS IN CLOUD
00773           
00774           if (debug_) std::cout<<"Finished loop over layers in cloud.  Trajectory now has " <<newTrajectory.measurements().size()
00775                                << " hits. " << std::endl;
00776           if (debug_) std::cout<<"The trajectory has " << newTrajectory.measurements().size() <<" hits on " << final_layers.size()
00777                                << " layers from a cloud of " << RecHitsByLayer.size() 
00778                                << " layers and a chunk of " << (RecHitsByLayer.size() - ilayer0) << " layers " << std::endl;
00779           if (newTrajectory.measurements().size() >= theNumHitCut)
00780             extendedChunks.insert(extendedChunks.end(), newTrajectory);
00781           if (final_layers.size() >= (RecHitsByLayer.size() - ilayer0)){
00782             if (debug_) std::cout<<"All layers of the chunk have been used..." << std::endl;
00783             all_layers_used = true;
00784           }
00785         }  // END ELSE TO RECOVER SKIPPED LAYERS
00786         if (all_layers_used) {
00787           if (debug_) std::cout << "All layers were used, so break....." << std::endl;
00788           all_chunk_layers_used = true;
00789           break;
00790         }
00791         if (debug_) std::cout<<"Going to next clean chunk....." << std::endl;
00792       } // END LOOP OVER CLEAN CHUNKS
00793       if (debug_) std::cout<< "Now Clean the " << extendedChunks.size() << " extended chunks " <<std::endl;
00794       if (extendedChunks.size() > 1) theTrajectoryCleaner->clean(extendedChunks);
00795       for (std::vector<Trajectory>::const_iterator it = extendedChunks.begin();
00796            it != extendedChunks.end(); it++) {
00797         if (it->isValid()) CloudTrajectories.push_back(*it);
00798       }
00799       if (all_chunk_layers_used) break;
00800     }
00801     
00802     // ********************* END NEW ADDITION
00803     
00804     if (debug_) std::cout<< "Finished with the cloud, so clean the " 
00805                          << CloudTrajectories.size() << " cloud trajectories "<<std::endl ;
00806     theTrajectoryCleaner->clean(CloudTrajectories);
00807     for (std::vector<Trajectory>::const_iterator it = CloudTrajectories.begin();
00808          it != CloudTrajectories.end(); it++) {
00809       if (it->isValid()) FinalTrajectories.push_back(*it);
00810     }
00811     
00812     if (debug_) std::cout<<" Final trajectories now has size " << FinalTrajectories.size()<<std::endl ;
00813     
00814   } // End loop over Cloud Collection
00815 
00816 
00817   if (debug_) std::cout<< " Finished loop over all clouds " <<std::endl;
00818 
00819   output = PrepareTrackCandidates(FinalTrajectories);
00820 
00821   delete theAloPropagator;
00822   delete theRevPropagator; 
00823   delete theAnalyticalPropagator;
00824   delete theHitMatcher;
00825   delete theSmoother;
00826   
00827   if (debug_ || debugCosmics_) std::cout<< "Found " << output.size() << " track candidate(s)."<<std::endl;
00828 
00829 }

edm::ParameterSet RoadSearchTrackCandidateMakerAlgorithm::conf_ [private]

Definition at line 113 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by RoadSearchTrackCandidateMakerAlgorithm(), and run().

bool RoadSearchTrackCandidateMakerAlgorithm::CosmicReco_ [private]

Definition at line 117 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by chooseStartingLayers(), initialTrajectory(), initialTrajectoryFromTriplet(), RoadSearchTrackCandidateMakerAlgorithm(), and run().

double RoadSearchTrackCandidateMakerAlgorithm::cosmicSeedPt_ [private]

Definition at line 125 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by initialTrajectory(), and RoadSearchTrackCandidateMakerAlgorithm().

bool RoadSearchTrackCandidateMakerAlgorithm::CosmicTrackMerging_ [private]

Definition at line 119 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by PrepareTrackCandidates(), and RoadSearchTrackCandidateMakerAlgorithm().

bool RoadSearchTrackCandidateMakerAlgorithm::debug_ [private]

Definition at line 131 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by createSeedTrajectory(), extrapolateTrajectory(), FindBestHit(), FindBestHits(), initialTrajectory(), initialTrajectoryFromTriplet(), PrepareTrackCandidates(), RoadSearchTrackCandidateMakerAlgorithm(), and run().

bool RoadSearchTrackCandidateMakerAlgorithm::debugCosmics_ [private]

Definition at line 132 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by PrepareTrackCandidates(), RoadSearchTrackCandidateMakerAlgorithm(), and run().

double RoadSearchTrackCandidateMakerAlgorithm::initialVertexErrorXY_ [private]

Definition at line 123 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by initialTrajectory(), initialTrajectoryFromTriplet(), and RoadSearchTrackCandidateMakerAlgorithm().

bool RoadSearchTrackCandidateMakerAlgorithm::lstereo[128] [private]

Definition at line 151 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by chooseStartingLayers(), and run().

const MagneticField* RoadSearchTrackCandidateMakerAlgorithm::magField [private]

Definition at line 137 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by initialTrajectory(), initialTrajectoryFromTriplet(), and run().

double RoadSearchTrackCandidateMakerAlgorithm::maxPropagationDistance [private]

Definition at line 127 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by extrapolateTrajectory(), FindBestHit(), FindBestHits(), FindBestHitsByDet(), PrepareTrackCandidates(), RoadSearchTrackCandidateMakerAlgorithm(), and run().

std::string RoadSearchTrackCandidateMakerAlgorithm::measurementTrackerName_ [private]

Definition at line 129 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by RoadSearchTrackCandidateMakerAlgorithm(), and run().

int RoadSearchTrackCandidateMakerAlgorithm::MinChunkLength_ [private]

Definition at line 120 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by RoadSearchTrackCandidateMakerAlgorithm(), and run().

int RoadSearchTrackCandidateMakerAlgorithm::nFoundMin_ [private]

Definition at line 121 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by chooseStartingLayers(), RoadSearchTrackCandidateMakerAlgorithm(), and run().

bool RoadSearchTrackCandidateMakerAlgorithm::NoFieldCosmic_ [private]

Definition at line 118 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by initialTrajectory(), and run().

bool RoadSearchTrackCandidateMakerAlgorithm::splitMatchedHits_ [private]

Definition at line 124 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by PrepareTrackCandidates(), and RoadSearchTrackCandidateMakerAlgorithm().

PropagatorWithMaterial* RoadSearchTrackCandidateMakerAlgorithm::theAloPropagator [private]

Definition at line 140 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by initialTrajectory(), initialTrajectoryFromTriplet(), PrepareTrackCandidates(), and run().

AnalyticalPropagator* RoadSearchTrackCandidateMakerAlgorithm::theAnalyticalPropagator [private]

Definition at line 142 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by PrepareTrackCandidates(), and run().

double RoadSearchTrackCandidateMakerAlgorithm::theChi2Cut [private]

Definition at line 116 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by RoadSearchTrackCandidateMakerAlgorithm().

MeasurementEstimator* RoadSearchTrackCandidateMakerAlgorithm::theEstimator [private]

Definition at line 144 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by createSeedTrajectory(), extrapolateTrajectory(), FindBestHit(), FindBestHits(), FindBestHitsByDet(), PrepareTrackCandidates(), RoadSearchTrackCandidateMakerAlgorithm(), run(), and ~RoadSearchTrackCandidateMakerAlgorithm().

SiStripRecHitMatcher* RoadSearchTrackCandidateMakerAlgorithm::theHitMatcher [private]

Definition at line 145 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by createSeedTrajectory(), extrapolateTrajectory(), and run().

const MeasurementTracker* RoadSearchTrackCandidateMakerAlgorithm::theMeasurementTracker [private]

Definition at line 134 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by extrapolateTrajectory(), FindBestHit(), FindBestHits(), FindBestHitsByDet(), PrepareTrackCandidates(), and run().

unsigned int RoadSearchTrackCandidateMakerAlgorithm::theNumHitCut [private]

Definition at line 115 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by RoadSearchTrackCandidateMakerAlgorithm(), and run().

PropagatorWithMaterial* RoadSearchTrackCandidateMakerAlgorithm::thePropagator [private]

Definition at line 139 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by createSeedTrajectory(), extrapolateTrajectory(), FindBestHit(), FindBestHits(), FindBestHitsByDet(), initialTrajectory(), initialTrajectoryFromTriplet(), and run().

PropagatorWithMaterial* RoadSearchTrackCandidateMakerAlgorithm::theRevPropagator [private]

Definition at line 141 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by initialTrajectory(), initialTrajectoryFromTriplet(), and run().

KFTrajectorySmoother* RoadSearchTrackCandidateMakerAlgorithm::theSmoother [private]

Definition at line 146 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by PrepareTrackCandidates(), and run().

TrajectoryCleanerBySharedHits* RoadSearchTrackCandidateMakerAlgorithm::theTrajectoryCleaner [private]

Definition at line 149 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by PrepareTrackCandidates(), RoadSearchTrackCandidateMakerAlgorithm(), run(), and ~RoadSearchTrackCandidateMakerAlgorithm().

TrajectoryStateTransform* RoadSearchTrackCandidateMakerAlgorithm::theTransformer [private]

Definition at line 148 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by createSeedTrajectory(), PrepareTrackCandidates(), RoadSearchTrackCandidateMakerAlgorithm(), and ~RoadSearchTrackCandidateMakerAlgorithm().

TrajectoryStateUpdator* RoadSearchTrackCandidateMakerAlgorithm::theUpdator [private]

Definition at line 143 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by createSeedTrajectory(), extrapolateTrajectory(), FindBestHit(), FindBestHits(), FindBestHitsByDet(), PrepareTrackCandidates(), RoadSearchTrackCandidateMakerAlgorithm(), run(), and ~RoadSearchTrackCandidateMakerAlgorithm().

const TrackerGeometry* RoadSearchTrackCandidateMakerAlgorithm::trackerGeom [private]

Definition at line 135 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by createSeedTrajectory(), extrapolateTrajectory(), FindBestHit(), FindBestHits(), FindBestHitsByDet(), initialTrajectory(), initialTrajectoryFromTriplet(), PrepareTrackCandidates(), and run().

const TransientTrackingRecHitBuilder* RoadSearchTrackCandidateMakerAlgorithm::ttrhBuilder [private]

Definition at line 136 of file RoadSearchTrackCandidateMakerAlgorithm.h.

Referenced by createSeedTrajectory(), extrapolateTrajectory(), FindBestHit(), FindBestHits(), FindBestHitsByDet(), and run().

RoadSearchTrackCandidateMakerAlgorithm Class Reference

Public Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


Public Member Functions
bool	chooseStartingLayers (std::vector< std::pair< const DetLayer , RoadSearchCloud::RecHitVector > > &RecHitsByLayer, std::vector< std::pair< const DetLayer , RoadSearchCloud::RecHitVector > >::iterator ilyr0, const std::multimap< int, const DetLayer * > &layer_map, std::set< const DetLayer * > &good_layers, std::vector< const DetLayer * > &middle_layers, RoadSearchCloud::RecHitVector &recHits_middle)
Trajectory	createSeedTrajectory (FreeTrajectoryState &fts, const TrackingRecHit InnerHit, const DetLayer innerHitLayer)
std::vector< Trajectory >	extrapolateTrajectory (const Trajectory &inputTrajectory, RoadSearchCloud::RecHitVector &theLayerHits, const DetLayer innerHitLayer, const TrackingRecHit outerHit, const DetLayer *outerHitLayer)
std::vector < TrajectoryMeasurement >	FindBestHit (const TrajectoryStateOnSurface &tsosBefore, const std::set< const GeomDet * > &theDets, edm::OwnVector< TrackingRecHit > &theHits)
std::vector < TrajectoryMeasurement >	FindBestHits (const TrajectoryStateOnSurface &tsosBefore, const std::set< const GeomDet * > &theDets, const SiStripRecHitMatcher *theHitMatcher, edm::OwnVector< TrackingRecHit > &theHits)
std::vector < TrajectoryMeasurement >	FindBestHitsByDet (const TrajectoryStateOnSurface &tsosBefore, const std::set< const GeomDet * > &theDets, edm::OwnVector< TrackingRecHit > &theHits)
FreeTrajectoryState	initialTrajectory (const edm::EventSetup &es, const TrackingRecHit InnerHit, const TrackingRecHit OuterHit)
FreeTrajectoryState	initialTrajectoryFromTriplet (const edm::EventSetup &es, const TrackingRecHit InnerHit, const TrackingRecHit MiddleHit, const TrackingRecHit *OuterHit)
TrackCandidateCollection	PrepareTrackCandidates (std::vector< Trajectory > &theTrajectories)
	RoadSearchTrackCandidateMakerAlgorithm (const edm::ParameterSet &conf)
void	run (const RoadSearchCloudCollection *input, const edm::Event &e, const edm::EventSetup &es, TrackCandidateCollection &output)
	Runs the algorithm.
	~RoadSearchTrackCandidateMakerAlgorithm ()
Private Attributes
edm::ParameterSet	conf_
bool	CosmicReco_
double	cosmicSeedPt_
bool	CosmicTrackMerging_
bool	debug_
bool	debugCosmics_
double	initialVertexErrorXY_
bool	lstereo [128]
const MagneticField *	magField
double	maxPropagationDistance
std::string	measurementTrackerName_
int	MinChunkLength_
int	nFoundMin_
bool	NoFieldCosmic_
bool	splitMatchedHits_
PropagatorWithMaterial *	theAloPropagator
AnalyticalPropagator *	theAnalyticalPropagator
double	theChi2Cut
MeasurementEstimator *	theEstimator
SiStripRecHitMatcher *	theHitMatcher
const MeasurementTracker *	theMeasurementTracker
unsigned int	theNumHitCut
PropagatorWithMaterial *	thePropagator
PropagatorWithMaterial *	theRevPropagator
KFTrajectorySmoother *	theSmoother
TrajectoryCleanerBySharedHits *	theTrajectoryCleaner
TrajectoryStateTransform *	theTransformer
TrajectoryStateUpdator *	theUpdator
const TrackerGeometry *	trackerGeom
const TransientTrackingRecHitBuilder *	ttrhBuilder