文|梁佐佐

應唐光毅博士/后之約，對于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模擬效率提升。

總結

以上是生活随笔為你收集整理的geant4运行例子_Geant4--一次编译，运行多个Run，极大提升模拟效率的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。