#include <BesSensitiveManager.hh>

Public Member Functions
	BesSensitiveManager ()

	~BesSensitiveManager ()

void	AddSensitiveDetector (BesSensitiveDetector *detector)

G4int	GetCurrentTrackIndex () const

std::vector< BesTruthTrack * > *	GetTrackList ()

std::vector< BesTruthVertex * > *	GetVertexList ()

void	BeginOfTruthEvent (const G4Event *)

void	EndOfTruthEvent (const G4Event *)

G4int	CheckType (const HepMC::GenEvent *hepmcevt)

void	SaveParticlesFromGenerator ()

void	SetVertex0 (const G4Event *)

void	UpdatePrimaryTrack (const G4Track *)

G4bool	CheckDecayTrack (const G4Track *)

void	UpdateVertex (BesTStats, const G4Track *)

G4bool	MatchDaughterTrack (const G4Track *)

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

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

void	BeginOfTrack (const G4Track *track)

void	EndOfTrack (const G4Track track, G4TrackingManager )

G4int	GetNumberTracks () const

G4int	GetNumberVertices () const

void	ClearEvent ()

void	SetLogLevel (G4int level)

Static Public Member Functions
static BesSensitiveManager *	GetSensitiveManager ()

Protected Member Functions
void	MakeNewTrack (BesTStats &stat, const G4Track *track)

BesTStats	FollowTrack (const G4Track *track)

Protected Attributes
std::vector< BesSensitiveDetector * >::iterator	iter

std::vector< BesSensitiveDetector * >::iterator	iter_end

G4int	m_trackIndex

G4int	m_trackFlag

std::vector< BesTruthTrack * > *	m_trackList

std::vector< BesTruthVertex * > *	m_vertexList

std::vector< BesSensitiveDetector * >	clients

std::vector< BesTStats >	chain

G4int	m_count

G4ThreeVector	m_pos0

G4double	m_t0

G4int	m_logLevel

HepMC::GenEvent *	m_hepmcevt

Static Protected Attributes
static BesSensitiveManager *	m_sensitiveManager =0

Detailed Description

Definition at line 73 of file BesSensitiveManager.hh.

Constructor & Destructor Documentation

◆ BesSensitiveManager()

BesSensitiveManager::BesSensitiveManager ( )

inline

Definition at line 75 of file BesSensitiveManager.hh.

                       :m_trackIndex(-1),m_trackFlag(0),m_count(0) 
  {
    m_hepmcevt = 0; 
    if(m_sensitiveManager) 
      G4Exception("BesSensitiveManager constructed twice.");
    m_sensitiveManager = this ;
  }

◆ ~BesSensitiveManager()

BesSensitiveManager::~BesSensitiveManager ( )

inline

Definition at line 83 of file BesSensitiveManager.hh.

83{;}

Member Function Documentation

◆ AddSensitiveDetector()

void BesSensitiveManager::AddSensitiveDetector ( BesSensitiveDetector * detector )

inline

Definition at line 88 of file BesSensitiveManager.hh.

  {
    clients.push_back( detector );
  }

Referenced by BesSensitiveDetector::BesSensitiveDetector().

◆ BeginOfTrack()

void BesSensitiveManager::BeginOfTrack ( const G4Track * track )

Definition at line 252 of file BesSensitiveManager.cc.

{
  m_trackFlag = 0; //if m_trackFlag=1, this track is from generator
 
  // Append our track in the parentage history
  chain.push_back(FollowTrack(track));
 
  // Get TruthTrack index
  //int i = chain.size();
 
  /*for(;;) 
  {
    if(--i<0) G4cerr<<"::parents are corrupted"<<G4endl;
    m_trackIndex = chain[i].trackIndex;
    if (m_trackIndex != -1) break;
  }*/
  
  // Invoke clients
  //
  iter=clients.begin();
  iter_end=clients.end();
  while(iter!=iter_end) 
  {
    (*iter)->BeginOfTrack( track );
    iter++;
  }
}

Referenced by BesTrackingAction::PreUserTrackingAction().

◆ BeginOfTruthEvent()

void BesSensitiveManager::BeginOfTruthEvent ( const G4Event * evt )

Definition at line 43 of file BesSensitiveManager.cc.

{
  // Allocate our lists. We'll own these until EndOfEvent
 
  m_trackList  = new std::vector<BesTruthTrack*>;
  m_vertexList = new std::vector<BesTruthVertex*>;
  
  SetVertex0(evt);
 
  m_count = 0;
  SaveParticlesFromGenerator();
  m_count += m_trackList->size();
 
  // Tell clients
  iter=clients.begin();
  iter_end=clients.end();
  while( iter != iter_end ) 
  {
    (*iter)->BeginOfTruthEvent(evt);
    iter++;
  }
}

Referenced by BesEventAction::BeginOfEventAction().

◆ CheckDecayTrack()

G4bool BesSensitiveManager::CheckDecayTrack ( const G4Track * track )

Definition at line 590 of file BesSensitiveManager.cc.

