Geant4模拟程序基本框架

Geant4入门文档:https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/GettingStarted/gettingStarted.html
Example B1是Geant4最简单的示例,描述了粒子在生物体液、骨骼、组织环境中的运动。
Example B1文档:https://geant4-userdoc.web.cern.ch/Doxygen/examples_doc/html/ExampleB1.html
Example B1位置:/path/to/geant4-v11.0.1-install/share/Geant4-11.0.1/examples/basic/B1

本文通过梳理ExampleB1的代码逻辑,快速介绍一个基本的Geant4模拟程序的各部分组成。

ExampleB1程序清单:
| | Example B1 |
|————————|——————|
|Description |Simple application for accounting dose in a selected volume|
|Geometry | solids: box, cons, trd, simple placements with translation|
|Physics | Geant4 physics list: QBBC|
|Primary generator| Particle gun|
|Scoring | User action classes|
|Vis/GUI |显示探测器和粒子轨迹|
|Stacking |无|
|Analysis |无||

main()

main()函数结构介绍:https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/GettingStarted/mainProgram.html
示例B1的main()函数在exampleB1.cc中。
main()函数中,首先建立一个G4RunManager,需要引入头文件G4RunManagerFactory.hh

1
2
auto* runManager =
G4RunManagerFactory::CreateRunManager(G4RunManagerType::Default);

接下来在runManager中加入(1) DetectorConstruction、(2) PhysicsList和(3) ActionInitialization,其中DetectorConstruction是自定义的,写在src/DetectorConstruction.cc中,头文件是include/DetectorConstruction.hh,ActionInitialization也是自定义的,代码也在对应的位置。PhysicsList选择QBBC,头文件是Geant4的QBBC.h,包含强和电磁相互作用,比较适合描述1GeV以下薄层靶实验。这里创建了QBBC的physicsList,还设置了参数。这三部分是必要的。

1
2
3
4
5
6
7
8
9
10
11
12
// Set mandatory initialization classes
//
// Detector construction
runManager->SetUserInitialization(new DetectorConstruction());

// Physics list
G4VModularPhysicsList* physicsList = new QBBC;
physicsList->SetVerboseLevel(1);
runManager->SetUserInitialization(physicsList);

// User action initialization
runManager->SetUserInitialization(new ActionInitialization());

接下来开启可视化。如果没有这部分,通过后面UImanager运行init_vis.mac打开的窗口是没有探测器和事例显示的。

1
2
3
4
5
6
// Initialize visualization
//
G4VisManager* visManager = new G4VisExecutive;
// G4VisExecutive can take a verbosity argument - see /vis/verbose guidance.
// G4VisManager* visManager = new G4VisExecutive("Quiet");
visManager->Initialize();

设置用户接口G4UImanager,用来接受并运行用户对程序的指令。需要引入头文件G4UImanager.hh。在交互模式下开启一个session,如果没有SessionStart(),程序不会等待用户输入指令,也不会显示可视化窗口。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
// Get the pointer to the User Interface manager
G4UImanager* UImanager = G4UImanager::GetUIpointer();

// Process macro or start UI session
//
if ( ! ui ) {
// batch mode
G4String command = "/control/execute ";
G4String fileName = argv[1];
UImanager->ApplyCommand(command+fileName);
}
else {
// interactive mode
UImanager->ApplyCommand("/control/execute init_vis.mac");
ui->SessionStart();
delete ui;
}

这个地方用到了main()开头我们跳过部分的前两行,用来根据运行参数确定是交互模式(无参数)还是批处理模式(参数为宏文件),如果是交互模式,新建一个G4UIExecutive,用到了头文件G4UIExecutive.hh

1
2
3
4
// Detect interactive mode (if no arguments) and define UI session
//
G4UIExecutive* ui = nullptr;
if ( argc == 1 ) { ui = new G4UIExecutive(argc, argv); }

此外,交互模式需要一个init_vis.mac宏文件,以及这个文件调用的vis.mac文件。这两个文件结构比较简单,在多数情况下可以直接调用,也可以在此基础上修改。

最后,批处理或者交互模式结束后,删除visManager和runManager。

1
2
3
4
5
6
7
// Job termination
// Free the store: user actions, physics_list and detector_description are
// owned and deleted by the run manager, so they should not be deleted
// in the main() program !

delete visManager;
delete runManager;

探测器构造概述

探测器几何的说明:https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/GettingStarted/geometryDef.html
探测器材料的说明:https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/GettingStarted/materialDef.html

Geant4的物体层级

Geant4中,探测器的几何是由各种volume组成的。最大的volume称为World,容纳所有其他容器,容器中还可以容纳其他容器。

