文|梁佐佐
应唐光毅博士/后之约，对于Geant4模拟，我们看是否能解决这么一个问题：我现在想模拟探测器不同角度下的响应，每次模拟需要/run/beamOn 100， 可是我真的不想一遍一遍的去DetectorConstruction.cc中修改几何放置角度，然后编译完怒敲exampleB1 run1.mac；或者，我想只编译运行一次G4就可以跑几百次/run/beamOn 100 且需要每次Run的时候射线源的出射位置、能量等参数不同？
这么机智的事情，有助于解决做模拟会哭的问题。让我们开始吧！
以G4中的basic/B5 例子为基础，我们现在要模拟第一个场景：
a. 设置一个探测器，绕Y轴可设置不同的旋转角度θ，θ范围为0°-45°，分别 间隔5°采样一次；
b. 射线源在每个角度下/run/beamOn 100;
c. 要求得到每个角度下探测器探测到的计数，可以认为此目的是对比探测器在不同射线入射角度下的探测效率；
d. 总共10个角度，定义一个输出文件，总共输出10个数值，代表不同角度下的测得计数。
以G4中的basic/B5 例子为基础，我们可以分以下几步实现上述场景：
1. 定义宏命令/B5/detector/armAngle X deg 用以在*.mac文件中设置探测器角度,B5中，这是现成的！
2. 关键点——定义一个loop.mac 和一个angle.mac
2.1 loop.mac
/run/initialize
/gun/particle gamma
/gun/energy 611.7 keV
/control/loop angle.mac angle 0.0 45.1 5.0
## 0.045.15.0表示从0.0°开始，每间隔5.0°赋予一次数值，到小于45.1为止，和for循环很像
2.2 angle.mac
/B5/detector/armAngle {angle} deg
/run/beamOn 100
3.在SteppingAciton.cc中累积每个Event的能量沉积(需要在B5例子中添加SteppingAciton函数，可仿照B1)，在B5EventAction.cc中抽取计数信息并输出文件。
4.运行exampleB5 loop.mac 大功告成！
那么Geant4中具体应该怎样实现？以B5例子为依托，上代码！
第1步：
B5DetectorConstruction.cc中给出了怎样通过UI命令调整探测器臂角度的事例：
B5DetectorConstruction::B5DetectorConstruction(): G4VUserDetectorConstruction(),
fMessenger(nullptr),
fHodoscope1Logical(nullptr), fHodoscope2Logical(nullptr),
fWirePlane1Logical(nullptr), fWirePlane2Logical(nullptr),
fCellLogical(nullptr), fHadCalScintiLogical(nullptr),
fMagneticLogical(nullptr),
fVisAttributes(),
fArmAngle(30.*deg), fArmRotation(nullptr),
fSecondArmPhys(nullptr)
// fArmAngle参数就是要改变的角度θ
{
fArmRotation = new G4RotationMatrix();
fArmRotation->rotateY(fArmAngle);
// fArmRotation为改变探测器角度的实际“参数”
//define commands for this class
DefineCommands();
// DefineCommands()这个函数放在初始化列表中就是为了通过UI命令直接定义初始化
//B5DetectorConstruction.cc，毫无疑问，DefineCommands()中一定是
//定义了有关fArmAngle怎样改变探测器角度的方法
}
//~~~~~~~~~~~~~~~~~~~~//
//second arm
auto secondArmSolid =new G4Box("secondArmBox",2.*m,2.*m,3.5*m);
auto secondArmLogical =new G4LogicalVolume(secondArmSolid,air,"secondArmLogical");
auto x = -5.*m * std::sin(fArmAngle);
auto z = 5.*m * std::cos(fArmAngle);
fSecondArmPhys =new G4PVPlacement(fArmRotation,G4ThreeVector(x,0.,z),
secondArmLogical,"fSecondArmPhys",worldLogical,false,0,checkOverlaps);
// fSecondArmPhys就是我们要放置的“探测器”，它的旋转角度以及三维位置都是由
//fArmAngle直接决定
//~~~~~~~~~~~~~~~~~~~~~~~~~//
void B5DetectorConstruction::SetArmAngle(G4double val)
{
if(!fSecondArmPhys) {
G4cerr << "Detector has not yet been constructed."<< G4endl;
return;
}
fArmAngle = val;
*fArmRotation = G4RotationMatrix();
// make it unit vector
fArmRotation->rotateY(fArmAngle);
//fArmRotation为探测器的角度设置
auto x = -5.*m * std::sin(fArmAngle);
auto z = 5.*m * std::cos(fArmAngle);
fSecondArmPhys->SetTranslation(G4ThreeVector(x,0.,z));
//fSecondArmPhys->SetTranslation(…)为探测器的位置设置
//tell G4RunManager that we change the geometry
G4RunManager::GetRunManager()->GeometryHasBeenModified();
}
//SetArmAngle()这个函数确实也应该包含在DefineCommands()中，借以实际改变探测角度及位置
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
void B5DetectorConstruction::DefineCommands()
{
//Define /B5/detector command directory using generic messenger class
fMessenger = new G4GenericMessenger(this,"/B5/detector/","Detector control");
//armAngle command
auto& armAngleCmd=
fMessenger->DeclareMethodWithUnit("armAngle","deg",
&B5DetectorConstruction::SetArmAngle,
"Set rotation angle of the second arm.");
armAngleCmd.SetParameterName("angle", true);
armAngleCmd.SetRange("angle>=0. && angle<180.");
armAngleCmd.SetDefaultValue("30.");
}
// DefineCommands()作为调整探测器角度位置的核心函数，定义了一个可以在*.mac中调用的命令，
// /B5/detector/armAngle 10 deg 表示将fArmAngle设置为10°，相应的探测角度和位置也与之对应改变
第2步：
设置loop.mac 和angle.mac，略
第3步：
1）SteppingAction.cc中设置
//判断当前step位于探测器几何中。。。。。。
G4double edep = aStep->GetTotalEnergyDeposit();
theEvent->AddEdep(edep);
//将edep累加给fEdep，fEdep为每个Event沉积的总能量，在EventAction.cc中需要放在BeginOfEventAction(const G4Event* aEvent)中初始化。
2) EventAction.hh和EventAction.cc设置
在EventAction.hh中设置私有变量 realcounts=0 和tempcouts=0。注意这两个变量不能放在BeginOfEventAction()中初始化。具体的用法如下：
void EventAction::EndOfEventAction(const G4Event* aEvent)
{
if(fEdep>0.0) realcounts++;
G4long event_id = aEvent->GetEventID()+1;
fstream datafile;if((event_id) % 100 == 0) {
G4cout<<"Event"<<event_id<<" is over"<<G4endl;
datafile.open("outputcounts.xls",ios::out|ios::app);
G4int outcounts=realcounts-tempcounts;
datafile <<outcounts<<G4endl;
tempcounts = realcounts;
datafile.close();
}
//关键部分，每跑100个粒子输出一次探测器计数
第4步：
make - -> exampleB5 loop.mac 即可。
总结：
通过 /control/loop 配合UI改变角度参数进而一次性运行多次Run，每次Run对应的角度参数不同，在EventAction中设置输出参数，realcounts=0 和tempcouts=0需要放置在EventAction.hh中初始化，tempcouts总是等于上一次Run之后的realcounts数值，巧妙利用EventID识别第几次Run完结，作为输出计数和文件的节点。
第二个场景：
跑几百次Run，每次Run的射线源位置或者属性不同。
第1步：
void MYPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
const G4Run *nowrun=G4RunManager::GetRunManager()->GetCurrentRun();
G4int runid=nowrun->GetRunID();
//此时，runid就是一个反映当前第几个Run的变量，以每次Run100个粒子为例，也可通过设置int(eventID/100)来替代runid，二者等价
G4int posx,posy,posz;
posx=int(floor(runid/64));
posy=int(floor((runid%64)/8));
posz=int(floor((runid%64)%8));
particleGun->SetParticlePosition(G4ThreeVector((posx-3.5)*3.2*mm,(posy-3.5)*3.2*mm,(posz-3.5)*3.2*mm));
}
第2步：
定义一个loop.mac
/run/initialize
/control/loop rungun.mac 0 1 512
##总共跑512次，0-511
定义一个rungun.mac
/run/beamOn 100
第3步：
同场景1一样，略。
第4步：
Make - -> exampleMY loop.mac 大功告成！
大总结：
/control/loop用好这个命令，助力G4模拟效率提升。
网友评论