{ 
  G4int flag = 0;
  G4int trackID = track->GetTrackID();
  G4int parentID = track->GetParentID();
  G4int pdgcode = track->GetDefinition()->GetPDGEncoding();
  G4ThreeVector p3 = track->GetMomentum();
  if(m_logLevel<=4)
    G4cout<<"CheckDecayTrack p3: "<<p3<<"  "<<track->GetTotalEnergy()<<G4endl;
 
  //HepLorentzVector  p4(track->GetMomentum(), track->GetTotalEnergy());
  
  BesTruthTrack* truthTrack=0;
  G4int nTrack = m_trackList->size();
  G4int parentTrackIndex=-99;
  G4int terminalVertexIndex=-99;
  BesTruthTrack* pTrack;
  for(int i=0;i<nTrack;i++)
  {
    truthTrack = (*m_trackList)[i];
    if(truthTrack->GetG4TrackId() == parentID)
    {
      //parentTrackIndex = i;
      parentTrackIndex = truthTrack->GetIndex();
      if(truthTrack->GetTerminalVertex())
      {
        terminalVertexIndex = truthTrack->GetTerminalVertex()->GetIndex();
        pTrack = truthTrack;
        if(m_logLevel<=4)
        {
          G4cout<<"parentTrackIndex:"<<parentTrackIndex<<" "
                <<"parent terminate at vertex: "<<terminalVertexIndex<<G4endl;
          if(parentTrackIndex != i)
            G4cout<<"i: "<<i<<std::endl;
        }
        break;
      }
    }
  }
  if(parentTrackIndex==-99 || terminalVertexIndex==-99)
  {
    G4cout<<"MatchDecayError!"<<G4endl;
    //return false;
    exit(1);
  }
 
  //match decayed track with p3
  /*G4bool matchDauFlag=false; 
  if(terminalVertexIndex>0)
    matchDauFlag = MatchDaughterTrack(track);
  if(matchDauFlag)
    return true;*/
 
  //match decayed track with vertex index
  G4double minDiffP=1e99;
  G4int truthI=-9;
  /*for(int i=0;i<nTrack;i++)
  {
    truthTrack = (*m_trackList)[i];
    G4String source = truthTrack->GetSource();
    G4int pdgcode_tru = truthTrack->GetPDGCode();
    if(pdgcode_tru==-22)  pdgcode_tru = 22;
    if(truthTrack->GetVertex()->GetIndex() == terminalVertexIndex &&
       pdgcode_tru == pdgcode &&source=="FromGenerator"&&truthTrack->NotFound())
    {
      HepLorentzVector tp4 = truthTrack->GetP4();
      G4ThreeVector tp3(tp4.x(), tp4.y(), tp4.z());
      G4double diffP = (p3-tp3).mag2();
      if(m_logLevel<=4)
        G4cout<<"diffP: "<<diffP<<G4endl;
      if(diffP<minDiffP )
      {  minDiffP = diffP;  truthI = i; }
    }
  }*/
  
  /*if(truthI!=-9)
  {
    if(m_logLevel<=4)
      G4cout<<"MatchDecayedTrackWithVertexIndex, trackIndex:"<<truthI<<G4endl;
    m_trackIndex = truthI;
    truthTrack=(*m_trackList)[truthI];
    //truthTrack->SetP4(p4);
    truthTrack->SetPDGCharge(track->GetDefinition()->GetPDGCharge());
    truthTrack->SetParticleName(track->GetDefinition()->GetParticleName());
    truthTrack->SetG4TrackId(track->GetTrackID());
    truthTrack->Found();
    G4int size = truthTrack->GetDaughterIndexes().size();
    if(size>0)
    {
      G4int minDau = (truthTrack->GetDaughterIndexes())[0];
      G4int maxDau = (truthTrack->GetDaughterIndexes())[size-1];
      if(m_logLevel<=4)
        G4cout<<"daughters: "<<minDau<<" "<<maxDau<<G4endl;
    }
    else
    {
      if(m_logLevel<=4)
      G4cout<<G4endl;
    }
    return true; 
  }*/ //match decay track with vertex index
 
  //match decay track
  
  if(m_logLevel<=4)
    std::cout<<"chain.back().vertices.size(): "<<chain.back().vertices.size()<<std::endl;
  //if particle is from generator, chain.back().vertices.size()=0;
  //if particle is not from generator, but decayed during flight in Geant4, chain.back().vertices.size()=1;
  if(chain.back().vertices.size()<1)
  {
    if(m_logLevel<=4)
      std::cout<<"trackList size: "<<m_trackList->size()<<std::endl;
    std::vector<int>* vDau = new std::vector<int>;
    GetDaughterVertexes(pTrack, vDau);
    G4int sizev = vDau->size();
    if(sizev>0)
    {
      for(int i=0;i<nTrack;i++)
      {
        truthTrack = (*m_trackList)[i];
        G4int vIndex = truthTrack->GetVertex()->GetIndex();
        G4int pdgT = truthTrack->GetPDGCode();
        if(pdgT==-22)  pdgT = 22;
        if(pdgT == pdgcode)
        {
          G4bool matchFlag = MatchVertex(vIndex, vDau);
          if(matchFlag)
          {
            HepLorentzVector tp4 = truthTrack->GetP4();
            G4ThreeVector tp3(tp4.x(), tp4.y(), tp4.z());
            G4double diffP = (p3-tp3).mag2();
            if(m_logLevel<=4)
            {
              G4cout<<"index: "<<truthTrack->GetIndex()<<"tp3: "<<tp3<<G4endl;
              G4cout<<"diffP: "<<diffP<<G4endl;
            }
            if(diffP<minDiffP)
            {  
              minDiffP = diffP;  truthI = i; 
              if(m_logLevel<=4)
                G4cout<<"truthI: "<<i<<G4endl;
            }
          }
        }
      }
      if(vDau)
        delete vDau;
    }
  }
 
  if(truthI!=-9)
  {
    m_trackIndex = truthI;
    if(m_logLevel<=4)
    {
      G4cout<<"MatchDecayedTrackWithTrueMother"<<G4endl;
      G4cout<<"MatchWithTrueMother m_trackIndex: "<<m_trackIndex<<G4endl;
      if(minDiffP>1e-9)
        G4cout<<"large minDiffP: "<<minDiffP<<G4endl;
    }
    truthTrack=(*m_trackList)[truthI];
    G4double mass = truthTrack->GetP4().m();
    G4double E = sqrt(mass*mass+p3.x()*p3.x()+p3.y()*p3.y()+p3.z()*p3.z());    
    truthTrack->GetP4().setX(p3.x());
    truthTrack->GetP4().setY(p3.y());
    truthTrack->GetP4().setZ(p3.z());
    truthTrack->GetP4().setE(E);
    HepLorentzVector  p4 = truthTrack->GetP4();
    if(m_logLevel<=4)
    {
      G4cout<<"primary p: "<<p4.x()<<" "<<p4.y()<<" "<<p4.z()<<" "<<p4.m()<<G4endl;
    }
    truthTrack->SetPDGCharge(track->GetDefinition()->GetPDGCharge());
    truthTrack->SetParticleName(track->GetDefinition()->GetParticleName());
    truthTrack->SetG4TrackId(track->GetTrackID());
    truthTrack->Found();
    //truthTrack->GetVertex()->SetPosition(track->GetPosition());
    //truthTrack->GetVertex()->SetTime(track->GetGlobalTime());
    return true;
  }
  return false;   
}

