核物理实验数据处理 Experimental Data Analysis in Nuclear Physics

5.1 γ-γ 符合矩阵与开窗谱

本节从上一节生成的 eurica_addback.root 文件出发，构造 γ-γ 符合矩阵。

上一节已经完成事件重构、时间修正和 addback。本节不再重复这些处理过程，而是直接使用 addback 后的 γ hit，讨论如何从事件中的多条 γ 构造二维符合矩阵，并利用开窗得到开窗谱。

从探测阵列到 γ-γ 符合矩阵

在多探测器 γ 探测阵列中，多个探测器从不同方向同时记录 γ 的能量和时间。一次核退激过程可能连续发射多条 γ。如果这些 γ 在同一个重构事件中被探测到，就可能形成 γ-γ 符合。

本节将各探测器的能量换算到共同的 keV 标尺，合并不同探测器对的能量关联，不保留角度坐标。对强度作定量解释时，仍需考虑各探测器对的效率及角关联；合并并不意味着所有探测器响应相同。

假设要检查 $121\ {\rm keV}$ γ 是否与 $244\ {\rm keV}$ γ 符合。直接做法是：先在某个探测器或某一组探测器中寻找 $121\ {\rm keV}$ γ，然后在同一个事件中遍历其他探测器，检查是否同时存在 $244\ {\rm keV}$ γ。对于一对指定的 γ，这种检查是可以完成的。

但是在实际 γ 能谱中，往往存在许多 γ 峰。我们并不总是事先知道哪些 γ 彼此符合。如果对每一条 γ 都逐一设置条件，再遍历所有探测器和所有事件去寻找它的符合 γ，分析过程会非常繁琐，也不利于整体观察能级级联关系。

为了避免上述繁琐过程，γ 谱学中通常先构造 γ-γ 符合矩阵。其基本做法是：对每一个重构事件，将其中所有 γ hit 两两配对，并把每一对 γ 的能量填入二维直方图：

$$ E_{\gamma 1}, \qquad E_{\gamma 2}. $$

这样，所有可能的 γ-γ 符合关系都被预先整理到同一个二维矩阵中。之后，如果想检查某一条 γ 与哪些 γ 符合，只需要在矩阵的一个能量轴上设置开窗，再将对应区域投影到另一个能量轴上。得到的开窗谱就显示了与该 γ 符合的其他 γ。

为了构造对称 γ-γ 符合矩阵，通常对同一个 γ-γ pair 同时填入

$$ (E_i,E_j), \qquad (E_j,E_i). $$

对称填充来自交换同一 pair 的两个坐标，不要求各探测器效率和分辨率完全相同。两轴采用相同的刻度单位与选择条件，就可把它们作为无序能量对来显示。镜像的两个计数来自同一 pair，不能当成两次独立测量。

---

γ-γ 符合矩阵与 开窗谱（gated spectrum）

γ-γ 符合矩阵的构造可以先从一个简化的事件列表理解。

图中左侧给出了几个事件。每个事件包含若干条 γ hit。例如，一个事件中包含 $422\ {\rm keV}$ 和 $885\ {\rm keV}$ 两条 γ；另一个事件包含 $277\ {\rm keV}$、$615\ {\rm keV}$ 和 $339\ {\rm keV}$ 三条 γ。

需要把其中所有 γ hit 两两配对。如果一个事件中有两条 γ，能量分别为 $E_i$ 和 $E_j$，则在二维直方图中填入这个能量对。对于对称 γ-γ 符合矩阵，通常同时填入 $$ (E_i,E_j), \qquad (E_j,E_i). $$ 例如，包含 $422\ {\rm keV}$ 和 $885\ {\rm keV}$ 的事件会贡献两个矩阵点：

$$ (422,885), \qquad (885,422). $$

如果一个事件中有三条 γ：

$$ E_1,\quad E_2,\quad E_3, $$

则会贡献六个有序 pair：

$$ (E_1,E_2),\quad (E_2,E_1), $$

$$ (E_1,E_3),\quad (E_3,E_1), $$

$$ (E_2,E_3),\quad (E_3,E_2). $$

这就是 γ-γ 符合矩阵关于对角线对称的原因。

矩阵中的一个点表示：在同一个事件中观测到了对应能量的两条 γ。不过，这个点并不能自动解释为真实的核能级级联关系。偶然符合、Compton 本底、峰下本底等都可能对矩阵产生贡献。因此，矩阵中的结构还需要通过时间条件、开窗谱和本底扣除进一步判断。

