McTkrHitValsTool class Reference

00108                                : ValBase( type, name, parent ),
00109                                  m_partToTrajTab(0),
00110                                  m_pointToPosHitTab(0),
00111                                  m_clusToPosHitTab(0)
00112 {    
00113     // Declare additional interface
00114     declareInterface<IValsTool>(this); 
00115 }

00180 {
00181     StatusCode sc = StatusCode::SUCCESS;
00182     MsgStream  log( msgSvc(), name() );
00183 
00184     // Grab the event header to get the event time for printing a message in debug mode
00185     SmartDataPtr<Event::EventHeader>   header(m_pEventSvc,    EventModel::EventHeader);
00186     double t = header->time();
00187     log << MSG::DEBUG << "Event time: " << t << endreq;;
00188 
00189     // First task is to recover the McParticle Collection and categorize the event
00190     SmartDataPtr<Event::McParticleCol> particleCol(m_pEventSvc, EventModel::MC::McParticleCol);
00191 
00192     // No collection then we do not have MC! 
00193     if (particleCol == 0) return sc;
00194 
00195     // Recover Relational tables that we will need...
00196     SmartDataPtr<Event::McPartToTrajectoryTabList> 
00197         mcPartToTraj(m_pEventSvc, "/Event/MC/McPartToTrajectory");
00198     m_partToTrajTab = new Event::McPartToTrajectoryTab(mcPartToTraj);
00199 
00200     // Recover the McPositionHits to trajectory points relational tables
00201     SmartDataPtr<Event::McPointToPosHitTabList> 
00202         mcPointToPosHitList(m_pEventSvc, "/Event/MC/McPointToPosHit");
00203     m_pointToPosHitTab = new Event::McPointToPosHitTab(mcPointToPosHitList);
00204 
00205     // Recover the McIntegratingHits to trajectory points relational tables
00206     //SmartDataPtr<Event::McPointToIntHitTabList> 
00207     //    mcPointToIntHitList(m_pEventSvc, "/Event/MC/McPointToIntHit");
00208     //m_pointToIntHitTab = new Event::McPointToIntHitTab(mcPointToIntHitList);
00209 
00210     // Retrieve the table giving us TkrDigi to McPositionHit relations
00211     SmartDataPtr<Event::ClusMcPosHitTabList> 
00212         clusMcPosHitList(m_pEventSvc, EventModel::Digi::TkrClusterHitTab);
00213 
00214     // It can happen that there are no clusters...
00215     if (clusMcPosHitList != 0) m_clusToPosHitTab = new Event::ClusMcPosHitTab(clusMcPosHitList);
00216 
00217     // Our first task is to identify the primary in the event
00218     Event::McParticle* primary    = 0;
00219     std::string        primString = "primary";
00220 
00221     Event::McParticleCol::iterator mcPartIter = particleCol->begin();
00222 
00223     for(; mcPartIter != particleCol->end(); mcPartIter++)
00224     {
00225         Event::McParticle* mcPart = *mcPartIter;
00226 
00227         // Is it the primary?
00228         if (mcPart->getProcess() == primString) 
00229         {
00230             primary = mcPart;
00231             break;
00232         }
00233     }
00234 
00235     // Everything follows if we have found a primary
00236     if (primary)
00237     {
00238         // We now count the total number of hits in the Tracker and the number
00239         // due to the primary and its products
00240         CntTotalHits(primary);
00241 
00242         // Count the number of McPositionHits due to primary, electrons and positrons
00243         CntMcPosHits(primary);
00244 
00245         // Next up is to get the number of hits for each of McParticles in the list
00246         // Set up the map to hold the results
00247         std::map<const Event::McParticle*,int> mcPartHitMap;
00248         mcPartHitMap.clear();
00249 
00250         // Reset the McParticle iterator
00251         mcPartIter = particleCol->begin();
00252 
00253         for(++mcPartIter; mcPartIter != particleCol->end(); mcPartIter++)
00254         {
00255             Event::McParticle* mcPart = *mcPartIter;
00256 
00257             // Get the number of hits from this track 
00258             int numHits = GetTrackHits(mcPart);
00259 
00260             mcPartHitMap[mcPart] = numHits;
00261         }
00262 
00263         // Recover the primary's particle identification and number of hits
00264         m_primType    = primary->particleProperty();
00265         m_primNumHits = mcPartHitMap[primary];
00266 
00267         // Get some of the primary's characteristics
00268 //        ParticleProperty* ppty     = m_ppsvc->findByStdHepID(m_primType);
00269 //        std::string       partName = ppty->particle(); 
00270 
00271         // We will need pointers to the "two" daughters to determine shared hits
00272         const Event::McParticle* daughter1 = 0;
00273         const Event::McParticle* daughter2 = 0;
00274 
00275         // Loop through the daughters keeping track of the number of hits
00276         const SmartRefVector<Event::McParticle>& daughterList = primary->daughterList();
00277 
00278         for(SmartRefVector<Event::McParticle>::const_iterator daughterIter = daughterList.begin();
00279             daughterIter != daughterList.end(); daughterIter++)
00280         {
00281             const Event::McParticle* daughter = *daughterIter;
00282 
00283             int numHits = mcPartHitMap[daughter];
00284 
00285             // Check to see if this exceeds the current daughter1 count
00286             if (numHits > m_dght1NumHits || m_dght1Type == 0)
00287             {
00288                 // Check to see if we need to switch down to second daughter
00289                 if (m_dght1NumHits > m_dght2NumHits || m_dght2Type == 0)
00290                 {
00291                     m_dght2NumHits = m_dght1NumHits;
00292                     m_dght2Type    = m_dght2Type;
00293                     daughter2      = daughter1;
00294                 }
00295 
00296                 m_dght1NumHits = numHits;
00297                 m_dght1Type    = daughter->particleProperty();
00298                 daughter1      = daughter;
00299             }
00300             // Otherwise, check to see if this is the hit count for the next particle
00301             else if (numHits > m_dght2NumHits || m_dght2Type == 0)
00302             {
00303                 m_dght2NumHits = numHits;
00304                 m_dght2Type    = daughter->particleProperty();
00305                 daughter2      = daughter;
00306             }
00307         }
00308 
00309         // Finally... get the number of hits these two daughters share
00310         m_dghtSharedHits = GetSharedHits(daughter1, daughter2);
00311     }
00312 
00313     // Clean up after ourselves
00314     if (m_partToTrajTab)    {delete m_partToTrajTab;    m_partToTrajTab    = 0;}
00315     if (m_pointToPosHitTab) {delete m_pointToPosHitTab; m_pointToPosHitTab = 0;}
00316     //if (m_pointToIntHitTab) {delete m_pointToIntHitTab; m_pointToIntHitTab = 0;}
00317     if (m_clusToPosHitTab)  {delete m_clusToPosHitTab;  m_clusToPosHitTab  = 0;}
00318     
00319     log << MSG::DEBUG << " returning. " << endreq;
00320 
00321     return sc;
00322 }

00545 {
00546     // Define codes for electrons and positrons
00547     static Event::McParticle::StdHepId gamma    =  22;
00548     static Event::McParticle::StdHepId electron =  11;
00549     static Event::McParticle::StdHepId positron = -11;
00550 
00551     // Clear counters
00552     m_posHitPrimary   = 0;
00553     m_posHitGammas    = 0;
00554     m_posHitElectrons = 0;
00555     m_posHitPositrons = 0;
00556     m_posHitOthers    = 0;
00557 
00558     // First task is to recover the McParticle Collection and categorize the event
00559     SmartDataPtr<Event::McPositionHitCol> positionHitCol(m_pEventSvc, EventModel::MC::McPositionHitCol);
00560 
00561     // If collection then proceed
00562     if (positionHitCol)
00563     {
00564         // Loop over the collection of McPositionHits
00565         for(Event::McPositionHitCol::iterator posHitIter = positionHitCol->begin();
00566             posHitIter != positionHitCol->end(); posHitIter++)
00567         {
00568             Event::McPositionHit* posHit = *posHitIter;
00569 
00570             // Check the volume identifier to cut on tracker hits
00571             idents::VolumeIdentifier volId = posHit->volumeID();
00572 
00573             // Do not count ACD hits
00574             if (volId[0] != 0) continue;
00575 
00576             // Get parent information
00577             const Event::McParticle* particle = posHit->mcParticle();
00578             const Event::McParticle* mother   = 0;
00579             
00580             if (particle) mother = particle->getMother();
00581 
00582             // Check to see if hit is associated to primary or immediate daughters
00583             if      (particle                  == primary)  m_posHitPrimary++;
00584             else if (mother                    == primary)  m_posHitPrimary++;
00585 
00586             // What kind of particle made the hit?
00587             if      (posHit->getMcParticleId() == gamma)    m_posHitGammas++;
00588             else if (posHit->getMcParticleId() == electron) m_posHitElectrons++;
00589             else if (posHit->getMcParticleId() == positron) m_posHitPositrons++;
00590             else                                            m_posHitOthers++;
00591         }
00592     }
00593 
00594 
00595     return;
00596 }

00326 {
00327     m_totalHits     = 0;
00328     m_totalPrimHits = 0;
00329 
00330     // If there is no cluster table then there are no hits (by definition)
00331     if (m_clusToPosHitTab)
00332     {
00333         // Retrieve pointer to McPositionHit Collection in the TDS
00334         SmartDataPtr<Event::McPositionHitCol> posHitCol(m_pEventSvc, EventModel::MC::McPositionHitCol);
00335 
00336         // Define the results containers
00337         std::map<int, std::set<Event::TkrCluster*> > idToHitMap;    // For each track, number of actual hits made
00338         std::set<int>                                primTrkIdSet;  // For primary, number of tracks resulting
00339         std::set<Event::TkrCluster*>                 primClusSet;   // Number of clusters by tracks from primary
00340 
00341         idToHitMap.clear();
00342         primTrkIdSet.clear();
00343         primClusSet.clear();
00344 
00345         // Loop through the McPositionHit collection
00346         for (Event::McPositionHitCol::iterator posHitIter = posHitCol->begin(); posHitIter != posHitCol->end(); posHitIter++)
00347         {
00348             Event::McPositionHit* posHit = *posHitIter;
00349 
00350             // Is this hit related to a Cluster (ie did this make a hit in the Tracker?)
00351             Event::ClusMcPosHitVec clusHitVec = m_clusToPosHitTab->getRelBySecond(posHit);
00352 
00353             // Do we have anything? (note that this will eliminate McPositionHits from the ACD)
00354            
00355             if (!clusHitVec.empty())
00356             {
00357                 // An McPositionHit belongs to one cluster (but a cluster can be made up of
00358                 // several McPositionHits...)
00359                 Event::TkrCluster* cluster = (*(clusHitVec.begin()))->getFirst();
00360 
00361                 // Check our current map for this hit, start by getting the track id for this hit
00362                 int trackId = posHit->getPackedFlags(); // Track ID is hiding in here...
00363 
00364                 idToHitMap[trackId].insert(cluster);
00365 
00366                 // Does this particle descend from the primary?
00367                 //if (posHit->originMcParticle() == primary) // ok, this turns out to not be valid for gammas... 
00368                 if (posHit->originMcParticle() != 0)         // temporary kludge...
00369                 {
00370                     primClusSet.insert(cluster);
00371                     primTrkIdSet.insert(trackId);
00372                 }
00373             }
00374         }
00375 
00376         // Now pass through our map to count total number of hits from each track 
00377         for (std::map<int,std::set<Event::TkrCluster*> >::iterator idToHitIter = idToHitMap.begin(); 
00378              idToHitIter != idToHitMap.end(); idToHitIter++)
00379         {
00380             m_totalHits += idToHitIter->second.size();
00381         }
00382 
00383         // Go through and count total number of hits from the primary and its descendants
00384         for (std::set<int>::iterator trkIdIter = primTrkIdSet.begin(); trkIdIter != primTrkIdSet.end(); trkIdIter++)
00385         {
00386             m_totalPrimHits += idToHitMap[*trkIdIter].size();
00387         }
00388 
00389         // And now the total number of clusters from the primary and its descendants
00390         m_totalPrimClusHits = primClusSet.size();
00391     }
00392 
00393     return;
00394 }

00457 {
00458     int numShared = 0;
00459 
00460     // Plan? 
00461     // loop through hits on daughter1 track and get clusters. 
00462     //      use table to get mcposhits from cluster
00463     //      loop through mcposhits and look for match to daughter2
00464 
00465     // If there is no cluster table then there are no hits, also protect that both daughters exist
00466     if (m_clusToPosHitTab && daughter1 && daughter2)
00467     {
00468         // Look up trajectory for this particle 
00469         // (use the second daughter because it is shorter)
00470         Event::McPartToTrajectoryVec trajVec = m_partToTrajTab->getRelByFirst(daughter2);
00471         Event::McTrajectory* mcTraj          = (*(trajVec.begin()))->getSecond();
00472 
00473         // Set up to loop through Trajectory points
00474         std::vector<Event::McTrajectoryPoint*> points = mcTraj->getPoints();
00475         std::vector<Event::McTrajectoryPoint*>::const_iterator pointIter;
00476 
00477         // We need a list of clusters resulting from the daughter track, we will then try to relate
00478         // those clusters to those used by the other daughter track 
00479         std::set<Event::TkrCluster*> trackClusters;
00480         trackClusters.clear();
00481 
00482         // Loop through all the McTrajectoryPoints on this track and store away the clusters
00483         for(pointIter = points.begin(); pointIter != points.end(); pointIter++) 
00484         {
00485             // De-reference the McTrajectoryPoint pointer
00486             Event::McTrajectoryPoint* mcPoint = *pointIter;
00487 
00488             // Look for McTrajectoryPoint <--> McPositionHit
00489             Event::McPointToPosHitVec posRelVec = m_pointToPosHitTab->getRelByFirst(mcPoint);
00490 
00491             // Did this trajectory point leave a mark?
00492             if (!posRelVec.empty())
00493             {
00494                 // Can only be one McPositionHit per McTrajectoryPoint... by definition
00495                 Event::McPointToPosHitRel* hitRel = *(posRelVec.begin());
00496 
00497                 Event::McPositionHit* mcPosHit = hitRel->getSecond();
00498 
00499                 // Get any TkrClusters associated to this McPositionHit (note: this eliminates ACD)
00500                 Event::ClusMcPosHitVec clusHitVec = m_clusToPosHitTab->getRelBySecond(mcPosHit);
00501 
00502                 // Do we have anything?
00503                 if (!clusHitVec.empty())
00504                 {
00505                     // An McPositionHit belongs to one cluster (but a cluster can be made up of
00506                     // several McPositionHits...)
00507                     Event::TkrCluster* cluster = (*(clusHitVec.begin()))->getFirst();
00508 
00509                     // "cluster" will be unique at the end 
00510                     trackClusters.insert(cluster);
00511                 }
00512             }
00513         }
00514 
00515         // Go through the set of clusters to look up any relation they might have to the other track
00516         for (std::set<Event::TkrCluster*>::iterator clusIter = trackClusters.begin(); clusIter != trackClusters.end(); clusIter++)
00517         {
00518             Event::TkrCluster* cluster = *clusIter;
00519 
00520             // Look up all McPositionHits associated with this cluster
00521             Event::ClusMcPosHitVec clusHitVec = m_clusToPosHitTab->getRelByFirst(cluster);
00522 
00523             // Proceed if it might be shared
00524             if (clusHitVec.size() > 1)
00525             {
00526                 for(Event::ClusMcPosHitVec::iterator relHitIter = clusHitVec.begin(); relHitIter != clusHitVec.end(); relHitIter++)
00527                 {
00528                     Event::McPositionHit* relPosHit = (*relHitIter)->getSecond();
00529 
00530                     // Check that a valid pointer exists and see if it matches the first daughter
00531                     if (relPosHit && relPosHit->mcParticle() == daughter1)
00532                     {
00533                         numShared++;
00534                         break;
00535                     }
00536                 }
00537             }
00538         }
00539     }
00540 
00541     return numShared;
00542 }

00397 {
00398     int nTrackHits = 0;
00399 
00400     // If there is no cluster table then there are no hits
00401     if (m_clusToPosHitTab)
00402     {
00403         // Look up trajectory for this particle
00404         Event::McPartToTrajectoryVec trajVec = m_partToTrajTab->getRelByFirst(particle);
00405         Event::McTrajectory* mcTraj          = (*(trajVec.begin()))->getSecond();
00406 
00407         // Set up to loop through Trajectory points
00408         std::vector<Event::McTrajectoryPoint*> points = mcTraj->getPoints();
00409         std::vector<Event::McTrajectoryPoint*>::const_iterator pointIter;
00410 
00411         // Use a set to keep track of unique clusters (since a cluster may be associated
00412         // to several McPositionHits). We'll use this to count Tracker hits
00413         std::set<Event::TkrCluster*> trackClusters;
00414         trackClusters.clear();
00415 
00416         // Loop through all the McTrajectoryPoints on this track
00417         for(pointIter = points.begin(); pointIter != points.end(); pointIter++) 
00418         {
00419             // De-reference the McTrajectoryPoint pointer
00420             Event::McTrajectoryPoint* mcPoint = *pointIter;
00421 
00422             // Look for McTrajectoryPoint <--> McPositionHit
00423             Event::McPointToPosHitVec posRelVec = m_pointToPosHitTab->getRelByFirst(mcPoint);
00424 
00425             // Did this trajectory point leave a mark?
00426             if (!posRelVec.empty())
00427             {
00428                 // Can only be one McPositionHit per McTrajectoryPoint... by definition
00429                 Event::McPointToPosHitRel* hitRel = *(posRelVec.begin());
00430 
00431                 Event::McPositionHit* mcPosHit = hitRel->getSecond();
00432 
00433                 // Get any TkrClusters associated to this McPositionHit (note: this eliminates ACD)
00434                 Event::ClusMcPosHitVec clusHitVec = m_clusToPosHitTab->getRelBySecond(mcPosHit);
00435 
00436                 // Do we have anything?
00437                 if (!clusHitVec.empty())
00438                 {
00439                     // An McPositionHit belongs to one cluster (but a cluster can be made up of
00440                     // several McPositionHits...)
00441                     Event::TkrCluster* cluster = (*(clusHitVec.begin()))->getFirst();
00442 
00443                     // "cluster" will be unique at the end 
00444                     trackClusters.insert(cluster);
00445                 }
00446             }
00447         }
00448 
00449         // ok, the size of trackClusters is, then, the number of hits this track left
00450         nTrackHits = trackClusters.size();
00451     }
00452 
00453     return nTrackHits;
00454 }