Referenced by FollowTrack().

◆ CheckType()

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

Definition at line 79 of file BesSensitiveManager.cc.

{
  G4int flag =0;
  for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
          vitr != hepmcevt->vertices_end(); ++vitr ) { // loop for vertex ...
    for (HepMC::GenVertex::particle_iterator
             pitr= (*vitr)->particles_begin(HepMC::children);
                pitr != (*vitr)->particles_end(HepMC::children); ++pitr) {
      if((*pitr)->end_vertex())
        {flag = 1; break;}
    }
    if(flag) break;
  }
  if(m_logLevel <= 4)
  {
    if(flag == 0)
      G4cout<<G4endl<<"generator is GENBES!"<<G4endl;
    else
      G4cout<<G4endl<<"generator is not GENBES!"<<G4endl;
  }
  return flag;
}

Referenced by SaveParticlesFromGenerator().

◆ ClearEvent()

void BesSensitiveManager::ClearEvent ( )

Definition at line 234 of file BesSensitiveManager.cc.

{
  if(m_trackList)
  {
    for(size_t i=0;i<m_trackList->size();i++)
      delete (*m_trackList)[i];
    m_trackList->clear();
    delete m_trackList;
  }
  if(m_vertexList)
  {
    for(size_t i=0;i<m_vertexList->size();i++)
      delete (*m_vertexList)[i];
    m_vertexList->clear();
    delete m_vertexList;
  }
}

Referenced by BesEventAction::EndOfEventAction().

◆ EndOfTrack()

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

Definition at line 279 of file BesSensitiveManager.cc.

{
  if(chain.back().savedByDefault==true)
  { 
    BesTStats &stat = chain.back();
    G4int endVtxFlag = 0;
    if(m_trackFlag==1)
    {
      BesTruthTrack* truthTrack = (*m_trackList)[m_trackIndex];
      if(truthTrack->GetTerminalVertex())
      {
        UpdateVertex(stat,track); endVtxFlag = 1;}
    } 
    if(endVtxFlag == 0)
    {
      // Tracks saved to BesTruthTrack require additional work
      G4StepPoint *finalStep = track->GetStep()->GetPostStepPoint();
      G4StepStatus stepStatus = finalStep->GetStepStatus();
      if (track->GetNextVolume() != 0 ||
          (stepStatus != fGeomBoundary && stepStatus != fWorldBoundary) )
      {
        if(m_logLevel<=4)
          G4cout<<"Particle died. make a terminal vertex: ";
        // Particle died. We always make a terminal vertex
        const G4ThreeVector pos(finalStep->GetPosition());
        G4int vindex = m_vertexList->size();
        if(m_logLevel<=4)
          G4cout<<vindex<<G4endl;
        stat.vertices.push_back(vindex);
        BesTruthVertex *newVertex = new BesTruthVertex();
        newVertex->SetPosition(pos);
        newVertex->SetTime( finalStep->GetGlobalTime());
        if(m_trackFlag==1)
          newVertex->SetParentTrack( (*m_trackList)[m_trackIndex] );
        else
          newVertex->SetParentTrack( m_trackList->back() );
        newVertex->SetTerminal( true );
        newVertex->SetIndex( vindex );
        
        //set minDaughter index equal to m_count
        newVertex->SetMinDau(m_count);  
        
        //Get how many decayed daughters this track has, if no, 0 is set.
        G4TrackVector* trackVector = trackingManager->GimmeSecondaries();
        G4int nSecon = trackVector->size();
        G4Track* seconTrack;
        G4int nDau=0;
        if(nSecon>0)
        {
          for(G4int i=0;i<nSecon;i++)
          {
            seconTrack = (*trackVector)[i];
            G4String processName = seconTrack->GetCreatorProcess()->GetProcessName();
            if(processName == "Decay")
              nDau++;
          }
        }
        m_vertexList->push_back( newVertex );
  
        m_count += nDau;
        
        // Give this terminal vertex to the track
        //(*m_trackList)[stat.trackIndex]->SetTerminalVertex( m_vertexList->back() );
        // for particles from generator, 
        // m_trackList->back() may not be current track.
 
        if(m_trackFlag==1)
          (*m_trackList)[m_trackIndex]->SetTerminalVertex( m_vertexList->back() );
        else
          m_trackList->back()->SetTerminalVertex( m_vertexList->back() );
      }
    }
  } 
  // Invoke clients
  iter=clients.begin();
  iter_end=clients.end();
  while( iter!=iter_end) {
    (*iter)->EndOfTrack( track );
    iter++;
  }
}

Referenced by BesTrackingAction::PostUserTrackingAction().

◆ EndOfTruthEvent()

void BesSensitiveManager::EndOfTruthEvent ( const G4Event * evt )

Definition at line 222 of file BesSensitiveManager.cc.

{
  // Tell clients
  iter=clients.begin();
  iter_end=clients.end();
  while( iter != iter_end ) 
  {
    (*iter)->EndOfTruthEvent(evt);
    iter++;
  }
}

Referenced by BesEventAction::EndOfEventAction().

◆ FollowTrack()

BesTStats BesSensitiveManager::FollowTrack ( const G4Track * track )

protected

Definition at line 505 of file BesSensitiveManager.cc.

