McValsTool class Reference

00129                        : ValBase( type, name, parent )
00130 {    
00131     // Declare additional interface
00132     declareInterface<IValsTool>(this); 
00133 }

00255 {
00256     StatusCode sc = StatusCode::SUCCESS;
00257     
00258     // Recover Track associated info. 
00259     SmartDataPtr<Event::TkrTrackCol>  pTracks(m_pEventSvc,EventModel::TkrRecon::TkrTrackCol); 
00260     SmartDataPtr<Event::TkrVertexCol>       pVerts(m_pEventSvc,EventModel::TkrRecon::TkrVertexCol);
00261     // Recover MC Pointer
00262     SmartDataPtr<Event::McParticleCol> pMcParticle(m_pEventSvc, EventModel::MC::McParticleCol);
00263     // this is avoid creating the object as a side effect!!
00264     SmartDataLocator<Event::MCEvent> pMcEvent(m_pEventSvc, EventModel::MC::Event);
00265 
00266     if(!pMcParticle) return sc;
00267 
00268     if(pMcEvent) {
00269         char temp[2] = "_";
00270         MC_SourceId = pMcEvent->getSourceId();
00271         strncpy(MC_SourceName, pMcEvent->getSourceName().c_str(),80);
00272         if (MC_SourceName=="") strncpy(MC_SourceName, temp, 80);
00273     }
00274     
00275     MC_Energy = -1;
00276     
00277     if (pMcParticle) {
00278         
00279         Event::McParticleCol::const_iterator pMCPrimary = pMcParticle->begin();
00280         // Skip the first particle... it's for bookkeeping.
00281         // The second particle is the first real propagating particle.
00282         // except in the case of beamtest data!!!
00283         MC_NumIncident = (*pMCPrimary)->daughterList().size();
00284 
00285         // if there are no incident particles
00286         if (MC_NumIncident==0) return sc;
00287 
00288         // ok, go ahead and call the ACD stuff now
00289         getAcdReconVars();
00290 
00291         // if there is one incident particle, it's okay to use it as before
00292         // if there are more than one, I don't know what to do, for now
00293         //     use the mother particle. That way, at least the energy will
00294         //     be correct.
00295 
00296         if(MC_NumIncident == 1) {
00297 
00298             Event::McParticle::StdHepId hepid= (*pMCPrimary)->particleProperty();
00299             MC_Id = (double)hepid;
00300             ParticleProperty* ppty = m_ppsvc->findByStdHepID( hepid );
00301             if (ppty) {
00302                 std::string name = ppty->particle(); 
00303                 MC_Charge = ppty->charge();          
00304             }
00305 
00306             pMCPrimary++;
00307         }
00308         
00309         HepPoint3D Mc_x0;
00310         // launch point for charged particle; conversion point for neutral
00311         Mc_x0 = (MC_Charge==0 ? (*pMCPrimary)->finalPosition() : (*pMCPrimary)->initialPosition());
00312         CLHEP::HepLorentzVector Mc_p0 = (*pMCPrimary)->initialFourMomentum();
00313         // there's a method v.m(), but it does something tricky if m2<0
00314         double mass = sqrt(std::max(Mc_p0.m2(),0.0));
00315         
00316         Vector Mc_t0 = Vector(Mc_p0.x(),Mc_p0.y(), Mc_p0.z()).unit();
00317         
00318         //Pure MC Tuple Items
00319         MC_Energy = std::max(Mc_p0.t() - mass, 0.0);
00320         MC_LogEnergy = log10(MC_Energy);
00321         
00322         MC_x0     = Mc_x0.x();
00323         MC_y0     = Mc_x0.y();
00324         MC_z0     = Mc_x0.z();
00325         
00326         MC_xdir   = Mc_t0.x();
00327         MC_ydir   = Mc_t0.y();
00328         MC_zdir   = Mc_t0.z();
00329 
00330         // convert to (ra, dec)
00331         // moved to McCoordsAlg to accomodate Interleave
00332  
00333         //Attempt to estimate energy exiting the tracker
00334         MC_TkrExitEne = getEnergyExitingTkr(*pMCPrimary);
00335 
00336         if((*pMCPrimary)->daughterList().size() > 0) {
00337             SmartRefVector<Event::McParticle> daughters = (*pMCPrimary)->daughterList();
00338             SmartRef<Event::McParticle> pp1 = daughters[0]; 
00339             std::string interaction = pp1->getProcess();
00340             if(interaction == "conv") { // Its a photon conversion; For comptons "compt" or brems "brem"  
00341                 CLHEP::HepLorentzVector Mc_p1 = pp1->initialFourMomentum();
00342                 SmartRef<Event::McParticle> pp2 = daughters[1];
00343                 CLHEP::HepLorentzVector Mc_p2 = pp2->initialFourMomentum();
00344                 double e1 = Mc_p1.t();
00345                 double e2 = Mc_p2.t();
00346                 MC_EFrac = e1/MC_Energy; 
00347                 if(e1 < e2) MC_EFrac = e2/MC_Energy;
00348                 Vector Mc_t1 = Vector(Mc_p1.x(),Mc_p1.y(), Mc_p1.z()).unit();
00349                 Vector Mc_t2 = Vector(Mc_p2.x(),Mc_p2.y(), Mc_p2.z()).unit();
00350                 double dot_prod = Mc_t1*Mc_t2;
00351                 if(dot_prod > 1.) dot_prod = 1.;
00352                 MC_OpenAngle = acos(dot_prod);
00353             }  
00354         }
00355 
00356         // This should really be a test on vertices, not tracks
00357         if(!pTracks) return sc; 
00358         int num_tracks = pTracks->size(); 
00359         if(num_tracks <= 0 ) return sc;
00360         
00361         // Get track energies and event energy
00362         Event::TkrTrackColConPtr pTrack1 = pTracks->begin();
00363         const Event::TkrTrack*   track_1 = *pTrack1;
00364         
00365         double e1 = track_1->getInitialEnergy();
00366         double gamEne = e1; 
00367         double e2 = 0.; 
00368         if(num_tracks > 2) {
00369             pTrack1++;
00370             const Event::TkrTrack* track_2 = *pTrack1;
00371             e2 = track_2->getInitialEnergy();
00372             gamEne += e2;
00373         }
00374         
00375         //Make sure we have valid reconstructed data
00376         if (pVerts) {
00377             
00378             // Get the first Vertex - First track of first vertex = Best Track
00379             Event::TkrVertexConPtr pVtxr = pVerts->begin(); 
00380             Event::TkrVertex*   gamma = *pVtxr++; 
00381             Point  x0 = gamma->getPosition();
00382             Vector t0 = gamma->getDirection();
00383 
00384             // Reference position errors at the start of recon track(s)
00385             double arc_len = (x0.z()-Mc_x0.z())/Mc_t0.z();
00386             HepPoint3D x_start = Mc_x0 + arc_len*Mc_t0;
00387             MC_x_err  = x0.x()-x_start.x(); 
00388             MC_y_err  = x0.y()-x_start.y();
00389             // except for z, use the difference between MC conversion point and start of vertex track
00390             MC_z_err  = x0.z()-Mc_x0.z();
00391             
00392             MC_xdir_err = t0.x()-Mc_t0.x(); 
00393             MC_ydir_err = t0.y()-Mc_t0.y();
00394             MC_zdir_err = t0.z()-Mc_t0.z();
00395 
00396                         bool VTX_set = false;
00397                         for(;pVtxr != pVerts->end(); pVtxr++) {
00398                 Event::TkrVertex* vtxN = *pVtxr; 
00399                                 if(vtxN->getStatusBits()& Event::TkrVertex::NEUTRALVTX) {
00400                                         Vector tN = vtxN->getDirection();
00401                     double acostNtMC = acos(tN*Mc_t0);
00402                                         if(!(VTX_set)) {
00403                                                 MC_dir_errN  = acostNtMC;
00404                                                 VTX_set = true;
00405                                         }
00406                                         MC_dir_errN1 = acostNtMC; // Assumes last VTX is 1Tkr Neutral Vtx
00407                         }   }   
00408 
00409             
00410             double cost0tMC = t0*Mc_t0;
00411             
00412             MC_dir_err  = acos(cost0tMC);
00413             
00414             int nParticles = gamma->getNumTracks(); 
00415             
00416             SmartRefVector<Event::TkrTrack>::const_iterator pTrack1 = gamma->getTrackIterBegin();  
00417             const Event::TkrTrack* track_1 = *pTrack1;
00418             
00419             Point  x1 = track_1->front()->getPoint(Event::TkrTrackHit::SMOOTHED);
00420             Vector t1 = track_1->front()->getDirection(Event::TkrTrackHit::SMOOTHED);
00421             double cost1tMC = t1*Mc_t0;
00422             
00423             MC_TKR1_dir_err  = acos(cost1tMC);
00424             
00425             if(nParticles > 1) {
00426                 pTrack1++;
00427                 const Event::TkrTrack* track_2 = *pTrack1;
00428                 Point  x2 = track_2->front()->getPoint(Event::TkrTrackHit::SMOOTHED);
00429                 Vector t2 = track_2->front()->getDirection(Event::TkrTrackHit::SMOOTHED);
00430                 double cost2tMC = t2*Mc_t0;
00431                 
00432                 MC_TKR2_dir_err  = acos(cost2tMC);
00433             }
00434         } 
00435     }
00436 
00437     return sc;
00438 }

