Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

BesSensitiveManager Class Reference

#include <BesSensitiveManager.hh>

List of all members.

Public Member Functions

void AddSensitiveDetector (BesSensitiveDetector *detector)
void AddSensitiveDetector (BesSensitiveDetector *detector)
void BeginOfTrack (const G4Track *track)
void BeginOfTrack (const G4Track *track)
void BeginOfTruthEvent (const G4Event *)
void BeginOfTruthEvent (const G4Event *)
 BesSensitiveManager ()
 BesSensitiveManager ()
G4bool CheckDecayTrack (const G4Track *)
G4bool CheckDecayTrack (const G4Track *)
G4int CheckType (const HepMC::GenEvent *hepmcevt)
G4int CheckType (const HepMC::GenEvent *hepmcevt)
void ClearEvent ()
void ClearEvent ()
void EndOfTrack (const G4Track *track, G4TrackingManager *)
void EndOfTrack (const G4Track *track, G4TrackingManager *)
void EndOfTruthEvent (const G4Event *)
void EndOfTruthEvent (const G4Event *)
G4int GetCurrentTrackIndex () const
G4int GetCurrentTrackIndex () const
void GetDaughterVertexes (BesTruthTrack *pTrack, std::vector< int > *vDau)
void GetDaughterVertexes (BesTruthTrack *pTrack, std::vector< int > *vDau)
G4int GetNumberTracks () const
G4int GetNumberTracks () const
G4int GetNumberVertices () const
G4int GetNumberVertices () const
std::vector< BesTruthTrack * > * GetTrackList ()
std::vector< BesTruthTrack * > * GetTrackList ()
std::vector< BesTruthVertex * > * GetVertexList ()
std::vector< BesTruthVertex * > * GetVertexList ()
G4bool MatchDaughterTrack (const G4Track *)
G4bool MatchDaughterTrack (const G4Track *)
G4bool MatchVertex (G4int vIndex, std::vector< int > *vDau)
G4bool MatchVertex (G4int vIndex, std::vector< int > *vDau)
void SaveParticlesFromGenerator ()
void SaveParticlesFromGenerator ()
void SetLogLevel (G4int level)
void SetLogLevel (G4int level)
void SetVertex0 (const G4Event *)
void SetVertex0 (const G4Event *)
void UpdatePrimaryTrack (const G4Track *)
void UpdatePrimaryTrack (const G4Track *)
void UpdateVertex (BesTStats, const G4Track *)
void UpdateVertex (BesTStats, const G4Track *)
 ~BesSensitiveManager ()
 ~BesSensitiveManager ()

Static Public Member Functions

BesSensitiveManager * GetSensitiveManager ()
BesSensitiveManager * GetSensitiveManager ()

Protected Member Functions

BesTStats FollowTrack (const G4Track *track)
BesTStats FollowTrack (const G4Track *track)
void MakeNewTrack (BesTStats &stat, const G4Track *track)
void MakeNewTrack (BesTStats &stat, const G4Track *track)

Protected Attributes

std::vector< BesTStats > chain
std::vector< BesTStats > chain
std::vector< BesSensitiveDetector * > clients
std::vector< BesSensitiveDetector * > clients
std::vector< BesSensitiveDetector
* >::iterator 		iter
std::vector< BesSensitiveDetector
* >::iterator 		iter
std::vector< BesSensitiveDetector
* >::iterator 		iter_end
std::vector< BesSensitiveDetector
* >::iterator 		iter_end
G4int m_count
HepMC::GenEvent * m_hepmcevt
HepMC::GenEvent * m_hepmcevt
G4int m_logLevel
G4ThreeVector m_pos0
G4double m_t0
G4int m_trackFlag
G4int m_trackIndex
std::vector< BesTruthTrack * > * m_trackList
std::vector< BesTruthTrack * > * m_trackList
std::vector< BesTruthVertex * > * m_vertexList
std::vector< BesTruthVertex * > * m_vertexList

Static Protected Attributes

BesSensitiveManager * m_sensitiveManager
BesSensitiveManager * m_sensitiveManager = 0

Constructor & Destructor Documentation

BesSensitiveManager::BesSensitiveManager (   )  [inline]
  

00075                        :m_trackIndex(-1),m_trackFlag(0),m_count(0) 
00076   {
00077     m_hepmcevt = 0; 
00078     if(m_sensitiveManager) 
00079       G4Exception("BesSensitiveManager constructed twice.");
00080     m_sensitiveManager = this ;
00081   }

BesSensitiveManager::~BesSensitiveManager (   )  [inline]
  

00083 {;}

BesSensitiveManager::BesSensitiveManager (   )  [inline]
  

00075                        :m_trackIndex(-1),m_trackFlag(0),m_count(0) 
00076   {
00077     m_hepmcevt = 0; 
00078     if(m_sensitiveManager) 
00079       G4Exception("BesSensitiveManager constructed twice.");
00080     m_sensitiveManager = this ;
00081   }

BesSensitiveManager::~BesSensitiveManager (   )  [inline]
  

00083 {;}

Member Function Documentation

void BesSensitiveManager::AddSensitiveDetector (  BesSensitiveDetector *  detector  )  [inline]
  

00089   {
00090     clients.push_back( detector );
00091   }

void BesSensitiveManager::AddSensitiveDetector (  BesSensitiveDetector *  detector  )  [inline]
  

00089   {
00090     clients.push_back( detector );
00091   }

void BesSensitiveManager::BeginOfTrack (  const G4Track *  track  ) 
 

void BesSensitiveManager::BeginOfTrack (  const G4Track *  track  ) 
  