{
  // Some default stats for this track
  BesTStats stat( track->GetTrackID(),
          HepLorentzVector( track->GetMomentum(), track->GetTotalEnergy() ),
          track->GetGlobalTime());
 
  // Check immediate parent
  G4int parentId = track->GetParentID();
  G4int pdg = track->GetDefinition()->GetPDGEncoding();
  G4String particleName = track->GetDefinition()->GetParticleName();
  G4ThreeVector trackPos = track->GetPosition();
  G4double diffT = m_t0-track->GetGlobalTime();
  G4ThreeVector diffPos = (m_pos0-trackPos).mag2();
 
  if (parentId == 0) 
  {
    if(m_logLevel<=4)
    {
      G4cout<<G4endl;
      G4cout<<"trackId:"<<track->GetTrackID()<<" ";
      G4cout<<"parentId:"<<parentId<<" ";
      G4cout<<"pdg:"<<pdg<<" ";
      G4cout<<"name:"<<particleName<<" ";
      G4cout<<"timeDiff:"<<diffT<<" ";
      G4cout<<"posDiff:"<<diffPos<<G4endl;
    }
 
    // Primary particle: wipe decay chain
    chain.clear();
    // Always save
    stat.savedByDefault = true;
    // Make TruthTrack
    if(m_hepmcevt==0)
      MakeNewTrack( stat, track );  
    else
    {
      UpdatePrimaryTrack(track);
      m_trackFlag = 1;
    }
    return stat;
  }
    
  // Not primary particle: trim down chain until
  // we uncover the parent. The parent must be there!
  for(;;) 
  {
    if (chain.back().G4index == parentId) break;
    chain.pop_back();
  }
 
  // This track is a candidate for saving by default
  // only if its parent was saved by default
  if (chain.back().savedByDefault) 
  {
    // There are three ways particles are saved by default:
    //   1. if they are a primary particle
    //   2. if they are the result of a decay
    //
    G4String processName = track->GetCreatorProcess()->GetProcessName();
    if (processName=="Decay") 
    {
      if(m_logLevel<=4)
      {
        G4cout<<G4endl;
        G4cout<<"trackId: "<<track->GetTrackID()<<" ";
        G4cout<<"parentId: "<<parentId<<" ";
        G4cout<<"pdg: "<<pdg<<" ";
        G4cout<<"pname: "<<particleName<<G4endl;
      }
      stat.savedByDefault = true;
 
      if(CheckDecayTrack(track))
        m_trackFlag = 1;
      else
        MakeNewTrack( stat, track);
      return stat;
    }
  }
  
  //now this track will not be saved as BesTruthTrack
  stat.savedByDefault = false;
  return stat;
}

Referenced by BeginOfTrack().

◆ GetCurrentTrackIndex()

G4int BesSensitiveManager::GetCurrentTrackIndex ( ) const

inline

Definition at line 93 of file BesSensitiveManager.hh.

93{ return m_trackIndex; }

Referenced by BesSensitiveDetector::GetCurrentTrackIndex().

◆ GetDaughterVertexes()

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

Definition at line 773 of file BesSensitiveManager.cc.

{
  G4int size = pTrack->GetDaughterIndexes().size();
  if(size>0)
  {
    G4int minDau = (pTrack->GetDaughterIndexes())[0];
    G4int maxDau = (pTrack->GetDaughterIndexes())[size-1];
    //note! here, dau<=maxDau, not dau<maxDau
    for(G4int dau=minDau;dau<=maxDau;dau++)
    {
      if(dau<m_trackList->size())
      {
        BesTruthTrack* truthTrack = (*m_trackList)[dau];
        if(truthTrack->GetVertex())
        {
          vDau->push_back(truthTrack->GetVertex()->GetIndex());
          GetDaughterVertexes(truthTrack,vDau);
        }
      } 
    }
  }
}

Referenced by CheckDecayTrack(), and GetDaughterVertexes().

◆ GetNumberTracks()

G4int BesSensitiveManager::GetNumberTracks ( ) const

inline

Definition at line 119 of file BesSensitiveManager.hh.

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

◆ GetNumberVertices()

G4int BesSensitiveManager::GetNumberVertices ( ) const

inline

Definition at line 120 of file BesSensitiveManager.hh.

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

◆ GetSensitiveManager()

static BesSensitiveManager * BesSensitiveManager::GetSensitiveManager ( )

inlinestatic

Definition at line 85 of file BesSensitiveManager.hh.

85{return m_sensitiveManager; }

Referenced by BesEventAction::BeginOfEventAction(), BesSensitiveDetector::BesSensitiveDetector(), BesEventAction::EndOfEventAction(), BesSensitiveDetector::GetCurrentTrackIndex(), BesTrackingAction::PostUserTrackingAction(), BesTrackingAction::PreUserTrackingAction(), BesEmcSD::ProcessHits(), BesMcTruthWriter::SaveMcParticle(), BesAsciiIO::SaveTrackTruth(), and BesAsciiIO::SaveVertexTruth().

◆ GetTrackList()

std::vector< BesTruthTrack * > * BesSensitiveManager::GetTrackList ( )

inline

Definition at line 94 of file BesSensitiveManager.hh.

94{return m_trackList;}

Referenced by BesSensitiveDetector::GetCurrentTrackIndex(), BesEmcSD::ProcessHits(), BesMcTruthWriter::SaveMcParticle(), and BesAsciiIO::SaveTrackTruth().

◆ GetVertexList()

std::vector< BesTruthVertex * > * BesSensitiveManager::GetVertexList ( )

inline

Definition at line 95 of file BesSensitiveManager.hh.

95{return m_vertexList;}

Referenced by BesMcTruthWriter::SaveMcParticle(), and BesAsciiIO::SaveVertexTruth().

◆ MakeNewTrack()

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

protected

Definition at line 381 of file BesSensitiveManager.cc.