00507                                  {
00508 
00509   MsgStream log(msgSvc(), name());  
00510 
00511   double bestActDist(-2000.);
00512   idents::AcdId bestId;
00513   std::map<idents::AcdId, double> energyIdMap;
00514 
00515   SmartDataPtr<Event::AcdRecon>           pACD(m_pEventSvc,EventModel::AcdRecon::Event);
00516   if (pACD) {
00517     // Make a map relating AcdId to energy in the tile
00518     const std::vector<idents::AcdId>& tileIds = pACD->getIdCol();
00519     const std::vector<double>& tileEnergies   = pACD->getEnergyCol();
00520     
00521     int maxTiles = tileIds.size();
00522     for (int i = 0; i<maxTiles; i++) {
00523       energyIdMap[tileIds[i]] = tileEnergies[i];
00524     }
00525   }
00526 
00527   // Here we will loop over calculated poca's to find best (largest) active distance
00528   SmartDataPtr<Event::AcdTkrHitPocaCol> acdTkrHits(m_pEventSvc,EventModel::MC::McAcdTkrHitPocaCol); 
00529   if ( ! acdTkrHits ) {
00530     // no poca's found, set guard values.
00531     log << "no AcdTkrHitPocas found on TDS" << endreq;
00532     MC_AcdActiveDist3D = bestActDist;
00533     MC_AcdActDistTileId = bestId.id();
00534     MC_AcdActDistTileEnergy = -1.;
00535   } else {  
00536     for ( Event::AcdTkrHitPocaCol::const_iterator itr = acdTkrHits->begin();
00537           itr != acdTkrHits->end(); itr++ ) {
00538       const Event::AcdTkrHitPoca* aHitPoca = *itr;
00539       if ( aHitPoca->getDoca() > bestActDist ) {
00540         bestId = aHitPoca->getId();
00541         bestActDist = aHitPoca->getDoca();
00542       }
00543     }  
00544     // latch values
00545     MC_AcdActiveDist3D = bestActDist;
00546     MC_AcdActDistTileId = bestId.id();
00547     MC_AcdActDistTileEnergy = energyIdMap[bestId];
00548   }
00549 
00550   // Here we will get the point the MC parent enters the ACD volume
00551   SmartDataPtr<Event::AcdTkrPointCol> acdPoints(m_pEventSvc,EventModel::MC::McAcdTkrPointCol);
00552   if ( ! acdPoints || acdPoints->size() < 1 ) {
00553     // no point's found, set guard values.
00554     log << "no AcdTkrPoints found on TDS" << endreq;
00555     MC_AcdXEnter = MC_AcdYEnter = MC_AcdZEnter = -2000.;
00556   } else {
00557     const Event::AcdTkrPoint* entryPoint = (*acdPoints)[0];
00558     const Point& thePoint = entryPoint->point();
00559     
00560     // latch values
00561     MC_AcdXEnter = thePoint.x(); 
00562     MC_AcdYEnter = thePoint.y(); 
00563     MC_AcdZEnter = thePoint.z(); 
00564   }
00565 
00566 }

