$\gamma-\gamma$ 对称矩阵的整体减本底方法：Radware approach¶

本节讨论一种二维对称矩阵的整体减本底方法，它的基本思路是：

  1. 由二维 $\gamma-\gamma$ 矩阵得到 total projection；
  2. 将 total projection 近似分解为 peak component 和 background component；
  3. 根据这两个一维成分构造二维 background matrix；
  4. 从原始二维矩阵中扣除 background matrix；
  5. 在扣本底后的矩阵上进行开窗(gate)，检查符合关系。

ROOT文件¶

本节继续使用 5.0 中得到的 addback 后数据文件：

eurica_addback.root

构建二维 $\gamma-\gamma$ 对称矩阵¶

二维 $\gamma-\gamma$ 矩阵的一个 bin $(E_i,E_j)$ 表示在同一个 prompt coincidence 中，两个 addback hit 的能量分别为 $E_i$ 和 $E_j$。

对于同一事件中的两个 addback hit，先计算时间差

$$ \Delta t = t_i - t_j. $$

若

$$ |\Delta t| < t_{\mathrm{prompt}}, $$

则认为这两个 hit 属于 prompt coincidence，并将其能量填入二维矩阵。

为了构造对称矩阵，每一个 $\gamma-\gamma$ pair 同时填入

$$ (E_i,E_j) $$

和

$$ (E_j,E_i). $$

这样得到的矩阵满足

$$ M_{ij}=M_{ji}. $$

下面的代码直接从 eurica_addback.root 构造二维对称矩阵。为了检查时间窗，同时保存一个对称填充的时间差谱 hdt_sym。这里 hdt_sym 只用于观察时间结构，不用于归一化。

In [2]:
%%cpp -d

#include <iostream>
#include <cmath>
#include <algorithm>

const int MAXHIT = 1024;

void make_symmetric_gg_matrix(const char *infile  = "eurica_addback.root",
                              const char *outfile = "eurica_ggmatrix.root",
                              int nbins = 1500,
                              double emax = 1530.0,
                              double emin = 30.0,
                              double dtPrompt = 180.0)
{
  TFile *fin = new TFile(infile);
  if (!fin || fin->IsZombie()) {
    std::cout << "cannot open " << infile << std::endl;
    return;
  }

  TTree *tree = (TTree*)fin->Get("tree");
  if (!tree) {
    std::cout << "cannot find tree" << std::endl;
    fin->Close();
    return;
  }

  int ahit;
  int aid[MAXHIT];
  double ae[MAXHIT];
  double at[MAXHIT];

  tree->SetBranchAddress("ahit", &ahit);
  tree->SetBranchAddress("aid", aid);
  tree->SetBranchAddress("ae", ae);
  tree->SetBranchAddress("at", at);

  TH1::SetDefaultSumw2(kTRUE);

  TH2D *hgg = new TH2D("hgg",
                       "#gamma-#gamma symmetric matrix;E_{x} (keV);E_{y} (keV)",
                       nbins, 0, emax,
                       nbins, 0, emax);

  TH1D *hdt_sym = new TH1D("hdt_sym",
                           "symmetric time difference spectrum;#Deltat (ns);counts",
                           800, -4000, 4000);

  Long64_t nentries = tree->GetEntries();

  for (Long64_t ientry = 0; ientry < nentries; ientry++) {

    tree->GetEntry(ientry);

    int nhit = std::min(ahit, MAXHIT);

    for (int i = 0; i < nhit; i++) {
      if (ae[i] < emin) continue;
      for (int j = i + 1; j < nhit; j++) {
        if (ae[j] < emin) continue;
        if (aid[i] == aid[j]) continue;

        double dt = at[i] - at[j];
        double adt = std::fabs(dt);

        hdt_sym->Fill(dt);
        hdt_sym->Fill(-dt);

        if (adt < dtPrompt) {

          hgg->Fill(ae[i], ae[j]);
          hgg->Fill(ae[j], ae[i]);

        }
      }
    }

    if (ientry % 100000 == 0) {
      std::cout << "processing " << ientry << " / " << nentries << "\r";
      std::cout.flush();
    }
  }

  std::cout << std::endl;

  TFile *fout = new TFile(outfile, "recreate");
  hdt_sym->Write();
  hgg->Write();
  fout->Close();

  fin->Close();

  std::cout << "symmetric gamma-gamma matrix saved to "
            << outfile << std::endl;
}
In [3]:
make_symmetric_gg_matrix();