{
  if (stat.originVertex < 0 && chain.size() > 0) 
  {
    if(m_logLevel<=4)
      G4cout<<"MakeNewTrack for decayed particles,";
    //for decayed particles, there may already be a BesTruthVertex that is suitable for this track.
    //If so, it's always the terminal vertex of its immediate parent( chain.back() ).
 
    G4int parentVstat = chain.back().vertices.size();
    while( --parentVstat >= 0)
    {
      G4int vindex = chain.back().vertices[parentVstat];
 
      G4ThreeVector diff((*m_vertexList)[vindex]->GetPosition());
      diff = diff - track->GetPosition();
      if (diff.mag2() < 1E-9) 
      {
          stat.originVertex = vindex;
        if(m_logLevel<=4)
          G4cout<<"find a vertex:"<<vindex;
        (*m_vertexList)[vindex]->AddCurrentDau();
        break;
      }
    }
    
  }
  
  if(stat.originVertex < 0 && chain.size() == 0)
  {
    //for primary track, check if there is already a vertex suitable for it.
    G4int nVertex = m_vertexList->size();     
    for(G4int i=0;i<nVertex;i++)
    {
      G4ThreeVector diff((*m_vertexList)[i]->GetPosition());
      diff = diff -track->GetPosition();
      if(diff.mag2() < 1E-9)
      {
        stat.originVertex = i; 
        if(m_logLevel<=4)
          G4cout<<"MakeNewTrack for primary particles,find a vertex:"
                <<i;
        break;
      }
    }
  }
  
  if (stat.originVertex < 0)
  {
    if(m_logLevel<=4) 
      G4cout<<"MakeNewTrack for primary particles,make new vertex:"
            <<m_vertexList->size();
    // it's a primary track(there's no vertex suitable for it)
    // Make a new BesTruthVertex
    const G4ThreeVector v(track->GetPosition());
        
    stat.originVertex = m_vertexList->size();
 
    BesTruthVertex *newVertex = new BesTruthVertex();
    newVertex->SetPosition(v);
    newVertex->SetTime( track->GetGlobalTime());
    
    //if (chain.size() > 0) {
    //  G4int trackIndex = -1;
    //  G4int i = chain.size();
    //  while(--i>=0 && trackIndex < 0) trackIndex = chain[i].trackIndex;
    //  newVertex->SetParentTrack( trackIndex < 0 ? 0 : (*m_trackList)[trackIndex] );
    //}
    
    newVertex->SetParentTrack( 0 );
    newVertex->SetTerminal( false );
    newVertex->SetIndex( stat.originVertex );
    if(track->GetCreatorProcess())
      newVertex->SetProcessName(track->GetCreatorProcess()->GetProcessName());
    else
      newVertex->SetProcessName("generator");    
    
    m_vertexList->push_back( newVertex );
  }
    
  // Now, finally make our new BesTruthTrack
  // We'll leave the assignment of terminalVertex until
  // we know what happens to this track     
  BesTruthTrack *newTrack = new BesTruthTrack();
  newTrack->SetP4( stat.p4 );
  newTrack->SetPDGCode(track->GetDefinition()->GetPDGEncoding());
  newTrack->SetPDGCharge(track->GetDefinition()->GetPDGCharge());
  newTrack->SetParticleName(track->GetDefinition()->GetParticleName());
  newTrack->SetSource("FromG4");
  BesTruthVertex* vertex = (*m_vertexList)[stat.originVertex];
  newTrack->SetVertex(vertex);
  
  if(track->GetParentID()==0)
  {
    stat.trackIndex = m_count;
    m_count++;
  }
  else
    stat.trackIndex = vertex->GetMinDau() + vertex->GetCurrentDau()-1;
  
  newTrack->SetIndex( stat.trackIndex );
  m_trackIndex = stat.trackIndex;
  if(m_logLevel<=4)
    G4cout<<"  m_trackIndex:"<<m_trackIndex<<G4endl;
  
  newTrack->SetG4TrackId( track->GetTrackID()); 
  m_trackList->push_back( newTrack );
  //tell the parent track its daughter track index
  BesTruthTrack* parent = newTrack->GetVertex()->GetParentTrack();
  if(parent)
  {
    parent->AddDaughterIndex(stat.trackIndex);
    if(m_logLevel<=4)
      G4cout<<"add this daughter to:"<<parent->GetIndex()<<G4endl;
  }
}

Referenced by FollowTrack().

◆ MatchDaughterTrack()

G4bool BesSensitiveManager::MatchDaughterTrack ( const G4Track * track )

Definition at line 810 of file BesSensitiveManager.cc.

{
  /*G4int flag=0;
  G4int pdgcode = track->GetDefinition()->GetPDGEncoding();
  G4ThreeVector p = track->GetMomentum();
  HepLorentzVector  p4(track->GetMomentum(), track->GetTotalEnergy());
  G4ThreeVector pos = track->GetPosition();
  G4double time = track->GetGlobalTime();
  BesTruthTrack* truthTrack=0;
  G4int size = m_trackList->size();
  if(size>0)
  { 
    for(G4int i=0;i<size;i++)
    {
      truthTrack=(*m_trackList)[i];
      G4String source = truthTrack->GetSource(); 
      G4int pdgcode_tru = truthTrack->GetPDGCode();
      if(pdgcode_tru==-22)  pdgcode_tru = 22;
      HepLorentzVector tp4 = truthTrack->GetP4();
      G4ThreeVector tp3(tp4.x(), tp4.y(), tp4.z());
      G4double diffP = (p-tp3).mag2();
      if(source=="FromGenerator"&&pdgcode_tru==pdgcode&&diffP<1e-9&&truthTrack->NotFound())
      {
        m_trackIndex = i;
        //truthTrack->SetP4(p4);
        truthTrack->SetPDGCharge(track->GetDefinition()->GetPDGCharge());
        truthTrack->SetParticleName(track->GetDefinition()->GetParticleName());
        truthTrack->SetG4TrackId(track->GetTrackID());
        truthTrack->Found();
        //truthTrack->GetVertex()->SetPosition(pos);
        //truthTrack->GetVertex()->SetTime(time);
        flag=1;
        if(m_logLevel<=4)
        {
          G4cout<<"MatchDecayedTrackWithP3!"<<"trackIndex:"<<m_trackIndex
                <<" pdgcode:"<<pdgcode<<G4endl;
          G4cout<<"parentIndex:"<<truthTrack->GetVertex()->GetParentTrack()->GetIndex()
                <<"  PDGCode:"<<truthTrack->GetVertex()->GetParentTrack()->GetPDGCode()<<G4endl; 
        }
        return true;
      }
    }
  }
  if(flag!=1)
    return false;*/
  return true;
}

