Hands-on 5 Storing hits

Important:

Every time you open a new shell (tcsh or bash) you need to setup your environment by running setup-geant4.csh or setup-geant4-sh.

You also have to make sure that all other environment variables are set (G4EXE, ...)

Introduction

Before starting the exercise it would be a good idea to check the errata page for any errors or updates that may have been found recently.

To start this exercise, unpack HandsOn_5.tgz to your working directory.

Compile/link it using following commands.

$ cd HandsOn_5

$ make

You may see two warning messages because the example code is a skeleton, please do not mind.

Once you have successfully compiled and linked the example, try it to see how it runs.

$ ./bin/$G4SYSTEM/beamTest

After initialization, a prompt Idle> should be appear. At this point, enter the following UI command for executing the provided macro file "run.mac".

Idle> /control/execute run.mac

Idle> exit

You will not see any output from scorers. In this Hands-on, users are required to implement a Sensitive detector on the scoring sphere. The scoring result is given as a "Hit" which is a snapshot of each step rather than the accumulated physical quantity.

Defining Hit class

First of all, you have to implement your own sensitive detector class and hit class. Let's start from implementing a Hit class. The header file and skeleton of the implementation is provided for storing information about a strip number (copy No.) of the geometry, position and momentum of the track.

Let's look at BeamTestHit.hh, and try to add new information about kinetic energy and particle definition of the track as well as its access methods.

A snippet of BeamTestHit.hh

#ifndef BeamTestHit_h

#define BeamTestHit_h 1

#include "G4VHit.hh"

#include "G4THitsCollection.hh"

#include "G4Allocator.hh"

#include "G4ThreeVector.hh"

#include "G4ParticleDefinition.hh"

class BeamTestHit : public G4VHit

{

public:

BeamTestHit();

~BeamTestHit();

G4int operator==(const BeamTestHit &right) const;

inline void *operator new(size_t);

inline void operator delete(void *aHit);

void Draw();

void Print();

private:

G4int stripNo;

G4ThreeVector position;

G4ThreeVector momentum;

G4double energy;

G4ParticleDefinition* particle;

public:

inline void SetStripNo(G4int strip)

{ stripNo=strip; }

inline G4int GetStripNo()

{ return stripNo; }

inline void SetPosition(G4ThreeVector pos)

{ position=pos; }

inline G4ThreeVector GetPosition()

{ return position; }

inline void SetMomentum(G4ThreeVector mom)

{ momentum = mom; }

inline G4ThreeVector GetMomentum()

{ return momentum; }

inline void SetEnergy(G4double ene)

{ energy = ene; }

inline G4double GetEnergy()

{ return energy; }

inline void SetParticle(G4ParticleDefinition* pdef)

{ particle = pdef; }

inline G4ParticleDefinition* GetParticle()

{ return particle; }

};

typedef G4THitsCollection<BeamTestHit> BeamTestHitsCollection;

--- snipped ---

As you see in the code above, BeamTestHit class is derived from G4VHit base class. It can store the strip number and the position of a Hit, snapshot each step. BeamTestHitsCollection is an instantiation of the G4THitsCollection template class, which will store BeamTestHit class objects.

Defining Sensitive Detector class

Taking the above Hit class into account, let's start implementing each method of BeamTestSD.cc which is a derived class from G4VSensitiveDetector base class.

G4VSensitiveDetector class object maintains "collectionName" vector which is the name of hits collection defined in the sensitive detector object.

In the constructor, the name of the hits collection must be defined. (HCID is an integer variable for hits collection ID, and should be initialized -1.)

Constructor of BeamTestSD.cc

#include "BeamTestSD.hh"

#include "BeamTestHit.hh"

#include "G4Step.hh"

#include "G4Track.hh"

#include "G4HCofThisEvent.hh"

#include "G4TouchableHistory.hh"

BeamTestSD::BeamTestSD(G4String name)

:G4VSensitiveDetector(name)

{

G4String HCname;

collectionName.insert(HCname="hitsCollection");

HCID = -1;

}

Initialize() method is invoked at the beginning of each event. Here you must instantiate a hits collection object and set it to the G4HCofThisEvent object. In this example, the variable "hitsCollection" is defined in the header file, and an empty collection without any hit is instantiated. The collection ID is stable for all runs in this example. Thus, the collection ID is gotten only once.

