Day 13 root 相关说明--以 ANAEX01 为实例

Day 13 root 相关说明–以 ANAEX01 为实例

一、项目，分为4个部分：

1. 建模：

• ​ 反应模型 DetectorConstruction，模型 + 模型材料

• ​ 粒子模型 PrimaryGeneratorAction

2. Physics list ：

粒子反应的物理过程，一般使用G4已经写好的

3. 收集物理信息，四种方式：

• ​ 通过动作类； 通过 SD + hit； G4Step，G4Track； …

4. 把收集到的信息进行输出：

写进输出文件（root）里面输出出来，方便之后再对这些信息进行再加工，比如画图、能谱分布、位置信息、能量信息等，另外的还有 HDF5、AIDA XML、CSV，这里主说 root 文件。储存（数据）信息的两种形式：Ntuple、histogram。root 是一个工具包，用来做数据分析；

二、Root 在 geant4 中的使用分为两个部分：

• 第一个部分是 root 文件的生成（写入），也就是在模拟过程中把数据写进去，然后模拟结束后生成一个 .root 文件；
• 第二部分是 root 文件中数据的读取（重新加工），即 root 文件处理。

1. 怎么生成一个 .root 文件：

首先 cmake 编译，然后执行一个 .mac 文件；

2. 查看命令：

在终端输入 root，然后 TBrowser 或者 new TBrowser ，会出来操作框，然后双击 root 文件，然后点击需要查看的数据就可以看见图片，退出命令是 .q 。

通常在模拟时，将想要做的命令都写在 run.mac 文件中，通常一个 .mac文件就一行 /run/beamOn …（发射一次run），正常来说一个模拟实验就是一个 /run/beamOn …，但是有例外，看自己的需求，G4会根据我们写的 /run/beamOn …的个数，自动编码，从0开始，一个一个的进行。

3. .root 文件

正常来说，一个 /run/beamOn 就只会生成一个 .root 文件（单线程）。

4. 什么是多线程

于单线程而言，先跑一个 run ，run 内部有很多个 event，去循环每一个event（一个粒子一个event，n 个粒子 n 个 event），每个 event 有一个 eventAction，有一个 beginOfEventAction 和 endOfEventAction，然后，有runAction。

假如有100个粒子，我想要分成两个线程，那么自动把100个粒子分为两个部分，两个部分粒子数和为100，每个部分都有一个单独的 /run/beamOn …，每个部分走各自的线程，每一个子线程都有自己的runAction、eventAction、TrackingAction、SteppingAction…，每一个线程都会有一个自己的 .root 文件，会生成多个 .root 文件，但最后都会有一个 merge 操作，将几个 .root 文件进行汇总，实现一个 run 只有一个 root，这个汇总的 .root 文件才是我们需要的文件。

// 针对使用的Ntuple元组，使用这一句话进行merge
analysisManager -> SetNtupleMerging(true);
// 针对使用的histogram直方图，不需要使用这一句话进行merge，自动实现merge

有时候一个 run 会生成很多个 .root 文件，因为使用了Ntuple但缺少了一个命令 Nutuple 相关语句，导致 root 没有写好，出现问题，本来应该存储在一个 .root 文件中的，但是被分散在了多个 .root 文件中，进而没有 merge 在一个文件中。

5. 什么时候打开文件、关闭文件：

root 文件的打开关闭通常都是放在runAction文件中进行的（多线程），但是程序模拟的操作都需要在主函数 main中进行。单线程的话，在 main 函数中就可以进行操作打开关闭。但是一般都写在runAction中的。

写入（固定的套路写入）

三、Histogram + Ntuple

1. Histogram（柱状图、直方图）

存在一维、二维、三维的直方图（ histogram ），但是最常用的是一维；运用实例：统计粒子的能量分布、统计 event 的能量分布。

使用 histogram 必须先设置 bin ，bin 的选取直接决定了最后的直方图效果。

• 统计粒子的能量分布：

  是从粒子进入世界就开始进行统计，于 histogram 而言，首先需要设计 bin，假如指定有 100 个 bin，粒子能量指定总共 2Mev，那么每个 bin 的横间距就是 2/100

