How can I customize trajectory?


You can customize the trajectory and the trajectory point classes for

Base classes of trajectory and trajectory point are G4VTrajectory and G4VTrajectoryPoint, respectively, and Geant4 provides concrete classes, G4Trajectory and G4TrajectoryPoint, which could be utilized for your example.

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 "G4ios.hh" 
#include "g4std/vector"
#include "globals.hh" 
#include "G4ParticleDefinition.hh" 
#include "G4TrajectoryPoint.hh"
#include "G4Track.hh"
#include "G4Step.hh"

typedef G4std::vector T01TrajectoryPointContainer;

class G4Polyline;                   // Forward declaration.

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);} 

// Get/Set functions 
     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; }

// Other member functions
     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;

   public:
     virtual int GetPointEntries() const
     { return positionRecord->size(); }
     virtual G4VTrajectoryPoint* GetPoint(G4int i) const 
     { return (*positionRecord)[i]; }
};

extern G4Allocator 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 "G4TrajectoryPoint.hh"
#include "G4ParticleTable.hh"
#include "G4ParticleTypes.hh"
#include "G4ThreeVector.hh"
#include "G4Polyline.hh"
#include "G4Circle.hh"
#include "G4Colour.hh"
#include "G4VisAttributes.hh"
#include "G4VVisManager.hh"

G4Allocator myTrajectoryAllocator;

T01Trajectory::T01Trajectory()
{
   fpParticleDefinition = 0;
   ParticleName = "";
   PDGCharge = 0;
   PDGEncoding = 0;
   fTrackID = 0;
   fParentID = 0;
   positionRecord = 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()));
}

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;isize();i++)
  {
    G4TrajectoryPoint* rightPoint = (G4TrajectoryPoint*)((*(right.positionRecord))[i]);
    positionRecord->push_back(new G4TrajectoryPoint(*rightPoint));
  }
}

T01Trajectory::~T01Trajectory()
{
  size_t i;
  for(i=0;isize();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 << "  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;ipush_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 and set it to G4TrackingManager is G4UserTrackingAction. Here is the sample implementation.

Source file (T01TrackingAction.cc)

#include "T01TrackingAction.hh"
#include "G4TrackingManager.hh"
#include "G4Track.hh"
#include "G4TrackVector.hh"
#include "T01TrackInformation.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)
{
}


Makoto Asai