---

总投影谱（Total projection spectrum）

将 γ-γ 符合矩阵沿一个能量轴投影，得到的能谱称为总投影谱。它不能简单等同于普通的一维能谱。

在普通一维能谱中，一条 γ 被探测到一次，就计数一次。而在 γ-γ 符合矩阵的投影中，一条 γ 会按照与其符合的 γ 数目重复计数。因此，如果一个事件中的 γ 重数为 $M$，那么其中每条 γ 在矩阵投影中会贡献 $M-1$ 次计数。

数学上，沿 $E_x$ 轴的投影可写为

$$ P(E_x)=\sum_{E_y} M(E_x,E_y), $$

其中 $M(E_x,E_y)$ 是 γ-γ 符合矩阵。

因此，$P(E_x)$ 表示的是与 $E_x$ 发生符合的 pair 数目，而不是 $E_x$ 在原始一维 γ 能谱中出现的次数。

总投影谱的作用是帮助选择可能的开窗 γ 峰。那些经常与其他 γ 发生符合的跃迁，在总投影谱中会更加明显。

---

开窗（gate） 与 开窗谱（gated spectrum）

开窗是在 γ-γ 符合矩阵的一个能量轴上选择某一条 γ 峰。例如，在 $615\ {\rm keV}$ 上开窗，就是在一个能量轴上选择 $615\ {\rm keV}$ 附近的计数，并把这部分计数投影到另一个能量轴上。得到的能谱称为 $615\ {\rm keV}$ 开窗谱，表示与 $615\ {\rm keV}$ γ 符合的其他 γ 的能量分布。

在理想情况下，可以把开窗理解为选择 $E_x=E_g$。但在实际能谱中，γ 峰具有有限能量分辨率，峰面积分布在若干个 channel 中。因此，实际开窗必须选取一个有限能量范围。对于中心能量为 $E_g$、半宽为 $\Delta E$ 的开窗，开窗谱可以写成

$$ G(E_y \mid E_g)=\sum_{E_x \in [E_g-\Delta E,\ E_g+\Delta E]} M(E_x,E_y). $$

这个能谱回答的是一个条件问题：如果一个事件中观测到能量位于 $E_g$ 附近的 γ，那么另一个与它符合的 γ 的能量分布是什么？

开窗宽度 $\Delta E$ 一般根据总投影谱中该 γ 峰的实际峰形和能量分辨率选择。

实际选择时通常遵循以下原则：

1. 开窗应覆盖该 γ 峰的主要全能峰区域；
2. 开窗边界应尽量避开相邻 γ 峰；
3. 若峰附近本底较高，开窗不宜过宽；
4. 若后续要做峰旁窗本底扣除，峰旁窗的宽度应与 peak 窗保持一致或明确归一化；
5. 不同能量处的探测器分辨率不同，因此高能 γ 峰可能需要比低能 γ 峰稍宽的开窗。

在示意图中，对 $615\ {\rm keV}$ 开窗后，开窗谱中出现 $186\ {\rm keV}$、$277\ {\rm keV}$、$339\ {\rm keV}$ 和 $761\ {\rm keV}$ 的峰。这说明在事件列表中，$615\ {\rm keV}$ 与这些 γ 出现在同一个事件中。由于 $615\ {\rm keV}$ 在这个例子中同时出现在多个级联分支中，所以 $615\ {\rm keV}$ 开窗谱会同时包含来自不同级联分支的 γ。

相比之下，对 $186\ {\rm keV}$ 开窗后，只得到 $615\ {\rm keV}$ 和 $761\ {\rm keV}$。这个开窗更干净，因为 $186\ {\rm keV}$ 在示意图中只出现在一个级联分支中。

这说明开窗谱的质量不仅取决于开窗宽度，也取决于所选 γ 峰本身是否干净。一个孤立、强度适中、峰下本底较低的 γ 峰，通常比一个被多个级联分支共享或被邻近峰污染的 γ 峰更适合作为开窗对象。

%jsroot on

TCanvas *c1=new TCanvas("c1","c1");
TFile *fin = new TFile("eurica_addback.root");
TTree *tree = (TTree*)fin->Get("tree");
TH1::SetDefaultSumw2(); // 本底相减之前保留各 bin 的计数方差

