前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >Geant4--是怎样使用的?--(1.信息抽取)

Geant4--是怎样使用的?--(1.信息抽取)

作者头像
梁佐佐
发布2020-09-04 17:24:29
1.9K0
发布2020-09-04 17:24:29
举报
文章被收录于专栏:人芳觅人芳觅

对于Geant4模拟,我们关心它到底是怎样使用的,到底是怎样获取我们想要的信息,即信息抽取。前面几篇入门教程有提到,Geant4的模拟流程中从信息流的整合来看,物理过程框架可从大到小分为Run、Event、(Track)、Step(对应各自的*Action.cc)。

让我们再来看一下4者之间的组合关系:一个Run包含多个Event,每个Event包含多个Track,每个Track包含多个Step。如下图所示:

图1. Geant4关键函数调用关系

从B1例子来看,如果我们想知道每个Event总共沉积多少能量,只需要在SteppingAction中调用一个Event的储值变量,将该Event下的每个Step沉积能量累加到Event的储值变量中便可。

每个Event由很多Track及Step组成,其中,Track用来表征一个不变的粒子,并描述其径迹,Step是Geant4中最基本的蒙卡抽样概念,用以“试探”当前粒子的物理过程,当发生能量交换时,便会衍生出新的粒子即径迹,比如,一个gamma射线在某一Step发生康普顿效应时,便产生一个新粒子- ->一个吸收了能量的新电子,然后再各自描述其径迹。三者间的联系见下图:

图2. 一个Event,其中每段小节都是一个Step,Parent ID为当前粒子所继承上者的TrackID。

了解了这些基本概念,我们才能更加清晰地拿G4来干活,这些基础知识,需要多看教材、反复去琢磨才能更加透彻。继续以B1例子为基础讲解,怎样抽取你想要的物理过程信息。

进入/B1/,打开run1.mac,修改/tracking/verbose 1为2,数字越大,屏幕打印信息越丰富。编译B1,执行run1.mac。(敲黑板!/tracking/verbose2 真是个好东西呀,自己多看看打印信息,了解G4模拟时都能给出哪些信息,简单易读!

在exampleB1.cc中设置线程为1,便于我们分析。

代码语言:javascript
复制
#ifdef G4MULTITHREADED
 G4MTRunManager* runManager = new G4MTRunManager;
 runManager->SetNumberOfThreads(1);
#else
 G4RunManager* runManager = new G4RunManager;
#endif

运行后,得到结果如下图:

我们继续分析,发射第一个Event粒子为伽玛,在第2个Step时(Step#1)发生康普顿效应,产生次级粒子即电子,然后我们转到第二个Track信息——电子e-,从母粒子伽玛衍生出来,在Step#1时,发生电离(ProcessName “eIoni”),沉积能量(dE)为1.23MeV,Step#1的初始位置即Step#0的位置(X,Y,Z),Step#1结束时的动能(KinE)为2.51MeV。

重点开始了,以代码事例来看,怎样抽取信息。事例1:获取一个探测器“/B1/Shape1”的能谱;事例2:获得多个探测器“/B1/Shape1”的计数分布。

事例1:

1)B1SteppingAction.cc中需要收集每个Event中“Shape1”中的沉积能量和;

代码语言:javascript
复制
voidB1SteppingAction::UserSteppingAction(const G4Step* step)
{
 
  //get volume of the current step
G4LogicalVolume*volume= step->GetPreStepPoint()->GetTouchableHandle()
->GetVolume()->GetLogicalVolume();
G4Stringvolname=volume->GetName();
// check if weare in scoring volume
 
if (volname =="Shape1")
{
  //collect energy deposited in this step
 G4double edepStep = step->GetTotalEnergyDeposit();
fEventAction->AddEdep(edepStep);
}
}

2)B1EventAction.hh中定义AddEdep()和fEdep;

代码语言:javascript
复制
public:
voidAddEdep(G4double edep) { fEdep += edep; } 
private:
   G4double     fEdep;

3)B1EventAction.cc中在BeginOfEventAction(constG4Event*)初始化fEdep=0,在EndOfEventAction(const G4Event*aEvent)中输出每个Event沉积的能量。

代码语言:javascript
复制
#include <fstream>
void B1EventAction::BeginOfEventAction(constG4Event*)
{
 fEdep1 = 0.;
}
void B1EventAction::EndOfEventAction(constG4Event*)
{
fstream dataFile1;
dataFile1.open("outspectrum.dat",ios::app|ios::out);
dataFile1<< fEdep1 <<endl;

//在每个Event结束时,将沉积的总能量输出至文件,最后用Matlab等数据处理软件做直方图统计便可
G4int eventID= aEvent->GetEventID();
if ((eventID+1)%100000==0)G4cout<<" Now run out of "<< eventID<<endl;
}

