粒子物理蒙特卡罗模拟库Geant4开发之一个实例

项目内容

9 MeV光子束在均匀介质中的能量沉积的模拟;介质是一个圆柱体，半径为5 cm，充满C6F6。粒子列表只包含γ、电子、正电子。物理事件列表包含“标准”电磁过程(光电效应、Compton散射、伽玛转换等等)。

架构逻辑是标准的Geant4程序架构，详见上一篇博客的“最基础的主函数”。

---

源码内涵

• ActionInitialization.cc : 设置RunAction、PrimaryGeneratorAction、EventAction等等。
• PhysicsList.cc : 粒子的构造、粒子作用过程的注册汇总、阶段距离设定
• SteppingAction.cc :G4VUserSteppingAction的实例。在每个 step后定义具体事件。本例中它统计了每一步的总能量沉积。
• DetectorConstruction.cc :定义了探测器材料尺寸以及“world”场景。
• PrimaryGeneratorAction.cc : 用户描述初级事件的初始状态。本例中定义了粒子枪以及角度、方向。
• SteppingVerbose.cc :输出信息之类的
• EventAction.cc :G4UserEventAction的实例。有beginOfEventAction()和EndOfEventAction()两个关键方法：一般为一
  个特殊的事件初始化和后处理。本例中是在初始化总能量沉积和最后统计。
• RunAction.cc :G4UserRunAction的实例。有BeginOfRunAction()和EndOfRunAction()两个虚方法。分别在 beamOn() 方法开始和结束后调用。还有一个RunAction()方法。本例用来初始化G4AnalysisManager及数据写出到文件。

---

探测器及场景构建

G4VPhysicalVolume* DetectorConstruction::Construct()
{
  //
  // define a material from its elements.   case 1: chemical molecule
  // 
  G4double a, z;
  G4double density;  
  G4int ncomponents, natoms;

  G4Element* C = new G4Element("Carbon"  ,"C" , z= 6., a= 12.01*g/mole);
  G4Element* F = new G4Element("Fluorine","N" , z= 9., a= 18.99*g/mole);

  G4Material* C6F6 = 
  new G4Material("FluorCarbonate", density= 1.61*g/cm3, ncomponents=2);
  C6F6->AddElement(C, natoms=6);
  C6F6->AddElement(F, natoms=6);

  G4cout << C6F6 << G4endl;

  //     
  // Container
  //  
  G4double Rmin=0., Rmax=5*cm, deltaZ= 5*cm, Phimin=0., deltaPhi=360*degree;

  G4Tubs*  
  solidWorld = new G4Tubs("C6F6",                        //its name
                   Rmin,Rmax,deltaZ,Phimin,deltaPhi);        //its size

  G4LogicalVolume*                         
  logicWorld = new G4LogicalVolume(solidWorld,                //its solid
                                   C6F6,                //its material
                                   "C6F6");                //its name
  G4VPhysicalVolume*                                   
  physiWorld = new G4PVPlacement(0,                        //no rotation
                                   G4ThreeVector(),        //at (0,0,0)
                                 logicWorld,                //its logical volume
                                 "C6F6",                //its name
                                 0,                        //its mother  volume
                                 false,                        //no boolean operation
                                 0);                        //copy number

  //
  //always return the physical World
  //  
  return physiWorld;
}

---

粒子及粒子反应构建

void PhysicsList::ConstructEM()
{
  G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();

  auto particleIterator=GetParticleIterator();
  particleIterator->reset();
  while( (*particleIterator)() ){
    G4ParticleDefinition* particle = particleIterator->value();
    G4String particleName = particle->GetParticleName();

    if (particleName == "gamma") {

      ph->RegisterProcess(new G4PhotoElectricEffect, particle);
      ph->RegisterProcess(new G4ComptonScattering,   particle);
      ph->RegisterProcess(new G4GammaConversion,     particle);

    } else if (particleName == "e-") {

      ph->RegisterProcess(new G4eMultipleScattering, particle);
      ph->RegisterProcess(new G4eIonisation,         particle);
      ph->RegisterProcess(new G4eBremsstrahlung,     particle);      

    } else if (particleName == "e+") {

      ph->RegisterProcess(new G4eMultipleScattering, particle);
      ph->RegisterProcess(new G4eIonisation,         particle);
      ph->RegisterProcess(new G4eBremsstrahlung,     particle);
      ph->RegisterProcess(new G4eplusAnnihilation,   particle);      
    }
  }
}

---

粒子枪构建

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

  // default particle kinematic

  G4ParticleDefinition* particle
           = G4ParticleTable::GetParticleTable()->FindParticle("gamma");
  fParticleGun->SetParticleDefinition(particle);
  fParticleGun->SetParticleMomentumDirection(G4ThreeVector(1.,0.,0.));
  fParticleGun->SetParticleEnergy(9*MeV);
  fParticleGun->SetParticlePosition(G4ThreeVector(0.*cm,0.*cm,0.*cm));

}

void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
  //this function is called at the begining of event
  //
  //distribution uniform in solid angle
  //
  G4double cosTheta = 2*G4UniformRand() - 1., phi = twopi*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));

  fParticleGun->GeneratePrimaryVertex(anEvent);
}

---

编译运行

root@master:# cmake -DGeant4_DIR=/download/InstallForGeant/lib/Geant4-10.0.1 /download/geant4-master/examples/extended/electromagnetic/
root@master:# make -j2
root@master:# ./TestEm4
............
PhysicsList::SetCuts:CutLength : 1 (mm)
/control/cout/ignoreThreadsExcept 0
/run/verbose 2
#
/run/printProgress 10000
#
/run/beamOn 100000

phot:   for  gamma    SubType= 12  BuildTable= 0
      LambdaPrime table from 200 keV to 100 TeV in 61 bins 
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
       PhotoElectric :  Emin=        0 eV    Emax=      100 TeV   AngularGenSauterGavrila

compt:   for  gamma    SubType= 13  BuildTable= 1
      Lambda table from 100 eV  to 1 MeV, 7 bins per decade, spline: 1
      LambdaPrime table from 1 MeV to 100 TeV in 56 bins 
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
       Klein-Nishina :  Emin=        0 eV    Emax=      100 TeV

conv:   for  gamma    SubType= 14  BuildTable= 1
      Lambda table from 1.022 MeV to 100 TeV, 18 bins per decade, spline: 1
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
        BetheHeitler :  Emin=        0 eV    Emax=       80 GeV
     BetheHeitlerLPM :  Emin=       80 GeV   Emax=      100 TeV
............

---

图像交互：显示探测器

Idle> /vis/sceneHandler/create OGLIX
Idle> /vis/scene/create
Idle> /vis/specify C6F6
Idle> /vis/viewer/flush

未完待续