New Track Event Model: TrueTracksCreator Class Reference

00040   : GaudiAlgorithm( name, pSvcLocator )
00041   , m_otTim2MCHit(0)
00042   , m_itClus2MCP(0)
00043   , m_veloClus2MCP(0)
00044   , m_otTracker(0)
00045   , m_itTracker(0)
00046   , m_velo(0)
00047   , m_trackSelector(0)
00048   , m_stateCreator(0)
00049 //  , m_tracksFitter(0)
00050 //  , m_veloTracksFitter(0)
00051 //  , m_seedTracksFitter(0)
00052 {
00054   declareProperty( "AddOTTimes",      m_addOTTimes = true );
00055   declareProperty( "AddITClusters",   m_addITClusters = true );
00056   declareProperty( "AddVeloClusters", m_addVeloClusters = true );
00057   declareProperty( "InitState",       m_initState = true );
00058   declareProperty( "FitTracks",       m_fitTracks = true );
00059   declareProperty( "FitUpstream",     m_upstream = true );
00060   declareProperty( "TrueStatesAtMeasZPos", m_trueStatesAtMeas = false );
00061   declareProperty( "TracksTESPath",
00062                    m_tracksTESPath = "Rec/Track/Ideal" );
00063   declareProperty( "RelationTable",
00064                    m_relationTable = "Rec/Relations/IdealTr2MCP" );
00065   declareProperty( "OTGeometryPath",
00066                    m_otTrackerPath = DeOTDetectorLocation::Default );
00067   declareProperty( "ITGeometryPath",
00068                    m_itTrackerPath = DeSTDetectorLocation::Default );
00069   declareProperty( "VeloGeometryPath",
00070                    m_veloPath = "/dd/Structure/LHCb/Velo" );
00071   declareProperty( "MinNHits", m_minNHits = 6 );
00072   declareProperty( "ErrorX2",  m_errorX2  = 4.0 );
00073   declareProperty( "ErrorY2",  m_errorY2  = 400.0 );
00074   declareProperty( "ErrorTx2", m_errorTx2 = 6.e-5 );
00075   declareProperty( "ErrorTy2", m_errorTy2 = 1.e-4 );
00076   declareProperty( "ErrorP",   m_errorP   = 0.15 );
00077 }

00442 {
00443   unsigned int nITMeas = 0;
00444 
00446   ITCluster2MCParticleAsct::FromRange range = m_itClus2MCP -> rangeTo( mcPart );
00447   ITCluster2MCParticleAsct::FromIterator iClus;
00448   for ( iClus = range.begin(); iClus != range.end(); ++iClus) {
00449     ITCluster* aCluster = iClus->to();
00450     ITMeasurement meas = ITMeasurement( *aCluster, *m_itTracker );
00451     track -> addToMeasurements( meas );
00452     nITMeas++;
00453   }
00454 
00455   debug() << "- " << nITMeas << " ITMeasurements added" << endreq;
00456 
00457   return StatusCode::SUCCESS;
00458 }

00390 {
00391   unsigned int nOTMeas = 0;
00392 
00393   // Loop over outer tracker clusters
00394   OTTimes::const_iterator itOTTime = times -> begin();
00395   while ( itOTTime != times->end() ) {
00396     OTTime* aTime = *itOTTime;
00397 
00398     // retrieve MCHit associated with this cluster
00399     OTTime2MCHitAsct::ToRange range = m_otTim2MCHit -> rangeFrom(aTime);
00400     OTTime2MCHitAsct::ToIterator iTim;
00401     for ( iTim = range.begin(); iTim != range.end(); ++iTim ) {
00402       MCHit* aMCHit = iTim -> to();
00403       MCParticle* aParticle = aMCHit -> mcParticle();
00404       if ( aParticle == mcPart ) {
00405         // get the ambiguity from the MCHit
00406         const OTChannelID channel = aTime -> channel();
00407         DeOTModule* module = m_otTracker -> module( channel );
00408         const HepPoint3D entryP = aMCHit -> entry();
00409         const HepPoint3D exitP = aMCHit -> exit();
00410         double deltaZ = exitP.z() - entryP.z();
00411         if (0.0 == deltaZ) continue; // curling track inside layer
00412         const double tx = (exitP.x() - entryP.x()) / deltaZ;
00413         const double ty = (exitP.y() - entryP.y()) / deltaZ;
00414         double mcDist = module -> distanceToWire( channel.straw(),
00415                                                   entryP, tx, ty );
00416         int ambiguity = 1;
00417         if ( mcDist < 0.0 ) ambiguity = -1;
00418         // Get the tu from the MCHit
00419         double angle = module -> stereoAngle();
00420         double tu = tx * cos(angle) + ty * sin(angle);
00421         // Make the TrMeasurement
00422         OTMeasurement otTim = OTMeasurement( *aTime, *m_otTracker,
00423                                              ambiguity, tu );
00424         track -> addToMeasurements( otTim );
00425         nOTMeas++;
00426       }
00427     }
00428     // next cluster
00429     ++itOTTime;
00430   } // loop over outer tracker clusters
00431 
00432   debug() << "- " << nOTMeas << " OTMeasurements added" << endreq;
00433 
00434   return StatusCode::SUCCESS;
00435 }