00250 {
00251   m_trackFlag = 0; //if m_trackFlag=1, this track is from generator
00252 
00253   // Append our track in the parentage history
00254   chain.push_back(FollowTrack(track));
00255 
00256   // Get TruthTrack index
00257   //int i = chain.size();
00258 
00259   /*for(;;) 
00260   {
00261     if(--i<0) G4cerr<<"::parents are corrupted"<<G4endl;
00262     m_trackIndex = chain[i].trackIndex;
00263     if (m_trackIndex != -1) break;
00264   }*/
00265   
00266   // Invoke clients
00267   //
00268   iter=clients.begin();
00269   iter_end=clients.end();
00270   while(iter!=iter_end) 
00271   {
00272     (*iter)->BeginOfTrack( track );
00273     iter++;
00274   }
00275 }

void BesSensitiveManager::BeginOfTruthEvent (  const G4Event *   ) 
 

void BesSensitiveManager::BeginOfTruthEvent (  const G4Event *   ) 
  

00044 {
00045   // Allocate our lists. We'll own these until EndOfEvent
00046 
00047   m_trackList  = new std::vector<BesTruthTrack*>;
00048   m_vertexList = new std::vector<BesTruthVertex*>;
00049   
00050   SetVertex0(evt);
00051 
00052   m_count = 0;
00053   SaveParticlesFromGenerator();
00054   m_count += m_trackList->size();
00055 
00056   // Tell clients
00057   iter=clients.begin();
00058   iter_end=clients.end();
00059   while( iter != iter_end ) 
00060   {
00061     (*iter)->BeginOfTruthEvent(evt);
00062     iter++;
00063   }
00064 }

G4bool BesSensitiveManager::CheckDecayTrack (  const G4Track *   ) 
 

G4bool BesSensitiveManager::CheckDecayTrack (  const G4Track *   ) 
  

00587 { 
00588   G4int flag = 0;
00589   G4int trackID = track->GetTrackID();
00590   G4int parentID = track->GetParentID();
00591   G4int pdgcode = track->GetDefinition()->GetPDGEncoding();
00592   G4ThreeVector p3 = track->GetMomentum();
00593   HepLorentzVector  p4(track->GetMomentum(), track->GetTotalEnergy());
00594   
00595   BesTruthTrack* truthTrack=0;
00596   G4int nTrack = m_trackList->size();
00597   G4int parentTrackIndex=-99;
00598   G4int terminalVertexIndex=-99;
00599   for(int i=0;i<nTrack;i++)
00600   {
00601     truthTrack = (*m_trackList)[i];
00602     if(truthTrack->GetG4TrackId() == parentID)
00603     {
00604       //parentTrackIndex = i;
00605       parentTrackIndex = truthTrack->GetIndex();
00606       if(truthTrack->GetTerminalVertex())
00607       {
00608         terminalVertexIndex = truthTrack->GetTerminalVertex()->GetIndex();
00609         if(m_logLevel<=4)
00610         {
00611           G4cout<<"parentTrackIndex:"<<parentTrackIndex<<" "
00612                 <<"parent terminate at vertex: "<<terminalVertexIndex<<G4endl;
00613         }
00614         break;
00615       }
00616     }
00617   }
00618   if(parentTrackIndex==-99 || terminalVertexIndex==-99)
00619   {
00620     G4cout<<"MatchDecayError!"<<G4endl;
00621     return false;
00622   }
00623  
00624   //match decayed track with p3
00625   G4bool matchDauFlag=false; 
00626   if(terminalVertexIndex>0)
00627     matchDauFlag = MatchDaughterTrack(track);
00628   if(matchDauFlag)
00629     return true;
00630 
00631   //match decayed track with vertex index
00632   G4ThreeVector minDiffP(1e99,1e99,1e99);
00633   G4int truthI=-9;
00634   for(int i=0;i<nTrack;i++)
00635   {
00636     truthTrack = (*m_trackList)[i];
00637     G4String source = truthTrack->GetSource();
00638     G4int pdgcode_tru = truthTrack->GetPDGCode();
00639     if(pdgcode_tru==-22)  pdgcode_tru = 22;
00640     if(truthTrack->GetVertex()->GetIndex() == terminalVertexIndex &&
00641        pdgcode_tru == pdgcode &&source=="FromGenerator"&&truthTrack->NotFound())
00642     {
00643       HepLorentzVector tp4 = truthTrack->GetP4();
00644       G4ThreeVector tp3(tp4.x(), tp4.y(), tp4.z());
00645       G4ThreeVector diffP = (p3-tp3).mag2();
00646       if(diffP<minDiffP)
00647       {  minDiffP = diffP;  truthI = i; }
00648     }
00649   }
00650   
00651   if(truthI!=-9)
00652   {
00653     if(m_logLevel<=4)
00654       G4cout<<"MatchDecayedTrackWithVertexIndex, trackIndex:"<<truthI<<" ";
00655     m_trackIndex = truthI;
00656     truthTrack=(*m_trackList)[truthI];
00657     //truthTrack->SetP4(p4);
00658     truthTrack->SetPDGCharge(track->GetDefinition()->GetPDGCharge());
00659     truthTrack->SetParticleName(track->GetDefinition()->GetParticleName());
00660     truthTrack->SetG4TrackId(track->GetTrackID());
00661     truthTrack->Found();
00662     G4int size = truthTrack->GetDaughterIndexes().size();
00663     if(size>0)
00664     {
00665       G4int minDau = (truthTrack->GetDaughterIndexes())[0];
00666       G4int maxDau = (truthTrack->GetDaughterIndexes())[size-1];
00667       if(m_logLevel<=4)
00668         G4cout<<"daughters: "<<minDau<<" "<<maxDau<<G4endl;
00669     }
00670     else
00671     {
00672       if(m_logLevel<=4)
00673       G4cout<<G4endl;
00674     }
00675     return true; 
00676   }//match decay track with vertex index
00677 
00678   //match decay track, if its mother is not defined in geant4
00679   if(chain.back().vertices.size()<1)
00680   {
00681     BesTruthTrack* pTrack = (*m_trackList)[parentTrackIndex];
00682     std::vector<int>* vDau = new std::vector<int>;
00683     GetDaughterVertexes(pTrack, vDau);
00684     G4int sizev = vDau->size();
00685     if(sizev>0)
00686     {
00687       G4ThreeVector minDiffP(1e99,1e99,1e99);
00688       for(int i=0;i<nTrack;i++)
00689       {
00690         truthTrack = (*m_trackList)[i];
00691         G4int vIndex = truthTrack->GetVertex()->GetIndex();
00692         G4int pdgT = truthTrack->GetPDGCode();
00693         if(pdgT==-22)  pdgT = 22;
00694         if(pdgT == pdgcode)
00695         {
00696           G4bool matchFlag = MatchVertex(vIndex, vDau);
00697           if(matchFlag)
00698           {
00699             HepLorentzVector tp4 = truthTrack->GetP4();
00700             G4ThreeVector tp3(tp4.x(), tp4.y(), tp4.z());
00701             G4ThreeVector diffP = (p3-tp3).mag2();
00702             if(diffP<minDiffP)
00703             {  minDiffP = diffP; truthI = i;}
00704           }
00705         }
00706       }
00707     }
00708     if(vDau)
00709       delete vDau;
00710   }
00711   if(truthI!=-9)
00712   {
00713     if(m_logLevel<=4)
00714       G4cout<<"MatchDecayedTrackWithTrueMother"<<G4endl;
00715     m_trackIndex = truthI;
00716     truthTrack=(*m_trackList)[truthI];
00717     //truthTrack->SetP4(p4);
00718     truthTrack->SetPDGCharge(track->GetDefinition()->GetPDGCharge());
00719     truthTrack->SetParticleName(track->GetDefinition()->GetParticleName());
00720     truthTrack->SetG4TrackId(track->GetTrackID());
00721     truthTrack->Found();
00722     truthTrack->GetVertex()->SetPosition(track->GetPosition());
00723     truthTrack->GetVertex()->SetTime(track->GetGlobalTime());
00724     return true;
00725   }
00726   return false;   
00727 }