事例2:

图3. 在B1例子上修改,放置10个改变体积后的Shape1,取消Shape2

1)B1SteppingAction.cc中在确定Step位于“Shape1”时,给出“Shape1”的CopyNumber;

代码语言:javascript
复制
if (volname == "Shape1")
{
G4StepPoint*      point2 =step->GetPostStepPoint();                      
G4TouchableHandletouch2 = point2 ->GetTouchableHandle();  
G4int copyno=touch2->GetCopyNumber();               
G4cout<<"copynumis "<<copyno<<G4endl;
  // collect energy deposited in this step
G4double edepStep = step->GetTotalEnergyDeposit();
fEventAction->AddEdep(copyno,edepStep);
//此时要修改AddEdep,将沉积能量累加到对应的逻辑体中
}

2)假如我们放置10个“Shape1”,copynumber从0-9,来看一下探测器构建、AddEdep()和一维数组fEdep的定义;

在B1DetectorConstruction.cc中修改:

代码语言:javascript
复制
for (int i=0;i<10;i++)
{
 G4ThreeVector pos1 = G4ThreeVector(0, 0, -i*2*cm);
  newG4PVPlacement(0,                      //no rotation
                    pos1,                    //at position
                    logicShape1,             //its logical volume
                    "Shape1",                //its name
                    logicEnv,                //its mother  volume
                    false,                   //no boolean operation
                    i,                       //copy number
                    checkOverlaps);          //overlaps checking
}
代码语言:javascript
复制
class B1EventAction : publicG4UserEventAction
{ 
 public:
   B1EventAction(B1RunAction* runAction);
   virtual ~B1EventAction();
    
   virtual void BeginOfEventAction(const G4Event* event);
   virtual void EndOfEventAction(const G4Event* event);
    
   void AddEdep(G4int copyno,G4double edep) { fEdep[copyno] += edep; }
  
 private:
   B1RunAction* fRunAction;
   G4double     fEdep[10];
};#endif
代码语言:javascript
复制
voidB1EventAction::BeginOfEventAction(const G4Event*)
{    
memset(fEdep, 0, sizeof(fEdep));
}

3)在B1RunAction中设立一个一维计数统计数组Counts[10],在每个Event结束的时候判断,如果fEdep[2]>0,则Counts[2]++,在Run结束时输出数组Counts[10]。

代码语言:javascript
复制
void B1EventAction::EndOfEventAction(constG4Event*)
{   
for (int i=0;i<10;i++)
{
if(fEdep[i]>0.) {fRunAction->Counts[i]++;}
} 
}
代码语言:javascript
复制
class B1RunAction : public G4UserRunAction
{ 
 public:
   B1RunAction();
   virtual ~B1RunAction();
   // virtual G4Run* GenerateRun();
   virtual void BeginOfRunAction(const G4Run*);
   virtual void   EndOfRunAction(const G4Run*);
    void AddEdep (G4double edep);
Counts[10]; //添加为共有变量,可以在B1EventAcion.cc中调用。
  private:
   G4Accumulable<G4double> fEdep;
   G4Accumulable<G4double> fEdep2;
}; 
#endif
 
代码语言:javascript
复制
void B1RunAction::BeginOfRunAction(constG4Run*)
{
 memset(Counts, 0, sizeof(Counts));
}
void B1RunAction::EndOfRunAction(constG4Run* run)
{

for (int i=0;i<10;i++)
{
G4cout<<Counts[i]<<"\t";//当然可以输入到一个文件中
}
G4cout<<G4endl;
}

这两个事例都是基于获取能量的探测器计数,当然在SteppingAction中灵活性可以更大,考虑到step所在的Track获取当前Step的粒子名称,以及Step所在的位置,还是非常灵活的。诸君加油,虽然不难,但写这些东西真是要麻烦哭了==

https://www.youtube.com/watch?v=nz31wvR83hI这是个介绍Geant4-QT界面的视频演示,很到位。

图4. 视频教程演示QT操作

当然了,大家可能还需要某种能使用谷歌学术啦、谷歌啦之类的福利,关注公众号,后台回复“学术”即可给链接,简单易用靠谱!感谢读者竟然能读Geant4教程到本文结尾!!

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

本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2019-04-29,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 人芳觅 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体分享计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档