00465 {
00466   unsigned int nVeloRMeas   = 0;
00467   unsigned int nVeloPhiMeas = 0;
00468 
00470   VeloCluster2MCParticleAsct::FromRange range = m_veloClus2MCP->rangeTo(mcPart);
00471   VeloCluster2MCParticleAsct::FromIterator iClus;
00472   for ( iClus = range.begin(); iClus != range.end(); ++iClus) {
00473     VeloCluster* aCluster = iClus->to();
00474     // Check if velo cluster is r or phi clusters
00475     if ( m_velo -> isRSensor( aCluster->sensor() ) ) {
00476 
00477       double z = m_velo -> zSensor( aCluster -> sensor() );
00478       double phi = 999.0;
00479       State* tempState;
00480       StatusCode sc = m_stateCreator -> createState( mcPart, z, tempState );
00481       if ( sc.isSuccess() ) {
00482         HepVector vec = tempState -> stateVector();
00483         phi = atan2( vec[1], vec[0] );
00484       }
00485       delete tempState;
00486 
00487       VeloRMeasurement meas = VeloRMeasurement( *aCluster, *m_velo, phi );
00488       track -> addToMeasurements( meas );
00489       nVeloRMeas++;
00490     }
00491     else {
00492       VeloPhiMeasurement meas = VeloPhiMeasurement( *aCluster, *m_velo );
00493       track -> addToMeasurements( meas );
00494       nVeloPhiMeas++;
00495     }
00496   }
00497 
00498   debug() << "- " << nVeloRMeas << " / " << nVeloPhiMeas
00499           << " Velo R/Phi Measurements added" << endreq;
00500 
00501   return StatusCode::SUCCESS;
00502 }

00558 {
00559   for ( std::vector<State*>::iterator it  = track->states().begin();
00560                                       it != track->states().end(); it++) {
00561     track -> removeFromStates( (*it) );
00562   }
00563   
00564   return StatusCode::SUCCESS;
00565 }