00441 {
00442     //This function attempts to estimate the energy which exist the tracker during 
00443     //an event. The idea is to recursively traverse the McParticle tree and look for
00444     //particles which start in the tracker (so, for example, the electron and positron
00445     //from a gamma conversion) and then exit, keep track of the energy of these particles.
00446     //Due to the vagaries of our simulation, several interesting problems can arise and
00447     //an attempt is made to deal with these (see notes below).
00448     //A flaw here is that NO attempt is made to account for continuous energy loss of 
00449     //charged particles in the tracker. Hopefully, a more sophisticated routine will 
00450     //appear soon.
00451     double tkrExitEne = 0.;
00452     double partEnergy = 0.;
00453     double partLostE  = 0.;
00454 
00455     //idents::VolumeIdentifier initial = mcPart->getInitialId();
00456     idents::VolumeIdentifier final   = mcPart->getFinalId();
00457 
00458     //This should check that the initial point of the particle is in the tracker
00459     //and that its final point is outside. So, it started in the tracker and 
00460     //carried its energy outside (could be anywhere)
00463     if (final.size() == 9 && (final[0] != 0 || final[3] != 1) )
00464     {
00465         //Set total initial energy of this particle...
00466         partEnergy = mcPart->initialFourMomentum().t();
00467     }
00468 
00469     //Next step is to go through the daughter list to find the tracker exiting
00470     //energy due to them
00471     SmartRefVector<Event::McParticle> daughters = mcPart->daughterList();
00472 
00473     for(SmartRefVector<Event::McParticle>::iterator partIter = daughters.begin(); 
00474                                                     partIter < daughters.end(); 
00475                                                     partIter++)
00476     {
00477         Event::McParticle* daughter = *partIter;
00478 
00479         //First we check to see if the daughter was created in the tracker. 
00480         idents::VolumeIdentifier daughterInitial = daughter->getInitialId();
00481 
00482         if (daughterInitial.size() == 9 && daughterInitial[0] == 0 && daughterInitial[3] ==1)
00483         {
00484             //Keep track of the energy this daughter takes from its parent
00485             partLostE += daughter->initialFourMomentum().t();
00486 
00487             //Now get the energy exiting the tracker due to this particle
00488             tkrExitEne += getEnergyExitingTkr(daughter);
00489         }
00490     }
00491 
00492     //Ok, now correct the mcPart energy for that it lost traversing the tracker by 
00493     //creating daughter particles
00494     partEnergy -= partLostE;
00495 
00496     if (partEnergy < 0.) partEnergy = 0.;
00497 
00498     //Finally, add this to the energy exiting the tracker
00499     tkrExitEne += partEnergy;
00500 
00501     //This should now be the energy exiting the tracker due to those particles created 
00502     //in the tracker. 
00503     return tkrExitEne;
00504 }

