Geant4 从入门到精通

目录

一、通过 Geant4 的简单例子入门

1 搭建 main() 函数框架

2 创建探测器

2.1 探测器的形状

2.2 探测器的材料

2.3 探测器的位置

3 创建粒子

3.1 具体实现

3.2 设置截断

4 指定物理过程

---

常用函数：

initialize()

建立探测器结构，创建物理过程，计算截面并且建立 run。

ApplyCommand()

发命令给应用程序让其打印信息，例如打印运行 run 、事件 event 和粒子跟踪 tracking 的信息。

BeamOn()

启动事件的一个 run。

G4cout()

由 Geant4 定义的 iostream 对象，被 G4UImanager 处理，输出需要的信息。

G4cerr()

由 Geant4 定义的 iostream 对象，被 G4UImanager 处理，输出错误信息。

一、通过 Geant4 的简单例子入门

1 搭建 main() 函数框架

G4中 main() 函数必须自己构造，以下为一个最简单的例子，之后可根据具体情况开发。 

#include "G4RunManager.hh"
#include "G4UImanager.hh"

#include "ExN01DetectorConstruction.hh"
#include "ExN01PhysicsList.hh"
#include "ExN01PrimaryGeneratorAction.hh"

int main()
{
    // 构造缺省 run manager
    G4RunManager* runManager = new G4RunManager;
    // 设置必须的初始化类
    runManager->SetUserInitialization(new ExN01DetectorConstruction);
    runManager->SetUserInitialization(new ExN01PhysicsList);
    // 设置必须的用户行为类
    runManager->SetUserAction(new ExN01PrimaryGeneratorAction);
    // 初始化 G4 内核
    runManager->initialize();
    // 获取指向 UI manager 的指针并设置 verbosities
    G4UImanager* UI = G4UImanager::GetUIpointer();
    UI->ApplyCommand("/run/verbose 1");
    UI->ApplyCommand("/event/verbose 1");
    UI->ApplyCommand("/tracking/verbose 1");
    // 开始运行三个按顺序处理的事件，启动一个 run
    int numberOfEvent = 3;
    runManager->BeamOn(numberOfEvent);
    // 任务结束
    delete runManager;
    return 0;
}

运行管理类：

G4RunManager 类是 G4 内核中唯一的一个运行管理类。需要传递所有必须的信息给运行控制进程 G4RunManager，包括：探测器将如何构建、被模拟的所有粒子和所有物理过程、在一个事件中的初级粒子将如何产生和其他模拟必须的信息。

其它管理类在运行管理类创建的时候创建。其中一个是用户接口管理类 G4UImanager，在main() 函数中，必须获取指向用户接口管理进程的指针。

用户初始化类：

ExN01DetectorConstruction 由 G4VUserDetectorConstruction 类派生 

可定义：探测器的几何形状，在探测器中使用的材料，探测器的敏感区域和这些敏感区域的读出方式。

ExN01PhysicsList 由 G4VUserPhysicsList 类派生

可定义：在模拟中将被使用的粒子，这些粒子的截断范围和所有将被模拟的物理过程。

用户行为类：

ExN01PrimaryGeneratorAction 由 G4VUserPrimaryGeneratorAction 派生

可定义：初级事件的初始状态。

G4 提供了 5 个可选用户行为类：

G4UserRunAction，G4UserEventAction，G4UserStackingAction，G4UserTrackingAction ，G4UserSteppingAction 

以下为一个使用交互式终端和可视化的 main() 函数例子。

#include "G4RunManager.hh"
#include "G4UImanager.hh"
#include "G4UIterminal.hh"

#include "N02VisManager.hh"
#include "N02DetectorConstruction.hh"
#include "N02PhysicsList.hh"
#include "N02PrimaryGeneratorAction.hh"
#include "N02RunAction.hh"
#include "N02EventAction.hh"
#include "N02SteppingAction.hh"

#include "g4templates.hh"

