文|梁佐佐
應唐光毅博士/后之約，對于Geant4模擬，我們看是否能解決這么一個問題：我現在想模擬探測器不同角度下的響應，每次模擬需要/run/beamOn 100， 可是我真的不想一遍一遍的去http://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.0 45.1 5.0 表示從0.0°開始，每間隔5.0°賦予一次數值，到小于45.1為止，和for循環很像
2.2 angle.mac
/B5/detector/armAngle {angle} deg
/run/beamOn 100
3.在http://SteppingAciton.cc中累積每個Event的能量沉積(需要在B5例子中添加SteppingAciton函數，可仿照B1)，在http://B5EventAction.cc中抽取計數信息并輸出文件。
4.運行exampleB5 loop.mac 大功告成！
那么Geant4中具體應該怎樣實現？以B5例子為依托，上代碼！
第1步：
http://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 armauto 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 geometryG4RunManager::GetRunManager()->GeometryHasBeenModified();
}
//SetArmAngle()這個函數確實也應該包含在DefineCommands()中，借以實際改變探測角度及位置//~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
void B5DetectorConstruction::DefineCommands()
{//Define /B5/detector command directory using generic messenger classfMessenger = new G4GenericMessenger(this,"/B5/detector/","Detector control");//armAngle commandauto& 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）http://SteppingAction.cc中設置
//判斷當前step位于探測器幾何中。。。。。。
G4double edep = aStep->GetTotalEnergyDeposit();
theEvent->AddEdep(edep);
//將edep累加給fEdep，fEdep為每個Event沉積的總能量，在http://EventAction.cc中需要放在BeginOfEventAction(const G4Event* aEvent)中初始化。
2) EventAction.hh和http://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模擬效率提升。
設計講解)