G4int BesSensitiveManager::CheckType (  const HepMC::GenEvent *  hepmcevt  ) 
 

G4int BesSensitiveManager::CheckType (  const HepMC::GenEvent *  hepmcevt  ) 
  

00080 {
00081   G4int flag =0;
00082   for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
00083           vitr != hepmcevt->vertices_end(); ++vitr ) { // loop for vertex ...
00084     for (HepMC::GenVertex::particle_iterator
00085              pitr= (*vitr)->particles_begin(HepMC::children);
00086                 pitr != (*vitr)->particles_end(HepMC::children); ++pitr) {
00087       if((*pitr)->end_vertex())
00088         {flag = 1; break;}
00089     }
00090     if(flag) break;
00091   }
00092   if(m_logLevel <= 4)
00093   {
00094     if(flag == 0)
00095       G4cout<<G4endl<<"generator is GENBES!"<<G4endl;
00096     else
00097       G4cout<<G4endl<<"generator is not GENBES!"<<G4endl;
00098   }
00099   return flag;
00100 }

void BesSensitiveManager::ClearEvent (   ) 
 

void BesSensitiveManager::ClearEvent (   ) 
  

00232 {
00233   if(m_trackList)
00234   {
00235     for(size_t i=0;i<m_trackList->size();i++)
00236       delete (*m_trackList)[i];
00237     m_trackList->clear();
00238     delete m_trackList;
00239   }
00240   if(m_vertexList)
00241   {
00242     for(size_t i=0;i<m_vertexList->size();i++)
00243       delete (*m_vertexList)[i];
00244     m_vertexList->clear();
00245     delete m_vertexList;
00246   }
00247 }

void BesSensitiveManager::EndOfTrack (  const G4Track *  track,
G4TrackingManager * 
) 
 

void BesSensitiveManager::EndOfTrack (  const G4Track *  track,
G4TrackingManager * 
) 
  

00277 {
00278   if(chain.back().savedByDefault==true)
00279   { 
00280     BesTStats &stat = chain.back();
00281     G4int endVtxFlag = 0;
00282     if(m_trackFlag==1)
00283     {
00284       BesTruthTrack* truthTrack = (*m_trackList)[m_trackIndex];
00285       if(truthTrack->GetTerminalVertex())
00286         { UpdateVertex(stat,track); endVtxFlag = 1;}
00287     } 
00288     if(endVtxFlag == 0)
00289     {
00290       // Tracks saved to BesTruthTrack require additional work
00291       G4StepPoint *finalStep = track->GetStep()->GetPostStepPoint();
00292       G4StepStatus stepStatus = finalStep->GetStepStatus();
00293       if (track->GetNextVolume() != 0 ||
00294               (stepStatus != fGeomBoundary && stepStatus != fWorldBoundary) )
00295       {
00296         if(m_logLevel<=4)
00297           G4cout<<"Particle died. make a terminal vertex: ";
00298         // Particle died. We always make a terminal vertex
00299         const G4ThreeVector pos(finalStep->GetPosition());
00300         G4int vindex = m_vertexList->size();
00301         if(m_logLevel<=4)
00302           G4cout<<vindex<<G4endl;
00303         stat.vertices.push_back(vindex);
00304         BesTruthVertex *newVertex = new BesTruthVertex();
00305         newVertex->SetPosition(pos);
00306         newVertex->SetTime( finalStep->GetGlobalTime());
00307         if(m_trackFlag==1)
00308           newVertex->SetParentTrack( (*m_trackList)[m_trackIndex] );
00309         else
00310           newVertex->SetParentTrack( m_trackList->back() );
00311         newVertex->SetTerminal( true );
00312         newVertex->SetIndex( vindex );
00313         
00314         //set minDaughter index equal to m_count
00315         newVertex->SetMinDau(m_count);  
00316         
00317         //Get how many decayed daughters this track has, if no, 0 is set.
00318         G4TrackVector* trackVector = trackingManager->GimmeSecondaries();
00319         G4int nSecon = trackVector->size();
00320         G4Track* seconTrack;
00321         G4int nDau=0;
00322         if(nSecon>0)
00323         {
00324           for(G4int i=0;i<nSecon;i++)
00325           {
00326             seconTrack = (*trackVector)[i];
00327             G4String processName = seconTrack->GetCreatorProcess()->GetProcessName();
00328             if(processName == "Decay")
00329               nDau++;
00330           }
00331         }
00332         m_vertexList->push_back( newVertex );
00333   
00334         m_count += nDau;
00335         
00336         // Give this terminal vertex to the track
00337         //(*m_trackList)[stat.trackIndex]->SetTerminalVertex( m_vertexList->back() );
00338         // for particles from generator, 
00339         // m_trackList->back() may not be current track.
00340 
00341         if(m_trackFlag==1)
00342           (*m_trackList)[m_trackIndex]->SetTerminalVertex( m_vertexList->back() );
00343         else
00344           m_trackList->back()->SetTerminalVertex( m_vertexList->back() );
00345       }
00346     }
00347   } 
00348   // Invoke clients
00349   iter=clients.begin();
00350   iter_end=clients.end();
00351   while( iter!=iter_end) {
00352     (*iter)->EndOfTrack( track );
00353     iter++;
00354   }
00355 }

