McAnalValsTool class Reference

00159                                : ValBase( type, name, parent )
00160 {    
00161     // Declare additional interface
00162     declareInterface<IValsTool>(this); 
00163 }

00572 {
00573     CLHEP::Hep3Vector       mainVec    = m_mcEvent->getTrackDirection(mcPart);
00574     CLHEP::HepLorentzVector main4mom   = mcPart->initialFourMomentum();
00575     double                  mainMom    = main4mom.vect().mag();
00576     CLHEP::HepLorentzVector main4cls(mainMom*mainVec,main4mom.e());
00577     double                  cosAng2Gam = main4mom.vect().unit().dot(mcGamma->initialFourMomentum().vect().unit());
00578 
00579     if (ang2Gam < -1.) 
00580     {
00581         ang2Gam = -1.;
00582     }
00583     if (ang2Gam >  1.)
00584     {
00585         ang2Gam = 1.;
00586     }
00587     ang2Gam     = acos(cosAng2Gam);
00588 
00589     cosAng2Gam  = main4cls.vect().unit().dot(mcPart->initialFourMomentum().vect().unit());
00590     cls2Gam     = acos(cosAng2Gam);
00591 
00592     mcTrkRms    = m_mcEvent->getTrackStraightness(mcPart);
00593 
00594     prtType     = mcPart->particleProperty();
00595     prtEnergy   = main4mom.e();
00596     prtCosDirX  = main4mom.vect().unit().x();
00597     prtCosDirY  = main4mom.vect().unit().y();
00598     prtCosDirZ  = main4mom.vect().unit().z();
00599 
00600     prtNHits    = (m_mcEvent->getMcPartTrack(mcPart)).size();
00601     prtNClstrs  = m_mcEvent->getNumClusterHits(mcPart);
00602     prtNGaps    = m_mcEvent->getNumGaps(mcPart);
00603 
00604     //If track gaps, find size and check to see if first one shortens the track
00605     if (prtNGaps > 0)
00606     {
00607         prt1stGapSz = m_mcEvent->getGapSize(mcPart, 0);
00608         prtNHits2Gp = m_mcEvent->getGapStartHitNo(mcPart, 0) - 1;
00609     }
00610 
00611     return;
00612 }

00617 {
00618     int nBremsInt,nDeltasInt,nDeltaHtInt;
00619 
00620     lastHitE = m_mcEvent->getTrackELastHit(mcPart);
00621     radLossE = m_mcEvent->getTrackBremELoss(mcPart, nBremsInt);
00622     dltaRayE = m_mcEvent->getTrackDeltaELoss(mcPart, nDeltasInt, nDeltaHtInt);
00623     nDeltas  = m_mcEvent->getTrackDeltaRange(mcPart, aveRange, maxRange);
00624 
00625     nBrems   = nBremsInt;
00626     nDeltas  = nDeltasInt;
00627     nDeltaHt = nDeltaHtInt;
00628 
00629     return;
00630 }