addback能谱

.L peaks.C

tree->Draw("ae>>hae(1500,0,1500)", "", "hist");
peaks("hae");
c1->SetLogy();
c1->Draw();

从 addback 后 hit 构造 γ-γ pair

对每个事件，将所有 addback 后 γ hit 两两配对。同一个 hit 不能与自身配对。下面的 macro 生成一个新的 tree，其中保存每个 γ-γ pair 的两个能量、两个时间和两个探测器编号。

本例要求两条 addback hit 的能量均不低于 30 keV，代表晶体编号不同，避免同一晶体重复脉冲的组合。循环保留两个坐标顺序；若全部 pair 都通过选择，一个 M 重事件贡献 M(M−1) 行，而不是 M 行。

%%cpp -d
void make_ggpair()
{
  TFile *fin = new TFile("eurica_addback.root");
  TTree *tree = (TTree*)fin->Get("tree");

  const int maxhit = 200;

  int ahit;
  int aid[maxhit];
  double ae[maxhit];
  double at[maxhit];
  if (tree->GetMaximum("ahit")>maxhit) throw std::runtime_error("Increase maxhit before reading");
  tree->SetBranchAddress("ahit",&ahit);
  tree->SetBranchAddress("aid",aid);
  tree->SetBranchAddress("ae",ae);
  tree->SetBranchAddress("at",at);

  TFile *fout = new TFile("eurica_ggpair.root","recreate");
  TTree *tout = new TTree("tree","gamma-gamma pairs");

  double ex, ey;
  double tx, ty;
  int idx, idy;

  tout->Branch("ex",&ex,"ex/D");
  tout->Branch("ey",&ey,"ey/D");
  tout->Branch("tx",&tx,"tx/D");
  tout->Branch("ty",&ty,"ty/D");
  tout->Branch("idx",&idx,"idx/I");
  tout->Branch("idy",&idy,"idy/I");

  Long64_t nentries = tree->GetEntries();
  for (Long64_t entry = 0; entry < nentries; entry++) {
    tree->GetEntry(entry);
    for (int i = 0; i < ahit; i++) {
      if (ae[i]<30) continue; // i 的条件每轮外循环只判断一次
      for (int j = 0; j < ahit; j++) {
        if (i == j || aid[i] == aid[j]) continue;
        if (ae[j]<30) continue;

        ex  = ae[i];
        ey  = ae[j];
        tx  = at[i];
        ty  = at[j];
        idx = aid[i];
        idy = aid[j];

        tout->Fill();
      }
    }
  }

  std::cout << "Ordered pair rows = " << tout->GetEntries() << std::endl;
  fout->Write();
  fout->Close();
}

make_ggpair();

Ordered pair rows = 7135608

TFile *fin = new TFile("eurica_ggpair.root");
TTree *tree = (TTree*)fin->Get("tree");

时间差与 prompt 条件

γ-γ pair tree 中同时包含 prompt 符合和 random 符合。为了区分它们，首先检查 pair 中两条 γ 的时间差：

$$ \Delta t = t_x - t_y. $$

c1->Clear();
gPad->SetLogy(0);
tree->Draw("tx-ty>>hdt(400,-1000,1000)", "", "hist");
c1->Draw();

根据上一节中时间修正后得到的时间分辨，本节采用

$$ |\Delta t| < 180\ {\rm ns} $$

作为 prompt 条件, 同时也可以定义一个远离 prompt 峰的 random 时间窗。 prompt 条件选择的是时间上接近的 γ-γ pair，因此它们是真实级联符合的候选。random 时间窗选择的是远离 prompt 峰的 pair，主要代表偶然符合。

prompt 与 random 矩阵在这里用于对照结构，尚未相减。若随机时间差密度与接受度在两个窗口中相同，才可按总宽度使用 $\alpha=360/5200$；本例 3 μs 分组使大时间差的 pair 接受度改变，定量相减应先确定有效暴露比例。避免把 prompt 尾部或真实 delayed 级联放入 random 窗。

TCut prompt_cut = "abs(tx-ty)<180";
TCut random_cut = "abs(tx-ty)>400 && abs(tx-ty)<3000";

TCanvas *c2 = new TCanvas("c2","random gamma-gamma matrix",1000,500);
c2->Divide(2,1);
c2->cd(1);
tree->Draw("ey:ex>>hgg_prompt(1500,0,1500,1500,0,1500)",prompt_cut,"colz");
gPad->SetLogz();
c2->cd(2);
tree->Draw("ey:ex>>hgg_random(1500,0,1500,1500,0,1500)",random_cut,"colz");
gPad->SetLogz();
c2->Draw();

在 prompt 矩阵中，真实 γ 级联关系会表现为对应 γ 能量交叉处的局域符合点。random 矩阵不代表真实核能级级联。但强 γ 峰仍然可能形成明显的横线、竖线或交叉结构，因为强峰有较高概率与其他无关 γ 发生偶然符合。

因此，prompt γ-γ 符合矩阵中的一个点不能自动解释为真实级联。它还需要结合时间条件、random 本底、开窗谱以及已知能级方案进行判断。

总投影谱

总投影谱中152Eu和133Ba的一系列$\gamma$ 峰清晰可见。

c1->Clear();
TH2 *hgg = (TH2*)gDirectory->Get("hgg_prompt");
TH1D *hproj = hgg->ProjectionX("hproj",1,hgg->GetNbinsY());
peaks("hproj");
c1->Draw();

对 $121\ {\rm keV}$ γ 开窗

下面将开窗的概念应用到真实 γ 峰。这里以 $121\ {\rm keV}$ γ 为例。

开窗宽度 $\Delta E$ 的选择需要在“保留峰面积”和“减少本底污染”之间折中。开窗太窄，会丢失全能峰两侧的有效计数，使开窗谱统计降低；开窗太宽，则会引入更多 Compton 本底、峰下本底以及邻近峰污染。 实际选择时通常遵循以下原则：

1. 开窗应覆盖该 γ 峰的主要全能峰区域；
2. 开窗边界应尽量避开相邻 γ 峰；
3. 不同能量处的探测器分辨率不同，因此高能 γ 峰可能需要比低能 γ 峰稍宽的开窗。

TCut gate121 = "abs(ex-121)<2";

这表示在 $E_x$ 轴上选择

119 keV<$E_x$<123 keV.

然后将这个区域投影到 $E_y$ 轴，得到 $121\ {\rm keV}$ 开窗谱：

c1->Clear();
tree->Draw("ey>>hgate121(1000,0,1000)",prompt_cut && gate121,"hist");
peaks("hgate121");
c1->Draw();

如果某个 γ 峰在 $121\ {\rm keV}$ 开窗谱中明显增强，说明这条 γ 可能与 $121\ {\rm keV}$ γ 处于同一级联过程中。

但是，仅根据开窗谱判断符合关系是不够的。$121\ {\rm keV}$ 开窗中不仅包含全能峰的计数，也包含 peak 下方连续本底的计数。这些本底计数同样会产生自己的符合谱。因此，需要回到二维 γ-γ 符合矩阵中观察 peak 窗和峰旁窗的结构。

在二维矩阵中观察 $121\ {\rm keV}$ peak 窗和side-gate spectrum（左右侧本底窗）

为了理解 $121\ {\rm keV}$ 开窗谱中的成分，可以在 prompt γ-γ 符合矩阵上同时标出 peak 窗和峰side-gate spectrum。

TCanvas *c3 = new TCanvas("c3","121-keV peak and side gates",1000,700);
  tree->Draw("ey:ex>>hgg_gate_view(300,0,300,1000,0,1000)",prompt_cut,"colz");

  // peak gate: 121 +- 2 keV
  TLine *p1 = new TLine(119,0,119,1000);
  TLine *p2 = new TLine(123,0,123,1000);
  p1->SetLineColor(kRed);
  p2->SetLineColor(kRed);
  p1->SetLineWidth(2);
  p2->SetLineWidth(2);
  p1->Draw("same");
  p2->Draw("same");

  // left side gate: 116 +- 2 keV
  TLine *l1 = new TLine(114,0,114,1000);
  TLine *l2 = new TLine(118,0,118,1000);
  l1->SetLineColor(kOrange+1);
  l2->SetLineColor(kOrange+1);
  l1->Draw("same");
  l2->Draw("same");

  // right side gate: 126 +- 2 keV
  TLine *r1 = new TLine(124,0,124,1000);
  TLine *r2 = new TLine(128,0,128,1000);
  r1->SetLineColor(kOrange+1);
  r2->SetLineColor(kOrange+1);
  r1->Draw("same");
  r2->Draw("same");
  gPad->SetLogz();
  c3->Draw();

二维图的三个区域分别对应：

TCut gate121       = "abs(ex-121)<2";
TCut gate121_left  = "abs(ex-116)<2";
TCut gate121_right = "abs(ex-126)<2";

在二维矩阵中可以这样理解：

• 红色 peak 窗选择 $121\ {\rm keV}$ peak 附近的 γ-γ pair；
• left/right side gate 选择 $121\ {\rm keV}$ peak 两侧的连续本底区域；
• peak 窗的投影包含真实 $121\ {\rm keV}$ 符合，也包含 peak 下方本底的符合；
• side gate 的投影用于近似估计 peak 下方本底产生的符合谱。

一般地，本底扣除可写为

$$ G_{\rm sub}(E) = G_{\rm peak}(E) - \alpha G_{\rm side}(E), $$

其中 $G_{\rm peak}(E)$ 是 peak 开窗谱，$G_{\rm side}(E)$ 是side-gate spectrum估计的本底谱，$\alpha$ 是归一化因子。

这里三个窗等宽且中心对称，代码定义 $G_{\rm side}=(G_L+G_R)/2$，再做 $G_{\rm net}=G_P-(G_L+G_R)/2$。这一平均能估计局部线性本底在中心 peak 窗下的贡献；side gate 不能包含其他峰。若窗宽或位置不对称，应按宽度和本底斜率重新归一化。

一个实际做法是：在一维投影能谱中比较 peak 窗下方的本底计数与side-gate spectrum中的本底计数，用它们的比值确定 $\alpha$。

将这个本底谱从 peak 开窗谱中扣除后，剩余峰更接近于真正与 $121\ {\rm keV}$ 跃迁符合的 γ。

忽略不同窗之间的事件相关性时，逐 bin 有 $\sigma^2_{\rm net}=N_P+(N_L+N_R)/4$，ROOT 的 Sumw2、Scale 和 Add 会传播这些项。多重事件可同时贡献多个 pair，因此该式不是所有积分量的完整事件级误差；特别不能把矩阵两侧的镜像计数视作独立样本。

峰左右侧本底扣除

同时显示 peak 开窗谱、峰左右侧窗本底谱和扣除后的开窗谱：

tree->Draw("ey>>hgate121_peak(1000,0,1000)",prompt_cut && gate121, "goff");

tree->Draw("ey>>hgate121_left(1000,0,1000)",prompt_cut && gate121_left, "goff");

tree->Draw("ey>>hgate121_right(1000,0,1000)",prompt_cut && gate121_right, "goff");

  TH1 *hpeak  = (TH1*)gDirectory->Get("hgate121_peak");
  TH1 *hleft  = (TH1*)gDirectory->Get("hgate121_left");
  TH1 *hright = (TH1*)gDirectory->Get("hgate121_right");

  TH1 *hside = (TH1*)hleft->Clone("hside");
  hside->Add(hright);
  hside->Scale(0.5);

  TH1 *hsub = (TH1*)hpeak->Clone("hgate121_sub");
  hsub->Add(hside,-1.0);

  TCanvas *c4 = new TCanvas("c4","121-keV gate spectra",850,700);
  c4->Divide(1,2);
  c4->cd(1);
  hpeak->Draw("hist");
  hside->SetLineColor(kRed);
  hside->Draw("hist same");
  c4->cd(2);
  hsub->Draw("hist");
  c4->Draw();

总结

可以通过在峰两侧较平坦的本底区域设置side gate，估计并扣除 Compton 本底及峰下本底的贡献，从而得到更接近真实 γ-γ 符合关系的开窗谱。

这种side-gate method是二维对称 γ-γ 符合矩阵中最基本、最直观的本底扣除方法。在 γ 峰较孤立、峰两侧存在相对平坦本底区域时，推荐优先使用上述方法。

但是，当 γ 峰密度很高时，大多数 γ 峰两侧很难找到干净、平坦的本底区域。此时，这个方法不再可靠，也不便于逐个 γ 峰手动处理。

例如，在熔合蒸发反应中，复合核被布居到高自旋态后会通过级联 γ 退激，产生大量彼此接近的 γ 峰。此类 γ 能谱峰密度高、Compton 本底复杂，难以为每个 γ 峰选择合适的side gate。因此，需要寻找更系统、更方便的二维矩阵本底处理方法。