void BesSensitiveManager::EndOfTruthEvent (  const G4Event *   ) 
 

void BesSensitiveManager::EndOfTruthEvent (  const G4Event *   ) 
  

00220 {
00221   // Tell clients
00222   iter=clients.begin();
00223   iter_end=clients.end();
00224   while( iter != iter_end ) 
00225   {
00226     (*iter)->EndOfTruthEvent(evt);
00227     iter++;
00228   }
00229 }

BesTStats BesSensitiveManager::FollowTrack (  const G4Track *  track  )  [protected]
 

BesTStats BesSensitiveManager::FollowTrack (  const G4Track *  track  )  [protected]
  

00502 {
00503   // Some default stats for this track
00504   BesTStats stat( track->GetTrackID(),
00505                   HepLorentzVector( track->GetMomentum(), track->GetTotalEnergy() ),
00506                   track->GetGlobalTime());
00507 
00508   // Check immediate parent
00509   G4int parentId = track->GetParentID();
00510   G4int pdg = track->GetDefinition()->GetPDGEncoding();
00511   G4String particleName = track->GetDefinition()->GetParticleName();
00512   G4ThreeVector trackPos = track->GetPosition();
00513   G4double diffT = m_t0-track->GetGlobalTime();
00514   G4ThreeVector diffPos = (m_pos0-trackPos).mag2();
00515 
00516   if (parentId == 0) 
00517   {
00518     if(m_logLevel<=4)
00519     {
00520       G4cout<<G4endl;
00521       G4cout<<"trackId:"<<track->GetTrackID()<<" ";
00522       G4cout<<"parentId:"<<parentId<<" ";
00523       G4cout<<"pdg:"<<pdg<<" ";
00524       G4cout<<"name:"<<particleName<<" ";
00525       G4cout<<"timeDiff:"<<diffT<<" ";
00526       G4cout<<"posDiff:"<<diffPos<<G4endl;
00527     }
00528 
00529     // Primary particle: wipe decay chain
00530     chain.clear();
00531     // Always save
00532     stat.savedByDefault = true;
00533     // Make TruthTrack
00534     if(m_hepmcevt==0)
00535       MakeNewTrack( stat, track );  
00536     else
00537     {
00538       UpdatePrimaryTrack(track);
00539       m_trackFlag = 1;
00540     }
00541     return stat;
00542   }
00543         
00544   // Not primary particle: trim down chain until
00545   // we uncover the parent. The parent must be there!
00546   for(;;) 
00547   {
00548     if (chain.back().G4index == parentId) break;
00549     chain.pop_back();
00550   }
00551 
00552   // This track is a candidate for saving by default
00553   // only if its parent was saved by default
00554   if (chain.back().savedByDefault) 
00555   {
00556     // There are three ways particles are saved by default:
00557     //   1. if they are a primary particle
00558     //   2. if they are the result of a decay
00559     //
00560     G4String processName = track->GetCreatorProcess()->GetProcessName();
00561     if (processName=="Decay") 
00562     {
00563       if(m_logLevel<=4)
00564       {
00565         G4cout<<G4endl;
00566         G4cout<<"trackId: "<<track->GetTrackID()<<" ";
00567         G4cout<<"parentId: "<<parentId<<" ";
00568         G4cout<<"pdg: "<<pdg<<" ";
00569         G4cout<<"pname: "<<particleName<<G4endl;
00570       }
00571       stat.savedByDefault = true;
00572 
00573       if(CheckDecayTrack(track))
00574         m_trackFlag = 1;
00575       else
00576         MakeNewTrack( stat, track);
00577       return stat;
00578     }
00579   }
00580   
00581   //now this track will not be saved as BesTruthTrack
00582   stat.savedByDefault = false;
00583   return stat;
00584 }

G4int BesSensitiveManager::GetCurrentTrackIndex (   )  const [inline]
  

00093 { return m_trackIndex; }

G4int BesSensitiveManager::GetCurrentTrackIndex (   )  const [inline]
  

00093 { return m_trackIndex; }

void BesSensitiveManager::GetDaughterVertexes (  BesTruthTrack *  pTrack,
std::vector< int > *  vDau
) 
 

void BesSensitiveManager::GetDaughterVertexes (  BesTruthTrack *  pTrack,
std::vector< int > *  vDau
) 
  