◆ MatchVertex()

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

Definition at line 796 of file BesSensitiveManager.cc.

{
  G4int size = vDau->size();
  if(size>0)
  {
    for(G4int i=0;i<size;i++)
    {
      if(vIndex==(*vDau)[i])
        return true;
    }
  }
  return false;
}

Referenced by CheckDecayTrack().

◆ SaveParticlesFromGenerator()

void BesSensitiveManager::SaveParticlesFromGenerator ( )

Definition at line 102 of file BesSensitiveManager.cc.

{
  IDataProviderSvc* p_evtSvc=0;
  ISvcLocator* svcLocator = Gaudi::svcLocator();
  StatusCode sc=svcLocator->service("EventDataSvc", p_evtSvc);
  if (sc.isFailure())
      std::cout<<"BesHepMCInterface could not access EventDataSvc!!"<<std::endl;
 
  SmartDataPtr<McGenEventCol> mcCollptr( p_evtSvc, "/Event/Gen");
 
  if ( mcCollptr == 0 )
    std::cout << "no McGenEventCollection found."  << std::endl;
 
  else 
  {
    McGenEventCol::const_iterator it = mcCollptr->begin();
    McGenEvent* mcEvent = (McGenEvent* ) (*it);
    m_hepmcevt = mcEvent->getGenEvt();
    G4int typeFlag = CheckType(m_hepmcevt); 
    
 
    for(HepMC::GenEvent::vertex_const_iterator vitr= m_hepmcevt->vertices_begin();
              vitr != m_hepmcevt->vertices_end(); ++vitr ) 
    {
      G4int barcodeVtx = (*vitr)->barcode();
      if(m_logLevel<=4)
        G4cout<<G4endl<<"barcodeVtx:"<<barcodeVtx<<" ";
  
      G4LorentzVector xvtx((*vitr)->position().x(),(*vitr)->position().y(),(*vitr)->position().z(),(*vitr)->position().t());
      G4ThreeVector v(xvtx.x(), xvtx.y(), xvtx.z());
      BesTruthVertex* newVertex = new BesTruthVertex();
      newVertex->SetPosition(v);
      newVertex->SetTime(xvtx.t()*mm/c_light);
      if(m_logLevel<=4)
        G4cout<<"xyzt:"<<xvtx.x()<<" "<<xvtx.y()<<" "
            <<xvtx.z()<<" "<<xvtx.t()*mm/c_light<<" ";
      G4int nTrack = m_trackList->size();
      BesTruthTrack* parentTrack=0;
      G4int parentTrackIndex= -99;
      for(int i=0;i<nTrack;i++)
      {
        parentTrack = (*m_trackList)[i];
        if(parentTrack->GetBarcodeEndVtx() == barcodeVtx)
        {
          newVertex->SetParentTrack(parentTrack);
          parentTrackIndex = parentTrack->GetIndex();
          if(m_logLevel<=4)
            G4cout<<"parentTrack:"<<parentTrackIndex<<" ";
          parentTrack->SetTerminalVertex(newVertex);
          //break;
        }
      }
  
      G4int vtxFlag=0;
      G4int pOut = (*vitr)->particles_out_size();
      HepMC::GenVertex::particle_iterator pitr = (*vitr)->particles_begin(HepMC::children);
      G4int pdgcode= (*pitr)-> pdg_id();
      if(pOut == 1 && pdgcode == -11 && typeFlag==1)
        vtxFlag=1;
          
      if(vtxFlag==0)
      {
        m_vertexList->push_back(newVertex);
        newVertex->SetIndex(m_vertexList->size()-1);
        if(m_logLevel<=4)
          G4cout<<"vindex:"<<m_vertexList->size()-1<<G4endl;
  
        for (HepMC::GenVertex::particle_iterator pitr= (*vitr)->particles_begin(HepMC::children); pitr != (*vitr)->particles_end(HepMC::children); ++pitr)
        {
          G4LorentzVector p((*pitr)->momentum().x(), (*pitr)->momentum().y(),(*pitr)->momentum().z(),(*pitr)->momentum().t());
          G4LorentzVector pGeV(p.x()*GeV,p.y()*GeV,p.z()*GeV,p.t()*GeV);
          G4int pdgcode = (*pitr)->pdg_id();
  
          BesTruthTrack* newTrack = new BesTruthTrack;
          newTrack->SetP4(pGeV);
          newTrack->SetPDGCode(pdgcode);
          if(m_logLevel<=4)
          {  
            G4cout<<"pdg:"<<pdgcode<<" ";
            G4cout<<"p:"<<p.x()<<" "<<p.y()<<" "<<p.z()<<" "<<p.t()<<" ";
            G4double mass = sqrt(p.t()*p.t()-p.x()*p.x()-p.y()*p.y()-p.z()*p.z());
            G4cout<<mass<<" ";
          }
          newTrack->SetPDGCharge(99);
          newTrack->SetParticleName("unknown");
          newTrack->SetVertex(newVertex);
          newTrack->SetTerminalVertex(0);
          newTrack->SetSource("FromGenerator");
        
          G4int barcodeEndVtx=0;
          if((*pitr)->end_vertex())
          {
            barcodeEndVtx = (*pitr)->end_vertex()->barcode();
            if(m_logLevel<=4)
              G4cout<<"endVtx:"<<barcodeEndVtx<<" ";
          }
          newTrack->SetBarcodeEndVtx(barcodeEndVtx);
  
          m_trackList->push_back(newTrack);
          newTrack->SetIndex(m_trackList->size()-1);
          if(m_logLevel<=4)
            G4cout<<"index:"<<m_trackList->size()-1<<" ";
          //here, parentTrack->GetIndex can't be used,
          //use parentTrackIndex instead
          if(parentTrackIndex>=0)
          {
            if(m_logLevel<=4)
              G4cout<<"mother:"<<parentTrackIndex; 
            (*m_trackList)[parentTrackIndex]->AddDaughterIndex(m_trackList->size()-1);
          }
          if(m_logLevel<=4)
            G4cout<<G4endl;
        }  
      }
    } 
  } 
}