processing 10300000 / 10370481
symmetric gamma-gamma matrix saved to eurica_ggmatrix.root

检查时间差分布和二维矩阵¶

In [5]:
TCanvas *c1 = new TCanvas("c1", "c1", 800, 700);

TFile *fg = new TFile("eurica_ggmatrix.root");

TH1D *hdt_sym = (TH1D*)fg->Get("hdt_sym");
TH2D *hgg = (TH2D*)fg->Get("hgg");

二维 $\gamma-\gamma$ 对称矩阵¶

In [7]:
%jsroot off
hgg->Draw("colz");
gPad->SetLogz();
c1->Draw();
No description has been provided for this image

Total projection 与一维本底分解¶

设二维对称矩阵为 $M_{ij}$。沿一个方向投影得到 total projection：

$$ P_i=\sum_j M_{ij}. $$

其中 $P_i$ 表示二维矩阵中与第 $i$ 个能量 bin 相关的总符合计数。

In [9]:
.L peaks.C
In [10]:
TCanvas *c2 = new TCanvas("c2", "c2", 1000, 600);
TH1D *hproj = hgg->ProjectionX("hproj");

hproj->SetTitle("total projection spectrum;E_{#gamma} (keV);counts");
peaks("hproj");
gPad->SetLogy();
c2->Draw();
No description has been provided for this image

Total projection 中既包含 full-energy peaks，也包含 Compton continuum 和 quasi-continuum。Radware approach 首先将 total projection 近似分解为

$$ P_i=p_i+b_i. $$

其中：

  • $p_i$ 是第 $i$ 个 bin 中的 peak component；
  • $b_i$ 是第 $i$ 个 bin 中的 background component。

这里的 $p_i$ 和 $b_i$ 不是对每一个物理峰逐峰拟合得到的结果，而是从 total projection 中估计出的平均 peak / background 成分。

使用 TSpectrum::Background() 估计一维本底。

其中：

  • hproj 对应 $P_i$；
  • hbg1d 对应 $b_i$；
  • hpeak1d 对应 $p_i$。

TSpectrum::Background() 的参数需要通过图形检查。若本底线明显穿过强峰，说明本底估计过高；若本底线明显低于连续平台，说明本底估计过低。Radware approach 的效果首先取决于这一维分解是否合理。

In [12]:
TSpectrum *sp = new TSpectrum(500);

TH1D *hbg1d = (TH1D*)sp->Background(hproj, 6, "nosmoothing");
hbg1d->SetName("hbg1d");
hbg1d->SetTitle("1D background from total projection");

TH1D *hpeak1d = (TH1D*)hproj->Clone("hpeak1d");
hpeak1d->SetTitle("peak component of total projection");
hpeak1d->Add(hbg1d, -1.0);

TCanvas *c3 = new TCanvas("c3", "c3", 1100, 400);
c3->Divide(2,1);
c3->cd(1);
hproj->Draw();
hbg1d->SetLineColor(kRed);
hbg1d->Draw("same");
gPad->SetLogy();
c3->cd(2);
hpeak1d->SetLineColor(kBlack);
hpeak1d->Draw();
gPad->SetLogy();
c3->Draw();
No description has been provided for this image

Background matrix¶

在二维 $\gamma-\gamma$ 矩阵中，一个点 $(e_x,e_y)$ 表示在同一个符合事件中，两个探测器或探测单元分别测到了能量 $e_x$ 和 $e_y$。

设真实入射到两个探测器的 $\gamma$ 射线能量分别为 $E_{\gamma 1}$ 和 $E_{\gamma 2}$。由于探测器存在有限能量分辨、Compton 散射、阈值效应和 addback 处理，测得能量不一定等于入射 $\gamma$ 射线的全能量。

对于一个入射能量为 $E_\gamma$ 的 $\gamma$ 射线，探测器测得能量 $e$ 的条件概率写作

$$ p(e\mid E_\gamma). $$

该条件概率可以分成全能峰响应和非全能峰响应：

$$ p(e\mid E_\gamma) = p_{\mathrm{peak}}(e\mid E_\gamma) + p_{\mathrm{bg}}(e\mid E_\gamma). $$