Initialize() method of BeamTestSD.cc

void BeamTestSD::Initialize(G4HCofThisEvent* HCE)

{

hitsCollection = new BeamTestHitsCollection(SensitiveDetectorName,collectionName[0]);

if(HCID<0)

{ HCID = GetCollectionID(0); }

HCE->AddHitsCollection(HCID,hitsCollection);

}

The ProcessHits() is called each step in the scoring logical volume which the sensitive detector is set. In the ProcessHits() method, you should create a hit and add it to the hits collection. The hit is filled with the track information, e.g. the strip number of corresponding cell where the step takes place, position, momentum, energy, particle type of the track, in this example. The strip number of the cell is taken from the G4TouchableHistrory object.

(Note) The following example does not identify the geometry boundary as well as a surface of the geometry.

ProcessHits() method of BeamTestSD.cc

G4bool BeamTestSD::ProcessHits(G4Step* aStep, G4TouchableHistory*)

{

G4StepPoint* preStep = aStep->GetPreStepPoint();

G4TouchableHistory* touchable = (G4TouchableHistory*)(preStep->GetTouchable());

BeamTestHit* newHit = new BeamTestHit();

newHit->SetStripNo( touchable->GetReplicaNumber(0) );

newHit->SetPosition( aStep->GetPreStepPoint()->GetPosition() );

newHit->SetMomentum( aStep->GetPreStepPoint()->GetMomentum() );

newHit->SetEnergy( aStep->GetPreStepPoint()->GetTotalEnergy() );

newHit->SetParticle( aStep->GetTrack()->GetDefinition() );

hitsCollection->insert( newHit );

return true;

}

If you want to store only hits which come into the geometry, you may add following line in ProcessHits().

Note: this code still does not ensure a particular surface of the geometry.

ProcessHits() method of BeamTestSD.cc

G4bool BeamTestSD::ProcessHits(G4Step* aStep, G4TouchableHistory*)

{

G4StepPoint* preStep = aStep->GetPreStepPoint();

G4TouchableHistory* touchable = (G4TouchableHistory*)(preStep->GetTouchable());

// Ensure counting incoming tracks only.

if ( preStep->GetStepStatus() == fGeomBoundary ){

BeamTestHit* newHit = new BeamTestHit();

newHit->SetStripNo( touchable->GetReplicaNumber(0) );

newHit->SetPosition( aStep->GetPreStepPoint()->GetPosition() );

newHit->SetMomentum( aStep->GetPreStepPoint()->GetMomentum() );

newHit->SetEnergy( aStep->GetPreStepPoint()->GetTotalEnergy() );

newHit->SetParticle( aStep->GetTrack()->GetDefinition() );

hitsCollection->insert( newHit );

}

return true;

}

Setup Sensitive Detector

Now you have completed the implementation of BeamTestSD. The next step is to use BeamTestSD in BeamTestDetectorConstruction. This implementation proceeds in almost same as MultiFunctionalDetector in previous Hands-on.

A snippet of BeamTestDetectorConstruction

#include "BeamTestDetectorConstruction.hh"

--- snipped ----

include "BeamTestScoreParameterisation.hh"

#include "G4SDManager.hh"

#include "BeamTestSD.hh"

--- snipped ----

void BeamTestDetectorConstruction::SetupScorer(){

G4String SDname;

BeamTestSD* sd = new BeamTestSD(SDname = "BeamTestSD");

// The SD has to be registered to both SDManager and LogicalVolume.

G4SDManager* SDman = G4SDManager::GetSDMpointer();

SDman->AddNewDetector( sd );

fScoreLog->SetSensitiveDetector(sd);

}

Getting output of hits collection

Finally, modify BeamTestEventAction in order to getting output of the hits collection.

1) You need to know the hits collection ID only once in BeginOfEvent() method. Since GetCollectioID() method is a heavy operation, it should not be invoked for every event.

2) In EndOfEventAction() method, you have to get a pointer to G4HCofThisEvent (Hits Collection of This Event) object from G4Event object, and then get a pointer of your hits collection using hits collection ID.

BeamTestEventAction

