B1实例代码学习记录
三、粒子源定义(PrimaryGeneratorAction)
在这个 文件中我们需要完成以下内容的定义:
- 粒子的类型
- 粒子的能量
- 粒子发射的方向
- 粒子的产生的位置
接下来看B1是如何实现以上信息的定义的,仍然是首先看一下头文件:
class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
{
public:
PrimaryGeneratorAction();
~PrimaryGeneratorAction() override;
// method from the base class
void GeneratePrimaries(G4Event*) override;
// method to access particle gun
const G4ParticleGun* GetParticleGun() const { return fParticleGun; }
private:
G4ParticleGun* fParticleGun = nullptr; // pointer a to G4 gun class
G4Box* fEnvelopeBox = nullptr;
};
从头文件中可以看到,自定义的PrimaryGeneratorAction也是继承于以G4开头的G4VUserPrimaryGeneratorAction,除了构造函数与析构函数,其中实现粒子源定义的函数为GeneratePrimaries,剩下的GetParticleGun与粒子源定义无关,同时包含两个私有成员变量。
1.粒子的类型的定义
G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
G4String particleName;
G4ParticleDefinition* particle
= particleTable->FindParticle(particleName="gamma");
fParticleGun->SetParticleDefinition(particle);
粒子源类型的定义基本可以参照上面的固定格式,具体是如何实现的我认为不是初学的应该去关心的,应该先学会使用这个定义。具体除了particleName="gamma"还有哪些,我之后查一下再在这里列出来。
2.粒子的能量的定义
fParticleGun->SetParticleEnergy(6.*MeV);
利用fParticleGun的SetParticleEnergy方法即可实现对粒子源能量的定义,这里需要额外说的是6.*MeV,其中MeV也是一个double的数值1。
3.粒子发射方向的定义
fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));
这个方法可以定义粒子出射的方向,其中G4ThreeVector(0.,0.,1.)的三个坐标分别为X、Y、Z。需要注意的是SetParticleMomentumDirection这个方法的参数为G4ThreeVector。
上面这三个信息的定义都在PrimaryGeneratorAction的构造函数中定义了,但是我尝试了并不是一定要在构造函数中定义,也可以在接下来的GeneratePrimaries方法中定义。
4.粒子的产生的位置的定义
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);
}
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));
这里有很多方法值得我们参考,但是就像我说的作为初学者我们先学会使用,先知其然等精通之后再去究其所以然,这是我认为该有的学习方法,下面我们一起来看一下B1是如何实现的。
G4LogicalVolume* envLV= G4LogicalVolumeStore::GetInstance()->GetVolume("Envelope");
首先通过这行代码获取了图中蓝色体积逻辑体(LogicalVolume)的指针,并将其赋值给envLV。
if ( envLV ) fEnvelopeBox = dynamic_cast<G4Box*>(envLV->GetSolid());
这里稍微解释一下dynamic_cast<G4Box*>这是C++中基类指针向下(派生类)转换的命令,GetSolid方法的返回值为G4tgrSolid* GetSolid() const { return theSolid; }。
这条命令将Envelope的几何体的指针赋值给G4Box*类型的变量fEnvelopeBox
envSizeXY = fEnvelopeBox->GetXHalfLength()*2.;
envSizeZ = fEnvelopeBox->GetZHalfLength()*2.;
G4double size = 0.8;
G4double x0 = size * envSizeXY * (G4UniformRand()-0.5);
G4double y0 = size * envSizeXY * (G4UniformRand()-0.5);
G4double z0 = -0.5 * envSizeZ;
先看一下x、y的坐标获取,envSizeXY = fEnvelopeBox->GetXHalfLength()*2.;这条命令得到了Envelope在x、y轴上面的边长,然后G4double x0 = size * envSizeXY * (G4UniformRand()-0.5);这句话中有一个随机数G4UniformRand()它的范围是0-1,我们看到在几何体定义的时候用的是0.5*x是因为G4将几何体关于对称轴对称放置,那么(G4UniformRand()-0.5)即产生了(-0.5~0.5)的随机数,其结果是在x的范围内随机的产生一些点,另外的y也是一样。z0的定义 G4double z0 = -0.5 * envSizeZ;表明产生的点在其负半轴的那个面上,综合上述分析并结合G4示意图可以看出,这些代码在蓝色体积左侧表面0.8x*0.8y的范围内随机的产生一些点。
fParticleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));
最后将随机产生的点作为粒子源产生的位置。
经过以上方法即可自定义一个粒子源,并给出粒子源的类型、能量、出射位置与方向等信息。
5.扩展(自主学习或未来答疑)
- 这里使用的是固定能量的粒子,如果我们是根据能谱进行发射该如何更改代码?
- 除了ParticleGun这种类型的粒子源还有哪些粒子源?除了gamma粒子,还能发射哪些粒子?
- 前面有讲过PrimaryGeneratorAction这个类主要功能是定义上面提到的粒子源的信息,并不能实现粒子的发射,真正实现粒子发射的类是G4ParticleGun,也就是在程序中经常使用得到变量fParticleGun,所有定义的信息都要传送给fParticleGun,并最终通过下面这行代码发射初级粒子
fParticleGun->GeneratePrimaryVertex(anEvent);
- 所有的带有Action的类都称为动作类,都是给我们提供一个数据接口,方便我们在合适的位置对我们需要的信息进行提取。
到这一步基本上就完成了我们需要自己搭建的模型,当然这仅仅是在B1中。
四、物理过程的定义(PhysicsList)
auto physicsList = new QBBC;
PhysicsList在这里使用的是系统定义的QBBC,它包含了一个较全的物理过程,如果你在写自己的程序可以去自己研究一下如何去定义自己的PhysicsList,或者去看视频学习。
以上内容均为自学的学习笔记,如有错误希望帮忙指正。