构造容器:首先构造solid,指定外形;然后构造logical volume,指定其材质,添加其物理属性如磁场,添加传感器;最后构造physical volume,指定其旋转、位置坐标和父logical volume。

对于World,也需要构造physical volume,但是其位置需要在全局坐标系原点,不可旋转,且父logical volume为空指针。

在可视化窗口中观察各个物体的几何,可以通过窗口左侧的Scene Tree面板。物体在可视化窗口中的颜色和透明度可以通过用户指令设置,可以放在宏文件中。

DetectorConstruction类

ExampleB1的探测器构造在src/DetectorConstruction.cc中,显示了一个最简单的探测器几何构造案例。

ExampleB1的DetectorConstruction继承自Geant4的G4VUserDetectorConstruction类,需要头文件G4VUserDetectorConstruction.hh。构建探测器,需要做的是填写其Construct()函数。

搭建探测器时,如果还没有写好ActionInitialization,可以先注释掉main()函数中的
runManager->SetUserInitialization(new ActionInitialization());
编译后可以可视化查看探测器几何。

开工之前,获取一个NIST材料数据库:

1
2
// Get nist material manager
G4NistManager* nist = G4NistManager::Instance();

提前设定下Envelope的参数,方便引用:

1
2
3
4
// Envelope parameters
//
G4double env_sizeXY = 20*cm, env_sizeZ = 30*cm;
G4Material* env_mat = nist->FindOrBuildMaterial("G4_WATER");

打开Geant4的volume覆盖错误检查,在每次构造physical volume的时候用到,防止同一层次的volume在空间上形成覆盖导致程序错误:

1
2
3
// Option to switch on/off checking of volumes overlaps
//
G4bool checkOverlaps = true;

开工!首先构造World:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
//
// World
//
G4double world_sizeXY = 1.2*env_sizeXY;
G4double world_sizeZ = 1.2*env_sizeZ;
G4Material* world_mat = nist->FindOrBuildMaterial("G4_AIR");

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
world_mat, //its material
"World"); //its name

G4VPhysicalVolume* physWorld =
new G4PVPlacement(0, //no rotation
G4ThreeVector(), //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

如上,先构造了World的solid,是一个G4Box,然后构造其logical volume,材料是NIST数据库中的G4_AIR,最后将其放置在原点。注意:在vis.mac中,World被设为不可见:

1
/vis/geometry/set/visibility World 0 false

然后构造Envelope,是一个水质的长方体,比World略小:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
//
// Envelope
//
G4Box* solidEnv =
new G4Box("Envelope", //its name
0.5*env_sizeXY, 0.5*env_sizeXY, 0.5*env_sizeZ); //its size

G4LogicalVolume* logicEnv =
new G4LogicalVolume(solidEnv, //its solid
env_mat, //its material
"Envelope"); //its name

new G4PVPlacement(0, //no rotation
G4ThreeVector(), //at (0,0,0)
logicEnv, //its logical volume
"Envelope", //its name
logicWorld, //its mother volume
false, //no boolean operation
0, //copy number
checkOverlaps); //overlaps checking

然后构造Shape 1,几何为一个圆台,材质为生物组织:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
//
// Shape 1
//
G4Material* shape1_mat = nist->FindOrBuildMaterial("G4_A-150_TISSUE");
G4ThreeVector pos1 = G4ThreeVector(0, 2*cm, -7*cm);

// Conical section shape
G4double shape1_rmina = 0.*cm, shape1_rmaxa = 2.*cm;
G4double shape1_rminb = 0.*cm, shape1_rmaxb = 4.*cm;
G4double shape1_hz = 3.*cm;
G4double shape1_phimin = 0.*deg, shape1_phimax = 360.*deg;
G4Cons* solidShape1 =
new G4Cons("Shape1",
shape1_rmina, shape1_rmaxa, shape1_rminb, shape1_rmaxb, shape1_hz,
shape1_phimin, shape1_phimax);

G4LogicalVolume* logicShape1 =
new G4LogicalVolume(solidShape1, //its solid
shape1_mat, //its material
"Shape1"); //its name

new G4PVPlacement(0, //no rotation
pos1, //at position
logicShape1, //its logical volume
"Shape1", //its name
logicEnv, //its mother volume
false, //no boolean operation
0, //copy number
checkOverlaps); //overlaps checking

然后构造Shape 2,几何为一个四棱台,材质为骨骼:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
//
// Shape 2
//
G4Material* shape2_mat = nist->FindOrBuildMaterial("G4_BONE_COMPACT_ICRU");
G4ThreeVector pos2 = G4ThreeVector(0, -1*cm, 7*cm);

// Trapezoid shape
G4double shape2_dxa = 12*cm, shape2_dxb = 12*cm;
G4double shape2_dya = 10*cm, shape2_dyb = 16*cm;
G4double shape2_dz = 6*cm;
G4Trd* solidShape2 =
new G4Trd("Shape2", //its name
0.5*shape2_dxa, 0.5*shape2_dxb,
0.5*shape2_dya, 0.5*shape2_dyb, 0.5*shape2_dz); //its size

G4LogicalVolume* logicShape2 =
new G4LogicalVolume(solidShape2, //its solid
shape2_mat, //its material
"Shape2"); //its name

new G4PVPlacement(0, //no rotation
pos2, //at position
logicShape2, //its logical volume
"Shape2", //its name
logicEnv, //its mother volume
false, //no boolean operation
0, //copy number
checkOverlaps); //overlaps checking

fScoringVolume设为Shape 2的logical volume,将来用于判定入射粒子是否在Shape 2上沉积能量:

1
2
3
// Set Shape2 as scoring volume
//
fScoringVolume = logicShape2;

如果需要旋转几何体,可以在G4PVPlacement第一个参数中传入以下rotationMatrix:
G4RotationMatrix* rotationMatrix = new G4RotationMatrix(); rotationMatrix->rotateY(90.*deg);

最后返回World的physical volume,是Construct()函数的要求:

1
2
3
4
//
//always return the physical World
//
return physWorld;

设置粒子源

首先在src/ActionInitialization.ccBuild()函数中注册产生子行为:

1
SetUserAction(new PrimaryGeneratorAction);

然后在src/PrimaryGeneratorAction.cc中设置粒子源。在PrimaryGeneratorAction()中设置粒子的一次发射的数量(一般都是1)、类型、动量、能量:

1
2
3
4
5
6
7
8
9
10
11
G4int n_particle = 1;
fParticleGun = new G4ParticleGun(n_particle);

// default particle kinematic
G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
G4String particleName;
G4ParticleDefinition* particle
= particleTable->FindParticle(particleName="gamma");
fParticleGun->SetParticleDefinition(particle);
fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));
fParticleGun->SetParticleEnergy(6.*MeV);

GeneratePrimaries(G4Event* anEvent)中设置粒子产生的位置,先动态地获取探测器几何中的logical volume “Evelope”:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
//this function is called at the begining of each event
//

// In order to avoid dependence of PrimaryGeneratorAction
// on DetectorConstruction class we get Envelope volume
// from G4LogicalVolumeStore.

G4double envSizeXY = 0;
G4double envSizeZ = 0;

if (!fEnvelopeBox)
{
G4LogicalVolume* envLV
= G4LogicalVolumeStore::GetInstance()->GetVolume("Envelope");
if ( envLV ) fEnvelopeBox = dynamic_cast<G4Box*>(envLV->GetSolid());
}

if ( fEnvelopeBox ) {
envSizeXY = fEnvelopeBox->GetXHalfLength()*2.;
envSizeZ = fEnvelopeBox->GetZHalfLength()*2.;
}
else {
G4ExceptionDescription msg;
msg << "Envelope volume of box shape not found.\n";
msg << "Perhaps you have changed geometry.\n";
msg << "The gun will be place at the center.";
G4Exception("PrimaryGeneratorAction::GeneratePrimaries()",
"MyCode0002",JustWarning,msg);
}

这里用了
G4LogicalVolume* envLV = G4LogicalVolumeStore::GetInstance()->GetVolume("Envelope");
来动态地获取探测器的logic volume,这样可以获得其下的solid的各种几何参数。获取探测器部件的位置需要获取其physical volume,可以用类似的方法:
fSourcePV = G4PhysicalVolumeStore::GetInstance()->GetVolume("Source");
source_x = fSourcePV->GetTranslation().getX();

然后根据Envelope的尺寸设置粒子产生的坐标,在XY平面内随机分布:

1
2
3
4
5
6
G4double size = 0.8;
G4double x0 = size * envSizeXY * (G4UniformRand()-0.5);
G4double y0 = size * envSizeXY * (G4UniformRand()-0.5);
G4double z0 = -0.5 * envSizeZ;

fParticleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));

最后调用GeneratoPrimaryVertex()一次产生n_particle个完全相同的粒子(注意坐标都一样的。但是在后面传播中的随机过程会不一样):

1
fParticleGun->GeneratePrimaryVertex(anEvent);

设置完探测器几何和粒子属性,即使没有写好RunAction、EventAction和SteppingAction,也可以注释掉ActionInitialization::Build()中的对应的注册语句,进行编译并在可视化窗口的Session栏中尝试/run/BeamOn 100来观察粒子和探测器几何的交互效果。注意之前在main()函数中是注册过PhysicsList的。