#include "BeamTestEventAction.hh"

#include "G4Event.hh"

#include "BeamTestHit.hh"

#include "G4HCofThisEvent.hh"

#include "G4SDManager.hh"

BeamTestEventAction::BeamTestEventAction()

{

hitsCollID = -1;

}

BeamTestEventAction::~BeamTestEventAction()

{;}

void BeamTestEventAction::BeginOfEventAction(const G4Event*)

{

G4SDManager * SDman = G4SDManager::GetSDMpointer();

if(hitsCollID<0){

G4String colNam;

hitsCollID = SDman->GetCollectionID(colNam="hitsCollection");

}

void BeamTestEventAction::EndOfEventAction(const G4Event* evt)

{

G4cout << ">>> Event " << evt->GetEventID() << G4endl;

if(hitsCollID<0) return;

G4HCofThisEvent * HCE = evt->GetHCofThisEvent();

BeamTestHitsCollection* HC = 0;

if(HCE)

{

HC = (BeamTestHitsCollection*)(HCE->GetHC(hitsCollID));

}

if ( HC ) {

int n_hit = HC->entries();

for ( int i = 0 ; i < n_hit; i++){

G4int stripNo = (*HC)[i]->GetStripNo();

G4ThreeVector position = (*HC)[i]->GetPosition();

G4ThreeVector momentum = (*HC)[i]->GetMomentum();

G4double energy = (*HC)[i]->GetEnergy();

G4int pid = (*HC)[i]->GetParticle()->GetPDGEncoding();

G4cout << "---- Hit # " << i << G4endl;

G4cout << " Strip No " << stripNo << G4endl;

G4cout << " Position " << position/cm << " [cm] " <<G4endl;

G4cout << " Momentum " << momentum/keV << " [keV] " <<G4endl;

G4cout << " Energy " << energy/keV << " [keV] " <<G4endl;

G4cout << " PID " << pid << G4endl;

}

In this example, the implemented BeamTestEventAction class has already been registered with the G4RunManager in main program beamTest.cc.

You will find following lines in your main program beamTest.cc.

A snippet of beamTest.cc

---- snipped ------

#include "BeamTestPrimaryGeneratorAction.hh"

#include "BeamTestEventAction.hh"

#include "BeamTestPhysicsListLowEnergy.hh"

#include "QGSP_BERT.hh"

int main(int argc,char** argv) {

// Construct the default run manager

G4RunManager * runManager = new G4RunManager;

---- snipped -----

runManager->SetUserAction(new BeamTestPrimaryGeneratorAction());

runManager->SetUserAction(new BeamTestEventAction());

//Initialize G4 kernel

runManager->Initialize();

--- snipped ----

}

After compiling and linking, execute the program. After some initialization, a prompt Idle> should appear. Then execute "run.mac" for example.

You will see a following output on your display.

A snippet of output showing hits collection.

----- snipped ----

Start Run processing.

>>> Event 0

---- Hit # 0

Strip No 12

Position (-4.696553,21.802572,97.48123) [cm]

Momentum (-49.349844,213.24913,938.43824) [keV]

Energy 963.62696 [keV]

PID 22

HepRepFile writing to G4Data0.heprep

G4VisManager: Using G4TrajectoryDrawByCharge as default trajectory model.

See commands in /vis/modeling/trajectories/ for other options.

>>> Event 1

---- Hit # 0

Strip No 7

Position (-13.053411,1.6570497,99.130533) [cm]

Momentum (-106.45715,17.572421,809.98238) [keV]

Energy 817.1373 [keV]

PID 22

---- Hit # 1

Strip No 5

Position (-6.9892147,-6.8634566,99.519063) [cm]

Momentum (-313.49648,-292.12645,4518.7608) [keV]

Energy 4539.0326 [keV]

PID 22

>>> Event 2

>>> Event 3

>>> Event 4

---- Hit # 0

Strip No 48

Position (-72.023684,-19.070127,66.700219) [cm]

Momentum (-60.970418,-15.756816,54.45269) [keV]

Energy 83.251213 [keV]

---- snipped ---

You can download the complete source of this exercise from HandsOn_5.complete.tgz

Geant4 Tutorial Course

Tsukasa Aso

March 2006