文|梁佐佐
对于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界面的视频演示，很到位。
图4. 视频教程演示QT操作
当然了，大家可能还需要某种能使用谷歌学术啦、谷歌啦之类的福利，关注公众号，后台回复“学术”即可给链接，简单易用靠谱！感谢读者竟然能读Geant4教程到本文结尾！！
网友评论