00141 {
00142     StatusCode sc = StatusCode::SUCCESS;
00143     
00144     MsgStream log(msgSvc(), name());
00145 
00146     if( ValBase::initialize().isFailure()) return StatusCode::FAILURE;
00147   
00148     // get the services    
00149     if( serviceLocator() ) {
00150         if( service("ParticlePropertySvc", m_ppsvc, true).isFailure() ) {
00151             log << MSG::ERROR << "Service [ParticlePropertySvc] not found" << endreq;
00152         }
00153     } else {
00154         return StatusCode::FAILURE;
00155     }
00156  
00157         addItem("McTHPrimType",       &m_primType);
00158     addItem("McTHProcess",        &m_process);
00159         addItem("McTHPrimNumHits",    &m_primNumHits);
00160         addItem("McTHDght1Type",      &m_dght1Type);
00161         addItem("McTHDght1NumHits",   &m_dght1NumHits);
00162         addItem("McTHDght2Type",      &m_dght2Type);
00163         addItem("McTHDght2NumHits",   &m_dght2NumHits);
00164     addItem("McTHPosHitPrimary",  &m_posHitPrimary);
00165     addItem("McTHPosHitGamma",    &m_posHitGammas);
00166     addItem("McTHPosHitElectron", &m_posHitElectrons);
00167     addItem("McTHPosHitPositron", &m_posHitPositrons);
00168     addItem("McTHPosHitOthers",   &m_posHitOthers);
00169     addItem("McTHTotalHits",      &m_totalHits);
00170     addItem("McTHTotPrmHits",     &m_totalPrimHits);
00171     addItem("McTHTotPrmClusHits", &m_totalPrimClusHits);
00172     addItem("McTHDghtSharedHits", &m_dghtSharedHits);
00173 
00174     zeroVals();
00175     
00176     return sc;
00177 }

McTkrHitValsTool Class Reference

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


Public Member Functions
	McTkrHitValsTool (const std::string &type, const std::string &name, const IInterface *parent)
virtual	~McTkrHitValsTool ()
StatusCode	initialize ()
StatusCode	calculate ()
	calculate all values; implemented by each XxxValsTool
Private Member Functions
int	GetTrackHits (const Event::McParticle *paricle)
void	CntTotalHits (const Event::McParticle *primary)
int	GetSharedHits (const Event::McParticle daught1, const Event::McParticle daught2)
void	CntMcPosHits (const Event::McParticle *primary)
Private Attributes
int	m_primType
int	m_primNumHits
int	m_dght1Type
int	m_dght1NumHits
int	m_dght2Type
int	m_dght2NumHits
int	m_posHitPrimary
int	m_posHitGammas
int	m_posHitElectrons
int	m_posHitPositrons
int	m_posHitOthers
int	m_process
int	m_totalHits
int	m_totalPrimHits
int	m_totalPrimClusHits
int	m_dghtSharedHits
IParticlePropertySvc *	m_ppsvc
Event::McPartToTrajectoryTab *	m_partToTrajTab
Event::McPointToPosHitTab *	m_pointToPosHitTab
Event::ClusMcPosHitTab *	m_clusToPosHitTab