int main(int argc,char** argv)
{
    // 构造缺省的 run manager
    G4RunManager * runManager = new G4RunManager;
    // 设置必须的初始化类
    N02DetectorConstruction* detector = new N02DetectorConstruction;
    runManager->SetUserInitialization(detector);
    runManager->SetUserInitialization(new N02PhysicsList);
    // 可视化 manager
    G4VisManager* visManager = new N02VisManager;
    visManager->initialize();
    // 设置用户行为类
    runManager->SetUserAction(new N02PrimaryGeneratorAction(detector));
    runManager->SetUserAction(new N02RunAction);
    runManager->SetUserAction(new N02EventAction);
    runManager->SetUserAction(new N02SteppingAction);
    // 获取指向 UI manager 的指针
    G4UImanager* UI = G4UImanager::GetUIpointer();
    if(argc==1)
    // 定义交互式(G)UI 终端
    {
        G4UIsession * session = new G4UIterminal;
        UI->ApplyCommand("/control/execute prerun.g4mac");
        session->sessionStart();
        delete session;
    }
    else
    // 批处理模式
    {
        G4String command = "/control/execute ";
        G4String fileName = argv[1];
        UI->ApplyCommand(command+fileName);
    }
    // 任务结束
    delete visManager;
    delete runManager;
    return 0;
}

intercoms类：

Geant4 提供了一大类叫 intercoms 的类，G4UImanager 是这个大类中的管理类，可以通过调用这些类对象的 set 方法来使用这些类的功能。

2 创建探测器

在G4中最大的的几何体叫世界(World)，其它的几何体，都被创建并放置在世界(World)中。每个几何体都通过描述它的形状（由实体定义，如尺寸、形状）和物理特性（由逻辑几何体定义，如材料、磁场、是否包含探测器的敏感单元）来创建，然后放置到另一个几何体中。当一个几何体被放置到另一个中时，我们叫前者为子几何体，后者为母几何体，并用相对于母体坐标系的坐标指定子体放置的位置。

2.1 探测器的形状

想构建一个探测器，首先要创建一个实体。

G4double expHall_x = 3.0*m;
G4double expHall_y = 1.0*m;
G4double expHall_z = 1.0*m;
G4Box* experimentalHall_box = new G4Box("expHall_box",expHall_x,expHall_y,expHall_z);

G4double innerRadiusOfTheTube = 0.*cm;
G4double outerRadiusOfTheTube = 60.*cm;
G4double hightOfTheTube = 50.*cm;
G4double startAngleOfTheTube = 0.*deg;
G4double spanningAngleOfTheTube = 360.*deg;
G4Tubs* tracker_tube = new G4Tubs("tracker_tube",innerRadiusOfTheTube,outerRadiusOfTheTube,hightOfTheTube,startAngleOfTheTube,spanningAngleOfTheTube);

这段代码分别创建了一个名为 expHall_box 的盒子，大小为沿 X 轴方向从-3m 到 3m，沿 Y 轴方向从-1m 到 1m，沿 Z 轴方向从-1m 到 1m；和一个名为 tracker_tube 的圆柱体，半径为 60cm，高为 50cm。

2.2 探测器的材料

创建好实体后，通过定义实体的物理特性，便可以得到逻辑几何体。

首先定义几何体的材料，这里我们可以利用刚刚创建好的实体。

G4LogicalVolume* experimentalHall_log
= new G4LogicalVolume(experimentalHall_box,Ar,"expHall_log");

G4LogicalVolume* tracker_log
= new G4LogicalVolume(tracker_tube,Al,"tracker_log");

这段代码用氩气填充了盒子 expHall_box，生成逻辑体 expHall_log；用铝填充了圆柱体 tracker_tube，生成逻辑体 tracker_log。

实际上G4中材料的定义是非常精细的，为此设置了几个专门的类，每一个类都有一个表作为一个静态数据成员，用于跟踪这些类各自的实例。

G4Element 类 描述原子属性：原子序数，核子数 ，原子质量，壳层能量，原子截面等等。

G4Material类 描述物质的宏观属性：密度，状态，温度，压强，辐射长度，平均自由程，单位长度能损等等。

例如我们可以通过指定密度、摩尔质量和原子序数，创建名为液氩的材料，并将其重新填充给盒子 expHall_box。

G4double density = 1.390*g/cm3;
G4double a = 39.95*g/mole;
G4Material* lAr = new G4Material(name="liquidArgon", z=18., a, density);
G4LogicalVolume* experimentalHall_log
= new G4LogicalVolume(experimentalHall_box,lAr,"expHall_log");