00730 {
00731   G4int size = pTrack->GetDaughterIndexes().size();
00732   if(size>0)
00733   {
00734     G4int minDau = (pTrack->GetDaughterIndexes())[0];
00735     G4int maxDau = (pTrack->GetDaughterIndexes())[size-1];
00736     //note! here, dau<=maxDau, not dau<maxDau
00737     for(G4int dau=minDau;dau<=maxDau;dau++)
00738     {
00739       BesTruthTrack* truthTrack = (*m_trackList)[dau];
00740       if(truthTrack->GetTerminalVertex())
00741       {
00742         vDau->push_back(truthTrack->GetTerminalVertex()->GetIndex());
00743         GetDaughterVertexes(truthTrack,vDau);
00744       } 
00745     }
00746   }
00747 }

G4int BesSensitiveManager::GetNumberTracks (   )  const [inline]
  

00119 { return m_trackList ? m_trackList->size() : 0; }

G4int BesSensitiveManager::GetNumberTracks (   )  const [inline]
  

00119 { return m_trackList ? m_trackList->size() : 0; }

G4int BesSensitiveManager::GetNumberVertices (   )  const [inline]
  

00120 { return m_vertexList ? m_vertexList->size() : 0; }

G4int BesSensitiveManager::GetNumberVertices (   )  const [inline]
  

00120 { return m_vertexList ? m_vertexList->size() : 0; }

BesSensitiveManager* BesSensitiveManager::GetSensitiveManager (   )  [inline, static]
  

00085 {return m_sensitiveManager; }

BesSensitiveManager* BesSensitiveManager::GetSensitiveManager (   )  [inline, static]
  

00085 {return m_sensitiveManager; }

std::vector<BesTruthTrack*>* BesSensitiveManager::GetTrackList (   )  [inline]
  

00094 {return m_trackList;}

std::vector<BesTruthTrack*>* BesSensitiveManager::GetTrackList (   )  [inline]
  

00094 {return m_trackList;}

std::vector<BesTruthVertex*>* BesSensitiveManager::GetVertexList (   )  [inline]
  

00095 {return m_vertexList;}

std::vector<BesTruthVertex*>* BesSensitiveManager::GetVertexList (   )  [inline]
  

00095 {return m_vertexList;}

void BesSensitiveManager::MakeNewTrack (  BesTStats &  stat,
const G4Track *  track
)  [protected]
 

void BesSensitiveManager::MakeNewTrack (  BesTStats &  stat,
const G4Track *  track
)  [protected]
  

00378 {
00379   if (stat.originVertex < 0 && chain.size() > 0) 
00380   {
00381     if(m_logLevel<=4)
00382       G4cout<<"MakeNewTrack for decayed particles,";
00383     //for decayed particles, there may already be a BesTruthVertex that is suitable for this track.
00384     //If so, it's always the terminal vertex of its immediate parent( chain.back() ).
00385 
00386     G4int parentVstat = chain.back().vertices.size();
00387     while( --parentVstat >= 0)
00388     {
00389       G4int vindex = chain.back().vertices[parentVstat];
00390 
00391       G4ThreeVector diff((*m_vertexList)[vindex]->GetPosition());
00392       diff = diff - track->GetPosition();
00393       if (diff.mag2() < 1E-9) 
00394       {
00395               stat.originVertex = vindex;
00396         if(m_logLevel<=4)
00397           G4cout<<"find a vertex:"<<vindex;
00398         (*m_vertexList)[vindex]->AddCurrentDau();
00399         break;
00400       }
00401     }
00402     
00403   }
00404   
00405   if(stat.originVertex < 0 && chain.size() == 0)
00406   {
00407     //for primary track, check if there is already a vertex suitable for it.
00408     G4int nVertex = m_vertexList->size();         
00409     for(G4int i=0;i<nVertex;i++)
00410     {
00411       G4ThreeVector diff((*m_vertexList)[i]->GetPosition());
00412       diff = diff -track->GetPosition();
00413       if(diff.mag2() < 1E-9)
00414       {
00415         stat.originVertex = i; 
00416         if(m_logLevel<=4)
00417           G4cout<<"MakeNewTrack for primary particles,find a vertex:"
00418                 <<i;
00419         break;
00420       }
00421     }
00422   }
00423   
00424   if (stat.originVertex < 0)
00425   {
00426     if(m_logLevel<=4) 
00427       G4cout<<"MakeNewTrack for primary particles,make new vertex:"
00428             <<m_vertexList->size();
00429     // it's a primary track(there's no vertex suitable for it)
00430     // Make a new BesTruthVertex
00431     const G4ThreeVector v(track->GetPosition());
00432                 
00433     stat.originVertex = m_vertexList->size();
00434 
00435     BesTruthVertex *newVertex = new BesTruthVertex();
00436     newVertex->SetPosition(v);
00437     newVertex->SetTime( track->GetGlobalTime());
00438     
00439     //if (chain.size() > 0) {
00440     //  G4int trackIndex = -1;
00441     //  G4int i = chain.size();
00442     //  while(--i>=0 && trackIndex < 0) trackIndex = chain[i].trackIndex;
00443     //  newVertex->SetParentTrack( trackIndex < 0 ? 0 : (*m_trackList)[trackIndex] );
00444     //}
00445     
00446     newVertex->SetParentTrack( 0 );
00447     newVertex->SetTerminal( false );
00448     newVertex->SetIndex( stat.originVertex );
00449     if(track->GetCreatorProcess())
00450       newVertex->SetProcessName(track->GetCreatorProcess()->GetProcessName());
00451     else
00452       newVertex->SetProcessName("generator");    
00453     
00454     m_vertexList->push_back( newVertex );
00455   }
00456         
00457   // Now, finally make our new BesTruthTrack
00458   // We'll leave the assignment of terminalVertex until
00459   // we know what happens to this track         
00460   BesTruthTrack *newTrack = new BesTruthTrack();
00461   newTrack->SetP4( stat.p4 );
00462   newTrack->SetPDGCode(track->GetDefinition()->GetPDGEncoding());
00463   newTrack->SetPDGCharge(track->GetDefinition()->GetPDGCharge());
00464   newTrack->SetParticleName(track->GetDefinition()->GetParticleName());
00465   newTrack->SetSource("FromG4");
00466   BesTruthVertex* vertex = (*m_vertexList)[stat.originVertex];
00467   newTrack->SetVertex(vertex);
00468   
00469   if(track->GetParentID()==0)
00470   {
00471     stat.trackIndex = m_count;
00472     m_count++;
00473   }
00474   else
00475     stat.trackIndex = vertex->GetMinDau() + vertex->GetCurrentDau()-1;
00476   
00477   newTrack->SetIndex( stat.trackIndex );
00478   m_trackIndex = stat.trackIndex;
00479   if(m_logLevel<=4)
00480     G4cout<<"  m_trackIndex:"<<m_trackIndex<<G4endl;
00481   
00482   newTrack->SetG4TrackId( track->GetTrackID()); 
00483   m_trackList->push_back( newTrack );
00484   //tell the parent track its daughter track index
00485   BesTruthTrack* parent = newTrack->GetVertex()->GetParentTrack();
00486   if(parent)
00487   {
00488     parent->AddDaughterIndex(stat.trackIndex);
00489     if(m_logLevel<=4)
00490       G4cout<<"add this daughter to:"<<parent->GetIndex()<<G4endl;
00491   }
00492 }