00141 {
00142   debug() << "==> Execute" << endreq;
00143 
00144   // Event header info
00145   EventHeader* evtHdr = get<EventHeader>( EventHeaderLocation::Default );
00146   debug() << "Run " << evtHdr -> runNum()
00147           << "\tEvent " << evtHdr -> evtNum() << endreq;
00148 
00149   // Retrieve the MCParticle container
00150   MCParticles* particles = get<MCParticles>( MCParticleLocation::Default );
00151 
00153   OTTimes* otTimes           = get<OTTimes>( OTTimeLocation::Default );
00154   debug() << "- retrieved " <<  otTimes -> size() << " OTTimes." << endreq;
00155 
00156 //  ITClusters* itClusters     = get<ITClusters>( ITClusterLocation::Default );
00157 //  VeloClusters* veloClusters = get<VeloClusters>( VeloClusterLocation::Default );
00158 
00159   // Make container for tracks
00160   Tracks* tracksCont = new Tracks();
00161   //StatusCode sc = eventSvc() -> registerObject( m_tracksTESPath, tracksCont )
00162   //if ( sc.isFailure() ) {
00163     //error() << "Unable to store track container in TES: "
00164             //<< m_tracksTESPath << " (sc = " << sc.getCode() << " )." << endreq;
00165     //return sc;
00166   //}
00167 
00168   // create relation table and register it in the event transient store
00169   typedef RelationWeighted1D<Track,MCParticle,double>  Table;
00170   Table* table = new Table();
00171   StatusCode sc = eventSvc() -> registerObject( m_relationTable, table );
00172   if( sc.isFailure() ) {
00173     error() << "Unable to register the relation container = "
00174             << m_relationTable << " status = " << sc << endreq;
00175     return sc ;
00176   }
00177   else {
00178     verbose() << "Relations container " << m_relationTable
00179               << " registered" << endreq;
00180   }
00181 
00182   debug() << "Starting loop over the "
00183           << particles -> size() << " MCParticles ... " << endreq;
00184 
00185   // loop over MCParticles
00186   // =====================
00187   MCParticles::const_iterator iPart;
00188   for ( iPart = particles->begin(); particles->end() != iPart; ++iPart ) {
00189     MCParticle* mcParticle = *iPart;
00190     if ( m_trackSelector -> select( mcParticle ) ) {
00191       debug() << "- MCParticle of type "
00192               << m_trackSelector -> trackType( mcParticle ) << endreq
00193               << "    - momentum = " << mcParticle -> momentum() << " MeV" <<endreq
00194               << "    - P        = " << mcParticle -> momentum().vect().mag()
00195               << " MeV" <<endreq
00196               << "    - charge   = "
00197               << ( mcParticle -> particleID().threeCharge() / 3 ) << endreq;
00198 
00199 
00200       // Make a new track
00201       // ----------------
00202       Track* track = new Track();
00203 
00204       m_trackSelector -> setTrackType( mcParticle, track );
00205 
00206       // Check whether a Velo track is backward
00207       if ( TrackKeys::Velo == track -> type() ) {
00208         const double pz = mcParticle -> momentum().pz();
00209         if ( pz < 0.0 ) track -> setFlag( TrackKeys::Backward, true );
00210       }
00211 
00212       // Add Velo clusters
00213       // -----------------
00214       if ( m_addVeloClusters == true ) {
00215         sc = addVeloClusters( mcParticle, track );
00216         if ( sc.isFailure() ) {
00217           error() << "Unable to add velo R clusters" << endreq;
00218           return StatusCode::FAILURE;
00219         }
00220       }
00221 
00222       // Sort the measurements in z
00223       // --------------------------
00224       sortMeasurements( track );
00225 
00226       // Add IT clusters
00227       // ---------------
00228       if ( m_addITClusters == true ) {
00229         sc = addITClusters( mcParticle, track );
00230         if ( sc.isFailure() ) {
00231           error() << "Unable to add inner tracker clusters" << endreq;
00232           return StatusCode::FAILURE;
00233         }        
00234       }
00235       
00236       // Add OTTimes
00237       // -----------
00238       if ( m_addOTTimes == true && otTimes != 0 ) {
00239         sc = addOTTimes( otTimes, mcParticle, track );
00240         if ( sc.isFailure() ) {
00241           error() << "Unable to add outer tracker OTTimes" << endreq;
00242           return StatusCode::FAILURE;
00243         }
00244       }
00245 
00246       // Check if the track contains enough hits
00247       // ---------------------------------------
00248       if ( (int) track -> nMeasurements() < m_minNHits) {
00249         delete track;
00250         continue; // go to next track
00251         debug() << " -> track deleted. Had only " <<track -> nMeasurements()
00252                 << " hits" << endreq;
00253       }
00254 
00255       // Sort the measurements in z
00256       // --------------------------
00257       sortMeasurements( track );
00258 
00259       // Initialize a state
00260       // ------------------
00261       if ( m_initState ) {
00262         if ( m_upstream ) {
00263           std::vector<Measurement*>::const_reverse_iterator rbeginM =
00264             track -> measurements().rbegin();
00265           sc = this -> initializeState( (*rbeginM)->z(),
00266                                         track, mcParticle );
00267         }
00268         else {
00269           std::vector<Measurement*>::const_iterator beginM =
00270             track -> measurements().begin();
00271           sc = this -> initializeState( (*beginM)->z(),
00272                                         track, mcParticle );
00273         }
00274         if ( !sc.isSuccess() ) {
00275           delete track;
00276           continue; // go to next track
00277         }
00278       }
00279 
00280       // Set some of the track properties
00281       // --------------------------------
00282       track -> setStatus( TrackKeys::PatRec );
00283       track -> setFlag( TrackKeys::Valid, true );
00284       track -> setFlag( TrackKeys::Unique, true );
00285 
00286       // Fit the track
00287       //if ( m_initState && m_fitTracks ) {
00288         // select appropriate track fitter
00289         //ITrackFitter* fitter = m_tracksFitter;
00290         //if( track -> type() == TrackKeys::Velo
00291             //|| track -> checkFlag( TrackKeys::Backward ) )
00292           //fitter = m_veloTracksFitter;
00293         //if ( track->velo() || track->veloBack() ) fitter = m_veloTracksFitter;
00294         //if ( track -> type() == TrackKeys::Upstream ) fitter = m_veloTTTracksFitter;
00295         //if ( track->veloTT() ) fitter = m_veloTTTracksFitter;
00296         //if ( track -> type() == TrackKeys::Ttrack ) fitter = m_seedTracksFitter;
00297         //if ( track->seed() ) fitter = m_seedTracksFitter;
00298         // Fit the track 
00299         //if ( m_upstream ) { // do upstream fit
00300 //          sc = fitter->fitUpstream(track, track->rbeginM(), track->rendM() );
00301         //} else { // do downstream fit
00302 //          sc = fitter->fitDownstream(track, track->beginM(), track->endM() );
00303         //}
00304         // Set error flag if fit failed
00305         //if ( sc.isFailure() ) track -> setStatus( TrackKeys::Failed );
00306         //if ( sc.isFailure() ) track -> setErrorFlag(1);
00307       //}
00308 
00309       // Add true states at each measurement 
00310       // ===================================
00311       if ( m_trueStatesAtMeas ) {
00312         deleteStates( track );
00313         sc = StatusCode::SUCCESS;
00314         std::vector<Measurement*>::const_iterator iMeas =
00315           track -> measurements().begin();
00316         std::vector<Measurement*>::const_iterator endM =
00317           track -> measurements().end();
00318         //Track::const_measure_iterator iMeas = track->beginM();
00319         while ( sc.isSuccess() && iMeas != endM ) {
00320           State* tempState;
00321           sc = m_stateCreator -> createState( mcParticle,
00322                                               (*iMeas)->z(),
00323                                               tempState );
00324           if ( sc.isSuccess() ) track -> addToStates( *tempState );
00325           else delete tempState;
00326           ++iMeas;
00327         }
00328         if ( sc.isFailure() ) {
00329           delete track;
00330           continue;
00331         }        
00332       }
00333 
00334       // Add the track to the Tracks container and fill the relations table
00335       // ------------------------------------------------------------------
00336       tracksCont -> add( track );
00337       table -> relate( track, mcParticle, 1.0 );
00338 
00339     // debugging Track ...
00340     // -------------------
00341     debug()
00342       << "-> Track with key # " << track -> key() << endreq
00343       << "  * charge         = " << track -> charge() << endreq
00344       << "  * is Valid       = " << track -> checkFlag( TrackKeys::Valid ) << endreq
00345       << "  * is Unique      = " << track -> checkFlag( TrackKeys::Unique ) << endreq
00346       << "  * is of type     = " << track -> type() << endreq
00347       << "  * is Backward    = " << track -> checkFlag( TrackKeys::Backward ) << endreq
00348       << "  * # measurements = " << track -> nMeasurements() << endreq;
00349 
00350     } // is selected
00351   } // looping over MCParticles
00352 
00353   // Store the Tracks in the TES
00354   // ===========================
00355   sc = registerTracks( tracksCont );
00356   if( sc.isFailure() ) return sc;
00357 
00358   debug() << "Created " << tracksCont -> size() << " tracks." << endreq;
00359 
00360   return StatusCode::SUCCESS;
00361 }

