Introduction to GEANT4

From UCL HEP PBT Wiki

(Difference between revisions)
Jump to: navigation, search
m
m
Line 164: Line 164:
-
GEANT4 provides a variety of physics processes. These processes are decoupled from one another and you can select those which are relevant to your simulation. The processes are grouped in seven categories and their full list is available [http://geant4.cern.ch/support/proc_mod_catalog/processes/ here]:
+
GEANT4 provides a variety of physics processes. These processes are decoupled from one another and you can select those which are relevant to your simulation. The processes are grouped in seven categories and their list is available [http://geant4.cern.ch/support/proc_mod_catalog/processes/ here]:
<div style="column-count:2;-moz-column-count:2;-webkit-column-count:2">
<div style="column-count:2;-moz-column-count:2;-webkit-column-count:2">
Line 177: Line 177:
-
In method ConstructProcess() define your physics processes:
+
For each particle defined in ConstructParticle() assign all the physics processes that you want to
 +
consider in your simulation:  
    
    
Line 246: Line 247:
=== <span style="color:#000080"> Detailed physics lists </span> ===  
=== <span style="color:#000080"> Detailed physics lists </span> ===  
-
If you want to build more realistic physics list you have to use the class <span style="color:#ff0000"> G4VModularPhysicsList </span>. For example, the photon from the example above can undergo compton scattering apart from conversion. In <span style="color:#ff0000"> G4VModularPhysicsList </span> you can group the physics processes into separate modules: EM physics, hadronic physics, optical physics etc.  
+
If you want to build more realistic physics list you have to use the class <span style="color:#ff0000"> G4VModularPhysicsList </span>. For example, the photon from the example above can undergo compton scattering apart from conversion. In <span style="color:#ff0000"> G4VModularPhysicsList </span> you can group the physics processes into separate modules: EM physics, hadronic physics, decay physics etc.  
Line 252: Line 253:
-
<span style="color:#800000"> class MyPhysicsList:public G4VModularPhysicsList {
+
<span style="color:#800000"> class MyPhysicsList:public G4VModularPhysicsList { </span>
-
<span style="color:#800000"> public:
+
<span style="color:#800000"> public: </span>
-
<span style="color:#800000"> MyPhysicsList();
+
<span style="color:#800000"> MyPhysicsList(); </span>
-
<span style="color:#800000"> ~MyPhysicsList();
+
-
<span style="color:#800000"> virtual void ConstructParticle();
+
<span style="color:#800000"> ~MyPhysicsList(); </span>
 +
 
 +
<span style="color:#800000"> virtual void ConstructParticle(); </span>
      
      
-
<span style="color:#800000"> virtual void SetCuts();
+
<span style="color:#800000"> virtual void SetCuts(); </span>
          
          
-
<span style="color:#800000"> void AddPhysicsList(const G4String& name);
+
<span style="color:#800000"> void AddPhysicsList(const G4String& name); </span>
-
<span style="color:#800000"> virtual void ConstructProcess();
+
 
 +
<span style="color:#800000"> virtual void ConstructProcess(); </span>
        
        
-
<span style="color:#800000"> private:
+
<span style="color:#800000"> private: </span>
      
      
-
<span style="color:#800000"> G4String                            fEmName;
+
<span style="color:#800000"> G4String                            fEmName; </span>
-
<span style="color:#800000"> G4VPhysicsConstructor*              fEmPhysicsList;
+
-
<span style="color:#800000"> G4VPhysicsConstructor*              fDecPhysicsList;
+
-
<span style="color:#800000"> std::vector<G4VPhysicsConstructor*>  fHadronPhys; };
+
 +
<span style="color:#800000"> G4VPhysicsConstructor*              fEmPhysicsList; </span>
-
<span style="color:#800000"> MyPhysicsList::MyPhysicsList():G4ModularPhysicsList() { </span>
+
<span style="color:#800000"> G4VPhysicsConstructor*              fDecPhysicsList; </span>
-
<span style="color:#800000"> defaultCutValue = 1.0*mm; </span>
+
<span style="color:#800000"> std::vector<G4VPhysicsConstructor*>  fHadronPhysicsList; }; </span>
-
<span style="color:#800000"> RegisterPhysics(new ProtonPhysics()); // All physics processes having to do with protons </span>
 
-
<span style="color:#800000"> RegisterPhysics(new ElectronPhysics()); // All physics processes having to do with electrons </span>
+
Now we can build the physics lists:
-
<span style="color:#800000"> RegisterPhysics(new DecayPhysics()); // Physics of unstable particles </span>
 
 +
<span style="color:#800000"> void MyPhysicsList::AddPhysicsList(const G4String& name) { </span>
-
where for example class ProtonPhysics() is defined as:
+
<span style="color:#800000"> ... </span>
 +
<span style="color:#800000"> if (name == "emstandard_opt3") { </span>
 +
   
 +
<span style="color:#800000"> fEmName = name; </span> 
 +
   
 +
<span style="color:#800000"> delete fEmPhysicsList; </span>
 +
   
 +
<span style="color:#800000"> fEmPhysicsList = new G4EmStandardPhysics_option3(); </span> 
-
<span style="color:#800000"> class ProtonPhysics():public G4VPhysicsConstructor { </span>
+
<span style="color:#800000"> } else if (name == "emlivermore") { </span>
 +
 
 +
<span style="color:#800000"> fEmName = name; </span>
 +
   
 +
<span style="color:#800000"> delete fEmPhysicsList; </span>
 +
   
 +
<span style="color:#800000"> fEmPhysicsList = new G4EmLivermorePhysics(); </span>
-
<span style="color:#800000"> public: </span>
+
<span style="color:#800000"> } else if (name == "empenelope") { </span>
 +
   
 +
<span style="color:#800000"> fEmName = name; </span>
 +
   
 +
<span style="color:#800000"> delete fEmPhysicsList; </span>
 +
   
 +
<span style="color:#800000"> fEmPhysicsList = new G4EmPenelopePhysics(); </span>
-
<span style="color:#800000"> ProtonPhysics(const G4String& name="proton"); </span>
+
<span style="color:#800000"> } else if (name == "HElastic") { </span>
 +
   
 +
<span style="color:#800000"> fHadronPhysicsList.push_back( new G4HadronHElasticPhysics()); </span>
-
<span style="color:#800000"> virtual ~ProtonPhysics(); </span>
+
<span style="color:#800000"> } else if (name == "HInelastic") { </span>
 +
 
 +
<span style="color:#800000"> fHadronPhysicsList.push_back(new G4HadronInelasticQBBC()); </span>
 +
   
 +
<span style="color:#800000"> } ... } </span>
-
<span style="color:#800000"> virtual void ConstructParticle(); </span>
 
-
<span style="color:#800000"> virtual void ConstructProcess(); } </span>
+
and
-
Particle production thresholds are defined in method SetCuts():
+
<span style="color:#800000"> void MyPhysicsList::ConstructProcess() { </span>
-
 
+
-
 
+
<span style="color:#800000"> AddTransportation(); // transportation </span>
-
<span style="color:#800000"> void MyModPhysList::SetCuts() </span>
+
 
-
 
+
<span style="color:#800000"> fEmPhysicsList->ConstructProcess(); // electromagnetic physics list </span>
-
<span style="color:#800000"> {SetCutsWithDefault();} </span>
+
 +
<span style="color:#800000"> fDecPhysicsList->ConstructProcess(); // decay physics list </span>
 +
 
 +
<span style="color:#800000"> for(size_t i=0; i<fHadronPhys.size(); i++) { // hadronic physics lists </span>
 +
 
 +
<span style="color:#800000"> fHadronPhys[i]->ConstructProcess(); } } </span>
 +
 
=== <span style="color:#000080"> Pre-packaged physics lists </span> ===
=== <span style="color:#000080"> Pre-packaged physics lists </span> ===

Revision as of 10:02, 8 July 2014

Contents

Introduction

GEANT4 is a software toolkit based on C++. In your code you have to define:

  • Your experimental setup - geometry, materials and primary particles.
  • Which physics process you are interested in.
  • You may take actions during the simulation - inspect and store results.

The interaction with GEANT4 is done via base classes.

Mandatory classes
  • G4VUserDetectorConstruction : Describe the experimental setup, geometry and materials
  • G4VUserPhysicsList : Define particles, physics processes and range cuts
  • G4VUserPrimaryGeneratorAction : Describe particle source, source dimensions, initial position, energy spectrum, angular distributions
Optional classes
  • G4UserRunAction : Define and store histograms
  • G4UserEventAction : Event selection and analysis of simulation data
  • G4UserStackingAction : Customize priority of tracks
  • G4UserTrackingAction : Decide whether a trajectory should be stored or not
  • G4UserSteppingAction : Kill, suspend, postpone a track
Manager class
  • G4RunManager : Manages the simulation process

The function main()

The function main() defines the skeleton of your simulation code. Inside the function you instantiate G4RunManager and notify it of your mandatory and optional classes. This is example main() function, where MyDetectorConstruction, MyPhysicsList, MyPrimaryGeneratorAction, MyEventAction and MyRunAction are derived classes from the GEANT4 base classes:

{

...

// Run manager construction

G4RunManager* runManager = new G4RunManager;


// mandatory user initialization classes

runManager->SetUserInitialization(new MyDetectorConstruction);

runManager->SetUserInitialization(new MyPhysicsList);


// mandatory user action classes

runManager->SetUserAction(new MyPrimaryGeneratorAction);


// optional user action classes

runManager->SetUserAction(new MyEventAction);

runManager->SetUserAction(new MyRunAction);

...

}

Experimental setup

You can define your experiment by using three base classes:

  • Describing the shape and the size of your detector: G4VSolid
  • Adding properties - material and electromagnetic field: G4Logical Volume
  • Placing it in another volume - in one or many positions: G4VPhysical Volume


Simple example of G4VUserDetectorConstruction :


G4VSolid* pBoxSolid = new G4Box(“WaterBox”, 1.*m, 2.*m, 3.*m);

G4LogicalVolume* pBoxLog = new G4LogicalVolume( pBoxSolid, water, “WaterBox”);

G4VPhysicalVolume* aBoxPhys = new G4PVPlacement( pRotation, G4ThreeVector(posX, posY, posZ), pBoxLog, “WaterBox”, pWorldLog, false, copyNo);


Your detector is always placed in a mother volume called the world volume. The world volume is defined in a similar way:


G4VSolid* pWorld = new G4Box("World",5*m,5*m,5*m);

G4LogicalVolume* pWorldLog = new G4LogicalVolume(pWorld,vacuum, "World");

G4VPhysicalVolume* pWorldPhys = new G4PVPlacement(0,G4ThreeVector(),pWorldLog,"World",0,false,0);


The elements and materials used in the experiment are defined using classes G4Element and G4Material . For example water, hydrogen and oxygen are defined as:


G4Element* H = new G4Element("Hydrogen","H",z=1.,a= 1.01*g/mole);

G4Element* O = new G4Element("Oxygen","O",z=8.,a=16.00*g/mole);


density = 1.000*g/cm3;

G4Material* H2O = new G4Material("Water",density,ncomp=2);

H2O->AddElement(H, natoms=2);

H2O->AddElement(O, natoms=1);


Physics processes

You can build your own physics list or chose from already built physics lists. To build your own physics lists, you can use two base physics list classes: G4VUserPhysicsList and G4ModularPhysicsList . The class G4VUserPhysicsList is used for simple physics lists while G4ModularPhysicsList is used to build more complex physics lists. There exist also already built pre-packaged physics lists.


Simple physics lists

If the particles in your simulation undergo a descrete number of physics processes you can use the class G4VUserPhysicsList to define them. This class has three methods:

  • ConstructParticles(): define all necessary particles;
  • ConstructProcesses(): define all necessary processes and assign them to proper particles;
  • SetCuts(): define production thresholds in terms of range;


Simple example of G4VUserPhysicsList :


class MyPhysicsList:public G4VUserPhysicsList() {

public:

MyPhysicsList();

~MyPhysicsList();

void ConstructParticle();

void ConstructProcess();

void SetCuts(); }


Now implement the methods ConstructParticle(), ConstructProcess() and SetCuts():


void MyPhysicsList::ConstructParticle() {

// Define the particles

G4Electron::ElectronDefinition();

G4Positron::PositronDefinition();

G4Proton::ProtonDefinition();

G4Neutron::NeutronDefinition();

G4Gamma::GammaDefinition(); ... }


GEANT4 provides a variety of physics processes. These processes are decoupled from one another and you can select those which are relevant to your simulation. The processes are grouped in seven categories and their list is available here:

  • electromagnetic
  • hadronic
  • decay
  • photolepton-hadron
  • optical
  • parameterization
  • transportation


For each particle defined in ConstructParticle() assign all the physics processes that you want to consider in your simulation:


void MyPhysicsList::ConstructProcess() {

AddTransportation(); // Assign transportation process to all particles

ConstructEM(); // Electromagnetic processes

ConstructGeneral(); // Other processes }


In methods ConstructEM() and ConstructGeneral() assign the physics processes to the corresponding particles:


void MyPhysicsList::ConstructEM() {

aParticleIterator->reset();

while((*aParticleIterator)()){

G4ParticleDefinition* particle = aParticleIterator->value();

G4ProcessManager* pmanager = particle->GetProcessManager();

G4String particleName = particle->GetParticleName();

if (particleName == "gamma") {

pmanager->AddDiscreteProcess(new G4GammaConversion); ...}


void MyPhysicsList::ConstructGeneral() {

G4Decay* theDecayProcess = new G4Decay()

aParticleIterator->reset();

while((*aParticleIterator)()) {

G4ParticleDefinition* particle = aParticleIterator->value();

G4ProcessManager* pmanager = particle->GetProcessManager();

if theDecayProcess->IsApplicable(*particle)) {

pmanager->AddProcess(theDecayProcess);

pmanager->SetProcessOrdering(theDecayProcess,idxPostStep);

pmanager->SetProcessOrdering(theDecayProcess,idxAtRest); }}}


This is the full list of physics processes every particle can participate in. Finally, method SetCuts() is defined as:


void MyPhysicsList::SetCuts() {

defaultCutValue = 1.0*mm;

SetCutValue(defaultCutValue, "gamma");

SetCutValue(defaultCutValue, "e+");

SetCutValue(defaultCutValue, "e-");


Detailed physics lists

If you want to build more realistic physics list you have to use the class G4VModularPhysicsList . For example, the photon from the example above can undergo compton scattering apart from conversion. In G4VModularPhysicsList you can group the physics processes into separate modules: EM physics, hadronic physics, decay physics etc.


Simple example of G4VModularPhysicsList :


class MyPhysicsList:public G4VModularPhysicsList {

public:

MyPhysicsList();

~MyPhysicsList();

virtual void ConstructParticle();

virtual void SetCuts();

void AddPhysicsList(const G4String& name);

virtual void ConstructProcess();

private:

G4String fEmName;

G4VPhysicsConstructor* fEmPhysicsList;

G4VPhysicsConstructor* fDecPhysicsList;

std::vector<G4VPhysicsConstructor*> fHadronPhysicsList; };


Now we can build the physics lists:


void MyPhysicsList::AddPhysicsList(const G4String& name) {

...

if (name == "emstandard_opt3") {

fEmName = name;

delete fEmPhysicsList;

fEmPhysicsList = new G4EmStandardPhysics_option3();

} else if (name == "emlivermore") {

fEmName = name;

delete fEmPhysicsList;

fEmPhysicsList = new G4EmLivermorePhysics();

} else if (name == "empenelope") {

fEmName = name;

delete fEmPhysicsList;

fEmPhysicsList = new G4EmPenelopePhysics();

} else if (name == "HElastic") {

fHadronPhysicsList.push_back( new G4HadronHElasticPhysics());

} else if (name == "HInelastic") {

fHadronPhysicsList.push_back(new G4HadronInelasticQBBC());

} ... }


and


void MyPhysicsList::ConstructProcess() {

AddTransportation(); // transportation

fEmPhysicsList->ConstructProcess(); // electromagnetic physics list

fDecPhysicsList->ConstructProcess(); // decay physics list

for(size_t i=0; i<fHadronPhys.size(); i++) { // hadronic physics lists

fHadronPhys[i]->ConstructProcess(); } }


Pre-packaged physics lists

Some built in pre-packaged physics lists are available here. You can use them as a starting point of your simulation.


Simple example of pre-packaged physics list :


In function main() :

G4PhysListFactory factory* physListFactory = new G4PhysListFactory();

G4VUserPhysicsList* physicsList = physListFactory->GetReferencePhysList(“FTFP_BERT”);

runManager->SetUserInitialization(physicsList);


Most of the pre-packaged physics lists use "standard" electromagnetic physics processes. The "standard" EM processes are defined with classes: G4EmStandardPhysics, G4EmStandardPhysics_option1, G4EmStandardPhysics_option2 and G4EmStandardPhysics_option3. If you want to use "low energy" electromagnetic physics processes (defined with classes G4EmLivermorePhysics, G4EmLivermorePolarizedPhysics, G4EmPenelopePhysics and G4EmDNAPhysics) in the pre-packaged physics list you have to define them in your physics list class (see example ProtonPencilBeam).


For example, if you want to simulate clinical proton beam of energy 150 MeV you can use pre-packaged physics list e.g. QGSP_BIC, QGSP_BERT and FTFP_BERT. If you are interested in Bragg curve physics, use a physics list ending in "EMV" or "EMX" e.g. QGSP_BERT_EMV.


Generate primary particles

G4VUserPrimaryGeneratorAction is a mandatory user action class to control the generation of primary particles. The particle generation is done via classes G4ParticleGun and G4GeneralParticleSource .

  • G4ParticleGun is used to simulate a beam of particles. It shoots a primary particle of a certain energy and direction from a given point at a given time.
  • G4GeneralParticleSource simulates a beam of particles and the primary vertex is randomly chosen on surface of a given volume with pre-defined energy spectra, spatial and angular distribution.

Optional user classes

Personal tools