G4bool BesSensitiveManager::MatchDaughterTrack (  const G4Track *   ) 
 

G4bool BesSensitiveManager::MatchDaughterTrack (  const G4Track *   ) 
  

00764 {
00765   G4int flag=0;
00766   G4int pdgcode = track->GetDefinition()->GetPDGEncoding();
00767   G4ThreeVector p = track->GetMomentum();
00768   HepLorentzVector  p4(track->GetMomentum(), track->GetTotalEnergy());
00769   G4ThreeVector pos = track->GetPosition();
00770   G4double time = track->GetGlobalTime();
00771   BesTruthTrack* truthTrack=0;
00772   G4int size = m_trackList->size();
00773   if(size>0)
00774   { 
00775     for(G4int i=0;i<size;i++)
00776     {
00777       truthTrack=(*m_trackList)[i];
00778       G4String source = truthTrack->GetSource(); 
00779       G4int pdgcode_tru = truthTrack->GetPDGCode();
00780       if(pdgcode_tru==-22)  pdgcode_tru = 22;
00781       HepLorentzVector tp4 = truthTrack->GetP4();
00782       G4ThreeVector tp3(tp4.x(), tp4.y(), tp4.z());
00783       G4ThreeVector diffP = (p-tp3).mag2();
00784       if(source=="FromGenerator"&&pdgcode_tru==pdgcode&&diffP<1e-9&&truthTrack->NotFound())
00785       {
00786         m_trackIndex = i;
00787         //truthTrack->SetP4(p4);
00788         truthTrack->SetPDGCharge(track->GetDefinition()->GetPDGCharge());
00789         truthTrack->SetParticleName(track->GetDefinition()->GetParticleName());
00790         truthTrack->SetG4TrackId(track->GetTrackID());
00791         truthTrack->Found();
00792         //truthTrack->GetVertex()->SetPosition(pos);
00793         //truthTrack->GetVertex()->SetTime(time);
00794         flag=1;
00795         if(m_logLevel<=4)
00796         {
00797           G4cout<<"MatchDecayedTrackWithP3!"<<"trackIndex:"<<m_trackIndex
00798                 <<" pdgcode:"<<pdgcode<<G4endl;
00799           G4cout<<"parentIndex:"<<truthTrack->GetVertex()->GetParentTrack()->GetIndex()
00800                 <<"  PDGCode:"<<truthTrack->GetVertex()->GetParentTrack()->GetPDGCode()<<G4endl; 
00801         }
00802         return true;
00803       }
00804     }
00805   }
00806   if(flag!=1)
00807     return false;
00808 }

G4bool BesSensitiveManager::MatchVertex (  G4int  vIndex,
std::vector< int > *  vDau
) 
 

G4bool BesSensitiveManager::MatchVertex (  G4int  vIndex,
std::vector< int > *  vDau
) 
  

00750 {
00751   G4int size = vDau->size();
00752   if(size>0)
00753   {
00754     for(G4int i=0;i<size;i++)
00755     {
00756       if(vIndex==(*vDau)[i])
00757         return true;
00758     }
00759   }
00760   return false;
00761 }

void BesSensitiveManager::SaveParticlesFromGenerator (   ) 
 

void BesSensitiveManager::SaveParticlesFromGenerator (   ) 
  