00367 {
00368   debug() << "==> Finalize" << endreq;
00369 
00370   // Release all tools and services
00371   // ------------------------------
00372   if ( m_veloClus2MCP )     toolSvc()->releaseTool( m_veloClus2MCP );
00373   if ( m_itClus2MCP )       toolSvc()->releaseTool( m_itClus2MCP );
00374   if ( m_otTim2MCHit  )     toolSvc()->releaseTool( m_otTim2MCHit );
00375   if ( m_trackSelector )    toolSvc()->releaseTool( m_trackSelector );
00376 
00377   //if ( m_tracksFitter )     toolSvc()->releaseTool( m_tracksFitter );
00378   //if ( m_veloTracksFitter ) toolSvc()->releaseTool( m_veloTracksFitter );
00379   //if ( m_seedTracksFitter ) toolSvc()->releaseTool( m_seedTracksFitter );
00380 
00381   return GaudiAlgorithm::finalize();  // must be called after all other actions
00382 }

00088 {
00089   StatusCode sc = GaudiAlgorithm::initialize(); // must be executed first
00090 
00091   // Retrieve the MC associators
00092   // ---------------------------
00093   sc = toolSvc() -> retrieveTool( "OTTime2MCHitAsct", m_otTim2MCHit );
00094   if ( !sc ) { return sc; }
00095   sc = toolSvc() -> retrieveTool( "ITCluster2MCParticleAsct", m_itClus2MCP );
00096   if ( !sc ) { return sc; }
00097   sc = toolSvc() -> retrieveTool( "VeloCluster2MCParticleAsct", m_veloClus2MCP );
00098   if ( !sc ) { return sc; }
00099   debug() << "Associators retrieved." << endreq;
00100 
00101   // Load Geometry from XmlDDDB
00102   // --------------------------
00103    m_otTracker = getDet<DeOTDetector>( m_otTrackerPath );
00104 
00105    m_itTracker = getDet<DeSTDetector>( m_itTrackerPath );
00106 
00107    m_velo      = getDet<DeVelo>( m_veloPath );
00108    debug() << "Geometry read in." << endreq;
00109 
00110   // Retrieve track selection tool
00111   sc = toolSvc() -> retrieveTool("TrackSelector", "select", m_trackSelector,this);
00112   if ( !sc ) { return sc; }
00113 
00114   // Retrieve TrueStateCreator tool
00115   m_stateCreator = tool<IStateCreator>( "TrueStateCreator" );
00116 
00117   if ( m_fitTracks ) {
00118 //    sc = toolSvc()->retrieveTool( "TrFitAtAllPoints", "Fitter", 
00119 //                                  m_tracksFitter, this );
00120 //    if ( !sc ) { return sc; }
00121 //    sc = toolSvc()->retrieveTool( "TrFitAtAllPoints", "SeedFitter", 
00122 //                                  m_seedTracksFitter, this );
00123 //    if ( !sc ) { return sc; }
00124 //    sc = toolSvc()->retrieveTool( "TrFitAtAllPoints", "VeloFitter", 
00125 //                                  m_veloTracksFitter, this );
00126 //    if ( !sc ) { return sc; }
00127 //    sc = toolSvc()->retrieveTool( "TrFitAtAllPoints", "VeloTTFitter", 
00128 //                                  m_veloTTTracksFitter, this );
00129 //    if ( !sc ) { return sc; }
00130   }
00131 
00132   info() << "initialized succesfully" << endreq;
00133 
00134   return StatusCode::SUCCESS;
00135 }

