Hands on 3B: Build a voxelized phantom, retrieve dose results
Back to Agenda
In this third hands-on you will learn how to:
- Create a semi-realistic geometry
- Collect simulation output from a voxelized phantom
- Use the event user-action to dump event dose information from hits on screen
The code for this example can be found here, while the complete solution is here. Download the exercise code and unpack it.
Follow the instructions of Hands On
1 to configure with cmake
the code and build it. Try
out the application:
$ cd <tutorial>
#Copy the exercise tarball here
$ tar xf HandsOn3B.tar.gz
$ mkdir build-HandsOn3B
$ cd build-HandsOn3B
$ cmake -DGeant4_DIR=$G4CMAKEDIR <tutorial>/HandsOn3B
$ make [-jN]
$ ./Texastut
|
This geometry should be displayed:
The geometry is the same as in Hands On
2, as we will start from here to build a two-arms setup. The first arm
is already defined, and in the first exercise you will build
the second arm completed with a voxelized phantom:
- Each arm includes 5 drift-chamber planes to measure the
position of the passing particles (in green).
- Each arm includes a hodoscope made of scintillator plates to
measure the time-of-flight of the incoming particles (in red).
- A central magnetic system to deflect the charged particles
(white cylinder).
- Exercise: A voxelized phantom
(white, blue and red voxels in the picture).
The second arm can be rotated between runs. The
magnetic-field value can also be changed. User defined UI commands allow
to change arm rotation and magnetic field value at run time.
At the end of this hands on the complete geometry will look
like:
There are 5 steps involved in this exercise to build the geometry.
The
application will compile and work correctly only when the first 4
steps are completed (however it is a good idea to try to compile at
each step to fix early trivial errors).
The last step is optional because it has the goal
to change visualization attributes (colors of geometry elements) and
has no effect on simulation results.
Reminder on different ways to create a geometry setup:
- After creating solids and logical volumes you can place physical volumes
via
G4PVPlacement
(these have been already covered in Hands On 2).
- You you can place multiple copies of the same logical volumes via
multiple placements.
- Or you can use of
G4PVParametrised
to place multiple copies
of the same volume with dimensions/position parametrised
by the copy number.
- You can also use replicas to slice a larger volume in smaller
pieces.
Check the DetectorConstruction.hh
file, since many
variables you will need are already defined there.
Implement the second hodoscope.
The second hodoscope is composed of 25 planes of dimensions:
10x40x1 cm. The hodoscopes tiles are composed of scintillator material.
Instantiate
a single shape and a single logical volume. Place 25 physical volume placements
in the second arm mother volume (this mother volume is already created).
Each tile is positioned at Y=Z=0 with
respect to the mother volume, while the X coordinates depends on the
tile numnber.
Hint: Check what is done for the hodoscope of the first
arm. Remember dimensions passed to Geant4 solid classes are half
dimensions.
Solution
DetectorConstruction.cc File: |
// =============================================
// Exercise 1
// Complete the full geometry.
// Note that second arm, by default is rotated of
// 30 deg.
// Step 1: Add an hodoscope with dimensions (X,Y,Z):
// (10,40,1)cm made of scintillator.
// There are 25 planes placed at Y=Z=0 (w.r.t. mother volume)
// hodoscopes in second arm
G4VSolid* hodoscope2Solid = new G4Box("hodoscope2Box",5.*cm,20.*cm,0.5*cm);
fHodoscope2Logical = new G4LogicalVolume(hodoscope2Solid,scintillator,"hodoscope2Logical");
for (G4int i=0;i<25;i++)
{
G4double x2 = (i-12)*10.*cm;
new G4PVPlacement(0,G4ThreeVector(x2,0.,0.),fHodoscope2Logical,"hodoscope2Physical",secondArmLogical,false,i,checkOverlaps);
}
|
Build the drift chambers.
The second arm contains 5 drift chambers made of argon gas with
dimensions 300x60x2 cm. These are equally spaced inside the second arm
starting from -2.5 m to -0.5 m along the Z coordinate.
Hint: Use same methods used for step 1.
Solution
DetectorConstruction.cc File: |
// Step 2: Add 5 drift chambers made of argon, with dimensions (X,Y,Z):
// (300,60,2)cm
// These are placed equidistant inside the second arm at distances from -2.5m
// to -0.5m
// drift chambers in second arm
G4VSolid* chamber2Solid = new G4Box("chamber2Box",1.5*m,30.*cm,1.*cm);
G4LogicalVolume* chamber2Logical = new G4LogicalVolume(chamber2Solid,argonGas,"chamber2Logical");
for (G4int i=0;i<5;i++)
{
G4double z2 = (i-2)*0.5*m - 1.5*m;
new G4PVPlacement(0,G4ThreeVector(0.,0.,z2),chamber2Logical,"chamber2Physical",secondArmLogical,false,i,checkOverlaps);
}
|
Add a virtual wire plane in the drift chambers.
Add a plane of wires in the drift chambers of step
2. To simplify our problem we do not describe the single wires,
instead we add a new argon-filled volume of dimensions 300x60x0.02 cm
in the center of each of the five drift chambers.
This exercise is technically simple (a single placement), however it
shows a very useful concept: we create a single instance of this
volume and we place it once inside the mother logical volume (the
drift chamber logical volume), since the mother volume is repeated
five times, each chamber gets its own wire plane. We are
reducing the number of class instances needed for the description
of our geometry (and thus reducing the memory footprint of our
application, beside making the code more compact and readable).
Solution
DetectorConstruction.cc File: |
// Step 3: Add a virtual wire plane of (300,60,0.02)cm
// at (0,0,0) in the drift chamber
// virtual wire plane
G4VSolid* wirePlane2Solid = new G4Box("wirePlane2Box",1.5*m,30.*cm,0.1*mm);
fWirePlane2Logical = new G4LogicalVolume(wirePlane2Solid,argonGas,"wirePlane2Logical");
new G4PVPlacement(0,G4ThreeVector(0.,0.,0.),fWirePlane2Logical, "wirePlane2Physical",chamber2Logical,false,0,checkOverlaps);
|
Build a voxelized phantom.
A voxelized phantom has the goal to measure dose deposited per voxel.
Its dimensions are such that an electron or
gamma of the typical beam energy is fully absorbed, while hadrons
(such as protons), only leave a fraction of their
energy the phantom (because it is too
short).
Build a Phantom of (200,200,120)cm.
The Phantom is divided into 10 x 10 x 6 voxels.
We use Geant4's Nested Parameterization technique.
This allows us to have a large number of voxels,
each with its own material assignment,
without having to actually instantiate a separate G4Box for each voxel.
It gives both efficient memory use and fast navigation
even when there are tens of millions of voxels (a typical patient
scan has 512 x 512 x 256 = 70 million voxels).
This class parameterises the position of each cell depending
on its copy number.
The parametrization
class is already available for you in
PhantomParametrisation
. Check the method
PhantomParametrisation::ComputeTransformation(...)
to
understand how the phantom voxels are implemented.
The phantom should be placed at 2 m downstream along Z in the second arm
mother volume.
Solution
DetectorConstruction.cc File: |
// Step 4: Build a Phantom of (200,200,120)cm
// placed at (0,0,2)m in the second arm.
// The Phantom is divided into 10 x 10 x 6 voxels.
// We use Geant4's Nested Parameterization technique.
// This allows us to have a large number of voxels,
// each with its own material assignment,
// without having to actually instantiate a separate G4Box for each voxel.
// It gives both efficient memory use and fast navigation
// even when there are tens of millions of voxels (a typical patient
// scan has 512 x 512 x 256 = 70 million voxels).
// This class parameterises the position of each cell depending
// on its copy number.
// Set dimensions of phantom
G4ThreeVector phantomSize;
phantomSize.setX(200*cm);
phantomSize.setY(200*cm);
phantomSize.setZ(120*cm);
// Set number of voxels
fnxVoxels = 10;
fnyVoxels = 10;
fnzVoxels = 6;
// Calculate size of voxel
G4ThreeVector voxelSize;
voxelSize.setX(phantomSize.x()/(G4double)fnxVoxels);
voxelSize.setY(phantomSize.y()/(G4double)fnyVoxels);
voxelSize.setZ(phantomSize.z()/(G4double)fnzVoxels);
// Create overall phantom volume
G4VSolid* phantomSolid = new G4Box("PhantomBox",phantomSize.x()/2.,phantomSize.y()/2.,phantomSize.z()/2.);
fPhantomLogical = new G4LogicalVolume(phantomSolid,water,"PhantomLogical");
new G4PVPlacement(0,G4ThreeVector(0.,0.,2.*m),fPhantomLogical,
"PhantomPhysical",secondArmLogical,
false,0,checkOverlaps);
// Slice Phantom into Y slices using the Replica Volume technique.
G4String yRepName("RepY");
G4VSolid* solYRep = new G4Box(yRepName,phantomSize.x()/2.,voxelSize.y()/2.,phantomSize.z()/2.);
G4LogicalVolume* logYRep = new G4LogicalVolume(solYRep,water,yRepName);
new G4PVReplica(yRepName,logYRep,fPhantomLogical,kYAxis,fnyVoxels,voxelSize.y());
// Further slice these along X, again using Replica Volume techniqe.
G4String xRepName("RepX");
G4VSolid* solXRep = new G4Box(xRepName,voxelSize.x()/2.,voxelSize.y()/2.,phantomSize.z()/2.);
G4LogicalVolume* logXRep = new G4LogicalVolume(solXRep,water,xRepName);
new G4PVReplica(xRepName,logXRep,logYRep,kXAxis,fnxVoxels,voxelSize.x());
// If we now sliced these again into Z slices using the Replica Volume technique,
// we would have our voxels, but there would be no way to have each voxel get its own material.
// So instead we use Parameterization for this third axis of slicing.
// Define a single box the size of just one voxel.
// Note that while we assign the material of this voxel to water here,
// that material assignment will be overridden in the actual parameterization
G4String zVoxName("Voxel");
G4VSolid* solVoxel = new G4Box(zVoxName,voxelSize.x()/2.,voxelSize.y()/2.,voxelSize.z()/2.);
fVoxelLogical = new G4LogicalVolume(solVoxel,water,zVoxName);
// Create a vector containing the set of materials we will be using
std::vector phantomMat(3);
phantomMat[0]=air;
phantomMat[1]=water;
phantomMat[2]=bone;
// Define that parameterization that we will be using.
PhantomParameterization* paramPhantom = new PhantomParameterization(voxelSize/2.,fnzVoxels,phantomMat);
// Apply this parameterization.
// The parameterization will adjust the position and material of the voxel depending on which actual
// voxel it represents at any given placement.
new G4PVParameterised("PhantomVoxels",fVoxelLogical,logXRep,kUndefined,fnzVoxels,paramPhantom);
|
Provide visualization attributes for the second arm volumes.
Note
that overall box of the phantom, and the y and x divisions are made invisible
so that their colors do not override the colors of the individual voxels.
The colors of the individual voxels are set in PhantomParametrisation
depending on the material of the specific voxel.
Solution
DetectorConstruction File: |
// visualization attributes ------------------------------------------------
// Step 5: uncomment visualization attributes of the newly created volumes
G4VisAttributes* visAttributes = new G4VisAttributes(G4Colour(1.0,1.0,1.0));
visAttributes->SetVisibility(false);
worldLogical->SetVisAttributes(visAttributes);
fVisAttributes.push_back(visAttributes);
visAttributes = new G4VisAttributes(G4Colour(0.9,0.9,0.9)); // LightGray
fMagneticLogical->SetVisAttributes(visAttributes);
fVisAttributes.push_back(visAttributes);
visAttributes = new G4VisAttributes(G4Colour(1.0,1.0,1.0));
visAttributes->SetVisibility(false);
firstArmLogical->SetVisAttributes(visAttributes);
secondArmLogical->SetVisAttributes(visAttributes);
fVisAttributes.push_back(visAttributes);
visAttributes = new G4VisAttributes(G4Colour(0.8888,0.0,0.0));
fHodoscope1Logical->SetVisAttributes(visAttributes);
fHodoscope2Logical->SetVisAttributes(visAttributes);
fVisAttributes.push_back(visAttributes);
visAttributes = new G4VisAttributes(G4Colour(0.0,1.0,0.0));
chamber1Logical->SetVisAttributes(visAttributes);
chamber2Logical->SetVisAttributes(visAttributes);
fVisAttributes.push_back(visAttributes);
visAttributes = new G4VisAttributes(G4Colour(0.0,0.8888,0.0));
visAttributes->SetVisibility(false);
fWirePlane1Logical->SetVisAttributes(visAttributes);
fWirePlane2Logical->SetVisAttributes(visAttributes);
fVisAttributes.push_back(visAttributes);
visAttributes = new G4VisAttributes(G4Colour(0.8888,0.8888,0.0));
visAttributes->SetVisibility(false);
fPhantomLogical->SetVisAttributes(visAttributes);
fVisAttributes.push_back(visAttributes);
visAttributes = new G4VisAttributes(G4Colour(0.9,0.9,0.0));
visAttributes->SetVisibility(false);
logYRep->SetVisAttributes(visAttributes);
fVisAttributes.push_back(visAttributes);
visAttributes = new G4VisAttributes(G4Colour(0.9,0.9,0.0));
visAttributes->SetVisibility(false);
logXRep->SetVisAttributes(visAttributes);
fVisAttributes.push_back(visAttributes);
|
In this exercise we will cover basic aspects of retrieving
physics quantities from the simulation kernel.
The basic simulation output is called a
hit (a user-defined class inheriting from G4VHit
):
an energy deposit in space and time. Typically we are not
interested in hits in all detector elements, but instead we want to
retrieve information only for the relevant detector
components, to simulate the detector read-out (e.g. the voxels in the phantom).
In Geant4 this is achieved with the concepts of hits and
sensitive detectors (SD): you can attach a SD (a user class
inheriting from G4VSensitiveDetector
) to a logical
volume, in this way Geant4 will call your user-code when a particle is
tracked in this specific volume. Information can be retrieved from the
G4Step
(e.g. energy
deposited along the step) and a new hit is created (or an
existing hit is updated). Geant4 will keep track of all hits created
in the application. These can be retrieved at the end of the event for further
post-processing and writing to output.
We will show how to measure a quantity, for each event,
from the hodoscopes. The goal is to measure the accumulated dose in each voxel
that has any hits.
The exercise is divided in three parts, and you will have to modify
four files:
PhantomHit.hh
and PhantomHit.cc
files
implement the hit class for the Phantom.
PhantomSD.cc
implements the Phantom sensitive
detector.
DetectorConstruction.cc
instantiates the sensitive detector
and attaches it to the correct logical volume.
Create a hit class.
This concrete Hit class represents a data container for only two
quantities: an integer value, representing the index of the Phantom voxel
that fired; and a double value, representing the accumulated dose in the
Phantom voxel that was hit.
You will need to modify the PhantomHit
class. The class skeleton is
already prepared, you should add two data members that identify which Phantom
voxel has been hit and accumulate the dose of from the hit.
Note, that empty Constructor, the operators new and delete have been already
implemented. You should remove the empty implementation and
implement the correct methods. Implement/modify the Print
method to dump
the hit content.
Important note on operator new
and
operator delete
: hits can put some
pressure on CPU, because, for each event, many hits may be created
and deleted at the end of the event. Allocating on the heap is a
(relatively) CPU-intensive operation, thus the handling of hits may
cause some
performance degradation in a complex application.
To mitigate this we
use an allocator that allows for an efficient re-use memory and
avoid many calls to new/delete.
The first time a hit is created a memory pool is created that can hold
(like in an array) many hits. Each time a hit is created
with new operator we first look in this pool for an available
pre-allocated memory location. If an empty slot is available, we
re-use it, otherwise we grow the pool to contain more
hits.
With this technique we reduce substantially the new/delete cycles needed
for the simulation.
An additional complication is that in multi-threading
environments special attention is needed for the use of allocators.
We recognize this is an advanced topic that requires some
more advanced knowledge of C++. If you do not feel
comfortable with this discussion, you can remove from the PhantomHit.hh file
the lines defining the new and delete operators, the application will
work perfectly and since the hits are very simple and the simulation
program is not too complex you will not see any CPU penalty.
This exercise implements a single sensitive detector and one hit
type. In the other version of Hands On 3 and in
Hands On 4 additional
sensitive detectors are used with hits in the drift chambers and in
the calorimeters. You can study that code to see additional types
of hits (calorimeter hits are of some interest since accumulate energy
from several steps instead of creating a new hit at each step).
Solution
PhantomHit.hh file: |
class PhantomHit : public G4VHit
{
public:
PhantomHit(G4int i,G4double dose);
virtual ~PhantomHit() {}
inline void *operator new(size_t);
inline void operator delete(void*aHit);
void Print();
G4int GetID() const { return fId; }
void SetDose(G4double val) { fDose = val; }
G4double GetDose() const { return fDose; }
private:
G4int fId;
G4double G4double GetDose() const { return fDose; };
};
typedef G4THitsCollection<PhantomHit> PhantomHitsCollection;
extern G4ThreadLocal G4Allocator<PhantomHit>* PhantomHitAllocator;
inline void* PhantomHit::operator new(size_t)
{
if (!PhantomHitAllocator)
PhantomHitAllocator = new G4Allocator<PhantomHit>;
return (void*)PhantomHitAllocator->MallocSingle();
}
inline void PhantomHit::operator delete(void*aHit)
{
PhantomHitAllocator->FreeSingle((PhantomHit*) aHit);
}
|
PhantomHit.cc file: |
G4ThreadLocal G4Allocator<PhantomHit>* PhantomHitAllocator;
PhantomHit::PhantomHit(G4int i,G4double dose)
: G4VHit(), fId(i), fDose(dose)
{}
void PhantomHit::Print()
{
G4cout << " Phantom[" << fId << "] " << fDose << G4endl;
}
|
Create and manipulate Phantom hits.
For this exercise you will modify PhantomSD.cc
file.
Some part of the code is already implemented, in particular the
initialization of the hits collection, use this code as a reference
for your future applications: it is important to understand the details of
how the registering of hits with the Geant4 kernel works.
What you need to do for this exercise is to modify the method
ProcessHits
and implement the logic to extract voxel ID and
dose. This is the method that Genat4 kernel will call every time a
particle passes through the volume associated with this SD. The
G4Step
object encodes the information regarding the
simulation step in the geometry volume.
Hint 1: Given a G4Step
two points are defined
(G4StepPoint
) that delimit the step itself (pre- and post-).
From each
point you can retrieve which volume the step belongs to via the
touchable history:
G4TouchableHistory* touchable = static_cast<G4TouchableHistory*>(
stepPoint->GetTouchable() );
|
Since the sensitive volume is from a nested parameterization, there are
three different indices needed to uniquely identify the voxel.
These are the indices from the outermost replica volume (the Y replicas),
then the second replica volume (the X replicas),
and lastly the parameterization copy number (the parameterized volume).
From the Y, X and Z indices, we compute a single unified index we call copyNo.
G4int iY = touchable->GetReplicaNumber(2);
G4int iX = touchable->GetReplicaNumber(1);
G4int iZ = touchable->GetReplicaNumber(0);
G4int copyNo = iY * fnxVoxels * fnzVoxels + iX * fnzVoxels + iZ;
|
Hint 2: There are two G4StepPoint
defining a G4Step
, which one of the two should you use,
pre- or post- step point? Why? The answer to this question is one of
the most trickiest part of Geant4 for a new user, be sure to
understand the reason why the two points are not equivalent!
Hint 3: We are simulating a phantom that will
trigger only if some energy has been deposited (i.e. via ionization),
for example if a neutron passes through the phantom (without making
interactions) its passage should not be recorded. Check the energy
deposited in the step, if zero do not do anything.
Hint 4: More than one step can be done by the same particle in
a single volume (why?), in addition secondaries produced in the volume
will also make steps in the SD. This means that for a given primary particle we
can have more than one call to the ProcessHits
.
If the hit is the first one in a given voxel, we have to add a member to our
hits collection to hold this voxel's dose.
If the hit is not the first one in the given voxel, we already have an appropriate
member in our hits collection for that voxel, so we just accumulate the new dose
into the dose that has already been recorded in that voxel.
Solution
PhantomSD.cc file: |
G4bool PhantomSD::ProcessHits(G4Step* step, G4TouchableHistory*)
{
G4double edep = step->GetTotalEnergyDeposit();
if (edep==0.) return true;
G4StepPoint* preStepPoint = step->GetPreStepPoint();
G4TouchableHistory* touchable
= (G4TouchableHistory*)(preStepPoint->GetTouchable());
G4int iY = touchable->GetReplicaNumber(2);
G4int iX = touchable->GetReplicaNumber(1);
G4int iZ = touchable->GetReplicaNumber(0);
G4int copyNo = iY * fnxVoxels * fnzVoxels + iX * fnzVoxels + iZ;
G4cout << "iY: " << iY << ", iX:" << iX << ", iZ: " << iZ << ", copyNo: " << copyNo << G4endl;
G4double density = preStepPoint->GetMaterial()->GetDensity();
G4double dose = edep / ( density * touchable->GetSolid()->GetCubicVolume() );
// check if this voxel already has a hit
G4int ix = -1;
for (G4int i=0;ientries();i++)
{
if ((*fHitsCollection)[i]->GetID()==copyNo)
{
ix = i;
break;
}
}
if (ix>=0)
// if it has, then accumulate from previous dose
{
dose += (*fHitsCollection)[ix]->GetDose();
(*fHitsCollection)[ix]->SetDose(dose);
} else
// if not, create a new hit and set it to the collection
{
PhantomHit* hit = new PhantomHit(copyNo,dose);
fHitsCollection->insert(hit);
}
return true;
}
|
Construct the SD and attach it to the correct logical volume.
We can now create an instance of the PhantomSD and attach it
to the correct logical volume.
Give the name "/phantom" to this SD.
We are going to modify the method
ConstructSDandField
in the DetectorCostruction class.
If you are already a user of older version of Geant4
(up to version 9.6) this is one of the new
main features introduced in version 10.0 to be compatible with multi-threading.
To reduce memory consumption geometry is
shared among threads, but sensitive-detectors are not.
Solution
DetectorConstruction.cc file: |
void DetectorConstruction::ConstructSDandField()
{
// sensitive detector -----------------------------------------------------
G4SDManager* SDman = G4SDManager::GetSDMpointer();
G4String SDname;
G4VSensitiveDetector* phantomSD = new PhantomSD(SDname="/phantom", fnxVoxels, fnzVoxels);
SDman->AddNewDetector(phantomSD);
fVoxelLogical->SetSensitiveDetector(phantomSD);
// magnetic field ----------------------------------------------------------
fMagneticField = new MagneticField();
fFieldMgr = new G4FieldManager();
fFieldMgr->SetDetectorField(fMagneticField);
fFieldMgr->CreateChordFinder(fMagneticField);
G4bool forceToAllDaughters = true;
fMagneticLogical->SetFieldManager(fFieldMgr, forceToAllDaughters);
// Register the field and its manager for deleting
G4AutoDelete::Register(fMagneticField);
G4AutoDelete::Register(fFieldMgr);
}
|
In this exercise we modify one of the user-actions to print on
screen the information collected from Phantom at the end of each event.
User actions allow to interact with the simulation to retrieve and
control the simulation results at specific points during the simulation of a run.
Different user action provides
specific interfaces to control the different aspects of the simulation. For
example, the G4UserEventAction
class provides interfaces
to interact with Geant4 at the beginning and at the end of each
event. G4UserRunAction
allows for the creation of a user-custom
G4Run
object and it executes user-code at the beginning and at the
end of a run (this will be covered in the Hands On
4). G4VUserPrimaryGeneratorAction
controls the
creation of primaries, G4UserSteppingAction
allows to
retrieve information at each step (indipendentely of sensitive detectors), G4UserTrackingAction
allows for interaction with each G4Track
and finally
G4UserStackingAction
allows to control the urgency of
each new G4Track
(advanced).
Note for users of older versions of Geant4:
Multi-threading requires user actions to be thread-private (differently
from initialization classes that are shared among threads). A new user initialization class is available in
version 10: G4VUserActionInitialization
this provides a
method Build()
in which all user actions are instantiated
(this method is called by each worker thread). A second method
BuildForMaster
is called by the master thread. Among all
user actions the G4UserRunAction
is the only one that can
also be instantiated for the master
thread, this is to allow for reduction
of results from worker threads to master thread
(e.g. sum the partial results of each thread into a global
result). This will be covered in the Hands On 4.
Using a G4UserEventAction
print on screen the number
of hits and the time registered in the Phantom.
For this exercise you will need to modify in file EventAction.cc
the method EndOfEventAction
, this method is called by
Geant4 at the end of the simulation of each event. The pointer to the
current G4Event
is passed to the user-code. From this
object you will retrieve the hits collections for the
Phantom and dump to screen the collected information.
Part of the EventAction
code is already implemented.
In particular take a moment to study the method
BeginOfEventAction
: in this method we retrieve the IDs of
the two collections. Note the if
statement that allows
for an efficient search of the IDs, given the collection names, only
once. Searching with strings is a time consuming operation, this
method allows for reducing the CPU time, if many collections are
created this is an important optimization to consider.
Important: The code assumes you have called the two SDs:
"/phantom" and that they create a hit collection
called "hodosopeColl". Change these if you have modified the names.
The EventAction
is instantiated in the
ActionInitialization
class. Take a look at it and see how
the EventAction
is created.
The solution shows how to introduce some run-time checks of the
effective existence of the hits. While this is not necessary in this
simple code, this is a good code practice:
in large applications the presence of hits collections may be
decided at run time depending on the job configuration.
Solution
EventAction.cc file: |
void EventAction::EndOfEventAction(const G4Event* event)
{
// =============================================
// Exercise 3
// Print on screen the hits of the Phantom
// Step 1: Get the hits collection of this event
G4HCofThisEvent* hce = event->GetHCofThisEvent();
if (!hce)
{
G4ExceptionDescription msg;
msg << "No hits collection of this event found.\n";
G4Exception("EventAction::EndOfEventAction()","Code001", JustWarning, msg);
return;
}
// Step 2: Using the memorised ID get the collection
// corresponding to the Phantom
// Get hits collection
PhantomHitsCollection* hHC1 = static_cast<PhantomHitsCollection*>(hce->GetHC(fHHC1ID));
if ( (!hHC1) )
{
G4ExceptionDescription msg;
msg << "Some of hits collections of this event not found.\n";
G4Exception("EventAction::EndOfEventAction()","Code001", JustWarning, msg);
return;
}
//
// Print diagnostics
//
G4int printModulo = G4RunManager::GetRunManager()->GetPrintProgress();
if ( printModulo==0 || event->GetEventID() % printModulo != 0) return;
G4PrimaryParticle* primary = event->GetPrimaryVertex(0)->GetPrimary(0);
G4cout << G4endl
<< ">>> Event " << event->GetEventID() << " >>> Simulation truth : "
<< primary->GetG4code()->GetParticleName()
<< " " << primary->GetMomentum() << G4endl;
// Step 3: Dump on screen hits
G4int n_hit = hHC1->entries();
G4cout << "Phantom has scored dose in " << n_hit << " voxels." << G4endl;
for (G4int i=0;i<n_hit;i++)
{
PhantomHit* hit = (*hHC1)[i];
hit->Print();
}
}
|
Tutorial by:
Joseph Perl (perl AT slac DOT stanford DOT edu)
June 2018