• 统计 event 的能量沉积：

  每一个event（每一个初级粒子及其次级粒子）反应产生的次级粒子包括初级粒子本身在物体里面的能量沉积，并总结画出直方图；最终得到histogram时，会做一个归一化处理

2. Ntuple（table、元组）

同时统计很多个数据，例：同时统计粒子的能量以及粒子对应的射出靶物质的 X Y Z 坐标，并存放在一个 table 中，这些物理量都必须是来自同一个粒子的信息（histogram 只能一次统计 X 或 Y 或 Z），

// 针对使用的Ntuple元组，使用这一句话进行merge
analysisManager -> SetNtupleMerging(true);

最终还需要将多个里面的是数据读取出来再加工才可以用在文章中。会比较好看且直观。

3. 怎么样正式操作使用Histogram / Ntuple

HistoManager（约定俗成的类，里面包含如下）：

• 指明要生成什么类型的文件来储存数据，root 、xml、csv，在 HistoManager.hh 文件中需要包含对应的头文件；

  #include "g4root.hh"
  #include "g4xml.hh"
  #include "g4csv.hh"
  

• 打开文件：定义一个方法，里面包含打开这个文件、写数据、关闭这个文件…

  打开文件的方法：

  void HistoManager::Book()
  {G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();if (!fFactoryOn) {//有很多方法，根据需要选择analysisManager->SetVerboseLevel(1);analysisManager->SetNtupleMerging(true);//如果项目里面用到了Ntuple以及多线程，就一定要写这一句话analysisManager->SetHistoDirectoryName("histo");//取名字analysisManager->SetNtupleDirectoryName("ntuple");//取名字}// 打开文件的方法G4bool fileOpen = analysisManager->OpenFile("AnaEx01");//这里双引号内部填写的生成的.root文件的名字，analysisManager是上述所获取的对象if (!fileOpen) {G4cerr << "\n---> HistoManager::Book(): cannot open " << analysisManager->GetFileName()<< G4endl;return;}...
  }
  

• 需要在这个类里面指明要创建的Histogram / Ntuple：

  1. **Histogram：**选择自己要创建几维的直方图；注意每个 Histogram 的 ID 。

  // id = 0，五个参数分别为：名字，标题，bin的个数，最低能量，最高能量（这两个表示了x轴的取值范围）
  // CreateH1 表示这个 Histogram 图像是一维的
  analysisManager->CreateH1("EAbs", "Edep in absorber (MeV)", 100, 0., 800 * MeV);
  // id = 1
  analysisManager->CreateH1("EGap", "Edep in gap (MeV)", 100, 0., 100 * MeV);
  // id = 2
  analysisManager->CreateH1("LAbs", "trackL in absorber (mm)", 100, 0., 1 * m);
  // id = 3
  analysisManager->CreateH1("LGap", "trackL in gap (mm)", 100, 0., 50 * cm);
  //按顺序创建之后，id会默认编排
  

  本实例中共构建了4个 Histogram ，在最终的root窗口中histogram会有4个文件，分别为EAbs、EGap、LAbs、LGap。

  2. Ntuple：注意储存的数据类型（G4int、G4float、G4double、G4String）；最后不要忘记 FinishNtuple() ；注意每个 Ntuple 的 ID 。

  //创建第一个 ntuple (id = 0)，储存能量，分别存储gap、absorber
  analysisManager->CreateNtuple("Ntuple1", "Edep");
  //这个ntuple有两列，分别存absorber的和gap的
  analysisManager->CreateNtupleDColumn("Eabs");  // column Id = 0
  analysisManager->CreateNtupleDColumn("Egap");  // column Id = 1
  // 必须加上这一句
  analysisManager->FinishNtuple();// 创建第二个 ntuple (id = 1)，储存每一个event的总的tracklength
  analysisManager->CreateNtuple("Ntuple2", "TrackL");
  analysisManager->CreateNtupleDColumn("Labs");  // column Id = 0
  analysisManager->CreateNtupleDColumn("Lgap");  // column Id = 1
  analysisManager->FinishNtuple();
  

  CreateNtupleDColumn 里面的 D 表示存储的是 double 数据

  设置了两个 Ntuple ，所以分开设置，会在 root 里面出现对应的文件，统计的信息上面 histogram 一致，其中这两个 Nutple 可以写成一个 Ntuple

• **保存关闭文件：**创建好Histogram / Ntuple之后，数据需要进行filing填充，填充之后就要关闭文件 save() 方法

  void HistoManager::Save()
  {if (!fFactoryOn) {return;}
  // 固定的写法G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();analysisManager->Write();analysisManager->CloseFile();G4cout << "\n----> Histograms and ntuples are saved\n" << G4endl;
  }
  

4. 怎么样填充数据：

创建好Histogram / Ntuple之后，数据需要进行filing填充。一个项目进行数据信息收集，当一个event结束之后，会收集到一个event对应的总和的fEnergyAbs、fEnergyGap、fTrackLAbs、fTrackLGap，然后把这些数值存到 Histogram / Ntuple 中 ：

void EventAction::EndOfEventAction(const G4Event*)
{fRunAct->FillPerEvent(fEnergyAbs, fEnergyGap, fTrackLAbs, fTrackLGap);//调用方法进行填充//填充 histogramefHistoManager->FillHisto(0, fEnergyAbs);fHistoManager->FillHisto(1, fEnergyGap);fHistoManager->FillHisto(2, fTrackLAbs);fHistoManager->FillHisto(3, fTrackLGap);// 填充 ntuplefHistoManager->FillNtuple(fEnergyAbs, fEnergyGap, fTrackLAbs, fTrackLGap);
}

填充方法的设计：

//填充histograme的方法，参数分别为：id，对应的要存的信息，weight（默认等于1）
void HistoManager::FillHisto(G4int ih, G4double xbin, G4double weight)
{G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();analysisManager->FillH1(ih, xbin, weight);
}...//填充ntuple的方法
void HistoManager::FillNtuple(G4double energyAbs, G4double energyGap, G4double trackLAbs,G4double trackLGap)
{G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();// 第一个 ntuple ( id = 0)，参数：ntuple的参数，这个ntuple的第几列，这一列的填充什么数据analysisManager->FillNtupleDColumn(0, 0, energyAbs);analysisManager->FillNtupleDColumn(0, 1, energyGap);analysisManager->AddNtupleRow(0);// 第二个 ntuple ( id = 1)analysisManager->FillNtupleDColumn(1, 0, trackLAbs);analysisManager->FillNtupleDColumn(1, 1, trackLGap);analysisManager->AddNtupleRow(1);

细节：

• fill() 方法的调用位置：填充数据的方法；在每一个event结束后获取到对应的信息之后被调用。

  多线程：每一个线程结束后都会有一个对应的End of Run，

• book() 方法调用位置：book()方法包含指明头文件、打开文件、创建 histogram / ntuple；

• save() 方法调用位置：包含保存并关闭文件；

• HistoManager 初始化位置，初始化跟随 RunAction

五、补充：DetectorConstruction之G4VReplica

1. G4VReplica：

• 用来建立沿轴向的、重复的、没有间隙的 volume；
• 是 G4VPhysicalVolume 的之类，类似 G4PVPlacement；
• 一个实例 new G4PVReplica() 表示的就是重复之后的整体。

整体叫：CalorParameters

 fPhysiLayer = new G4PVReplica("Layer",  // its namefLogicLayer,  // its logical volumefLogicCalor,  // its motherkXAxis,  // axis of replicationfNbOfLayers,  // number of replicafLayerThickness);  // witdth of replica}

2. 为什么要使用 G4PReplica ：

建立一个 volume 分成三个步骤：首先根据模型尺寸建立 solid volume ，然后根据材料填充建立 logical volume ，最后判断物体是否需要旋转平移位置信息来建立 physics volume 。

还有很多其他的各种各样的类，实现重复建立 volume 的操作。

• G4PVDivision：有间隙的、重复的、沿轴向的；
• G4ReflectionFactory
• G4AssemblyVoulme
• …

于这个实例而言，一个 layer 由一个 absorber + gap 组成：

DetectorConstruction::DetectorConstruction()
{ComputeCalorParameters();//调用缩写的方法// 设置absorber和gap的材料为Pb和ArDefineMaterials();SetAbsorberMaterial("G4_Pb");SetGapMaterial("G4_lAr");fDetectorMessenger = new DetectorMessenger(this);
}
...G4VPhysicalVolume* DetectorConstruction::Construct()
{// 调用所写的方法return ConstructCalorimeter();
}void DetectorConstruction::DefineMaterials()
{// 定义真空的办法，下述两句G4NistManager* man = G4NistManager::Instance();fDefaultMaterial = man->FindOrBuildMaterial("G4_Galactic");// 设置材料man->FindOrBuildMaterial("G4_Pb");man->FindOrBuildMaterial("G4_lAr");G4cout << *(G4Material::GetMaterialTable()) << G4endl;
}G4VPhysicalVolume* DetectorConstruction::ConstructCalorimeter()
{G4GeometryManager::GetInstance()->OpenGeometry();G4PhysicalVolumeStore::GetInstance()->Clean();G4LogicalVolumeStore::GetInstance()->Clean();G4SolidStore::GetInstance()->Clean();ComputeCalorParameters();// 1.定义世界：solid、logical、physical volume...
// calorimeter相关设置fSolidCalor = nullptr;fLogicCalor = nullptr;fPhysiCalor = nullptr;fSolidLayer = nullptr;fLogicLayer = nullptr;fPhysiLayer = nullptr;if (fCalorThickness > 0.) {//定义box一定是定义半长fSolidCalor = new G4Box("Calorimeter",fCalorThickness / 2, fCalorSizeYZ / 2, fCalorSizeYZ / 2);fLogicCalor = new G4LogicalVolume(fSolidCalor, fDefaultMaterial, // fDefaultMaterial表示真空，想要layer去填充这个"Calorimeter");fPhysiCalor = new G4PVPlacement(nullptr,G4ThreeVector(),fLogicCalor, "Calorimeter", fLogicWorld, false, 0); 
// layer相关设置：三个volume的建立fSolidLayer = new G4Box("Layer",fLayerThickness / 2, fCalorSizeYZ / 2, fCalorSizeYZ / 2);  // sizefLogicLayer = new G4LogicalVolume(fSolidLayer, fDefaultMaterial,// fDefaultMaterial表示真空，因为要放入Absorber和Gap两个东西进行填充"Layer"); if (fNbOfLayers > 1) {fPhysiLayer = new G4PVReplica("Layer",fLogicLayer, fLogicCalor,//上一级物体，放在上一级物体之中kXAxis,     // 重复的方向，这里沿着x轴fNbOfLayers, // 重复多少个单元fLayerThickness); }else {fPhysiLayer = new G4PVPlacement(nullptr,G4ThreeVector(), fLogicLayer,  "Layer",  fLogicCalor,  false, 0); }}// Absorber相关设置fSolidAbsorber = nullptr;fLogicAbsorber = nullptr;fPhysiAbsorber = nullptr;
//Absorber三个volume的建立if (fAbsorberThickness > 0.) {fSolidAbsorber = new G4Box("Absorber",fAbsorberThickness / 2, fCalorSizeYZ / 2, fCalorSizeYZ / 2);fLogicAbsorber = new G4LogicalVolume(fSolidAbsorber,fAbsorberMaterial,fAbsorberMaterial->GetName()); fPhysiAbsorber = new G4PVPlacement(nullptr,G4ThreeVector(-fGapThickness / 2, 0., 0.), //位置设置，位置设置一般为指定物体的中心位置fLogicAbsorber,  fAbsorberMaterial->GetName(),  fLogicLayer, //上一级是LogicaLayer，因为要放在Layer中false,  0); }// Gap的相关设置fSolidGap = nullptr;fLogicGap = nullptr;fPhysiGap = nullptr;
//Gap三个volume的建立if (fGapThickness > 0.) {fSolidGap = new G4Box("Gap", fGapThickness / 2, fCalorSizeYZ / 2, fCalorSizeYZ / 2);fLogicGap = new G4LogicalVolume(fSolidGap, fGapMaterial, fGapMaterial->GetName());fPhysiGap = new G4PVPlacement(nullptr, G4ThreeVector(fAbsorberThickness / 2, 0., 0.),  // 位置fLogicGap, fGapMaterial->GetName(),  fLogicLayer,  //上一级是LogicaLayer，因为要放在Layer中false, 0); }...// 其他设置

**注：**头文件：.hh；

源文件 .cc。

获取每一个step的steplength

stepl = aStep -> GetStepLength();//固定的方法