00524 {
00525   State* state;
00526   StatusCode sc = m_stateCreator -> createState( mcPart, z, state );
00527   if ( sc.isSuccess() ) {
00528     // set covariance matrix to a somewhat larger value for the fit.
00529     HepSymMatrix& cov = state -> covariance();
00530     cov.fast(1,1) = m_errorX2;
00531     cov.fast(2,2) = m_errorY2;
00532     cov.fast(3,3) = m_errorTx2;
00533     cov.fast(4,4) = m_errorTy2;
00534     cov.fast(5,5) = pow( m_errorP * state->qOverP(), 2. );
00535 
00536     track -> addToStates( *state );
00537 
00538     double qP = state -> qOverP();
00539 
00540     debug() << "- State added:" << endreq
00541             << "  - position = " << state -> position() << " mm" <<endreq
00542             << "  - momentum = " << state -> momentum() << " MeV" <<endreq
00543             << "  - P        = " << state -> p() << " MeV" <<endreq
00544             << "  - charge   = " << ( qP != 0. ? int(fabs(qP)/qP) : 0  ) << endreq;
00545   }
00546   else {
00547     delete state;
00548     debug() << "-> unable to add a State!" << endreq;
00549   }
00550   
00551   return sc;
00552 }

00508 {
00509   StatusCode sc = put( tracksCont, m_tracksTESPath );
00510   if ( sc.isFailure() )
00511     error() << "Unable to register the output container at "
00512             << m_tracksTESPath << ". Status is " << sc << endreq;
00513   verbose() << "Tracks container stored in the TES" << endreq;
00514 
00515   return StatusCode::SUCCESS;
00516 }

