專案內容

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

未完待續