Referenced by BeginOfTruthEvent().

◆ SetLogLevel()

void BesSensitiveManager::SetLogLevel ( G4int level )

inline

Definition at line 124 of file BesSensitiveManager.hh.

124{m_logLevel = level;}

Referenced by BesSim::initialize().

◆ SetVertex0()

void BesSensitiveManager::SetVertex0 ( const G4Event * anEvent )

Definition at line 66 of file BesSensitiveManager.cc.

{
  //G4int n_vertex = anEvent->GetNumberOfPrimaryVertex();
  for( G4int i=0; i<1; i++ )
  {
    G4PrimaryVertex* primaryVertex = anEvent->GetPrimaryVertex(i);
    m_pos0 = primaryVertex->GetPosition();
    m_t0 = primaryVertex->GetT0();
  }
  if(m_logLevel<=4)
    G4cout<<"pos0:"<<m_pos0<<"   t0:"<<m_t0<<G4endl;
}

Referenced by BeginOfTruthEvent().

◆ UpdatePrimaryTrack()

void BesSensitiveManager::UpdatePrimaryTrack ( const G4Track * track )

Definition at line 858 of file BesSensitiveManager.cc.

{
  if(m_logLevel<=4)
  {
    G4cout<<"UpdatePrimaryTrack! ";
    G4cout<<"position: "<<track->GetPosition()<<G4endl;
    G4cout<<"time: "<<track->GetGlobalTime()<<G4endl;
    G4cout<<"momentum: "<<track->GetMomentum()<<"  "<<track->GetTotalEnergy()<<G4endl;
  }
  G4int pdgcode = track->GetDefinition()->GetPDGEncoding();
  G4ThreeVector p = track->GetMomentum();
  //HepLorentzVector  p4(track->GetMomentum(), track->GetTotalEnergy());
  G4int nTrack = m_trackList->size();
  BesTruthTrack* truthTrack;
  G4int flag = 0;
  G4double minDiffP=1e99;
  for(int i=0;i<nTrack;i++)
  {
    truthTrack = (*m_trackList)[i];
    HepLorentzVector tp4 = truthTrack->GetP4();
    G4ThreeVector tp3(tp4.x(), tp4.y(), tp4.z());
    G4double diffP=(p-tp3).mag2();
    G4int pdgT = truthTrack->GetPDGCode();
    if(pdgT==-22) pdgT=22;
    if(pdgcode == pdgT && diffP<minDiffP)
      minDiffP = diffP;
    if(pdgcode == pdgT && diffP < 1E-9 && 
       truthTrack->NotFound())
    {
      m_trackIndex = truthTrack->GetIndex();
      if(m_logLevel<=4)
        G4cout<<"m_trackIndex: "<<m_trackIndex<<" "<<G4endl;
      truthTrack->SetPDGCharge(track->GetDefinition()->GetPDGCharge());
      truthTrack->SetParticleName(track->GetDefinition()->GetParticleName());
      truthTrack->SetG4TrackId(track->GetTrackID());
      truthTrack->Found();
      G4double mass = truthTrack->GetP4().m();
      G4double E = sqrt(mass*mass+p.x()*p.x()+p.y()*p.y()+p.z()*p.z());
      truthTrack->GetP4().setX(p.x());
      truthTrack->GetP4().setY(p.y());
      truthTrack->GetP4().setZ(p.z());
      truthTrack->GetP4().setT(E);
      if(m_logLevel<=4)
      {
        HepLorentzVector  p4 = truthTrack->GetP4();
        G4cout<<"primary p: "<<p4.x()<<" "<<p4.y()<<" "<<p4.z()<<" "<<p4.m()<<G4endl;
      }
      G4int size = truthTrack->GetDaughterIndexes().size();
      if(size>0)
      {
        G4double minDau = (truthTrack->GetDaughterIndexes())[0];
        G4double maxDau = (truthTrack->GetDaughterIndexes())[size-1];
        if(m_logLevel<=4)
          G4cout<<"daughters: "<<minDau<<" "<<maxDau<<G4endl;
      }
      else
      {
        if(m_logLevel<=4)
          G4cout<<G4endl;
      }
      flag = 1; break;
    }
  }
  if(flag==0)
  {
    G4cout<<"MatchError! parentID=0, but match no track in trackList"<<G4endl;
    G4cout<<"pdgcode: "<<pdgcode<<" min p diff: "<<minDiffP<<G4endl;
    exit(1);
  }
}

Referenced by FollowTrack().

◆ UpdateVertex()

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

Definition at line 361 of file BesSensitiveManager.cc.

{
  BesTruthTrack* truthTrack = (*m_trackList)[m_trackIndex];
  G4int terminalIndex = truthTrack->GetTerminalVertex()->GetIndex();
  if(m_logLevel<=4)  
    G4cout<<"updateVertex:"<<terminalIndex<<" ";
  BesTruthVertex* vertex = (*m_vertexList)[terminalIndex];
  G4StepPoint *finalStep = track->GetStep()->GetPostStepPoint();
  const G4ThreeVector pos(finalStep->GetPosition());
  G4double time = finalStep->GetGlobalTime();
  if(m_logLevel<=4)
    G4cout<<"newPos:"<<pos<<" "<<"newTime:"<<time<<G4endl;
  vertex->SetPosition(pos);
  vertex->SetTime(time);
  vertex->SetTerminal( true );
  //G4int vindex = m_vertexList->size();
  //stat.vertices.push_back(terminalIndex);
 
}