其中 $p_{\mathrm{peak}}(e\mid E_\gamma)$ 描述全能峰，$p_{\mathrm{bg}}(e\mid E_\gamma)$ 描述 Compton continuum、本底散射和其他非全能峰响应。

如果两个探测器的响应近似独立，则在给定 $E_{\gamma 1}$ 和 $E_{\gamma 2}$ 时，同时测得 $(e_x,e_y)$ 的条件概率可以写成

$$ p(e_x,e_y\mid E_{\gamma 1},E_{\gamma 2}) \simeq p_x(e_x\mid E_{\gamma 1}) p_y(e_y\mid E_{\gamma 2}). $$

这一步是后面出现乘积形式的物理基础。这里相乘的是两个条件概率，而不是两个计数。

将两个方向的响应分解为 peak 和 background：

$$ p_x(e_x\mid E_{\gamma 1}) = p_{x,\mathrm{peak}}(e_x\mid E_{\gamma 1}) + p_{x,\mathrm{bg}}(e_x\mid E_{\gamma 1}), $$

$$ p_y(e_y\mid E_{\gamma 2}) = p_{y,\mathrm{peak}}(e_y\mid E_{\gamma 2}) + p_{y,\mathrm{bg}}(e_y\mid E_{\gamma 2}). $$

二维响应展开为

$$ \begin{aligned} p(e_x,e_y\mid E_{\gamma 1},E_{\gamma 2}) =&\ p_{x,\mathrm{peak}}(e_x\mid E_{\gamma 1}) p_{y,\mathrm{peak}}(e_y\mid E_{\gamma 2}) \\ &+ p_{x,\mathrm{bg}}(e_x\mid E_{\gamma 1}) p_{y,\mathrm{peak}}(e_y\mid E_{\gamma 2}) \\ &+ p_{x,\mathrm{peak}}(e_x\mid E_{\gamma 1}) p_{y,\mathrm{bg}}(e_y\mid E_{\gamma 2}) \\ &+ p_{x,\mathrm{bg}}(e_x\mid E_{\gamma 1}) p_{y,\mathrm{bg}}(e_y\mid E_{\gamma 2}). \end{aligned} $$

第一项是两个方向都测到全能峰的 peak-peak coincidence，是建立能级纲图时希望保留的部分。后三项分别对应 background-peak、peak-background 和 background-background，是二维矩阵中的连续本底贡献。

Radware 本底矩阵¶

  • D.C. Radford, Nucl. Instr. Meth. A361(1995)306

实际数据中，我们不知道每一个 bin $(i,j)$ 对应的真实入射能量 $E_{\gamma 1}$ 和 $E_{\gamma 2}$，也不能逐事件判断测量能量来自全能峰还是 Compton continuum。因此，Radware approach 用 total projection 来估计平均的 peak / background 概率。

Total projection 的总计数为

$$ T=\sum_i P_i. $$

第 $i$ 个能量 bin 在矩阵某一方向上出现的平均概率为

$$ q_i=\frac{P_i}{T}. $$

对应的 peak 和 background 平均概率为

$$ q_i^{(p)}=\frac{p_i}{T}, $$

$$ q_i^{(b)}=\frac{b_i}{T}. $$

并且

$$ q_i=q_i^{(p)}+q_i^{(b)}. $$

这里的 $q_i^{(p)}$ 和 $q_i^{(b)}$ 可以看作对 $p_{\mathrm{peak}}(e_i\mid E_\gamma)$ 和 $p_{\mathrm{bg}}(e_i\mid E_\gamma)$ 在所有真实 $\gamma$ 能量、所有探测器和所有符合关系上的平均。

在平均意义下，二维本底中三类贡献的期望计数为

$$ Tq_i^{(b)}q_j^{(p)}, $$

$$ Tq_i^{(p)}q_j^{(b)}, $$

$$ Tq_i^{(b)}q_j^{(b)}. $$

代入 $q_i^{(p)}=p_i/T$ 和 $q_i^{(b)}=b_i/T$，得到

$$ Tq_i^{(b)}q_j^{(p)} = \frac{b_ip_j}{T}, $$

$$ Tq_i^{(p)}q_j^{(b)} = \frac{p_ib_j}{T}, $$

$$ Tq_i^{(b)}q_j^{(b)} = \frac{b_ib_j}{T}. $$