我们也可以通过指定分子中各种原子的种类和数目，来创建水这种材料。

a = 1.01*g/mole;
G4Element* elH = new G4Element(name="Hydrogen",symbol="H" , z= 1., a);
a = 16.00*g/mole;
G4Element* elO = new G4Element(name="Oxygen" ,symbol="O" , z= 8., a);
density = 1.000*g/cm3;
G4Material* H2O = new G4Material(name="Water",density,ncomponents=2);
H2O->AddElement(elH, natoms=2);
H2O->AddElement(elO, natoms=1);

此外还可以通过指定各种成分的质量百分数，创建像空气这样的混合物。

a = 14.01*g/mole;
G4Element* elN = new G4Element(name="Nitrogen",symbol="N" , z= 7., a);
a = 16.00*g/mole;
G4Element* elO = new G4Element(name="Oxygen" ,symbol="O" , z= 8., a);
density = 1.290*mg/cm3;
G4Material* Air = new G4Material(name="Air ",density,ncomponents=2);
Air->AddElement(elN, fractionmass=70*perCent);
Air->AddElement(elO, fractionmass=30*perCent);

最后我们可以通过打印来查看已经创建的材料信息。

G4cout << H2O; \\ 打印水的材料信息
G4cout << *(G4Material::GetMaterialTable()); \\ 打印已创建的材料列表

2.3 探测器的位置

在创建好逻辑几何体后，我们必须定义这个几何体的中心位于另一个几何体的什么位置，以及如何旋转。需要说明的是，World作为包含所有其他几何体的最大的几何体，必须使用空母体指针通过 G4PVPlacement 来创建，它必须是不旋转的，必须被放置于全局坐标原点。这里我们用之前创建的 experimental hall 作为世界。

G4VPhysicalVolume* experimentalHall_phys
= new G4PVPlacement(0, //不旋转
G4ThreeVector(0.,0.,0.), //子体在母体中的坐标
experimentalHall_log, //子逻辑体指针
"expHall", //物理体名
0, //母逻辑体指针
false, //无布尔操作
0); //物理体的拷贝序号

一旦这些都确定了，那么就可以创建一个物理体了，物理体是一个已经确定位置的逻辑几何体实例。

G4double trackerPos_x = -1.0*meter;
G4double trackerPos_y = 0.0*meter;
G4double trackerPos_z = 0.0*meter;
G4VPhysicalVolume* tracker_phys
= new G4PVPlacement(0, // 不旋转
G4ThreeVector(trackerPos_x,trackerPos_y,trackerPos_z), // 子体在母体中的坐标
tracker_log, // 子逻辑体指针
"tracker", // 物理体名
experimentalHall_log, // 母逻辑体指针
false, // 无布尔操作
0); // 物理体的拷贝序号

这段代码将逻辑体 tracker_log 放置在母体 experimentalHall_log 的坐标原点，然后沿着X轴负方向移动1m，不旋转。生成的物理体叫 tracker ，它的拷贝序号为 0。

3 创建粒子

Geant4 为用户提供了各种类型的粒子，每个粒子都由各自的类来描述，设定好了名字，质量，电荷，自旋等。这些类都是由 G4ParticleDefinition 派生的，每一类粒子都被定义在 geant4/source/particles 的子目录中，分别有一个相应的库，属性为只读，用户无法修改。这些粒子主要分为 6 类：轻子、介子、重子、玻色子、短寿命粒子、离子，总数超过100个。

G4ParticleTable 类是一个粒子字典。它提供了各种实用方法，如：

FindParticle(G4String name):     //通过名字查找粒子

FindParticle(G4int PDGencoding):      //通过 PDG 编码查找粒子

3.1 具体实现

ConstructParticle() 是一个纯虚方法，所有在模拟中需要的粒子的静态成员函数应在这个方法中被调用，这保证了那些粒子的对象将被建立。注意，用户必须定义初级粒子和其它所有可能出现的次级粒子。

 例如，假定用户需要一个质子和一个 geantino，geantino 是一个虚粒子，它不与介质发生作用。ConstructParticle() 方法的实现如下： 

void ExN01PhysicsList::ConstructParticle()
{
    G4Proton::ProtonDefinition();
    G4Geantino::GeantinoDefinition();
}

