Geant4学习一：写一个简单程序

1、 基本介绍

        Geant4是由欧洲核子中心(CERN)主导开发的一款采用C++语言编写的蒙特卡洛通用粒子输运软件。Geant4能够模拟各种粒子与不同物质的相互作用以及输运过程。其开源性和C++优越的编程性能使得Geant4在核物理与辐射探测、放射性医学等领域都得到了广泛的应用，并且数据库非常全面，更新速度很快。Geant4拥有丰富、全面的物理模型以及适用于不同领域的物理模型列表可供用户选择，用户还可以根据需要自行定义。Geant4非常适合用来开发特定领域的专用模拟程序，如TOPAS、Gate都是基于Geant4开发的放射医学领域的专用程序。

        在Geant4中有三个逻辑层次，分别是Step、Event和Run，Run包含整个模拟的总信息，Event包含每个源粒子的全部输运过程，Step则表示某一个粒子的某一次相互作用，可以从中获得该粒子的能量（Energy）、动量（Momentum）、位置（Position）、输运时间（GlobleTime）、物理过程的名字（ProcessName）等信息。Step中记录的信息可以累积到Event当中处理，Event中累积的信息则可以导入到Run当中。

        Geant4可以通过代码编辑复杂的几何形状，可以使用布尔运算，支持从分子尺度到行星尺度的粒子输运计算。程序支持对各种同位素、混合物进行自定义，并提供了丰富的参考案例和模板，同时拥有强大的用户界面以及可视化引擎便于调试。

2、程序构成

        当我们要构建一个简单的Geant4程序来进行计算时需要在程序中设定几何结构、材料、粒子源等信息，并对感兴趣的物理信息进行获取和处理。以14MeV中子入射聚乙烯，记录出射的中子的能谱为例。

2.1、构建材料

        我们模拟的对象都由特定的材料构成，Geant4中提供了构建各种材料的多个类，同位素（G4Isotope）、元素（G4Element）、材料（G4Material）。其中元素可以由一个或者多个同位素构成，材料可以由一个或者多个元素构成，也可以是若干元素和若干材料合成。我们可以按照这一逻辑定义任意类型的材料。同时Geant4的源代码中已经提供了各种基本同位素、元素以及常用材料的定义可以通过指针调用。在中子入射聚乙烯这一模拟中我们需要构建聚乙烯这一材料，以及环境中的空气。聚乙烯由氢元素和氧元素按照2:1的数量比定义如下：

G4Element* elH = new G4Element("Hydrogen" ,"elH" , 1., 1.01  *g/mole);//氢元素 
G4Element* elC = new G4Element("carbon"   ,"elC" , 6., 12.01 *g/mole);//氧元素

G4Material* mPE = new G4Material("PE",0.95* g/cm3,2);//核数比定义聚乙烯 
  mPE->AddElement(elH, 2); 
  mPE->AddElement(elC, 1);

        空气我们可以直接调用Geant4提供的材料库

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

        G4NistManager中的常用材料可以在源代码中查看（geant4.10.05\source\materials\src\G4NistMaterialBuilder.cc）,其中对空气的定义如下：

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此时我们完成了两个材料的定义及空气（mAir）和聚乙烯（mPE）。随后我们将对定义好的几何结构添加这两种材料。

2.2、构建几何

        Geant4的几何设定有一个母体与子体的概念（mother&daughter），每一个程序都要由一个基本的计算空间，也就是最大的母体，命名为World，一般我们把它作为环境空间，材料通常是空气或者真空。子体是设置在母体内的，World没有母体，其他 体都必须设置在一个唯一的母体内。每个体的内部都可以设置0或者若干个子体。在这个简单程序里我们需要设置World和聚乙烯板这两个体，World作为母体，聚乙烯板是World的一个子体。Geant4中有多种几何体的类，常用的有立方体（G4Box）、球体（G4Sphere）、圆柱体（G4Tubs）等。这里我们将World和聚乙烯板都定义为立方体，聚乙烯板的尺寸为长宽20cm，后附3cm，世界的尺寸为长宽高均为24cm，下面为具体定义：

  // World
  G4double PE_sizeXY = 20*cm, PE_sizeZ = 3*cm;//定义尺寸参数
  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

//聚乙烯板
  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);//定义位置坐标
               
  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、构建粒子源

Geant4构建粒子源有两种方法，一个使用粒子枪（fParticleGun）在PrimaryGeneratorAction中直接定义，一个是使用宏命令（G4GeneralParticleSource，GPS），这里我们先使用fParticleGun，设置粒子源的形状、尺寸、位置、类型、能量以及发射方向，以下是一个各向同性的中子点源，能量为2.2MeV:

PrimaryGeneratorAction::PrimaryGeneratorAction()
: G4VUserPrimaryGeneratorAction(),fParticleGun(0)
{
  G4int n_particle = 1;
  fParticleGun  = new G4ParticleGun(n_particle);
//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&Gamma&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  G4ParticleDefinition* particle = G4ParticleTable::GetParticleTable()->FindParticle("neutron");
  fParticleGun->SetParticleDefinition(particle);
  
  G4double Energy=2.2 ;//gamma能量
  G4double Px=0,Py=0,Pz=-100;//位置参数

  fParticleGun->SetParticlePosition( G4ThreeVector(Px,Py,Pz) );//设定位置
  fParticleGun->SetParticleEnergy(Eg);//设定能量
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
PrimaryGeneratorAction::~PrimaryGeneratorAction()
{
  delete fParticleGun;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
  //设置各向同性源
  //生成随机数
  float Seed1=G4UniformRand();//生成随机数(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->SetParticleMomentumDirection(G4ThreeVector(ux,uy,uz));

  fParticleGun->GeneratePrimaryVertex(anEvent);//每个Event重新设定
}

2.4、获取信息

        这里需要说明的是Geant4中每个粒子从产生到消失（截断值控制），发生的每一次相互作用，以及每次相互作用的位置、时间、动量与能量变化均可以提取。从每个初始粒子的产生到该初始粒子以及其产生的次级粒子全部截止，为一个Event，各种粒子都有专属的代码来识别，比如22代表γ，11代表电子。初始粒子以及其产生的次级粒子所发生的每一个相互作用为一个Step，可以通过ProcessName来确定相互作用的类型。每一个粒子都有一条径迹，由一个专属的TrackID，初始粒子的TrackID为1。每个粒子的径迹由若干个Step的点串联而成，每一条径迹的Step都从1开始编号，程伟StepNumber。此外可以通过ParentID来确定产生该粒子的母粒子的TrackID。比如初始中子的TrackID为1，中子的代码为2112，该中子由粒子源产生之后第一次与氢原子核发生弹性碰撞，产生了一个反冲质子，质子的代码为2212，该质子的TrackID为2，这个质子的ParentID为1，说明这是由初始中子产生的。这一个Step的相互作用的ProcessName为“hadElastic”，StepNumber为1，这里这个StepNumber是属于初始中子的。此时还可以获取中子和质子的能量和动量。需要单独说明的的是，在Geant4中与一个强制的Step，其ProcessName为“Transportation”，在一个粒子的径迹穿过某个体的边界进入另一个体是就会产生这一过程，这便于我们分析粒子从实体中逃逸的情况。这里我们要记录穿过聚乙烯板的中子的能量，就需要用到这一过程。中子入射到聚乙烯板中会与氢核和碳核发生弹性散射、非弹性散射以及辐射俘获。一部分被吸收，一部分被散射后从其他方向出射，一部分不发生相互作用直接穿过。既有透射也会由背散射，我们想记录从聚乙烯板背面出射的中子的内能，除了要判断中子是从聚乙烯运动到外面的World，还要通过穿越点的坐标限定穿出的位置在背面：

 G4double Energy01;
 using std::ofstream;
 ofstream outfile1("./data/中子E.txt");
 
 if((aStep->GetTrack()->GetDefinition()->GetPDGEncoding()==2112)//判断为中子
     && aStep->GetPostStepPoint()->GetPhysicalVolume() != NULL)//Step未超出计算边界
	{
        if(aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() == "PEBox"//前一点在聚乙烯中
    	&& aStep->GetPostStepPoint()->GetPhysicalVolume()->GetName()== "World"//后一点在世界
        && aStep->GetTrack()->GetPosition0().z()>4.9)//Step的位置,边界z坐标为5
            {

                Energy01=(aStep->GetTrack()->GetKineticEnergy());
	            outfile1 << Energy01 << G4endl;
            }  
	}

2.5、程序运行

        在程序目录新建data文件夹用于存储数据。在init_vis.mac中设定运行参数，可以设定运行的粒子数然后直接运行，也可以调用可视化界面进行查看和操作。

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

#运行可视化界面
#/control/execute vis.mac
#直接运行run
/run/beamOn 20000000 

在终端输入make -j进行编译，然后运行程序即可。