Geant4--是怎样使用的？--（1.信息抽取）

梁佐佐

发布于 2020-09-04 17:24:29

2.2K0

发布于 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，便于我们分析。

#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”中的沉积能量和；

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；

public:
voidAddEdep(G4double edep) { fEdep += edep; } 
private:
   G4double     fEdep;

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

#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；

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中修改：

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
}

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

voidB1EventAction::BeginOfEventAction(const G4Event*)
{    
memset(fEdep, 0, sizeof(fEdep));
}

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

void B1EventAction::EndOfEventAction(constG4Event*)
{   
for (int i=0;i<10;i++)
{
if(fEdep[i]>0.) {fRunAction->Counts[i]++;}
} 
}

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

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界面的视频演示，很到位。