00103 {
00104   IDataProviderSvc* p_evtSvc=0;
00105   ISvcLocator* svcLocator = Gaudi::svcLocator();
00106   StatusCode sc=svcLocator->service("EventDataSvc", p_evtSvc);
00107   if (sc.isFailure())
00108       std::cout<<"BesHepMCInterface could not access EventDataSvc!!"<<std::endl;
00109 
00110   SmartDataPtr<McGenEventCol> mcCollptr( p_evtSvc, "/Event/Gen");
00111 
00112   if ( mcCollptr == 0 )
00113     std::cout << "no McGenEventCollection found."  << std::endl;
00114 
00115   else 
00116   {
00117     McGenEventCol::const_iterator it = mcCollptr->begin();
00118     McGenEvent* mcEvent = (McGenEvent* ) (*it);
00119     m_hepmcevt = mcEvent->getGenEvt();
00120     G4int typeFlag = CheckType(m_hepmcevt); 
00121     
00122     for(HepMC::GenEvent::vertex_const_iterator vitr= m_hepmcevt->vertices_begin();
00123               vitr != m_hepmcevt->vertices_end(); ++vitr ) 
00124     {
00125       G4int barcodeVtx = (*vitr)->barcode();
00126       if(m_logLevel<=4)
00127         G4cout<<G4endl<<"barcodeVtx:"<<barcodeVtx<<" ";
00128   
00129       G4LorentzVector xvtx= (*vitr)-> position();
00130       G4ThreeVector v(xvtx.x(), xvtx.y(), xvtx.z());
00131       BesTruthVertex* newVertex = new BesTruthVertex();
00132       newVertex->SetPosition(v);
00133       newVertex->SetTime(xvtx.t()*mm/c_light);
00134       if(m_logLevel<=4)
00135         G4cout<<"xyzt:"<<xvtx.x()<<" "<<xvtx.y()<<" "
00136             <<xvtx.z()<<" "<<xvtx.t()*mm/c_light<<" ";
00137       G4int nTrack = m_trackList->size();
00138       BesTruthTrack* parentTrack=0;
00139       G4int parentTrackIndex= -99;
00140       for(int i=0;i<nTrack;i++)
00141       {
00142         parentTrack = (*m_trackList)[i];
00143         if(parentTrack->GetBarcodeEndVtx() == barcodeVtx)
00144         {
00145           newVertex->SetParentTrack(parentTrack);
00146           parentTrackIndex = parentTrack->GetIndex();
00147           if(m_logLevel<=4)
00148             G4cout<<"parentTrack:"<<parentTrackIndex<<" ";
00149           parentTrack->SetTerminalVertex(newVertex);
00150           //break;
00151         }
00152       }
00153   
00154       G4int vtxFlag=0;
00155       G4int pOut = (*vitr)->particles_out_size();
00156       HepMC::GenVertex::particle_iterator pitr = (*vitr)->particles_begin(HepMC::children);
00157       G4int pdgcode= (*pitr)-> pdg_id();
00158       if(pOut == 1 && pdgcode == -11 && typeFlag==1)
00159         vtxFlag=1;
00160           
00161       if(vtxFlag==0)
00162       {
00163         m_vertexList->push_back(newVertex);
00164         newVertex->SetIndex(m_vertexList->size()-1);
00165         if(m_logLevel<=4)
00166           G4cout<<"vindex:"<<m_vertexList->size()-1<<G4endl;
00167   
00168         for (HepMC::GenVertex::particle_iterator pitr= (*vitr)->particles_begin(HepMC::children); pitr != (*vitr)->particles_end(HepMC::children); ++pitr)
00169         {
00170           G4LorentzVector p= (*pitr)-> momentum();
00171           G4LorentzVector pGeV(p.x()*GeV,p.y()*GeV,p.z()*GeV,p.t()*GeV);
00172           G4int pdgcode = (*pitr)->pdg_id();
00173   
00174           BesTruthTrack* newTrack = new BesTruthTrack;
00175           newTrack->SetP4(pGeV);
00176           newTrack->SetPDGCode(pdgcode);
00177           if(m_logLevel<=4)
00178           {  
00179             G4cout<<"pdg:"<<pdgcode<<" ";
00180             G4cout<<"p:"<<p.x()<<" "<<p.y()<<" "<<p.z()<<" "<<p.t()<<" ";
00181           }
00182           newTrack->SetPDGCharge(99);
00183           newTrack->SetParticleName("unknown");
00184           newTrack->SetVertex(newVertex);
00185           newTrack->SetTerminalVertex(0);
00186           newTrack->SetSource("FromGenerator");
00187         
00188           G4int barcodeEndVtx=0;
00189           if((*pitr)->end_vertex())
00190           {
00191             barcodeEndVtx = (*pitr)->end_vertex()->barcode();
00192             if(m_logLevel<=4)
00193               G4cout<<"endVtx:"<<barcodeEndVtx<<" ";
00194           }
00195           newTrack->SetBarcodeEndVtx(barcodeEndVtx);
00196   
00197           m_trackList->push_back(newTrack);
00198           newTrack->SetIndex(m_trackList->size()-1);
00199           if(m_logLevel<=4)
00200             G4cout<<"index:"<<m_trackList->size()-1<<" ";
00201           //here, parentTrack->GetIndex can't be used,
00202           //use parentTrackIndex instead
00203           if(parentTrackIndex>=0)
00204           {
00205             if(m_logLevel<=4)
00206               G4cout<<"mother:"<<parentTrackIndex; 
00207             (*m_trackList)[parentTrackIndex]->AddDaughterIndex(m_trackList->size()-1);
00208           }
00209           if(m_logLevel<=4)
00210             G4cout<<G4endl;
00211         }         
00212       }
00213     } 
00214   } 
00215 }

void BesSensitiveManager::SetLogLevel (  G4int  level  )  [inline]
  

00124 {m_logLevel = level;}

void BesSensitiveManager::SetLogLevel (  G4int  level  )  [inline]
  

00124 {m_logLevel = level;}

void BesSensitiveManager::SetVertex0 (  const G4Event *   ) 
 

void BesSensitiveManager::SetVertex0 (  const G4Event *   ) 
  

00067 {
00068   //G4int n_vertex = anEvent->GetNumberOfPrimaryVertex();
00069   for( G4int i=0; i<1; i++ )
00070   {
00071     G4PrimaryVertex* primaryVertex = anEvent->GetPrimaryVertex(i);
00072     m_pos0 = primaryVertex->GetPosition();
00073     m_t0 = primaryVertex->GetT0();
00074   }
00075   if(m_logLevel<=4)
00076     G4cout<<"pos0:"<<m_pos0<<"   t0:"<<m_t0<<G4endl;
00077 }

void BesSensitiveManager::UpdatePrimaryTrack (  const G4Track *   ) 
 

void BesSensitiveManager::UpdatePrimaryTrack (  const G4Track *   ) 
  