00076                                         {
00077     std::sort( track->measurements().begin(),
00078                track->measurements().end(),
00079                TrackFunctor::increasingByZ<Measurement>() );
00080   };


Public Member Functions
	TrueTracksCreator (const std::string &name, ISvcLocator *pSvcLocator)
	Standard constructor.
virtual	~TrueTracksCreator ()
	Destructor.
virtual StatusCode	initialize ()
	Algorithm initialization.
virtual StatusCode	execute ()
	Algorithm execution.
virtual StatusCode	finalize ()
	Algorithm finalization.
void	sortMeasurements (Track *track)
Private Member Functions
StatusCode	registerTracks (Tracks *tracksCont)
	Register the tracks in the Transient event store.
StatusCode	addOTTimes (OTTimes times, MCParticle mcPart, Track *track)
	Add outer tracker clusters.
StatusCode	addITClusters (MCParticle mcPart, Track track)
	Add inner tracker clusters.
StatusCode	addVeloClusters (MCParticle mcPart, Track track)
	Add velo r and phi clusters.
StatusCode	initializeState (double z, Track track, MCParticle mcPart)
	Initialize seed state.
StatusCode	deleteStates (Track *track)
	Delete all states on the track.
Private Attributes
OTTime2MCHitAsct::IAsct *	m_otTim2MCHit
ITCluster2MCParticleAsct::IAsct *	m_itClus2MCP
VeloCluster2MCParticleAsct::IAsct *	m_veloClus2MCP
DeOTDetector *	m_otTracker
	Pointer to OT detector element.
DeSTDetector *	m_itTracker
	Pointer to ST detector element.
DeVelo *	m_velo
	Pointer to VELO detector element.
std::string	m_otTrackerPath
	Name of the OT XML geom path.
std::string	m_itTrackerPath
	Name of the IT XML geom path.
std::string	m_veloPath
	Name of the Velo XML geom path.
ITrackSelector *	m_trackSelector
	Track selection tool.
IStateCreator *	m_stateCreator
	Create 'seed' state at last measurement.
bool	m_addOTTimes
	true if OT clusters should be put on track
bool	m_addITClusters
	true if IT clusters should be put on track
bool	m_addVeloClusters
	true if Velo R clusters should be put on track
bool	m_initState
	initialize seed state
bool	m_fitTracks
	if true then the tracks are fitted
bool	m_upstream
	upstream or downstream fit
bool	m_trueStatesAtMeas
	Store true states at each measurement position.
std::string	m_tracksTESPath
	TrTracks container path in EvDS.
std::string	m_relationTable
	Location of the associator table.
int	m_minNHits
	Minimum number of hits on the track.
double	m_errorX2
	Error^2 on x.
double	m_errorY2
	Error^2 on y.
double	m_errorTx2
	Error^2 on slope x.
double	m_errorTy2
	Error^2 on slope y.
double	m_errorP
	dp/p

TrueTracksCreator Class Reference

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation