数据提取¶

数据结构¶

模拟程序运行结束之后，会生成文件 outputdata_tXX.root, 文件中有 40 个 branch，

ana1

EventID 模拟事件标号
ParentID 当前粒子的父粒子编号（0表示粒子源发射的粒子）
TrackID 在一个事件中的粒子编号（从 1 开始）
CurrentStepNumber 当前step 在track中的编号（从1开始）
PName 粒子名称
TrackWeight 权重
CreatorProcess 相互作用类型，"##" 为程序设置
EDep 当前step的能损
TrackLength
StepLength
TrackStatus 状态
Mass 粒子的质量
Charge 粒子的电荷态
MagneticMoment
VolNamePre 所在逻辑体的名称。前面要求每个探测器初始化唯一的逻辑体，目的在于此处通过名字可以区分不同的探测器
VolNamePost
GlobalTimePre 全局时间，事件开始为时间起点
GlobalTimePost
LocalTimePre 当前 track 时间，track 开始为时间起点
LocalTimePost
ProperTimePre
ProperTimePost
StepStatusPre
StepStatusPost
EkPre 动能
EkPost
xPre 位置x
yPre 位置y
zPre 位置z
xPost
yPost
zPost
xMomentumDirectionPre 动量方向x分量
yMomentumDirectionPre
zMomentumDirectionPre
Velocitypre 速度
xMomentumDirectionPost
yMomentumDirectionPost
zMomentumDirectionPost
Velocitypost

t->Scan("EventID:ParentID:TrackID:CurrentStepNumber:PName:EDep:VolNamePre:VolNamePost")

ana2

ana3

t->MakeClass("dataana")

生成 dataana.h 和 dataana.C 两个文件。将 dataana.h 中 char 型数组的长度修改成一个较大的数值。（因为默认的长度是用来生成该代码的文件数据中的最大长度，但是其它文件可能有比它还大的）。如下修改为 128，也可以改成更大，如 1024 等

ana4

统计穿过 DSSD 的事件数¶

void dataana::Loop()
{

  int count = 0;

   if (fChain == 0) return;

   Long64_t nentries = fChain->GetEntriesFast();

   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=0; jentry<nentries; jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // if (Cut(ientry) < 0) continue;

       // 示例代码中，两个探测器名字分别为 Clover、DSSD
       // 如果名字第一个字符为 D
      if(VolNamePre[0] == 'D') count++;
   }

   std::cout<<count<<std::endl;
}

穿过 DSSD 的事件数除以模拟的总事件数，估计立体角。

材料设置为真空
粒子不与物质发生相互作用
4$\pi$ 均匀发射

统计在 Clover 上的能量沉积¶

void dataana::Loop()
{

  int ID = -1;
  TH1I *h = new TH1I("h", "", 2500, 0, 2.5);
  double ene = 0;

   if (fChain == 0) return;

   Long64_t nentries = fChain->GetEntriesFast();

   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=0; jentry<nentries; jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // if (Cut(ientry) < 0) continue;

      // 示例代码中，两个探测器名字分别为 Clover、DSSD
      if(VolNamePre[0] != 'C') continue;
      // 粒子有  e-、e+、gamma 
      if(PName[0] != 'e') continue;

      if(EventID != ID)//读取完一个 Event 的所有 entry 之后，进行操作
      {
      if(ID != -1)
        {
          if(ene > 0) h->Fill(ene);
        }

      ID = EventID;
      ene = 0;
      }

      //  一个 Event 的 entry 能量沉积相加
      ene += EDep;
   }

   TCanvas *c1 = new TCanvas("c1");
   h->Draw();

}

ana5

In [ ]:

一些技巧¶

for 循环快速构建探测器阵列

// wuDetectorConstruction.hh
  G4LogicalVolume* logicCrystal[10000];
  G4VPhysicalVolume* physCrystal[10000];
  G4LogicalVolume* logicPMT[10000];
  G4VPhysicalVolume* physPMT[10000];

// wuDetectorConstruction.cc
  char tempname[32];

  G4Tubs* solidCrystal =
    new G4Tubs("Crystal",
           0,//内半径
           crystalR,//外半径
           0.5*crystalZ,//Z轴方向的半长度
           0*degree,//圆周起始位置弧度值
           360*degree);//该实体的圆心角弧度值


  for (int i = 0; i < 10; ++i)
    for (int j = 0; j < 10; ++j)
    {
      int num = 10*i+j;
      sprintf(tempname,"Crystal%04d",num);

      logicCrystal[num] =
    new G4LogicalVolume(solidCrystal,            //its solid
                crystal_mat,             //its material
                tempname);         //its name

      physCrystal[num] =
    new G4PVPlacement(0,                       //no rotation set 0
              G4ThreeVector((31.86-7.08*i)*cm, (31.86-7.08*j)*cm, 20.0*cm+0.5*crystalZ), //at (0,0,0)
              logicCrystal[num],               //its logical volume
              tempname,              //its name
              logicWorld,              //its mother  volume
              false,                   //no boolean operation
              0,                       //copy number
              checkOverlaps);          //overlaps checking

    }

ana6

数据处理关键代码

#define N_CRYSTALS 10000
std::stringstream s2i;//sstream cstring
double energy[N_CRYSTALS];
short crystal;
int id;
int tmpi;
double de;

// 循环内
  if(id != EventID)
    {
      if(id != -1)
    {
      //fill
      for (int i = 0; i < N_CRYSTALS; ++i)
        {
          if(energy[i] > 0)
        {
          crystal = i;
          de = energy[i];
          treewrite->Fill();//loop
        }
        }
    }

      id = EventID;

      for (int i = 0; i < N_CRYSTALS; ++i)
    {
      energy[i] = 0;
    }
    }


  // 获得探测器编号
  std::string test(VolNamePre+7);//"Crystal%04d
  // std::cout<<test<<std::endl;
  s2i.clear();//重复使用前一定要清空
  s2i<<test;
  s2i>>tmpi;

  energy[tmpi] += EDep;

In [1]:

!jupyter nbconvert dataana.ipynb --to html

[NbConvertApp] Converting notebook dataana.ipynb to html
[NbConvertApp] Writing 589778 bytes to dataana.html