Learn 1: write a simple program

1. Basic introduction

        Geant4 It is a Monte Carlo general particle transport software written in C + + language led by European Nuclear Center (CERN). Geant4 can simulate the interaction and transport process of various particles with different substances. Its open source and excellent programming performance of C + + make Geant4 widely used in nuclear physics, radiation detection, radiation medicine and other fields, and the database is very comprehensive and updated quickly. Geant4 has rich and comprehensive physical models and a list of physical models applicable to different fields for users to choose from, and users can also define them according to their needs. Geant4 is very suitable for developing special simulation programs in specific fields. For example, TOPAS and Gate are special programs in the field of Radiology based on Geant4.

         In Geant4, there are three logical levels: Step, Event and run. Run contains the total information of the whole simulation, Event contains all the transport processes of each source particle, and Step represents a certain interaction of a particle, from which the Energy, Momentum, Position, globaltime and Energy of the particle can be obtained Information such as the name of the physical process (ProcessName). The information recorded in Step can be accumulated into Event for processing, and the information accumulated in Event can be imported into run.

         Geant4 can edit complex geometry through code and use Boolean operation to support particle transport calculation from molecular scale to planetary scale. The program supports customization of various isotopes and mixtures, provides rich reference cases and templates, and has a powerful user interface and visualization engine for debugging.

2. Program composition

         When we want to build a simple Geant4 program for calculation, we need to set the geometric structure, material, particle source and other information in the program, and obtain and process the physical information of interest. Taking the energy spectrum of emitted neutrons recorded by 14MeV neutrons incident on polyethylene as an example.

2.1 construction materials

        The objects we simulate are composed of specific materials. Geant4 provides multiple classes for constructing various materials, including isotopes (G4Isotope), elements (G4Element) and materials (G4Material). Among them, elements can be composed of one or more isotopes, materials can be composed of one or more elements, or several elements and several materials can be synthesized. We can define any type of material according to this logic. Meanwhile, the source code of Geant4 has provided definitions of various basic isotopes, elements and common materials, which can be called through pointers. In the first mock exam of polyethylene, we need to build polyethylene and the air in the environment. Polyethylene is defined by hydrogen and oxygen in the quantity ratio of 2:1 as follows:

G4Element* elH = new G4Element("Hydrogen" ,"elH" , 1., 1.01  *g/mole);//Hydrogen element 
G4Element* elC = new G4Element("carbon"   ,"elC" , 6., 12.01 *g/mole);//Oxygen element

G4Material* mPE = new G4Material("PE",0.95* g/cm3,2);//Definition of nuclear ratio polyethylene 
  mPE->AddElement(elH, 2); 
  mPE->AddElement(elC, 1);

        We can directly call the material library provided by Geant4

 G4NistManager* man = G4NistManager::Instance();
 G4Material *mAir  = man->FindOrBuildMaterial("G4_AIR");//air

         Common materials in G4NistManager can be viewed in the source code (geant4.10.05\source\materials\src\G4NistMaterialBuilder.cc). Air is defined as follows:

void G4NistMaterialBuilder::NistCompoundMaterials()
  AddMaterial("G4_AIR", 0.00120479, 0, 85.7, 4, kStateGas);
  AddElementByWeightFraction( 6, 0.000124);
  AddElementByWeightFraction( 7, 0.755267);
  AddElementByWeightFraction( 8, 0.231781);
  AddElementByWeightFraction(18, 0.012827);

        c at this time, we have completed the definition of two materials, air (mpair) and polyethylene (mPE). Then we will add these two materials to the defined geometry.

2.2 construction geometry

         The geometric setting of Geant4 has a concept of mother and daughter. Each program should be named world by a basic computing space, that is, the largest mother. Generally, we take it as the environmental space, and the material is usually air or vacuum. The descendant is set in the parent. There is no parent in the world. All other bodies must be set in a unique parent. 0 or several descendants can be set inside each body. In this simple program, we need to set world and polyethylene board. World is the parent, and polyethylene board is a child of world. There are many kinds of Geometry classes in Geant4, such as cube (G4Box), sphere (G4Sphere), cylinder (G4Tubs), etc. Here, we define world and polyethylene board as cubes. The size of polyethylene board is 20cm in length and width, attached is 3cm, and the size of the world is 24cm in length, width and height. The following is the specific definition:

  // World
  G4double PE_sizeXY = 20*cm, PE_sizeZ = 3*cm;//Define dimension parameters
  G4double world_sizeXY = 1.2*PE_sizeXY;
  G4double world_sizeZ  = 8*PE_sizeZ;

  G4Box* solidWorld =    
    new G4Box("World",                       //its name
       0.5*world_sizeXY, 0.5*world_sizeXY, 0.5*world_sizeZ);//its size
  G4LogicalVolume* logicWorld =                         
    new G4LogicalVolume(solidWorld,          //its solid
                        mAir,           //its material
                        "World");            //its name
  G4VPhysicalVolume* physWorld = 
    new G4PVPlacement(0,                     //no rotation
                      G4ThreeVector(),       //Position at (0,0,0)
                      logicWorld,            //its logical volume
                      "World",               //its name
                      0,                     //its mother  volume
                      false,                 //no boolean operation
                      0,                     //copy number
                      checkOverlaps);        //overlaps checking

//Polyethylene sheet
  G4Box* solidPE =    
    new G4Box("PEBox",                    //its name
        0.5*PE_sizeXY, 0.5*PE_sizeXY, 0.5*PE_sizeZ); //its size
  G4LogicalVolume* logicPE =                         
    new G4LogicalVolume(solidPE,            //its solid
                        mPE,             //its material
                        "PEBox");         //its name

G4ThreeVector pos1 = G4ThreeVector(0*cm, 0*cm, -1*cm);//Define location coordinates
  new G4PVPlacement(0,                //no rotation
                    pos1,         	  //Position
                    logicPE,          //its logical volume
                    "PEBox",          //its name
                    logicWorld,       //its mother  volume
                    false,            //no boolean operation
                    0,                //copy number
                    checkOverlaps);   //overlaps checking

2.3. Construction of particle source

There are two methods to build a particle source in Geant4. One is to use the particle gun (fParticleGun) to directly define it in the PrimaryGeneratorAction, and the other is to use the macro command (g4general particle source, GPS). Here, we first use fParticleGun to set the shape, size, position, type, energy and emission direction of the particle source. The following is an isotropic neutron point source, Energy is 2.2MeV:

: G4VUserPrimaryGeneratorAction(),fParticleGun(0)
  G4int n_particle = 1;
  fParticleGun  = new G4ParticleGun(n_particle);
  G4ParticleDefinition* particle = G4ParticleTable::GetParticleTable()->FindParticle("neutron");
  G4double Energy=2.2 ;//gamma energy
  G4double Px=0,Py=0,Pz=-100;//Position parameters

  fParticleGun->SetParticlePosition( G4ThreeVector(Px,Py,Pz) );//Set position
  fParticleGun->SetParticleEnergy(Eg);//Set energy
  delete fParticleGun;
void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
  //Set isotropic source
  //Generate random number
  float Seed1=G4UniformRand();//Generate random number (0,1)
  float Seed2=G4UniformRand();
  G4double cosTheta = 2* Seed1 - 1., phi = twopi* Seed2;
  G4double sinTheta = std::sqrt(1. - cosTheta*cosTheta);
  G4double ux = sinTheta*std::cos(phi),
           uy = sinTheta*std::sin(phi),
           uz = cosTheta;

  fParticleGun->GeneratePrimaryVertex(anEvent);//Reset each Event

2.4 obtaining information

        It should be noted here that each particle in Geant4 can be extracted from generation to disappearance (cut-off value control), each interaction, and the position, time, momentum and energy changes of each interaction. From the generation of each initial particle to the end of the initial particle and its secondary particles, it is an Event. All particles have their own codes to identify, such as 22 γ, 11 stands for electronics. Each interaction between the initial particle and its secondary particle is a Step, and the type of interaction can be determined through ProcessName. Each particle has a track with an exclusive TrackID, and the TrackID of the initial particle is 1. The track of each particle is composed of several Step points in series. The steps of each track are numbered from 1, Cheng Wei StepNumber. In addition, the TrackID of the parent particle generating the particle can be determined through the ParentID. For example, the TrackID of the initial neutron is 1 and the code of the neutron is 2112. After the neutron is generated by the particle source, it has an elastic collision with the hydrogen nucleus for the first time, resulting in a recoil proton. The code of the proton is 2212, the TrackID of the proton is 2 and the ParentID of the proton is 1, indicating that it is generated by the initial neutron. The ProcessName of the interaction of this Step is "hadElastic", and the StepNumber is 1. Here, this StepNumber belongs to the initial neutron. At this time, the energy and momentum of neutrons and protons can also be obtained. It should be noted separately that in Geant4, there is a mandatory Step with ProcessName of "Transportation". This process will occur when the track of a particle passes through the boundary of an individual and enters another individual, which is convenient for us to analyze the escape of particles from the entity. We need this process to record the energy of neutrons passing through the polyethylene plate. When neutrons are incident into polyethylene plate, they will have elastic scattering, inelastic scattering and radiation capture with hydrogen and carbon nuclei. Some are absorbed, some are scattered and emitted from other directions, and some pass through without interaction. Both transmission and backscattering will occur. We want to record the internal energy of neutrons emitted from the back of the polyethylene plate. In addition to judging that neutrons move from the polyethylene to the outside World, we also need to limit the penetration position on the back through the coordinates of the crossing point:

 G4double Energy01;
 using std::ofstream;
 ofstream outfile1("./data/neutron E.txt");
 if((aStep->GetTrack()->GetDefinition()->GetPDGEncoding()==2112)//Judged as neutron
     && aStep->GetPostStepPoint()->GetPhysicalVolume() != NULL)//Step does not exceed the calculation boundary
        if(aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() == "PEBox"//The previous point is in polyethylene
    	&& aStep->GetPostStepPoint()->GetPhysicalVolume()->GetName()== "World"//The latter point is in the world
        && aStep->GetTrack()->GetPosition0().z()>4.9)//The position of Step, and the boundary z coordinate is 5

	            outfile1 << Energy01 << G4endl;

2.5 program operation

        Create a new data folder in the program directory to store data. In init_ Set the running parameters in vis.mac, you can set the number of running particles and then run directly, or call the visual interface for viewing and operation.

# Macro file for the initialization of example TestEm4
# in interactive session
# Set some default verbose
/control/verbose 2
/run/verbose 2
# Change the default number of threads (in multi-threaded mode)
/run/numberOfThreads 1
# Initialize kernel
#Set truncation value
/run/setCut  0.1 mm
# Visualization setting

#Operation visualization interface
#/control/execute vis.mac
#run directly
/run/beamOn 20000000 

Input make -j in the terminal to compile, and then run the program.

Tags: C++ Ubuntu

Posted on Wed, 20 Oct 2021 14:00:47 -0400 by jmack159