由于有大量预定义的粒子，用这个方法列出所有的粒子非常麻烦。Geant4 有 6 个粒子大类，对应有 6 个实用类，它们可用于简化粒子构造的过程：

G4BosonConstructor、G4LeptonConstructor、G4MesonConstructor、G4BarionConstructor 、G4IonConstructor、G4ShortlivedConstructor

例如如果需要构造所有轻子，可以这样写：

void ExN05PhysicsList::ConstructLeptons()
{
    // Construct all leptons
    G4LeptonConstructor pConstructor;
    pConstructor.ConstructParticle();
}

3.2 设置截断

为了避免红外发散，一些电磁过程要求设定一个低阈，在阈值以下，将不产生任何次级粒子。因此gamma、电子、正电子要求用户设定阈值。阈值应定义为距离或者截断范围，它将自动转换为对应于不同介质的截断能量。这个阈值应在用户初始化部分用 G4VUserPhysicsList 的SetCuts() 方法定义。这里要注意顺序，粒子、材料、物理过程的构造应在引用 SetCuts()之前。

“唯一截断范围值”的思想是 Geant4 的一个重要特性，对于多数应用来说，用户只要确定一个截断范围值，这个值将以同样的方式用于 gamma，电子，正电子。在这种情况下，可以使用 SetCutsWithDefault() 方法，由基类 G4VuserPhysicsList 提供，它有一个 defaultCutValue 成员作为缺省截断范围值，SetCutsWithDefault() 使用这个值。defaultCutValue 在缺省情况下设置为 1.0 mm，你可以在你的 physics list 类 中设置新的缺省值。

ExN04PhysicsList::ExN04PhysicsList(): G4VUserPhysicsList()
{
    // 默认缺省截断范围1.0mm，此处修改为2.0mm
    defaultCutValue = 2.0*mm;
}

void ExN04PhysicsList::SetCuts()
{
    // G4VUserPhysicsList::SetCutsWithDefault() 方法为所有粒子设置统一的截断值
    SetCutsWithDefault();
}

我们也可以对 gamma、电子和正电子设定不同的截断范围值，并且对不同的几何区域设置不同的截断范围值。这里利用 SetCuts() 方法。

void ExN03PhysicsList::SetCuts()
{
    // 先为gamma设置截断值，再为电子设置，最后为正电子设置，因为一些正负电子的过程需要gamma的截断值
    SetCutValue(cutForGamma, "gamma");
    SetCutValue(cutForElectron, "e-");
    SetCutValue(cutForElectron, "e+");
}

4 指定物理过程

物理过程描述粒子如何与物质相互作用，Geant4 提供了 7 个大类的这些过程：电磁相互作用、强相互作用、输运、衰变、可见光、轻重子、参数化相互作用。

想要指定物理过程，用户必须构造一个从 G4VUserPhysicsList 类派生的类，并且实现它的纯虚方法 ConstructProcess()。例如，如果在模拟过程中只使用了 G4Geantino 粒子，那么只要注册输运过程就可以了。

void ExN01PhysicsList::ConstructProcess()
{
    // 注册输运过程
    AddTransportation();
}

 AddTransportation() 方法由G4VUserPhysicsList 类中提供，用于为所有粒子类注册G4Transportation 类。G4Transportation 类描述粒子在时空中的运动，它是实现粒子跟踪所必需的类。 

下面是为光子gamma注册电磁作用过程的示例。

void MyPhysicsList::ConstructProcess()
{
    // 注册输运过程
    AddTransportation();
    // 电磁过程
    ConstructEM();
}
void MyPhysicsList::ConstructEM()
{
    // 获得gamma的过程管理器指针
    G4ParticleDefinition* particle = G4Gamma::GammaDefinition();
    G4ProcessManager* pmanager = particle->GetProcessManager();
    // 构建gamma的物理过程
    G4PhotoElectricEffect * thePhotoElectricEffect = new G4PhotoElectricEffect();
    G4ComptonScattering * theComptonScattering = new G4ComptonScattering();
    G4GammaConversion* theGammaConversion = new G4GammaConversion();
    // 将物理过程注册到gamma的过程管理器中
    pmanager->AddDiscreteProcess(thePhotoElectricEffect);
    pmanager->AddDiscreteProcess(theComptonScattering);
    pmanager->AddDiscreteProcess(theGammaConversion);
}

