Geant4 PrimaryGeneratorAction 中的粒子源定义琐记

一般定义简单源

整体结构

PrimaryGeneratorAction::PrimaryGeneratorAction()
: G4VUserPrimaryGeneratorAction(),fParticleGun(0)
{
  G4int n_particle = 1;
  fParticleGun  = new G4ParticleGun(n_particle);

  // default particle kinematic
  
  G4ParticleDefinition* particle = G4ParticleTable::GetParticleTable()->FindParticle("proton");//gamma，neutron等
  fParticleGun->SetParticleDefinition(particle);
  fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1));
  fParticleGun->SetParticleEnergy(7.*MeV);
  fParticleGun->SetParticlePosition(G4ThreeVector(0.*cm,0.*cm,-1.*cm));
}

PrimaryGeneratorAction::~PrimaryGeneratorAction()
{
  delete fParticleGun;
}

• 此部分为PrimaryGenerationAction的构造函数和析构函数
• PrimaryGeneratorAction::PrimaryGeneratorAction(): G4VUserPrimaryGeneratorAction(),fParticleGun(0) 中需注意fParticleGun(0)部分
• 在.hh中声明 G4ParticleGun* GetParticleGun() {return fParticleGun;};，G4ParticleGun* fParticleGun; 。

GeneratePrimaries 函数

void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
  //this function is called at the begining of ecah event

//更详细的粒子源定义，如果不需要可以不填写

  fParticleGun->GeneratePrimaryVertex(anEvent);
}

• 此部分需注意在.hh中声明 virtual void GeneratePrimaries(G4Event*);
• 以下为可以在该函数中添加的源相关描述

粒子类型

// code to get pre-defined particle
G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
G4ParticleDefinition* fGamma = particleTable->FindParticle("gamma"),
	fElec = particleTable->FindParticle("e-"),
	fNeut = particleTable->FindParticle("neutron"),
	fProt = particleTable->FindParticle("proton");

粒子源其他的具体描述

随机数抽样

#include"Randomize.hh"
// THE MOST IMPORTANT ONE
// uniform random number in [0,1]
G4double u = G4UniformRand();

// other distributions
G4double uGauss = G4RandGauss::shoot(mu,sigma);
G4double uExpo = G4RandExponential::shoot(mean);
G4double uFlat = G4RandFlat::shoot(a, b);

• 符合各种标准分布的随机数抽样描述在Randomize.hh

示例

//圆形面源均匀分布。无分散，即方向沿z方向
  G4double size = 1 * mm;//源半径
  G4double r = size * (G4UniformRand());
  G4double theta = 2.*3.14159265*G4UniformRand();
  G4double x0 = r * cos(theta);
  G4double y0 = r * sin(theta);
  G4double z0 = -10.*mm;
  
  fParticleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));

//方形面源均匀分布。无分散，即方向沿z方向
  G4double size = 0.01 * mm;//源的边长
  G4double x0 = size * (G4UniformRand()-0.5);//因为size是边长，所以使用抽样函数时-0.5！
  G4double y0 = size * (G4UniformRand()-0.5);
  G4double z0 = -4000.*mm;
  
  fParticleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));

//定义发散角
  G4ThreeVector pDir(0, 0, 1);//向z正方向，定义一个单位长度的向量
  pDir.setTheta(0.5*degree*G4UniformRand());//球坐标的theta：[0,180]，也可以使用rad单位
  pDir.setPhi(360*degree*G4UniformRand());//球坐标的phi：[0,360]，也可以使用rad单位
  
  fParticleGun->SetParticleMomentumDirection(pDir);

//各向同性源
G4double cosTheta = 2*G4UniformRand() - 1., 
phi = 2.*3.14159265*G4UniformRand()
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));

//粒子位置高斯分布
G4double muP = 0.*mm;//中心值
G4double sigmaP = 2.*mm;//标准差
G4double x0 =  G4RandGauss::shoot(muP,sigmaP) ;//高斯抽样
G4double y0 = G4RandGauss::shoot(muP,sigmaP) ;
G4double z0 = -10.*mm;

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

//能量--高斯分布
G4double muE = 7.*MeV;//中心值
G4double sigmaE = 0.01395 *MeV;//标准差
G4double eng = G4RandGauss::shoot(muE,sigmaE);//高斯抽样

fParticleGun->SetParticleEnergy(eng);