00352 {
00353     StatusCode sc = StatusCode::SUCCESS;
00354     MsgStream  log( msgSvc(), name() );
00355 
00356     if (m_mcEvent==0 || m_mcTracks==0) return sc;
00357 
00358     // Retrieving pointers from the TDS 
00359     SmartDataPtr<Event::EventHeader>   header(m_pEventSvc,    EventModel::EventHeader);
00360     SmartDataPtr<Event::MCEvent>       mcheader(m_pEventSvc,  EventModel::MC::Event);
00361     SmartDataPtr<Event::McParticleCol> particles(m_pEventSvc, EventModel::MC::McParticleCol);
00362 
00363     double t = header->time();
00364     log << MSG::DEBUG << "Event time: " << t << endreq;;
00365 
00366     if (m_mcEvent)
00367     {
00368         //int numTracksTotal = m_mcEvent->getNumMcTracks();
00369         int classifyBits   = m_mcEvent->getClassificationBits();
00370 
00371         // Pointers to the primary particle and (eventually) its secondaries
00372         const Event::McParticle* mcPart   = m_mcEvent->getPrimaryParticle();
00373         const Event::McParticle* mcMain   = 0;
00374         const Event::McParticle* mcSecond = 0;
00375 
00376         idents::VolumeIdentifier curVolId    = const_cast<Event::McParticle*>(mcPart)->getFinalId();
00377         CLHEP::HepLorentzVector  primary4mom = mcPart->initialFourMomentum();
00378 
00379         // Retrieve information about the primary particle
00380         m_prmEnergy     = primary4mom.e();
00381         m_prmDecEne     = mcPart->finalFourMomentum().e();
00382         m_prmPosX       = mcPart->initialPosition().x();
00383         m_prmPosY       = mcPart->initialPosition().y();
00384         m_prmPosZ       = mcPart->initialPosition().z();
00385         m_prmCosDirX    = primary4mom.vect().unit().x();
00386         m_prmCosDirY    = primary4mom.vect().unit().y();
00387         m_prmCosDirZ    = primary4mom.vect().unit().z();
00388         m_prmDecPosX    = mcPart->finalPosition().x();
00389         m_prmDecPosY    = mcPart->finalPosition().y();
00390         m_prmDecPosZ    = mcPart->finalPosition().z();
00391         m_prmDecCosX    = mcPart->finalFourMomentum().vect().unit().x();
00392         m_prmDecCosY    = mcPart->finalFourMomentum().vect().unit().y();
00393         m_prmDecCosZ    = mcPart->finalFourMomentum().vect().unit().z();
00394         m_prmDecCode    = classifyBits;
00395         m_prmNDghtrs    = mcPart->daughterList().size();
00396         m_prmNumSecndry = m_mcEvent->getNumSecondaries();
00397         m_prmNumAsscted = m_mcEvent->getNumAssociated();
00398 
00399         // If charged particle then get base information from this path...
00400         if (classifyBits & Event::McEventStructure::CHARGED)
00401         {
00402             // Get the initial volume id
00403             curVolId = const_cast<Event::McParticle*>(mcPart)->getInitialId();
00404 
00405             // For charged tracks, the main secondary is the track itself
00406             mcMain   = mcPart;
00407         }
00408         // If gamma then this path
00409         else if (classifyBits & Event::McEventStructure::GAMMA)
00410         {
00411             primary4mom = mcPart->finalFourMomentum();
00412 
00413             // This should mean that the gamma converted in the tracker (or, unfortunately, the ACD)
00414             if (m_prmNumSecndry > 0)
00415             {
00416                 double mcMainE   = 0.;
00417                 double mcSecondE = 0.;
00418 
00419                 // Identify the "primary" and "secondary" tracks in the conversion
00420                 // Primary = highest energy particle
00421                 // Secondary = lower energy of pair
00422                 for(int idx = 0; idx < m_prmNumSecndry; idx++)
00423                 {
00424                     const Event::McParticle* mcPart1 = m_mcEvent->getSecondary(idx);
00425 
00426                     // We are only interested in the tracks resulting from gamma conversion
00427                     // Note, however, that if process was purely compton, etc., then we will pass this
00428                     if (classifyBits & Event::McEventStructure::TRKCONVERT && mcPart1->getProcess() != "conv") 
00429                         continue;
00430 
00431                     // If we have already found a candidate primary, check against secondary
00432                     if (mcMain)
00433                     {
00434                         mcSecond  = mcPart1;
00435                         mcSecondE = mcPart1->initialFourMomentum().e();
00436 
00437                         // Primary defined to be highest energy of pair
00438                         if (mcMainE < mcSecondE)
00439                         {
00440                             mcSecond  = mcMain;
00441                             mcSecondE = mcMainE;
00442                             mcMain    = mcPart1;
00443                             mcMainE   = mcPart1->initialFourMomentum().e();
00444                         }
00445                     }
00446                     else
00447                     {
00448                         mcMain  = mcPart1;
00449                         mcMainE = mcPart1->initialFourMomentum().e();
00450                     }
00451                 }
00452             }
00453         }
00454 
00455         // Make sure conversion is in a tracker tower
00456         if (curVolId[0] == 0 && curVolId[3] == 1 && curVolId.size() > 3)
00457         {
00458             m_prmDecLayer = 19;
00459 
00460             if (curVolId.size() > 6)
00461             {
00462                 int trayNum = curVolId[4];
00463                 int botTop  = curVolId[6];
00464 
00465                 m_prmDecInCnv = botTop;
00466                 m_prmDecLayer = trayNum;
00467             }
00468         }
00469 
00470         // Make sure we have found a primary particle
00471         if (mcMain)
00472         {
00473             // To get the sum of the angles
00474             CLHEP::Hep3Vector       mainVec    = m_mcEvent->getTrackDirection(mcMain);
00475             CLHEP::HepLorentzVector main4mom   = mcMain->initialFourMomentum();
00476             double                  mainMom    = main4mom.vect().mag();
00477             CLHEP::HepLorentzVector main4cls(mainMom*mainVec,main4mom.e());
00478 
00479             // Fill the angle related info
00480             calcMcAngleInfo(mcMain, mcPart, m_scd1Ang2Gam, m_scd1Cls2Gam, m_scd1McTrkRms, m_scd1Type,
00481                             m_scd1Energy, m_scd1CosDirX, m_scd1CosDirY, m_scd1CosDirZ, m_scd1NHits, m_scd1NClstrs,
00482                             m_scd1NGaps, m_scd11stGapSz, m_scd1NHits2Gp);
00483 
00484             // Look here at track energy loss (if possible)
00485             calcMcEnergyInfo(mcMain, m_scd1ELastHit, m_scd1RadELoss, m_scd1DeltaRay,
00486                              m_scd1NBrems, m_scd1NDeltas, m_scd1NDeltaHt, m_scd1AveRange, m_scd1MaxRange);
00487 
00488             // Want to compress into the double scdType
00489             //short int* typeWords  = reinterpret_cast<short int*>(&m_scd1Type);
00490 
00491             //typeWords[0] = partType;
00492             //typeWords[1] = 0; //unused
00493             //typeWords[2] = m_mcEvent->getTrackHitLayer(mcMain,   0);
00494             //typeWords[3] = m_mcEvent->getTrackHitLayer(mcMain, 100);
00495             m_scd1FirstLyr = m_mcEvent->getTrackHitLayer(mcMain,   0);
00496             m_scd1LastLyr  = m_mcEvent->getTrackHitLayer(mcMain, 100);
00497 
00498             if (m_scd1RadELoss > m_scd1Energy - m_scd1ELastHit)
00499             {
00500                 //int jj=0;
00501             }
00502 
00503             // Testing at this point
00504             //const Event::TkrPatCand* patCand = m_mcTracks->getBestTkrPatCand(mcMain);
00505             //if (patCand)
00506             //{
00507             //    int numHits = m_mcTracks->getNumMcHits(patCand,mcMain);
00508             //    int jjj=0;
00509             //}
00510 
00511             // If there are two gamma conversion tracks then do this part
00512             if (mcSecond)
00513             {
00514                 CLHEP::Hep3Vector       secondVec = m_mcEvent->getTrackDirection(mcSecond);
00515                 CLHEP::HepLorentzVector scnd4mom  = mcSecond->initialFourMomentum();
00516                 double                  scndMom   = scnd4mom.vect().mag();
00517                 CLHEP::HepLorentzVector scnd4cls(scndMom*secondVec,scnd4mom.e());
00518                 
00519                 calcMcAngleInfo(mcSecond, mcPart, m_scd2Ang2Gam, m_scd2Cls2Gam, m_scd2McTrkRms, m_scd2Type,
00520                                 m_scd2Energy, m_scd2CosDirX, m_scd2CosDirY, m_scd2CosDirZ, m_scd2NHits, m_scd2NClstrs,
00521                                 m_scd2NGaps, m_scd21stGapSz, m_scd2NHits2Gp);
00522                 
00523                 calcMcEnergyInfo(mcSecond, m_scd2ELastHit, m_scd2RadELoss, m_scd2DeltaRay,
00524                                  m_scd2NBrems, m_scd2NDeltas, m_scd2NDeltaHt, m_scd2AveRange, m_scd2MaxRange);
00525                 
00526                 m_scd2FirstLyr = m_mcEvent->getTrackHitLayer(mcSecond,   0);
00527                 m_scd2LastLyr  = m_mcEvent->getTrackHitLayer(mcSecond, 100);
00528 
00529                 // Check information on shared hits with other secondary tracks
00530                 unsigned int sharedInfo = m_mcEvent->getSharedHitInfo(mcMain,mcSecond);
00531                 //int          nShared    = sharedInfo & 0x7;
00532 
00533                 m_prmTrkPattern = sharedInfo >> 4;
00534 
00535                 // Update the MC direction from the four vector
00536                 main4mom += scnd4mom;
00537                 main4cls += scnd4cls;
00538             }
00539 
00540             // Look at reproducing the original gamma direction
00541             double cosMcGamAng  = main4mom.vect().unit().dot(primary4mom.vect().unit());
00542             double cosClsGamAng = main4cls.vect().unit().dot(primary4mom.vect().unit());
00543 
00544             if (cosMcGamAng > 1.0 ) cosMcGamAng = 1.0;
00545             if (cosClsGamAng > 1.0) cosClsGamAng = 1.0;
00546 
00547             m_prmMcAngle = acos(cosMcGamAng);
00548             m_prmClsAngle = acos(cosClsGamAng);
00549 
00550             if (m_prmMcAngle > 0.0001) 
00551             {
00552                 //int checkithere = 0;
00553             }
00554         }
00555         else
00556         {
00557             //int ijk = 0;
00558         }
00559     }
00560     
00561     log << MSG::DEBUG << " returning. " << endreq;
00562 
00563     return sc;
00564 }

00234 {
00235     StatusCode sc = StatusCode::SUCCESS;
00236     
00237     MsgStream log(msgSvc(), name());
00238 
00239     if( ValBase::initialize().isFailure()) return StatusCode::FAILURE;
00240   
00241     // get the services    
00242     if( serviceLocator() ) {
00243         if( service("ParticlePropertySvc", m_ppsvc, true).isFailure() ) {
00244             log << MSG::ERROR << "Service [ParticlePropertySvc] not found" << endreq;
00245         }
00246     } else {
00247         return StatusCode::FAILURE;
00248     }
00249 
00250     // TO DO here: gracefully return if tools not located, set up to NOT run the tool
00251     m_mcEvent = 0;
00252     sc = toolSvc()->retrieveTool("McGetEventInfoTool", m_mcEvent);
00253     if (sc.isFailure()) {
00254         log << MSG::INFO << " McGetEventInfoTool not found" << endreq;
00255         log << MSG::INFO << " Will not generate McAnalVals" << endreq;
00256         return StatusCode::SUCCESS;
00257     }
00258 
00259     m_mcTracks = 0;
00260     sc = toolSvc()->retrieveTool("McGetTrackInfoTool", m_mcTracks);
00261     if (sc.isFailure()) {
00262         log << MSG::INFO << " McGetTrackInfoTool not found!" << endreq;
00263         log << MSG::INFO << " Will not generate McAnalVals" << endreq;
00264         return StatusCode::SUCCESS;
00265     }
00266     
00267     // load up the map
00268 
00269         addItem("McaNumCalls",       &m_numCalls);
00270         addItem("McaPrmEnergy",      &m_prmEnergy);
00271         addItem("McaPrmDecEne",      &m_prmDecEne);
00272         addItem("McaDecPosX",        &m_prmPosX);
00273         addItem("McaDecPosY",        &m_prmPosY);
00274         addItem("McaDecPosZ",        &m_prmPosZ);
00275         addItem("McaPrmDecPosX",     &m_prmDecPosX);
00276         addItem("McaPrmDecPosY",     &m_prmDecPosY);
00277         addItem("McaPrmDecPosZ",     &m_prmDecPosZ);
00278         addItem("McaPrmCosDirX",     &m_prmCosDirX);
00279         addItem("McaPrmCosDirY",     &m_prmCosDirY);
00280         addItem("McaPrmCosDirZ",     &m_prmCosDirZ);
00281         addItem("McaPrmDecCosX",     &m_prmDecCosX);
00282         addItem("McaPrmDecCosY",     &m_prmDecCosY);
00283         addItem("McaPrmDecCosZ",     &m_prmDecCosZ);
00284         addItem("McaPrmNDghtrs",     &m_prmNDghtrs);
00285         addItem("McaPrmDecCode",     &m_prmDecCode);
00286 
00287 
00288     addItem("McaPrmNumSecndry",  &m_prmNumSecndry);
00289     addItem("McaPrmNumAsscted",  &m_prmNumAsscted);
00290     addItem("McaPrmTrkPattern",  &m_prmTrkPattern);
00291     addItem("McaPrmMcAngle",     &m_prmMcAngle);
00292     addItem("McaPrmClsAngle",    &m_prmClsAngle);
00293 
00294     addItem("McaPrmDecLayer",    &m_prmDecLayer);
00295     addItem("McaPrmDecInCnv",    &m_prmDecInCnv);
00296 
00297     addItem("McaScd1PartType",   &m_scd1Type);
00298     addItem("McaScd1FirstLyr",   &m_scd1FirstLyr);
00299     addItem("McaScd1LastLyr",    &m_scd1LastLyr);
00300     addItem("McaScd1Energy",     &m_scd1Energy);
00301     addItem("McaScd1CosDirX",    &m_scd1CosDirX);
00302     addItem("McaScd1CosDirY",    &m_scd1CosDirY);
00303     addItem("McaScd1CosDirZ",    &m_scd1CosDirZ);
00304     addItem("McaScd1NHits",      &m_scd1NHits);
00305     addItem("McaScd1NClstrs",    &m_scd1NClstrs);
00306     addItem("McaScd1NGaps",      &m_scd1NGaps);
00307     addItem("McaScd11stGapSz",   &m_scd11stGapSz);
00308     addItem("McaScd1NHits2Gp",   &m_scd1NHits2Gp);
00309     addItem("McaScd1McTrkRms",   &m_scd1McTrkRms);
00310     addItem("McaScd1Ang2Gam",    &m_scd1Ang2Gam);
00311     addItem("McaScd1Cls2Gam",    &m_scd1Cls2Gam);
00312     addItem("McaScd1ELastHit",   &m_scd1ELastHit);
00313     addItem("McaScd1RadEloss",   &m_scd1RadELoss);
00314     addItem("McaScd1DeltaRay",   &m_scd1DeltaRay);
00315     addItem("McaScd1NBrems",     &m_scd1NBrems);
00316     addItem("McaScd1NDeltas",    &m_scd1NDeltas);
00317     addItem("McaScd1NDeltaHt",   &m_scd1NDeltaHt);
00318     addItem("McaScd1AveRange",   &m_scd1AveRange);
00319     addItem("McaScd1MaxRange",   &m_scd1MaxRange);
00320 
00321     addItem("McaScd2PartType",   &m_scd2Type);
00322     addItem("McaScd2FirstLyr",   &m_scd2FirstLyr);
00323     addItem("McaScd2LastLyr",    &m_scd2LastLyr);
00324     addItem("McaScd2Energy",     &m_scd2Energy);
00325     addItem("McaScd2CosDirX",    &m_scd2CosDirX);
00326     addItem("McaScd2CosDirY",    &m_scd2CosDirY);
00327     addItem("McaScd2CosDirZ",    &m_scd2CosDirZ);
00328     addItem("McaScd2NHits",      &m_scd2NHits);
00329     addItem("McaScd2NClstrs",    &m_scd2NClstrs);
00330     addItem("McaScd2NGaps",      &m_scd2NGaps);
00331     addItem("McaScd21stGapSz",   &m_scd21stGapSz);
00332     addItem("McaScd2NHits2Gp",   &m_scd2NHits2Gp);
00333     addItem("McaScd2McTrkRms",   &m_scd2McTrkRms);
00334     addItem("McaScd2Ang2Gam",    &m_scd2Ang2Gam);
00335     addItem("McaScd2Cls2Gam",    &m_scd2Cls2Gam);
00336     addItem("McaScd2ELastHit",   &m_scd2ELastHit);
00337     addItem("McaScd2RadEloss",   &m_scd2RadELoss);
00338     addItem("McaScd2DeltaRay",   &m_scd2DeltaRay);
00339     addItem("McaScd2NBrems",     &m_scd2NBrems);
00340     addItem("McaScd2NDeltas",    &m_scd2NDeltas);
00341     addItem("McaScd2NDeltaHt",   &m_scd2NDeltaHt);
00342     addItem("McaScd2AveRange",   &m_scd2AveRange);
00343     addItem("McaScd2MaxRange",   &m_scd2MaxRange);
00344 
00345     zeroVals();
00346     
00347     return sc;
00348 }

McAnalValsTool Class Reference

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


Public Member Functions
	McAnalValsTool (const std::string &type, const std::string &name, const IInterface *parent)
virtual	~McAnalValsTool ()
StatusCode	initialize ()
StatusCode	calculate ()
	calculate all values; implemented by each XxxValsTool
Private Member Functions
void	calcMcAngleInfo (const Event::McParticle mcPart, const Event::McParticle mcGamma, double &ang2Gam, double &cls2Gam, double &mcTrkRms, double &prtType, double &prtEnergy, double &prtCosDirX, double &prtCosDirY, double &prtCosDirZ, double &prtNHits, double &prtNClstrs, double &prtNGaps, double &prt1stGapSz, double &prtNHits2Gp)
	For determing MC tuple values.
void	calcMcEnergyInfo (const Event::McParticle *mcPart, double &lastHitE, double &radLossE, double &dltaRayE, double &nBrems, double &nDeltas, double &nDeltaHt, double &aveRange, double &maxRange)
Private Attributes
double	m_numCalls
double	m_prmEnergy
double	m_prmPosX
double	m_prmPosY
double	m_prmPosZ
double	m_prmCosDirX
double	m_prmCosDirY
double	m_prmCosDirZ
double	m_prmDecEne
double	m_prmDecPosX
double	m_prmDecPosY
double	m_prmDecPosZ
double	m_prmDecCosX
double	m_prmDecCosY
double	m_prmDecCosZ
double	m_prmNDghtrs
double	m_prmDecCode
double	m_prmDecLayer
double	m_prmDecInCnv
double	m_prmNumSecndry
double	m_prmNumAsscted
double	m_prmTrkPattern
double	m_prmMcAngle
double	m_prmClsAngle
double	m_scd1Type
double	m_scd1FirstLyr
double	m_scd1LastLyr
double	m_scd1Energy
double	m_scd1CosDirX
double	m_scd1CosDirY
double	m_scd1CosDirZ
double	m_scd1NHits
double	m_scd1NClstrs
double	m_scd1NGaps
double	m_scd11stGapSz
double	m_scd1NHits2Gp
double	m_scd1McTrkRms
double	m_scd1Ang2Gam
double	m_scd1Cls2Gam
double	m_scd1ELastHit
double	m_scd1RadELoss
double	m_scd1DeltaRay
double	m_scd1NBrems
double	m_scd1NDeltas
double	m_scd1NDeltaHt
double	m_scd1AveRange
double	m_scd1MaxRange
double	m_scd2Type
double	m_scd2FirstLyr
double	m_scd2LastLyr
double	m_scd2Energy
double	m_scd2CosDirX
double	m_scd2CosDirY
double	m_scd2CosDirZ
double	m_scd2NHits
double	m_scd2NClstrs
double	m_scd2NGaps
double	m_scd21stGapSz
double	m_scd2NHits2Gp
double	m_scd2McTrkRms
double	m_scd2Ang2Gam
double	m_scd2Cls2Gam
double	m_scd2ELastHit
double	m_scd2RadELoss
double	m_scd2DeltaRay
double	m_scd2NBrems
double	m_scd2NDeltas
double	m_scd2NDeltaHt
double	m_scd2AveRange
double	m_scd2MaxRange
IParticlePropertySvc *	m_ppsvc
IMcGetEventInfoTool *	m_mcEvent
IMcGetTrackInfoTool *	m_mcTracks