00192 {
00193     StatusCode sc = StatusCode::SUCCESS;
00194     
00195     MsgStream log(msgSvc(), name());
00196 
00197     if( ValBase::initialize().isFailure()) return StatusCode::FAILURE;
00198      
00199     if( serviceLocator() ) {
00200         if( service("ParticlePropertySvc", m_ppsvc, true).isFailure() ) {
00201             log << MSG::ERROR << "Service [ParticlePropertySvc] not found" << endreq;
00202         }
00203     } else {
00204         return StatusCode::FAILURE;
00205     }
00206     
00207     // load up the map
00208 
00209     addItem("McSourceId",     &MC_SourceId);
00210     addItem("McSourceName",    MC_SourceName);
00211     addItem("McNumIncident",  &MC_NumIncident);
00212     addItem("McId",           &MC_Id);  
00213     addItem("McCharge",       &MC_Charge);
00214     addItem("McEnergy",       &MC_Energy);  
00215     addItem("McLogEnergy",    &MC_LogEnergy);
00216     addItem("McEFrac",        &MC_EFrac);
00217     addItem("McOpenAngle",    &MC_OpenAngle);
00218     addItem("McTkrExitEne",   &MC_TkrExitEne);
00219     addItem("McX0",           &MC_x0);           
00220     addItem("McY0",           &MC_y0);           
00221     addItem("McZ0",           &MC_z0);  
00222     
00223     addItem("McXDir",         &MC_xdir);         
00224     addItem("McYDir",         &MC_ydir);         
00225     addItem("McZDir",         &MC_zdir);         
00226     
00227     addItem("McXErr",         &MC_x_err);        
00228     addItem("McYErr",         &MC_y_err);        
00229     addItem("McZErr",         &MC_z_err);        
00230     
00231     addItem("McXDirErr",      &MC_xdir_err);     
00232     addItem("McYDirErr",      &MC_ydir_err);     
00233     addItem("McZDirErr",      &MC_zdir_err);     
00234     
00235     addItem("McDirErr",       &MC_dir_err);      
00236     addItem("McTkr1DirErr",   &MC_TKR1_dir_err); 
00237     addItem("McTkr2DirErr",   &MC_TKR2_dir_err); 
00238         addItem("McDirErrN",      &MC_dir_errN); 
00239         addItem("McDirErrN1",      &MC_dir_errN1); 
00240 
00241     addItem("McAcdXEnter",     &MC_AcdXEnter);
00242     addItem("McAcdYEnter",     &MC_AcdYEnter);    
00243     addItem("McAcdZEnter",     &MC_AcdZEnter);
00244 
00245     addItem("McAcdActiveDist3D", &MC_AcdActiveDist3D);
00246     addItem("McAcdActDistTileId", &MC_AcdActDistTileId);
00247     addItem("McAcdActDistTileEnergy", &MC_AcdActDistTileEnergy);
00248     
00249     zeroVals();
00250     
00251     return sc;
00252 }


