Geant4--G4ParticleGun定义射线源的发射能谱

“ G4ParticleGun作为Geant4模拟中常用的粒子产生器，本文代码讲解怎样模拟发射符合自定义能量分布的粒子。”

一、G4ParticleGun/G4GeneralParticleSource简介

二、场景一：粒子源的能量为多个分离值

三、场景二：粒子源的能量为多段能谱

—

G4ParticleGun/G4GeneralParticleSource简介

在Geant4模拟中，粒子源的设置与实现是关键部分，用以定义初级粒子的种类、能量、动量方向、产生位置、极化等特征。用户自己所定义的初级粒子的产生方法即类函数，基于G4VuserPrimaryGeneratorAction（G4VPrimaryGenerator）这个G4自带的粒子产生函数派生而来，通过generatePrimaryVertex(G4Event*)函数发射出想要的粒子。G4ParticleGun、G4GeneralParticleSource和G4HEPEvtInterface是Geant4提供的三种用户可直接调用的粒子发生器类，通常只用前两个即gun和gps。Gun和gps两者的应用区别主要在于：gun通常在*.cc源文件中通过函数方法自定义源的特征，而gps则在*.mac宏文件中通过宏命令/pgs/*来定义源。gps的宏命令更加丰富，比如对于复杂分布形状和能谱的源，gps均有现成的宏命令结构可以调用。图1给出两者的宏命令直观比较，参见UserGuide-2.7.3小节。

图1 gun和gps的宏命令比较

能谱是大多数模拟用户所关心的问题，下面将代码讲解G4ParticleGun如何通过随机种G4UniformRand()函数，实现粒子的抽样能量服从自定义的能量分布。

—

场景一：粒子源的能量为多个分离值

当需要模拟的射线源能量为多个不同的数值，比如模拟多个单能射线源时，下述方法适用。模拟场景为粒子源的能量为2、7、10MeV，丰度占比各为0.5、0.25、0.25。

//oooooooo From G4UniformRand() toEnergy-spectrum ooooooooooooo//
G4double outene; //定义变量值来获取最终抽取的能量
G4double randseed=G4UniformRand();
//数值变化为0~1的随机种
G4double enespec[]={2,7,10};
//定义三个独立的能量
G4double pdfspec[]={0.5,0.25,0.25};
//定义三种能量对应的丰度值
assert(sizeof(enespec) == sizeof(pdfspec));
//判断两个数组元素个数相同
const G4int pdfnEntries =sizeof(pdfspec)/sizeof(G4double);
//获取数组的元素个数
const G4int pdfintegralnEntries=pdfnEntries +1;
//定义概率密度的积分点
G4doublepdfintegral[pdfintegralnEntries]={0};
for (int i=0;i<pdfnEntries;i++)
{
pdfintegral[i+1]=pdfintegral[i]+pdfspec[i];
//计算每个能量对应的概率积分值
}

//G4double pdfintegral[]={0,0.5,0.75,1};
for (int i=0;i<pdfnEntries;i++)
{
       if(pdfintegral[i]<randseed&&pdfintegral[i+1]>randseed)
       {
                outene=enespec[i]*MeV;
  //判断当前随机种数值位于哪个能量对应的概率积分段，获取对应的能量
       }
}

‍

图2 射线源能量分布为多个单能数值

—

场景二：粒子源的能量为多段能谱

当模拟的射线源能量需要服从某种概率密度函数PDF时，可以认为是连续的能谱填入到了一定数量的bin当中，每个bin的y大小对应当前能量段的概率密度。下述模拟场景对应的能谱如图3所示。

//user-defined energy hist，连续分段能谱，概率密度函数
G4double enespec[]={1,2,7,10};
G4double pdfspec[]={0.5,0.25,0.25};
const G4int pdfnEntries =sizeof(pdfspec)/sizeof(G4double);
const G4int pdfintegralnEntries=pdfnEntries +1;
G4double pdfintegral[pdfintegralnEntries]={0};
for (int i=0;i<pdfnEntries;i++)
{
pdfintegral[i+1]=pdfintegral[i]+pdfspec[i];
}

assert(sizeof(enespec) ==sizeof(pdfintegral));
//G4double pdfintegral[]={0,0.5,0.75,1};
for (int i=0;i<pdfnEntries;i++)
{
       if(pdfintegral[i]<randseed&&pdfintegral[i+1]>randseed)
       {   
               outene=(enespec[i]+(enespec[i+1]-enespec[i])*(randseed-pdfintegral[i])/pdfspec[i])*MeV;                
       }
}

图3 服从自定义概率密度函数分布的能谱

对射线源的能量做统计检验：

/*
fstream datafile1;
datafile1.open("separateene.txt",ios::app|ios::out);
//datafile1.open("segmentene.txt",ios::ate|ios::out);
datafile1<<outene<<G4endl;
*/
fParticleGun->SetParticleEnergy(outene);
fParticleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));
fParticleGun->GeneratePrimaryVertex(anEvent);

喜欢的话，分享一下吧~^o^~

来留言呀~