Referenced by EndOfTrack().

Member Data Documentation

◆ chain

std::vector<BesTStats> BesSensitiveManager::chain

protected

Definition at line 148 of file BesSensitiveManager.hh.

Referenced by BeginOfTrack(), CheckDecayTrack(), EndOfTrack(), FollowTrack(), and MakeNewTrack().

◆ clients

std::vector<BesSensitiveDetector*> BesSensitiveManager::clients

protected

Definition at line 146 of file BesSensitiveManager.hh.

Referenced by AddSensitiveDetector(), BeginOfTrack(), BeginOfTruthEvent(), EndOfTrack(), and EndOfTruthEvent().

◆ iter

std::vector<BesSensitiveDetector*>::iterator BesSensitiveManager::iter

protected

Definition at line 130 of file BesSensitiveManager.hh.

Referenced by BeginOfTrack(), BeginOfTruthEvent(), EndOfTrack(), and EndOfTruthEvent().

◆ iter_end

std::vector<BesSensitiveDetector*>::iterator BesSensitiveManager::iter_end

protected

Definition at line 131 of file BesSensitiveManager.hh.

Referenced by BeginOfTrack(), BeginOfTruthEvent(), EndOfTrack(), and EndOfTruthEvent().

◆ m_count

G4int BesSensitiveManager::m_count

protected

Definition at line 158 of file BesSensitiveManager.hh.

Referenced by BeginOfTruthEvent(), EndOfTrack(), and MakeNewTrack().

◆ m_hepmcevt

HepMC::GenEvent* BesSensitiveManager::m_hepmcevt

protected

Definition at line 164 of file BesSensitiveManager.hh.

Referenced by BesSensitiveManager(), FollowTrack(), and SaveParticlesFromGenerator().

◆ m_logLevel

G4int BesSensitiveManager::m_logLevel

protected

Definition at line 163 of file BesSensitiveManager.hh.

Referenced by CheckDecayTrack(), CheckType(), EndOfTrack(), FollowTrack(), MakeNewTrack(), SaveParticlesFromGenerator(), SetLogLevel(), SetVertex0(), UpdatePrimaryTrack(), and UpdateVertex().

◆ m_pos0

G4ThreeVector BesSensitiveManager::m_pos0

protected

Definition at line 160 of file BesSensitiveManager.hh.

Referenced by FollowTrack(), and SetVertex0().

◆ m_sensitiveManager

BesSensitiveManager * BesSensitiveManager::m_sensitiveManager =0

staticprotected

Definition at line 128 of file BesSensitiveManager.hh.

Referenced by BesSensitiveManager(), and GetSensitiveManager().

◆ m_t0

G4double BesSensitiveManager::m_t0

protected

Definition at line 161 of file BesSensitiveManager.hh.

Referenced by FollowTrack(), and SetVertex0().

◆ m_trackFlag

G4int BesSensitiveManager::m_trackFlag

protected

Definition at line 135 of file BesSensitiveManager.hh.

Referenced by BeginOfTrack(), EndOfTrack(), and FollowTrack().

◆ m_trackIndex

G4int BesSensitiveManager::m_trackIndex

protected

Definition at line 134 of file BesSensitiveManager.hh.

Referenced by CheckDecayTrack(), EndOfTrack(), GetCurrentTrackIndex(), MakeNewTrack(), UpdatePrimaryTrack(), and UpdateVertex().

◆ m_trackList

std::vector<BesTruthTrack*>* BesSensitiveManager::m_trackList

protected

Definition at line 139 of file BesSensitiveManager.hh.

Referenced by BeginOfTruthEvent(), CheckDecayTrack(), ClearEvent(), EndOfTrack(), GetNumberTracks(), GetTrackList(), MakeNewTrack(), SaveParticlesFromGenerator(), and UpdatePrimaryTrack().

◆ m_vertexList

std::vector<BesTruthVertex*>* BesSensitiveManager::m_vertexList

protected

Definition at line 140 of file BesSensitiveManager.hh.

Referenced by BeginOfTruthEvent(), ClearEvent(), EndOfTrack(), GetNumberVertices(), GetVertexList(), MakeNewTrack(), and SaveParticlesFromGenerator().

The documentation for this class was generated from the following files:

Public Member Functions

Static Public Member Functions

Protected Member Functions

Protected Attributes

Static Protected Attributes

Detailed Description

Constructor & Destructor Documentation

◆ BesSensitiveManager()

◆ ~BesSensitiveManager()

Member Function Documentation

◆ AddSensitiveDetector()

◆ BeginOfTrack()

◆ BeginOfTruthEvent()

◆ CheckDecayTrack()

◆ CheckType()

◆ ClearEvent()

◆ EndOfTrack()

◆ EndOfTruthEvent()

◆ FollowTrack()

◆ GetCurrentTrackIndex()

◆ GetDaughterVertexes()

◆ GetNumberTracks()

◆ GetNumberVertices()

◆ GetSensitiveManager()

◆ GetTrackList()

◆ GetVertexList()

◆ MakeNewTrack()

◆ MatchDaughterTrack()

◆ MatchVertex()

◆ SaveParticlesFromGenerator()

◆ SetLogLevel()

◆ SetVertex0()

◆ UpdatePrimaryTrack()

◆ UpdateVertex()

Member Data Documentation

◆ chain

◆ clients

◆ iter

◆ iter_end

◆ m_count

◆ m_hepmcevt

◆ m_logLevel

◆ m_pos0

◆ m_sensitiveManager

◆ m_t0

◆ m_trackFlag

◆ m_trackIndex

◆ m_trackList

◆ m_vertexList