00811 {
00812   if(m_logLevel<=4)
00813     G4cout<<"UpdatePrimaryTrack! ";
00814   G4int pdgcode = track->GetDefinition()->GetPDGEncoding();
00815   G4ThreeVector p = track->GetMomentum();
00816   G4int nTrack = m_trackList->size();
00817   BesTruthTrack* truthTrack;
00818   G4int flag = 0;
00819   G4ThreeVector minDiffP(1e99,1e99,1e99);
00820   for(int i=0;i<nTrack;i++)
00821   {
00822     truthTrack = (*m_trackList)[i];
00823     HepLorentzVector tp4 = truthTrack->GetP4();
00824     G4ThreeVector tp3(tp4.x(), tp4.y(), tp4.z());
00825     G4ThreeVector diffP = (p-tp3).mag2();
00826     G4int pdgT = truthTrack->GetPDGCode();
00827     if(pdgT==-22) pdgT=22;
00828     if(pdgcode == pdgT && diffP<minDiffP)
00829       minDiffP = diffP;
00830     if(pdgcode == pdgT && diffP < 1E-9 && 
00831        truthTrack->NotFound())
00832     {
00833       m_trackIndex = truthTrack->GetIndex();
00834       if(m_logLevel<=4)
00835         G4cout<<"m_trackIndex: "<<m_trackIndex<<" "<<G4endl;
00836       truthTrack->SetPDGCharge(track->GetDefinition()->GetPDGCharge());
00837       truthTrack->SetParticleName(track->GetDefinition()->GetParticleName());
00838       truthTrack->SetG4TrackId(track->GetTrackID());
00839       truthTrack->Found();
00840       G4int size = truthTrack->GetDaughterIndexes().size();
00841       if(size>0)
00842       {
00843         G4double minDau = (truthTrack->GetDaughterIndexes())[0];
00844         G4double maxDau = (truthTrack->GetDaughterIndexes())[size-1];
00845         if(m_logLevel<=4)
00846           G4cout<<"daughters: "<<minDau<<" "<<maxDau<<G4endl;
00847       }
00848       else
00849       {
00850         if(m_logLevel<=4)
00851           G4cout<<G4endl;
00852       }
00853       flag = 1; break;
00854     }
00855   }
00856   if(flag==0)
00857   {
00858     G4cout<<"MatchError! parentID=0, but match no track in trackList"<<G4endl;
00859     G4cout<<"pdgcode: "<<pdgcode<<" min p diff: "<<minDiffP<<G4endl;
00860   }
00861 }

void BesSensitiveManager::UpdateVertex (  BesTStats  ,
const G4Track * 
) 
 

void BesSensitiveManager::UpdateVertex (  BesTStats  ,
const G4Track * 
) 
  

00358 {
00359   BesTruthTrack* truthTrack = (*m_trackList)[m_trackIndex];
00360   G4int terminalIndex = truthTrack->GetTerminalVertex()->GetIndex();
00361   if(m_logLevel<=4)  
00362     G4cout<<"updateVertex:"<<terminalIndex<<" ";
00363   BesTruthVertex* vertex = (*m_vertexList)[terminalIndex];
00364   G4StepPoint *finalStep = track->GetStep()->GetPostStepPoint();
00365   const G4ThreeVector pos(finalStep->GetPosition());
00366   G4double time = finalStep->GetGlobalTime();
00367   if(m_logLevel<=4)
00368     G4cout<<"newPos:"<<pos<<" "<<"newTime:"<<time<<G4endl;
00369   vertex->SetPosition(pos);
00370   vertex->SetTime(time);
00371   vertex->SetTerminal( true );
00372   //G4int vindex = m_vertexList->size();
00373   //stat.vertices.push_back(terminalIndex);
00374 
00375 }

Member Data Documentation

std::vector<BesTStats> BesSensitiveManager::chain [protected]
 

std::vector<BesTStats> BesSensitiveManager::chain [protected]
 

std::vector<BesSensitiveDetector*> BesSensitiveManager::clients [protected]
 

std::vector<BesSensitiveDetector*> BesSensitiveManager::clients [protected]
 

std::vector<BesSensitiveDetector*>::iterator BesSensitiveManager::iter [protected]
 

std::vector<BesSensitiveDetector*>::iterator BesSensitiveManager::iter [protected]
 

std::vector<BesSensitiveDetector*>::iterator BesSensitiveManager::iter_end [protected]
 

std::vector<BesSensitiveDetector*>::iterator BesSensitiveManager::iter_end [protected]
 

G4int BesSensitiveManager::m_count [protected]
 

HepMC::GenEvent* BesSensitiveManager::m_hepmcevt [protected]
 

HepMC::GenEvent* BesSensitiveManager::m_hepmcevt [protected]
 

G4int BesSensitiveManager::m_logLevel [protected]
 

G4ThreeVector BesSensitiveManager::m_pos0 [protected]
 

BesSensitiveManager* BesSensitiveManager::m_sensitiveManager [static, protected]
 

BesSensitiveManager * BesSensitiveManager::m_sensitiveManager = 0 [static, protected]
 

G4double BesSensitiveManager::m_t0 [protected]
 

G4int BesSensitiveManager::m_trackFlag [protected]
 

G4int BesSensitiveManager::m_trackIndex [protected]
 

std::vector<BesTruthTrack*>* BesSensitiveManager::m_trackList [protected]
 

std::vector<BesTruthTrack*>* BesSensitiveManager::m_trackList [protected]
 

std::vector<BesTruthVertex*>* BesSensitiveManager::m_vertexList [protected]
 

std::vector<BesTruthVertex*>* BesSensitiveManager::m_vertexList [protected]
 

The documentation for this class was generated from the following files:
Generated on Wed Feb 2 15:53:14 2011 for BOSS6.5.5 by  doxygen 1.3.9.1