因此，二维 background matrix 为

$$ B_{ij} = \frac{1}{T} \left( b_ip_j+p_ib_j+b_ib_j \right). $$

其中：

  • $b_ip_j/T$ 表示 $i$ 方向来自 background，$j$ 方向来自 peak；
  • $p_ib_j/T$ 表示 $i$ 方向来自 peak，$j$ 方向来自 background；
  • $b_ib_j/T$ 表示两个方向都来自 background。

由于

$$ P_i=p_i+b_i, $$

$$ P_j=p_j+b_j, $$

所以

$$ P_iP_j-p_ip_j = (p_i+b_i)(p_j+b_j)-p_ip_j = b_ip_j+p_ib_j+b_ib_j. $$

因此也可写成

$$ B_{ij} = \frac{P_iP_j-p_ip_j}{T}. $$

最后，从原始二维矩阵中扣除该本底矩阵：

$$ M_{ij}^{\mathrm{net}} = M_{ij}-B_{ij}. $$

这个公式的物理含义是：尽量保留两个方向都来自全能峰的 peak-peak coincidence，扣除至少有一个方向来自连续本底的平均贡献。

上述推导的近似条件¶

上述推导包含三个主要近似。

探测器响应独立近似¶

给定真实入射能量 $E_{\gamma 1}$ 和 $E_{\gamma 2}$ 后，两个方向的测量过程近似独立：

$$ p(e_x,e_y\mid E_{\gamma 1},E_{\gamma 2}) \simeq p_x(e_x\mid E_{\gamma 1}) p_y(e_y\mid E_{\gamma 2}). $$

如果存在明显 cross-talk、错误 addback、pile-up 或 summing，这个近似会变差。

用 total projection 代替真实响应函数¶

真实的 $p(e\mid E_\gamma)$ 依赖入射 $\gamma$ 能量、探测器几何、晶体材料、阈值和 addback 方法。Radware approach 不逐个 $E_\gamma$ 建立响应函数，而是用 $p_i$ 和 $b_i$ 估计平均 peak / background 概率。

因此，$p_i/T$ 和 $b_i/T$ 不是严格的单能 $\gamma$ 响应概率，而是整个数据集平均后的概率。

二维本底因子化近似¶

公式

$$ B_{ij} = \frac{P_iP_j-p_ip_j}{T} $$

假设二维本底的平均能量分布可以由两个方向的一维平均概率相乘得到。这个近似可以描述主要的 Compton continuum 和 quasi-continuum，但不能精确描述所有二维相关结构，例如多重散射、summing、强 cascade 分支和探测器间串扰。

  • 上述方法可以推广到 3-fold，4-fold的矩阵

其他方法：¶

  • G. Palameta and J.C. Waddington, Nucl. Instr. Meth. A234(1985)476
  • TSpectrum2的Background()方法
    • https://root.cern.ch/doc/master/classTSpectrum2.html

Radware 方法的 ROOT 实现¶

下面根据

$$ B_{ij} = \frac{P_iP_j-p_ip_j}{T} $$

构造 background matrix，并从原始矩阵中扣除。

In [15]:
%%cpp -d

#include <iostream>
#include <cmath>
#include <algorithm>

