Hands-On 4 Implement Scorers and Alternative Physics-Lists

Hands-on 4 Implementing scorers and alternative physics-lists

Introduction

Defining scorers

Accumulate event data into run data: Defining user-specific Run class

Output run data: Defining user-specific RunAction class

Making a plot of output data

Alternative physics-lists

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_4.tgz to your working directory.

Compile/link it using following commands.

$ cd HandsOn_4

$ make

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

HepRepFile will produce files called "G4Data0.heprep" and "G4Data1.heprep" which can be viewed by invoking the Wired visualization tool.

"G4Data1.heprep" containsinformation about trajectories. However, we need to calculate the photon yield at the spherical surface produced by electrons, as a function of solid angle and photon energy. In this Hands-on example, we will try to add scorers for counting photon yield at the spherical scoring geometry.

Defining scorers

In order to count photon yield at the spherical surface, as a function of solid angle and photon energy, we have to select an appropriate primitive scorer and SD-filter.

(Counting particle yield)

A primitive scorer G4PSSphereSurfaceCurrent class scores the number of tracks normalized by the area at the spherical surface.

(Filter photons and its energy interval)

A sensitive detector filter G4SDParticleWithEnergyFilter class filters only defined particles withina defined energy interval.

(Solid angle dependence)

We already have a spherical scoring geometry using a parameterized volume. The parameterization is defined in the BeamTestScoreParamterisation class, which stripes a spherical surface in 1 degree intervals of solid angle, i.e. the copy number of the parameterized volume indicates the solid angle. By attaching a scorer (sensitivedetector) to the corresponding logical volume, we will be able to get the angular dependence of the particle yield.

In order to count the photon yield as a function of solid angle and photon energy, you have to:

1) Include the appropriate scorer class header file

2) Setthe G4MultiFunctionalDetector class object to the corresponding logical volume and SD manager

3) Add a SD filter "G4SDParticleWithEnergyFilter" that specifies photons and the energy interval to be scored

to the primitive scorer "G4PSSphereSurfaceCurrent" , and then register the primitive scorer with the G4MultiFunctionalDetector class object.

Note:

* In this example, we will score photons with kinetic energy from 0 MeV to 12 MeV in 0.25 MeV energy bins. i.e a total of 48 bins.

* The name of the primitive scorers has to be unique in a G4MultiFunctionalDetector object, sincethat name is used for the corresponding hits collection name

A snippet of BeamTestDetectorConstruction.cc
#include "BeamTestDetectorConstruction.hh"
---- snipped ----
//
#include "BeamTestScoreParameterisation.hh"
//
// (Hands-on 4) Scorer headers
#include "G4MultiFunctionalDetector.hh"
#include "G4SDManager.hh"
#include "G4PSSphereSurfaceCurrent.hh"
#include "G4SDParticleWithEnergyFilter.hh"
//--
--- snipped ----
void BeamTestDetectorConstruction::SetupScorer(){
//
// ( Hands-On 4 ) Scorer Definition
//
// Attach SensitiveDetecotor (MultiFunctionalDetector)
// to the scoring sphere.
//
G4String MFDname = "MFDet";
G4String psName;
G4String fltName;
//
// MultiFunctionalDetector with name MFDname.
//
G4MultiFunctionalDetector* MFDet= new G4MultiFunctionalDetector(MFDname);
//
// The MFD has to be registered to both SDManager and LogicalVolume.
G4SDManager* SDman = G4SDManager::GetSDMpointer();
SDman->AddNewDetector( MFDet );
fScoreLog->SetSensitiveDetector(MFDet);
//
// Scorer for gamma with 0.25 MeV kinetic energy bins.
// This code created 48th G4PSSphereSurfaceCurrent Scorers for gamma with
// different kinetic energy bins. Those kinetic energy range and the particle
// are specified by G4SDParticleWithEnergyFilter.
// (1) Create PSScorer
// (2) Create SDFilter (optional)
// Attach SDFilter to Scorer
// (3) Register PSScorer to MultiFunctionalDetector
//
//
char name[10];
for ( G4int i = 0; i < 48; i++){
std::sprintf(name,"Current%03d",i);
G4String psName(name);
// Primitive Scorer
G4PSSphereSurfaceCurrent* scorer =
new G4PSSphereSurfaceCurrent(psName,fCurrent_In);
// Filter
G4double kmin = i*0.25*MeV;
G4double kmax = (i+1)*0.25*MeV;
G4SDParticleWithEnergyFilter* filter =
new G4SDParticleWithEnergyFilter(fltName="gamma filter");
filter->add("gamma");
filter->SetKineticEnergy(kmin,kmax);
//
// Attach the filter to primitive scorer.
scorer->SetFilter(filter);
// Attach the scorer to multi-functional-detector
MFDet->RegisterPrimitive(scorer);
}
}

After implementing the above, try to compile/link it. You should not see any difference on execution this time.

Accumulate event data into run data:Defining user-specific Run class

Accumulation of event data into run data should be done in a user-specific Run class object which should be a sub-class of G4Run class. The Hands-on example contains a BeamTestRun class for this purpose. The BeamTestRun class is generalized for accumulating event data stored by primitive scorers into run data which are also defined as G4THitsMap<double> objects.

BeamTestRun class automatically doesthe following by using the MultiFunctionalDetector names in the constructor:

1) Identify collection names

2) Obtain collection IDs for accessing Hit Collection of This Event.

3) Prepare G4THitsMap<double> objects for accumulating event data

4) Accumulate event data into run data

5) Provide access method for accumulated run data

In order to generate a user-specific run object for your application, you have to:

1) Generate a user-specific Run object in your specific RunAction class which should be a sub-class of G4UserRunAction.
Note: The Hands-on example also contains the BeamTestRunAction class for this purpose.
You need to edit BeamTestRunAction.cc to include BeamTestRun.hh, and also to generate a BeamTestRun object in the GenerateRun() method of BeamTestRunAction class.

2) Register the BeamTestRunAction class object with the RunManager.

The following shows how to generate a BeamTestRun class object.

A snippet of BeamTestRunAction.cc
#include "BeamTestRunAction.hh"
#include "G4Run.hh"
#include "G4RunManager.hh"
#include "G4UnitsTable.hh"
//
// (Hands-on 4) user-specific Run class
#include "BeamTestRun.hh"
--- snipped ----
G4Run* BeamTestRunAction::GenerateRun(){
//
// ------
// Hands-On 4 Generate User-specific Run
// Replace Run object to BeamTestRun object
// which is designed for MultiFunctionalDetector.
// MultiFunctionalDetctor names has to be given.
//
// return new G4Run()
std::vector<G4String> theSDNames(1,"MFDet");
return new BeamTestRun(theSDNames);
// ------
}

The following shows how to register the BeamTestRunAction object with the RunManager..

A snippet of beamTest.cc
--- snipped ---
#include "BeamTestDetectorConstruction.hh"
#include "BeamTestPhysicsList.hh"
#include "BeamTestPrimaryGeneratorAction.hh"
//
// (Hands-on 4) User-Specific RunAction
#include "BeamTestRunAction.hh"
//
--- snipped ----
#ifdef G4VIS_USE
// visualization manager
G4VisManager* visManager = new G4VisExecutive;
visManager->Initialize();
#endif
//
// set user action classes
//
runManager->SetUserAction(new BeamTestPrimaryGeneratorAction());
//
// Hands-on 4 user-specific RunAction class.
runManager->SetUserAction(new BeamTestRunAction);
//Initialize G4 kernel
runManager->Initialize();
// get the pointer to the User Interface manager
G4UImanager* UI = G4UImanager::GetUIpointer();
--- snipped ----

After implementing the above, try to compile/link it. You will see following output at execution.

A snippet of output
--- snipped ----
Start closing geometry.
G4GeometryManager::ReportVoxelStats -- Voxel Statistics
Total memory consumed for geometry optimisation: 1 kByte
Total CPU time elapsed for geometry optimisation: 0 seconds
Voxelisation: top CPU users:
Percent Total CPU System CPU Memory Volume
------
0.00 0.00 0.00 1k World
0.00 0.00 0.00 1k sphereLog
Voxelisation: top memory users:
Percent Memory Heads Nodes Pointers Total CPU Volume
------
55.53 0k 1 19 70 0.00 sphereLog
44.47 0k 3 14 36 0.00 World
++ MFDet/Current000 id 0
++ MFDet/Current001 id 1
++ MFDet/Current002 id 2
++ MFDet/Current003 id 3
++ MFDet/Current004 id 4
---- snipped ------
++ MFDet/Current011 id 45
++ MFDet/Current012 id 46
++ MFDet/Current013 id 47
Start Run processing.
HepRepFile writing to G4Data1.heprep
---- snipped -----

This list of collection names is printedout from the BeamTestRun class object.

You can confirm that the 48 scorers are registered for scoring and accumulating event data.

Output run data: Defining user-specific RunAction class

Now, let us try to get the final scoring result. This is done in the EndOfRunAction() method of the user-specific RunAction class object. The accumulated data for the entire run are obtained from the BeamTestRun class object.

In order to get the accumulated data, you have to:

1) Cast the G4Run object to the user-specific run object, i.e.the BeamTestRun object in the Hands-on example.

2) Access each G4THitsMap<double> object and dump the results.

Note: How to access each G4THitsMap<double> will depend on your user-specific Run class design and implementation.

BeamTestRun class has DumpAllScorer() method which dumps all results to standard output.

After implementing the above, try to compile/link it. You should see following output at execution.

A snippet of output
---- snipped ----
Run Summary
Number of events processed : 100
User=0.71s Real=1.69s Sys=0.99s
Number of Events Processed:100 events.
PrimitiveScorer RUN MFDet,Current000
Number of entries 29
copy no.: 2 Run Value : 0.00020905893
copy no.: 3 Run Value : 0.00014937331
copy no.: 4 Run Value : 0.00011622645
copy no.: 6 Run Value : 0.00024166334
copy no.: 7 Run Value : 0.00013972709
copy no.: 8 Run Value : 6.1694463e-05
copy no.: 9 Run Value : 5.5250861e-05
copy no.: 10 Run Value : 5.0039762e-05
copy no.: 12 Run Value : 8.4263893e-05
copy no.: 14 Run Value : 3.6420729e-05
copy no.: 20 Run Value : 5.2077842e-05
copy no.: 21 Run Value : 2.4881287e-05
---- snipped ----

In this example, the primitive name and copy number represent the corresponding photon energy intervals and solid angles. The Run Values are the number of photons normalized by unit area. The unit is in Geant4 default units. G4THitsMap<double> stores values only for the cells where tracks are injected , so you will not see a copy number printed out for some scorers..

(Output to file)

Let us try to write the results into files. We assume that the file names should be the name of primitive scorers with file extension of ".csv". The file format follows CSV (Comma Separated Variables).

What you have to do is:

1) Get the G4THitsMap<double> object from BeamTestRun class object

2) Create a file for the output and write the value of G4THitsMap<double> into it

Note: BeamTestRun class provides access methods. The Hands-on example uses the following methods.

G4int GetNumberOfHitsMap(); returns number of G4THitsMap<double> objects defined in BeamTestRun object.

G4THitsMap<G4double>* GetHitsMap(G4int i); returns a pointer of i-th G4THitsMap<double> object

Note: G4THitsMap has a [ ] operator taking the key value as an argument and returning the pointer of the value. If you get NULL from the [ ] operator, it does not mean the value is zero, but that the provided key does not exist. The value itself is accessible with an asterisk. It is advised to check the validity of the returned pointer before accessing the value.

A snippet of BeamTestRunAction.cc
void BeamTestRunAction::EndOfRunAction(const G4Run* aRun)
{
// Processed Number of Event
G4int NbOfEvents = aRun->GetNumberOfEvent();
G4cout < " Number of Events Processed:"
< NbOfEvents < " events. " <G4endl;
//
// ( Hands-On 4) Scorer output
//
BeamTestRun* theRun = (BeamTestRun*)aRun;
theRun->DumpAllScorer();
// Retriving Accumulatet data and output results.
//
// Dump result for all HitsMap collections.
//
G4int NhitsMap = theRun->GetNumberOfHitsMap();
for ( G4int im = 0; im < NhitsMap ; im++ ){
G4THitsMap<G4double>* RunMap =theRun->GetHitsMap(im);
if ( RunMap ) {
//
// Write into files. ( Or output for analysis tool?)
//
G4String filename = RunMap->GetName()+".csv";
std::ofstream out(filename.data());
for ( G4inti = 0; i < 180; i++){
G4double* unitAreaCurrent = (*RunMap)[i];
if( !unitAreaCurrent ) unitAreaCurrent = new G4double(0.0);
out < std::setw(5) < i < ",¥t"
< std::setw(15)< *unitAreaCurrent*mm*mm
< ",¥t"< " /mm2 "<G4endl;
}
out.close();
}
}
// ------
}

After implementing the above, try to compile/link it. You will see 48 files in your working directory.

Those file names will be Current000, Current001, ...., Current047.

Making a plot of output data

Try to execute with 100000 electrons, for example. And make a plot for that output data using OpenOffice.

Note: It is recommended to turn off visualization this case.

A snapshot of result

Alternative physics-lists

Physics-list is one of mandatory user classes where the user must define particles and their processes for his/her own simulation.

The Hands-on example uses the BeamTestPhysicsList class which defines only standard electro-magnetic processes and a decay process. As an alternative physics-list, the Hands-on example includes the BeamTestPhysicsListLowEnergy class which defines low energy electro-magnetic processes and a decay process.

Also, the Geant4 distribution contains physics-lists for hadrons in $G4INSTALL/physics_lists/hadronic. However, these physics lists are not complied automatically in Geant4 installation. User has to build library for these hadronic physics-lists.

Building hadronic physics list libraries:

Before you try to compile, make sure that environment variables are properly set:

G4INSTALL pointing to the installation directory of Geant4

G4SYSTEM should be the system name (e.g Linux-g++ , WIN32-VC, ... )

G4LISTS_BASE should be unset. If you have already set G4LISTS_BASE,

$ unset G4LISTS_BASE

Then, build libraries,

1) cd $G4INSTALL/physics_list/hadronic

2) make

The libraries will be created at $G4INSTALL/lib/plists/$G4SYSTEM/

Note: If you like to use these provided physics-lists for hadrons in your application, you may need to modify your GNUmakefile.

(GNUmakefile includes hadronic_list.gmk where that library for hadronic physics lists is defined.)

In this Hands-on example, you need to uncomment following two lines.

A snippet of hadronic_lists.gmk
---- snipped ----
#
# Select your physics lists to link against.
#
# EXTRALIBS += -lFTFC
# EXTRALIBS += -lFTFP
# EXTRALIBS += -lLBE
# EXTRALIBS += -lLHEP
# EXTRALIBS += -lLHEP_GN
# EXTRALIBS += -lLHEP_HP
# EXTRALIBS += -lLHEP_LEAD
# EXTRALIBS += -lLHEP_BERT_HP
# EXTRALIBS += -lLHEP_BIC_HP
# EXTRALIBS += -lLHEP_LEAD_HP
# EXTRALIBS += -lLHEP_PRECO
EXTRALIBS += -lQGSP_BERT
# EXTRALIBS += -lLHEP_PRECO_HP
# EXTRALIBS += -lQGSC
# EXTRALIBS += -lQGSC_LEAD
# EXTRALIBS += -lQGSC_LEAD_HP
# EXTRALIBS += -lQGSP
# EXTRALIBS += -lQGSP_GN
# EXTRALIBS += -lQGSP_HP
# EXTRALIBS += -lLHEP_BERT
# EXTRALIBS += -lLHEP_BIC
# EXTRALIBS += -lQGSP_BIC
EXTRALIBS += -lPackaging
---- snipped ----

If you would like to replace the physics-list with an alternative one, you have to include the header file of that physics-list class and register its object with the RunManager.. Note: The physics-list is an alternative; you should only register one physics list with the RunManger.

The following shows the implementation for alternative physics-lists.

A snippet of beamTest.cc
---- snipped ----
#include "BeamTestDetectorConstruction.hh"
#include "BeamTestPhysicsList.hh"
#include "BeamTestPrimaryGeneratorAction.hh"
//
// (Hands-on 4) User-Specific RunAction
#include "BeamTestRunAction.hh"
//
// (Hands-on) 4 Alternative Physics List
#include "BeamTestPhysicsListLowEnergy.hh"
#include "QGSP_BERT.hh"
int main(int argc,char** argv) {
//
// Construct the default run manager
//
G4RunManager * runManager = new G4RunManager;
//
// set mandatory initialization classes
//
BeamTestDetectorConstruction* detector = new BeamTestDetectorConstruction;
runManager->SetUserInitialization(detector);
//
// Hands-on 4 User has to choose one of these PhysicsLists.
// Select one of three physics lists.
//runManager->SetUserInitialization(new BeamTestPhysicsList);
//runManager->SetUserInitialization(new BeamTestPhysicsListLowEnergy);
runManager->SetUserInitialization(new QGSP_BERT);
G4UIsession* session=0;
---- snipped ----

After changing the physics-list, try to compile/link and execute the program.

Let us try to look difference between BeamTestPhysicsList and QGSP_BERT which is one of Geant4 hadronic physics-lists, by injecting 500 MeV proton as a primary particle.

The primary particle can be changed using the UI command as follows

$ ./bin/$G4SYSTEM/beamTest

Idle> /gun/particle proton

Idle> /gun/energy 500 MeV

Idle> /control/execute run.mac

Idle> exit

Because the BeamTestPhysicsList defines only electro-magnetic process and a decay process without hadronic interactions, if you shoot a proton as primary particle, the proton loses its energy only by the ionization process. There should not be hadronic interaction. Compared to the BeamTestPhysicsList, QGSP_BERT defines not only electro-magnetic process but also hadronic process. You should look for hadronic showers if you register QGSP_BERT to the RunManager.

Simulation with BeamTestPhysicsList by incident 500 MeV protons
Simulation with QGSP_BERT hadronic physics-list by incident 500 MeV protons

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

You can also try the advanced functionality of the example code: HandsOn.fullset.tgz.

Geant4 Tutorial Course

Gean4.v8.0p01

2006 March

Tsukasa Aso