专栏首页人芳觅Geant4--一次编译,运行多个Run,极大提升模拟效率

Geant4--一次编译,运行多个Run,极大提升模拟效率

应唐光毅博士/后之约,对于Geant4模拟,我们看是否能解决这么一个问题:我现在想模拟探测器不同角度下的响应,每次模拟需要/run/beamOn 100, 可是我真的不想一遍一遍的去DetectorConstruction.cc中修改几何放置角度,然后编译完怒敲exampleB1 run1.mac;或者,我想只编译运行一次G4就可以跑几百次/run/beamOn 100 且需要每次Run的时候射线源的出射位置、能量等参数不同?

这么机智的事情,有助于解决做模拟会哭的问题。让我们开始吧!

以G4中的basic/B5 例子为基础,我们现在要模拟第一个场景:

  1. 设置一个探测器,绕Y轴可设置不同的旋转角度θ,θ范围为0°-45°,分别 间隔5°采样一次;
  2. 射线源在每个角度下/run/beamOn 100;
  3. 要求得到每个角度下探测器探测到的计数,可以认为此目的是对比探测器在不同射线入射角度下的探测效率;
  4. 总共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.在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模拟效率提升。

喜欢的话,分享一下吧~^o^~

本文分享自微信公众号 - 人芳觅(lililiangzuozuo),作者:梁佐佐

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-04-23

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • Geant4的几何边界判断代码解析

    以伽玛在Boundary2处做讲解分析。首先判断PreStepPoint所在的physicalvolume是否为”phybgo”,然后通过以下三种方式给出Bou...

    梁佐佐
  • ROOT-数据读取-直方图-Roofit拟合基本流程-(入门实用)

    笔者最近在测核素能谱,测出的能谱需要分析,比如计算某全能峰的分辨率。用到的数据处理分析工具是ROOT(cern),整个能谱读取分析的流程可给各位看官当入门或干货...

    梁佐佐
  • “网上冲浪“--第一篇 巧用关键词搜索

    身处信息时代,每天带给我们更多满足感的,相比手机电脑本身,是极其精彩的网络世界。可是,“网上冲浪”、“科学上网”这些基本技能我们又了解多少。说起百度、谷歌、必应...

    梁佐佐
  • 聊聊ArrayList源码(基于JDK1.8)

    打个广告,楼主自己造的轮子,感兴趣的请点[github]: https://github.com/haifeiWu/lightconf

    haifeiWu
  • 【Java集合源码剖析】ArrayList源码剖析

    转载请注明出处:http://blog.csdn.net/ns_code/article/details/35568011

    bear_fish
  • Hive-1.2.1_03_DDL操作 1.1. Create Database1.2. Drop Database1.3. Use Database2.1. Cre

      注意:各个语句的版本时间,有的是在 hive-1.2.1 之后才有的,这些语句我们在hive-1.2.1中是不能使用的。

    踏歌行
  • Java集合(3)---Java集合ArrayList

    ArrayList 是一个数组队列,相当于 动态数组。与Java中的数组相比,它的容量能动态增长。它继承于AbstractList,实现了List, Rando...

    葆宁
  • Java集合源码剖析——ArrayList源码剖析

    ArrayList是基于数组实现的,是一个动态数组,其容量能自动增长,类似于C语言中的动态申请内存,动态增长内存。

    Java团长
  • Java集合深度解析之ArrayList

    KS Knowledge Sharing 知识分享 现在是资源共享的时代,同样也是知识分享的时代,如果你觉得本文能学到知识,请把知识与别人分享。 转自:...

    互扯程序
  • Java 集合系列03之 ArrayList详细介绍(源码解析)和使用示例

    概要 上一章,我们学习了Collection的架构。这一章开始,我们对Collection的具体实现类进行讲解;首先,讲解List,而List中ArrayLis...

    老白

扫码关注云+社区

领取腾讯云代金券