void make_radware_matrix(const char *infile  = "eurica_ggmatrix.root",
                         const char *matrixName = "hgg",
                         const char *outfile = "eurica_gg_radware.root",
                         int bgIter = 6)
{
  TFile *fin = new TFile(infile);
  if (!fin || fin->IsZombie()) {
    std::cout << "cannot open " << infile << std::endl;
    return;
  }

  TH2D *hgg = (TH2D*)fin->Get(matrixName);
  if (!hgg) {
    std::cout << "cannot find matrix " << matrixName << std::endl;
    fin->Close();
    return;
  }

  TH1::SetDefaultSumw2(kTRUE);

  TH1D *hproj = hgg->ProjectionX("hproj");
  hproj->SetTitle("total projection spectrum;E_{#gamma} (keV);counts");

  TSpectrum *sp = new TSpectrum(500);

  TH1D *hbg1d = (TH1D*)sp->Background(hproj, bgIter, "nosmoothing");
  hbg1d->SetName("hbg1d");
  hbg1d->SetTitle("1D background from total projection");

  TH1D *hpeak1d = (TH1D*)hproj->Clone("hpeak1d");
  hpeak1d->SetTitle("1D peak component from total projection");
  hpeak1d->Add(hbg1d, -1.0);

  int nx = hgg->GetNbinsX();
  int ny = hgg->GetNbinsY();

  TH2D *hggb = (TH2D*)hgg->Clone("hggb_radware");
  hggb->Reset();
  hggb->SetTitle("Radware background matrix;E_{x} (keV);E_{y} (keV)");

  TH2D *hggmat = (TH2D*)hgg->Clone("hgg_radware");
  hggmat->SetTitle("background-subtracted #gamma-#gamma matrix;E_{x} (keV);E_{y} (keV)");

  double T = hproj->Integral();

  for (int ix = 1; ix <= nx; ix++) {

    double Pi = hproj->GetBinContent(ix);
    double pi = hpeak1d->GetBinContent(ix);

    if (Pi < 0) Pi = 0;
    if (pi < 0) pi = 0;
    if (pi > Pi) pi = Pi;

    for (int iy = 1; iy <= ny; iy++) {

      double Pj = hproj->GetBinContent(iy);
      double pj = hpeak1d->GetBinContent(iy);

      if (Pj < 0) Pj = 0;
      if (pj < 0) pj = 0;
      if (pj > Pj) pj = Pj;

      double Bij = 0.0;

      if (T > 0) {
        Bij = (Pi * Pj - pi * pj) / T;
      }

      if (Bij < 0) Bij = 0;

      hggb->SetBinContent(ix, iy, Bij);
      hggb->SetBinError(ix, iy, std::sqrt(Bij));
    }
  }

  hggmat->Add(hggb, -1.0);

  TFile *fout = new TFile(outfile, "recreate");

  hgg->Write("hgg_input");
  hproj->Write();
  hbg1d->Write();
  hpeak1d->Write();
  hggb->Write();
  hggmat->Write();

  fout->Close();
  fin->Close();

  std::cout << "Radware background-subtracted matrix saved to "
            << outfile << std::endl;
}
In [16]:
make_radware_matrix();

Radware background-subtracted matrix saved to eurica_gg_radware.root

检查 background matrix 和扣本底矩阵¶

In [18]:
TFile *fr = new TFile("eurica_gg_radware.root");

TH2D *hgg_input = (TH2D*)fr->Get("hgg_input");
TH2D *hggb_radware = (TH2D*)fr->Get("hggb_radware");
TH2D *hgg_radware = (TH2D*)fr->Get("hgg_radware");
In [19]:
c1->Clear();
hggb_radware->Draw("colz");
gPad->SetLogy(0);
gPad->SetLogz();
c1->Draw();
No description has been provided for this image

Background matrix 应主要反映连续本底区域的整体强度分布，而不应产生尖锐孤立的 peak-peak 结构。

再看扣本底后的二维矩阵：

In [21]:
hgg_radware->Draw("colz");
gPad->SetLogz();
c1->Draw();
No description has been provided for this image

减本底后应观察到：

  • 连续平台降低；
  • 强的真实符合点保留；
  • 由 Compton continuum 产生的横向或纵向本底减弱；
  • 统计较弱的位置可能出现接近零或负计数的波动。

扣本底后的矩阵不能简单理解为“只剩真实符合”。它只是降低了连续本底的平均贡献。

比较扣本底前后的投影谱：

In [23]:
TCanvas *c2 = new TCanvas("c2","c2",1000,600);

Warning in <TCanvas::Constructor>: Deleting canvas with same name: c2
In [24]:
TH1D *hproj_raw = hgg_input->ProjectionX("hproj_raw");
TH1D *hproj_net = hgg_radware->ProjectionX("hproj_net");

hproj_raw->SetLineColor(kBlack);
hproj_net->SetLineColor(kRed);

hproj_raw->Draw("hist");
hproj_net->Draw("hist same");

gPad->SetLogy();
c2->Draw();
No description has been provided for this image

扣本底后，连续平台应降低，主要的全能峰应保留。若主要峰明显被削弱，说明一维本底估计可能过高；若连续平台几乎没有变化，说明本底估计可能过低。

对于扣本底后的 bin，

$$ N_{\mathrm{net}} = N_{\mathrm{raw}}-N_{\mathrm{bkg}}. $$

其统计不确定度近似为

$$ \sigma_{N_{\mathrm{net}}} \simeq \sqrt{ \sigma_{N_{\mathrm{raw}}}^2 + \sigma_{N_{\mathrm{bkg}}}^2 }. $$

因此，弱峰在扣本底后的相对不确定度可能变大，不能只凭二维图上的颜色判断其为真实符合。

只有尽量减少本底贡献，才能确认弱关联的gamma之间的符合关系

途径：每个探测器加入anti-compton shield；通过加入其它探测器的符合，压制本底...

In [26]:
void zoom2D(TH2D *h,
            double xmin, double xmax,
            double ymin, double ymax,
            const char *opt = "colz")
{
  h->GetXaxis()->SetRangeUser(xmin, xmax);
  h->GetYaxis()->SetRangeUser(ymin, ymax);
  h->Draw(opt);
  //gPad->SetLogz();
  gPad->Update();
}
In [27]:
TCanvas *c5 = new TCanvas("c5","c5",900,300);
c5->Divide(3,1);
c5->cd(1);
zoom2D(hgg_input, 50, 200, 50, 200);

c5->cd(2);
zoom2D(hggb_radware, 50, 200, 50, 200);

c5->cd(3);
zoom2D(hgg_radware, 50, 200, 50, 200);
c5->Draw();
No description has been provided for this image

减本底的矩阵中，来源于133Ba的(80,80)的符合清晰可见，而其他有偶然本底导致的交叉点(80,121),(121,80)消失。

在扣本底矩阵上进行开窗(gate)¶

定义一个通用 gate 函数。在二维对称矩阵中，如果在 $x$ 轴某个能量范围开 gate，则投影到 $y$ 轴。

In [30]:
TH1D* gate_gg(TH2D *hgg,
              double gateE,
              double gateWidth = 2.0,
              double xmax = 1300.0)
{
  static int igate = 0;
  igate++;

  hgg->GetXaxis()->SetRangeUser(0, xmax);
  hgg->GetYaxis()->SetRangeUser(0, xmax);

  int bin1 = hgg->GetXaxis()->FindBin(gateE - gateWidth);
  int bin2 = hgg->GetXaxis()->FindBin(gateE + gateWidth);

  TString hname;
  hname.Form("gate_%.1f_keV_%04d", gateE, igate);

  TH1D *hg = hgg->ProjectionY(hname, bin1, bin2);

  hg->SetTitle(Form("gated on %.1f keV;E_{#gamma} (keV);counts",gateE));

  hg->GetXaxis()->SetRangeUser(0, xmax);

  peaks(hname);

  return hg;
}

以 $121.8\,\mathrm{keV}$ 为例。先在原始矩阵中开 gate：

In [32]:
c2->Clear();
TH1D *h121_raw = gate_gg(hgg_input, 121.8, 2.0);
gPad->SetLogy(0);
c2->Draw();
No description has been provided for this image

再在 Radware 扣本底后的矩阵中开 gate：

In [34]:
c2->Clear();
TH1D *h121_net = gate_gg(hgg_radware, 121.8, 2.0);
gPad->SetLogy(0);
c2->Draw();
No description has been provided for this image

比较两张 gate 谱时，应看到原始 gate 谱中的连续平台在扣本底后降低，而与 $121.8\,\mathrm{keV}$ 有真实级联关系的峰仍然保留。弱峰是否可信，需要结合统计误差和反向 gate 检查。

用 $^{133}\mathrm{Ba}$ 检查符合关系¶

本实验源包含 $^{152}\mathrm{Eu}$ 和 $^{133}\mathrm{Ba}$。在扣本底后的二维矩阵上，可以利用 $^{133}\mathrm{Ba}$ 的已知级联关系检查矩阵是否合理。

常见的 $^{133}\mathrm{Ba}$ $\gamma$ 能量包括

$$ 53,\ 80,\ 160,\ 223,\ 276,\ 302,\ 356,\ 383\ \mathrm{keV}. $$

根据衰变纲图，可以预期不同 gate 中是否出现对应峰。下面用 O 表示预期有符合，用 x 表示预期没有明显符合。

gate \ peak $80$ $160$ $223$ $302$ $383$ $53$ $276$ $356$
$80$ O x O O x O O O
$223$ O O x x x O x x
$160$ x x O x x O O x
$302$ O x x x x O x x
$383$ x x x x x O x x
$53$ O O O O O x x x
$276$ O O x x x x x x
$356$ O x x x x x x x

逐个开 gate：

In [37]:
TH1D *g80 = gate_gg(hgg_radware, 80.0, 2.0, 400);
c2->Draw();
No description has been provided for this image
In [38]:
TH1D *g223 = gate_gg(hgg_radware, 223.0, 2.0, 400);
c2->Draw();
No description has been provided for this image
In [39]:
TH1D *g160 = gate_gg(hgg_radware, 160.0, 2.0, 400);
c2->Draw();
No description has been provided for this image

g160：80附近有很明显峰结构，而在g80：160附近没有峰结构。可以基本断定80和160没有符合关系。参照[68]的右侧二维图的(80,160)附近结构。

In [41]:
c2->Clear();
TH1D *g302 = gate_gg(hgg_radware, 302.0, 2.0, 400);
c2->Draw();
No description has been provided for this image
In [42]:
c2->Clear();
TH1D *g383 = gate_gg(hgg_radware, 383.0, 2.0, 400);
c2->Draw();
No description has been provided for this image
In [43]:
c2->Clear();
TH1D *g53 = gate_gg(hgg_radware, 53.0, 2.0, 400);
c2->Draw();
No description has been provided for this image
In [44]:
c2->Clear();
TH1D *g276 = gate_gg(hgg_radware, 276.0, 2.0, 400);
c2->Draw();
No description has been provided for this image
In [45]:
c2->Clear();
TH1D *g356 = gate_gg(hgg_radware, 356.0, 2.0, 400);
c2->Draw();
No description has been provided for this image

检查符合关系时不要只看一个方向。一个候选符合关系至少应通过反向 gate 检查。例如，若 gate(160) 中出现 $80\,\mathrm{keV}$ 峰，还应检查 gate(80) 中是否存在 $160\,\mathrm{keV}$ 峰。

同时还要考虑低能探测效率、中间分支、峰面积显著性、summing peak 和 X-ray 干扰。二维矩阵最直接回答的问题是：两条 $\gamma$ 射线是否属于同一次级联退激过程。完整能级纲图还需要结合能量和关系、相对强度、效率修正、内转换、角关联和偏振等信息。

关于偶然符合¶

Radware approach 主要处理 Compton continuum 和 quasi-continuum 形成的连续本底，不等同于时间偶然符合扣除。

在本节代码中，二维矩阵由 prompt 时间窗内的 $\gamma-\gamma$ pair 构成。如果实验计数率较低，偶然符合通常只是较小修正，可以先观察 Radware 减本底后的效果。

在高计数率束流实验中，低能区可能出现强 X-ray。例如 Au、Pb 等材料被束流激发后，会产生特征 X-ray。这些 X-ray 与反应产物的 $\gamma$ 射线通常没有真实级联关系，但可能由于计数率高而进入 prompt 窗，形成偶然符合。

此时，应进一步构造 accidental matrix。设：

  • PAM：prompt window 内的矩阵，包含 prompt 和 accidental；
  • AM：prompt window 外的 accidental matrix；
  • PM：扣除偶然符合后的 prompt matrix。

若 prompt 时间窗宽度为 $\Delta t_P$，accidental 时间窗总宽度为 $\Delta t_A$，则

$$ PM = PAM - \frac{\Delta t_P}{\Delta t_A}AM. $$

在实际分析中，通常先进行 accidental subtraction，再对扣除偶然符合后的 prompt matrix 做 Radware background subtraction。

原则上，偶然符合时间窗的选择应紧挨着符合时间窗的两个对称时间窗。图b和c展示了扣除偶然符合前后的投影谱对比图。黑色表示未扣除偶然符合的投影谱，而红色表示扣除偶然符合后的投影谱。可以看出，通过偶然符合的扣除，197Au的X射线影响大大减少。

作业：¶

  1. 参考 5.1 中的 $^{152}$Eu 简化纲图，对同一分支内较强的 $\gamma$ 跃迁开窗，验证它们之间的符合关系。所有符合关系应用反向开窗交叉检查。
  2. 验证左右两分支的gamma之间没有符合关系，并在的简化纲图基础上进行扩展。
In [ ]: