一个半径 5.0 cm 的 U 球（这里教学目的设置得比较小，如果设置更大，很可能一个事件得运行几个小时到几天），热中子从球中间发射，示例代码如下：
wuPrimaryGeneratorActionAll.cc
void wuPrimaryGeneratorActionAll::GeneratePrimaries(G4Event* anEvent)
{
G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
G4ParticleDefinition* pp = 0;
pp = particleTable->FindParticle("neutron");
if(pp)
particleGun->SetParticleDefinition(pp);
else
G4cout<<"##Null pp in wuPrimaryGeneratorAction::SetParticleGun()"<<G4endl;
//primary particle kinetic energy
particleGun->SetParticleEnergy(0.025*MeV);
// primary particle position
G4double x = 0;
G4double y = 0;
G4double z = 0.;
particleGun->SetParticlePosition(G4ThreeVector(x, y, z));
// primary particle moving direction
G4double DirectRotX = 0.0*pi/180.*rad;
G4double DirectRotY = 0.0*pi/180.*rad;
G4double DirectRotZ = 0.0*pi/180.*rad;
G4double theta;
G4double RangeThetaPri = 180*pi/180.*rad;
G4double phi = G4UniformRand()*2.0*pi;
G4double rg = 1.0-cos(RangeThetaPri);
theta = acos(1.0-G4UniformRand()*rg);
G4double cosPX = sin(theta)*cos(phi);
G4double cosPY = sin(theta)*sin(phi);
G4double cosPZ = cos(theta);
G4ThreeVector directPri(cosPX, cosPY, cosPZ);
directPri.rotateX(DirectRotX);
directPri.rotateY(DirectRotY);
directPri.rotateZ(DirectRotZ);
particleGun->SetParticleMomentumDirection(directPri);
particleGun->GeneratePrimaryVertex(anEvent);
}
wuDetectorConstruction.hh
// 添加
G4LogicalVolume* logicHEU;
G4VPhysicalVolume* physHEU;
wuDetectorConstruction.cc
void wuDetectorConstruction::DefineMaterials()
{
// 这里省略 ...
// 添加这部分代码
G4Isotope* U4 = new G4Isotope("U234",92,234,234.02*g/mole);
G4Isotope* U5 = new G4Isotope("U235",92,235,235.01*g/mole);
G4Isotope* U6 = new G4Isotope("U236",92,236,236.04*g/mole);
G4Isotope* U8 = new G4Isotope("U238",92,238,238.03*g/mole);
G4Element* HEU4568 = new G4Element("Highly-enriched Uranium 4568","HEU4568",4);
HEU4568->AddIsotope(U4,0.0097);
HEU4568->AddIsotope(U5,0.9315);
HEU4568->AddIsotope(U6,0.0024);
HEU4568->AddIsotope(U8,0.0564);
G4Material* matHEU4568 = new G4Material("HEU4568",18.75*g/cm3,1);
matHEU4568->AddElement(HEU4568, 1.0);
// 这里省略 ...
}
G4VPhysicalVolume*wuDetectorConstruction::DefineVolumes()
{
// 这里省略 ...
// 相应位置添加以下部分代码
G4Material* heu_mat = G4Material::GetMaterial("HEU4568");
G4double HEUR0 = 0.0*cm;
G4double HEUR1 = 5.0*cm;
G4Sphere* solidHEU =
new G4Sphere("HEU",
HEUR0,
HEUR1,
0*degree,
360*degree,
0*degree,
180*degree );
logicHEU =
new G4LogicalVolume(solidHEU,
heu_mat,
"HEU");
physHEU =
new G4PVPlacement(0,
G4ThreeVector(0.0,0.0,0.0),
logicHEU,
"HEU",
logicWorld,
false,
0,
checkOverlaps);
// 这里省略 ...
}
注意：本示例代码在不熟悉之前，一次不要运行太多事件。建议刚开始调试时候只运行 10 个事件！！！如果 10 个事件中没有运行出几万行的事件来，则再加 10 个事件，直到出来几万行 entry 的事件。
t->Scan("EventID:ParentID:TrackID:CurrentStepNumber:PName:CreatorProcess:EkPre:EDep")
#或者
t->Scan("EventID:ParentID:TrackID:CurrentStepNumber:PName:CreatorProcess:EkPre:EDep","","colsize=20")
t->Scan("EventID:ParentID:TrackID:CurrentStepNumber:PName:CreatorProcess:EkPre:EDep:GlobalTimePre","","colsize=18")
可以看到一个事件产生几万 entry，对一个事件来说：
从粒子可以看到 “neutron”， “gamma”， “e-”， “e+”， “U235”， “alpha”， “Th231[387.836]”， “Th231[185.718]”， “Th231”， “anti_nu_e”， “Pa231[101.409]”， “Pa231[84.215]”， “Pa231”， “Ac227[46.350]”， “Ac227[27.370]”， “Ac227”， “Th227”， “Ra223[61.424]”， “Ra223”， “Rn219[598.720]”， “Rn219”， “Po215”， “Pb211”， “Bi211[831.960]”， “Bi211”， “Tl207”， 等等。
从物理相互作用过程可以看到 “nFission”， “neutronInelastic”， “phot”， “compt”， “nCapture”， “conv”， “hadElastic”， “RadioactiveDecayBase”， “annihil”， 等等。
这几万 entry 数据中，很多是我们不感兴趣的。我们尝试在数据输出中把它屏蔽掉。
例如，想知道一个中子入射产生多少个次级中子，因此把非中子的 entry 屏蔽掉。
// 在数据输出前添加以下代码
if(PName.data()[0]!='n') return;
不做任何数据筛选之前，一个热中子裂变事件产生约 20000+ 多 entry 数据, 数据量约为 4 MB。筛选之后数据只剩 100+ entry，数据量约为 60 kB。此时发现同一个粒子，它的 tracking 中有多个 step, 我们只需要记录一个即可，再加以下限制条件，同一个事件只剩 50+ entry。
// 在数据输出前添加以下代码
if(CurrentStepNumber != 1) return;
这里需要提醒一点，“neutronInelastic” 过程中子的 Track ID 会发生改变，但是它不属于 “nFission” 产生的新粒子。
在进行条件筛选时，采用的条件必须尽可能简单，避免筛选条件不当导致模拟出现偏差。