Public Member Functions
	McValsTool (const std::string &type, const std::string &name, const IInterface *parent)
virtual	~McValsTool ()
StatusCode	initialize ()
StatusCode	calculate ()
	calculate all values; implemented by each XxxValsTool
Private Member Functions
double	getEnergyExitingTkr (Event::McParticle *mcPart)
void	getAcdReconVars ()
Private Attributes
float	MC_SourceId
char	MC_SourceName [80]
float	MC_NumIncident
float	MC_Id
float	MC_Charge
float	MC_Energy
float	MC_LogEnergy
float	MC_EFrac
float	MC_OpenAngle
float	MC_TkrExitEne
float	MC_x0
float	MC_y0
float	MC_z0
float	MC_xdir
float	MC_ydir
float	MC_zdir
float	MC_x_err
float	MC_y_err
float	MC_z_err
float	MC_xdir_err
float	MC_ydir_err
float	MC_zdir_err
float	MC_dir_err
float	MC_dir_errN
float	MC_dir_errN1
float	MC_TKR1_dir_err
float	MC_TKR2_dir_err
float	MC_AcdXEnter
float	MC_AcdYEnter
float	MC_AcdZEnter
float	MC_AcdActiveDist3D
float	MC_AcdActDistTileId
float	MC_AcdActDistTileEnergy
IParticlePropertySvc *	m_ppsvc

McValsTool Class Reference

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation