G4Track is a transient class and there is no given way to store the generation of the particles. On the other hand, G4VTrajectory objects are kept at the end of the event if they are generated and if they are not deleted by the tracking action. This tip is for how to get access to the trajectories of ancestor particles.
Please note that the default trajectory class Geant4 provides, i.e. G4Trajectory, does not store any physically meaningful information such as momentum of the track. Thus you must create your own trajectory class which stores all necessary information.
A sample implementation of G4VTrajectory is the following.
| Header file (T01Trajectory.hh) |
#ifndef T01Trajectory_h
#define T01Trajectory_h 1
#include "G4VTrajectory.hh"
#include "G4Allocator.hh"
#include <stdlib.h>
#include "G4ThreeVector.hh"
#include "G4ios.hh"
#include "g4std/vector"
#include "globals.hh"
#include "G4ParticleDefinition.hh"
#include "G4TrajectoryPoint.hh"
#include "G4Track.hh"
#include "G4Step.hh"
class G4Polyline;
typedef G4std::vector<G4VTrajectoryPoint*> T01TrajectoryPointContainer;
class T01Trajectory : public G4VTrajectory
{
public:
T01Trajectory();
T01Trajectory(const G4Track* aTrack);
T01Trajectory(T01Trajectory &);
virtual ~T01Trajectory();
inline void* operator new(size_t);
inline void operator delete(void*);
inline int operator == (const T01Trajectory& right) const
{return (this==&right);}
inline G4int GetTrackID() const
{ return fTrackID; }
inline G4int GetParentID() const
{ return fParentID; }
inline G4String GetParticleName() const
{ return ParticleName; }
inline G4double GetCharge() const
{ return PDGCharge; }
inline G4int GetPDGEncoding() const
{ return PDGEncoding; }
inline const G4ThreeVector& GetMomentum() const
{ return momentum; }
inline const G4ThreeVector& GetVertexPosition() const
{ return vertexPosition; }
inline G4double GetGlobalTime() const
{ return globalTime; }
virtual int GetPointEntries() const
{ return positionRecord->size(); }
virtual G4VTrajectoryPoint* GetPoint(G4int i) const
{ return (*positionRecord)[i]; }
virtual void ShowTrajectory() const;
virtual void DrawTrajectory(G4int i_mode=0) const;
virtual void AppendStep(const G4Step* aStep);
virtual void MergeTrajectory(G4VTrajectory* secondTrajectory);
G4ParticleDefinition* GetParticleDefinition();
private:
T01TrajectoryPointContainer* positionRecord;
G4int fTrackID;
G4int fParentID;
G4ParticleDefinition* fpParticleDefinition;
G4String ParticleName;
G4double PDGCharge;
G4int PDGEncoding;
G4ThreeVector momentum;
G4ThreeVector vertexPosition;
G4double globalTime;
};
extern G4Allocator<T01Trajectory> myTrajectoryAllocator;
inline void* T01Trajectory::operator new(size_t)
{
void* aTrajectory;
aTrajectory = (void*)myTrajectoryAllocator.MallocSingle();
return aTrajectory;
}
inline void T01Trajectory::operator delete(void* aTrajectory)
{
myTrajectoryAllocator.FreeSingle((T01Trajectory*)aTrajectory);
}
#endif
|
| Source file (T01Trajectory.cc) |
#include "T01Trajectory.hh"
#include "G4ParticleTable.hh"
#include "G4ParticleTypes.hh"
#include "G4Polyline.hh"
#include "G4Circle.hh"
#include "G4Colour.hh"
#include "G4VisAttributes.hh"
#include "G4VVisManager.hh"
#include "G4UnitsTable.hh"
G4Allocator<T01Trajectory> myTrajectoryAllocator;
T01Trajectory::T01Trajectory()
{
fpParticleDefinition = 0;
ParticleName = "";
PDGCharge = 0;
PDGEncoding = 0;
fTrackID = 0;
fParentID = 0;
positionRecord = 0;
momentum = G4ThreeVector(0.,0.,0.);
vertexPosition = G4ThreeVector(0.,0.,0.);
globalTime = 0.;
}
T01Trajectory::T01Trajectory(const G4Track* aTrack)
{
fpParticleDefinition = aTrack->GetDefinition();
ParticleName = fpParticleDefinition->GetParticleName();
PDGCharge = fpParticleDefinition->GetPDGCharge();
PDGEncoding = fpParticleDefinition->GetPDGEncoding();
fTrackID = aTrack->GetTrackID();
fParentID = aTrack->GetParentID();
positionRecord = new T01TrajectoryPointContainer();
positionRecord->push_back(new G4TrajectoryPoint(aTrack->GetPosition()));
momentum = aTrack->GetMomentum();
vertexPosition = aTrack->GetPosition();
globalTime = aTrack->GetGlobalTime();
}
T01Trajectory::T01Trajectory(T01Trajectory & right)
{
ParticleName = right.ParticleName;
fpParticleDefinition = right.fpParticleDefinition;
PDGCharge = right.PDGCharge;
PDGEncoding = right.PDGEncoding;
fTrackID = right.fTrackID;
fParentID = right.fParentID;
positionRecord = new T01TrajectoryPointContainer();
for(int i=0;i<right.positionRecord->size();i++)
{
G4TrajectoryPoint* rightPoint = (G4TrajectoryPoint*)((*(right.positionRecord))[i]);
positionRecord->push_back(new G4TrajectoryPoint(*rightPoint));
}
momentum = right.momentum;
vertexPosition = right.vertexPosition;
globalTime = right.globalTime;
}
T01Trajectory::~T01Trajectory()
{
size_t i;
for(i=0;i<positionRecord->size();i++){
delete (*positionRecord)[i];
}
positionRecord->clear();
delete positionRecord;
}
void T01Trajectory::ShowTrajectory() const
{
G4cout << G4endl << "TrackID =" << fTrackID
<< " : ParentID=" << fParentID << G4endl;
G4cout << "Particle name : " << ParticleName
<< " Charge : " << PDGCharge << G4endl;
G4cout << "Original momentum : " <<
G4BestUnit(momentum,"Energy") << G4endl;
G4cout << "Vertex : " << G4BestUnit(vertexPosition,"Length")
<< " Global time : " << G4BestUnit(globalTime,"Time") << G4endl;
G4cout << " Current trajectory has " << positionRecord->size()
<< " points." << G4endl;
for( size_t i=0 ; i < positionRecord->size() ; i++){
G4TrajectoryPoint* aTrajectoryPoint = (G4TrajectoryPoint*)((*positionRecord)[i]);
G4cout << "Point[" << i << "]"
<< " Position= " << aTrajectoryPoint->GetPosition() << G4endl;
}
}
void T01Trajectory::DrawTrajectory(G4int i_mode) const
{
G4VVisManager* pVVisManager = G4VVisManager::GetConcreteInstance();
G4ThreeVector pos;
G4Polyline pPolyline;
for (int i = 0; i < positionRecord->size() ; i++) {
G4TrajectoryPoint* aTrajectoryPoint = (G4TrajectoryPoint*)((*positionRecord)[i]);
pos = aTrajectoryPoint->GetPosition();
pPolyline.push_back( pos );
}
G4Colour colour(0.2,0.2,0.2);
if(fpParticleDefinition==G4Gamma::GammaDefinition())
colour = G4Colour(0.,0.,1.);
else if(fpParticleDefinition==G4Electron::ElectronDefinition()
||fpParticleDefinition==G4Positron::PositronDefinition())
colour = G4Colour(1.,1.,0.);
else if(fpParticleDefinition==G4MuonMinus::MuonMinusDefinition()
||fpParticleDefinition==G4MuonPlus::MuonPlusDefinition())
colour = G4Colour(0.,1.,0.);
else if(fpParticleDefinition->GetParticleType()=="meson")
{
if(PDGCharge!=0.)
colour = G4Colour(1.,0.,0.);
else
colour = G4Colour(0.5,0.,0.);
}
else if(fpParticleDefinition->GetParticleType()=="baryon")
{
if(PDGCharge!=0.)
colour = G4Colour(0.,1.,1.);
else
colour = G4Colour(0.,0.5,0.5);
}
G4VisAttributes attribs(colour);
pPolyline.SetVisAttributes(attribs);
if(pVVisManager) pVVisManager->Draw(pPolyline);
}
void T01Trajectory::AppendStep(const G4Step* aStep)
{
positionRecord->push_back( new G4TrajectoryPoint(aStep->GetPostStepPoint()->
GetPosition() ));
}
G4ParticleDefinition* T01Trajectory::GetParticleDefinition()
{
return (G4ParticleTable::GetParticleTable()->FindParticle(ParticleName));
}
void T01Trajectory::MergeTrajectory(G4VTrajectory* secondTrajectory)
{
if(!secondTrajectory) return;
T01Trajectory* seco = (T01Trajectory*)secondTrajectory;
G4int ent = seco->GetPointEntries();
for(int i=1;i<ent;i++) // initial point of the second trajectory should not be merged
{
positionRecord->push_back((*(seco->positionRecord))[i]);
}
delete (*seco->positionRecord)[0];
seco->positionRecord->clear();
}
|
Once you create above, the class you have to implement for creating this class objects is G4UserTrackingAction. Here is the sample implementation.
| Source file (T01TrackingAction.cc) |
#include "T01TrackingAction.hh"
#include "G4TrackingManager.hh"
#include "G4Track.hh"
#include "T01Trajectory.hh"
T01TrackingAction::T01TrackingAction()
{;}
T01TrackingAction::~T01TrackingAction()
{;}
void T01TrackingAction::PreUserTrackingAction(const G4Track* aTrack)
{
fpTrackingManager->SetStoreTrajectory(true);
fpTrackingManager->SetTrajectory(new T01Trajectory(aTrack));
}
void T01TrackingAction::PostUserTrackingAction(const G4Track* aTrack)
{;}
|
Now you can get the information with the access to the trajectories of ancestor particles. seconday G4Track object. Thus, for example in your sensitive detector implementation or in your G4UserSteppingAction, you can retrieve them as following.
| Source file for the case of a sensitive detector |
#include "T01Trajectory.hh"
#include "G4Track.hh"
#include "G4Step.hh"
#include "G4ios.hh"
#include "G4UnitsTable.hh"
#include "G4RunManager.hh"
#include "G4Event.hh"
#include "G4TrajectoryContainer.hh"
#include "T01Trajectory.hh"
G4bool T01EmCalorimeter::ProcessHits(G4Step*aStep,G4TouchableHistory*ROhist)
{
.....
G4int parentID = aStep->GetTrack()->GetParentID();
G4int trackID = aStep->GetTrack()->GetTrackID();
G4cout << G4endl << "Trace back track ID " << trackID << " -- "
<< aStep->GetTrack()->GetDefinition()->GetParticleName()
<< G4endl;
while(parentID!=0)
{
T01Trajectory* parentTrajectory = GetParentTrajectory(parentID);
if(parentTrajectory==0)
{
G4cout << ">> Trajectory trace back failed - no parent found." << G4endl;
break;
}
parentID = parentTrajectory->GetParentID();
parentTrajectory->ShowTrajectory();
}
.....
return true;
}
T01Trajectory* T01EmCalorimeter::GetParentTrajectory(G4int parentID)
{
G4TrajectoryContainer* container =
G4RunManager::GetRunManager()->GetCurrentEvent()->GetTrajectoryContainer();
if(container==0) return 0;
TrajectoryVector* vect = container->GetVector();
G4VTrajectory** tr = vect->begin();
while(tr!=vect->end())
{
T01Trajectory* tr1 = (T01Trajectory*)(*tr);
if(tr1->GetTrackID()==parentID) return tr1;
tr++;
}
return 0;
}
|
Please note that this tip is rather slow. If you concern only the primary particle, it is strongly recommened to use G4VTrackInformation, you can find the tip here. Also note that this tip is valid only for the case of storing all trajectories of an event. It is not recommended to store all trajectoris of an event for high energy EM shower because of memory consumption.
Makoto Asai
|
| |