G4ProcessManager 可以在一个 run 过程中，使用 ActivateProcess() 和 InActivateProcess()方法来打开或关闭一些物理过程。这些方法只能在物理过程注册之后使用，所以它们不能在用户预初始化部分使用。

5 产生初级事件

G4VuserPrimaryGeneratorAction 是一个用户必需的类，用户需要从它派生自己的类。在用户的 G4VUserPrimaryGeneratorAction 类的构造函数中，用户应该初始化初级粒子发生器。

G4VUserPrimaryGeneratorAction 类，有一个纯虚方法 generatePrimaries()，这个方法在每一个事件的开始被调用。在这个方法中，用户必须调用具体 G4VPrimaryGenerator 类的 generatePrimaryVertex() 方法。用户可以调用一个或多个粒子发生器，也可以多次调用同一个粒子发生器。混合使用几个发生器可以产生一个非常复杂的初级事件。

在用户派生的类中，必须指定如何产生一个初级事件。实际上，初级粒子的产生是由G4VPrimaryGenerator 类完成的，用户的 G4VUserPrimaryGeneratorAction 构造类只是描述了初级粒子的产生方式。Geant4 提供了两个 G4VPrimaryGenerator 类，一个是 G4ParticleGun 类，另一个是 G4HEPEvtInterface 类，这里我们先使用前者。

G4ParticleGun 是 Geant4 提供的一个粒子发生器类。这个类用一个给定的动量和位置产生初级粒子。它不提供任何的随机模式。G4ParticleGun 的构造函数使用一个整型参数用于指定产生多少个完全相同的粒子。通常，用户要求产生一个具有随机能量，动量，和位置的初级粒子。这种随机特性可以通过调用 G4ParticleGun 类提供的 set 方法来获得。引用这些方法应在用户的 G4VUserPrimaryGeneratorAction 类的 generatePrimaries()方法中实现，还应在调用 G4ParticleGun 的 generatePrimaryVertex() 方法之前。

以下为一个使用 G4ParticleGun 类的 G4VUserPrimaryGeneratorAction 例子。

#ifndef ExN01PrimaryGeneratorAction_h
#define ExN01PrimaryGeneratorAction_h 1

#include "G4VUserPrimaryGeneratorAction.hh"

class G4ParticleGun;
class G4Event;

class ExN01PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
{
    public:
    ExN01PrimaryGeneratorAction();
    ~ExN01PrimaryGeneratorAction();

    public:
    void generatePrimaries(G4Event* anEvent);

    private:
    G4ParticleGun* particleGun;
};

#endif

#include "ExN01PrimaryGeneratorAction.hh"
#include "G4Event.hh"
#include "G4ParticleGun.hh"
#include "G4ThreeVector.hh"
#include "G4Geantino.hh"
#include "globals.hh"

ExN01PrimaryGeneratorAction::ExN01PrimaryGeneratorAction()
{
    G4int n_particle = 1;
    particleGun = new G4ParticleGun(n_particle);
    particleGun->SetParticleDefinition(G4Geantino::GeantinoDefinition());
    particleGun->SetParticleEnergy(1.0*GeV);
    particleGun->SetParticlePosition(G4ThreeVector(-2.0*m,0.0*m,0.0*m));
}

ExN01PrimaryGeneratorAction::~ExN01PrimaryGeneratorAction()
{
    delete particleGun;
}

void ExN01PrimaryGeneratorAction::generatePrimaries(G4Event* anEvent)
{
    G4int i = anEvent->get_eventID() % 3;
    switch(i)
    {
        case 0:
            particleGun->SetParticleMomentumDirection(G4ThreeVector(1.0,0.0,0.0));
            break;
        case 1:
            particleGun->SetParticleMomentumDirection(G4ThreeVector(1.0,0.1,0.0));
            break;
        case 2:
            particleGun->SetParticleMomentumDirection(G4ThreeVector(1.0,0.0,0.1));
            break;
    }
    
    particleGun->generatePrimaryVertex(anEvent);
}

---

持续更新中，收藏不迷路