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

6.2 衰变实验：注入—衰变的位置和时间关联

学习目的：

1. 理解重离子注入事件与后续衰变事件之间的 position-time correlation 方法。
2. 掌握利用 DSSD 位置和 timestamp 建立注入-衰变关联的基本流程。
3. 理解随机关联本底的来源，并掌握正时间本底和负时间本底的基本扣除方法。
4. 通过 $^{211}\mathrm{Ac}$ 的 $\alpha$ 衰变，学习半衰期提取和两级 $\alpha$ 衰变关联的基本方法。

Data:

本节使用上一节生成的 sort1.root 文件。该文件已经完成了 DSSD 正背面匹配、相邻条能量共享修正以及 MWPC-DSSD prompt 符合判断。

TTree branch:

• timestamp：DSSD 事件 timestamp，单位为 ns。
• xstrip：DSSD 正面条号。
• ystrip：DSSD 背面条号。
• de：DSSD 能量。
• me：MWPC 能量。
  • me>0：该 DSSD 事件与 MWPC 有 prompt 符合，主要对应重离子注入事件。
  • me<0：该 DSSD 事件没有 MWPC prompt 符合，主要作为衰变候选事件。

上一节只是把原始单通道 hit 重构成 DSSD 事件，并利用 MWPC 信息初步区分注入和衰变候选。本节的核心问题是：如何判断某一个衰变事件是否来自某一个已经注入到 DSSD 中的重离子。

对于停止在 DSSD 中的重离子余核，其后续 $\alpha$ 衰变应发生在同一位置或非常邻近的位置，并且衰变时间应晚于注入时间。因此，注入-衰变关联需要同时使用位置信息和时间信息，这种方法称为 position-time correlation。

读取 sort1.root

%jsroot on

TCanvas *c1 = new TCanvas("c1","c1");

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

tree->Print();

ULong64_t timestamp;
Int_t xstrip, ystrip;
Float_t me, de;

tree->SetBranchAddress("timestamp", &timestamp);
tree->SetBranchAddress("xstrip", &xstrip);
tree->SetBranchAddress("ystrip", &ystrip);
tree->SetBranchAddress("de", &de);
tree->SetBranchAddress("me", &me);
TH1::SetDefaultSumw2();

******************************************************************************
*Tree    :tree      : sorted events                                          *
*Entries :   107133 : Total =         2580802 bytes  File  Size =    1454786 *
*        :          : Tree compression factor =   1.77                       *
******************************************************************************
*Br    0 :timestamp : timestamp/l                                            *
*Entries :   107133 : Total  Size=     859985 bytes  File Size  =     544208 *
*Baskets :       27 : Basket Size=      32000 bytes  Compression=   1.58     *
*............................................................................*
*Br    1 :me        : me/F                                                   *
*Entries :   107133 : Total  Size=     430071 bytes  File Size  =     246423 *
*Baskets :       14 : Basket Size=      32000 bytes  Compression=   1.74     *
*............................................................................*
*Br    2 :de        : de/F                                                   *
*Entries :   107133 : Total  Size=     430071 bytes  File Size  =     379987 *
*Baskets :       14 : Basket Size=      32000 bytes  Compression=   1.13     *
*............................................................................*
*Br    3 :xstrip    : xstrip/I                                               *
*Entries :   107133 : Total  Size=     430143 bytes  File Size  =     133098 *
*Baskets :       14 : Basket Size=      32000 bytes  Compression=   3.23     *
*............................................................................*
*Br    4 :ystrip    : ystrip/I                                               *
*Entries :   107133 : Total  Size=     430143 bytes  File Size  =     149602 *
*Baskets :       14 : Basket Size=      32000 bytes  Compression=   2.87     *
*............................................................................*

$\alpha$ 能谱

首先比较所有 DSSD 事件、与 MWPC 符合的 DSSD 事件以及无 MWPC 符合的 DSSD 事件的能量谱。

tree->Draw("de>>hde(3000,0,30000)", "", "hist");
tree->Draw("de>>hdem(3000,0,30000)", "me>0", "hist");
tree->Draw("de>>hdenm(3000,0,30000)", "me<0", "hist");

TH1F *hde   = (TH1F*)gROOT->FindObject("hde");
TH1F *hdem  = (TH1F*)gROOT->FindObject("hdem");
TH1F *hdenm = (TH1F*)gROOT->FindObject("hdenm");

hde->SetTitle("DSSD energy spectrum;Energy (keV);Counts");
hde->Draw("hist");

hdem->SetLineColor(kGreen);
hdenm->SetLineColor(kRed);

hdem->Draw("hist same");
hdenm->Draw("hist same");

c1->SetLogy(0);
c1->Draw();

其中：

• 黑色谱：所有 DSSD 事件。
• 绿色谱：与 MWPC 有 prompt 符合的 DSSD 事件，主要对应重离子注入。
• 红色谱：无 MWPC prompt 符合的 DSSD 事件，主要包含衰变候选事件。

明显的 $\alpha$ 衰变峰主要出现在 me<0 的谱中。这与物理图像一致：重离子注入时有 MWPC 信号，而 $\alpha$ 衰变发生在注入之后，通常不与 MWPC prompt 符合。

hdenm->GetXaxis()->SetRangeUser(5000,8000);
hdenm->Draw("hist");

c1->SetLogy(0);
c1->Draw();

DSSD 位置分布

在建立 position-time correlation 之前，先检查 DSSD 上所有事件的位置分布。

tree->Draw("ystrip:xstrip>>hxypos(48,-0.5,47.5,128,-0.5,127.5)", "", "colz");

c1->SetLogy(0);
c1->SetLogz(0);
c1->Draw();

DSSD 的位置分布反映了束流注入和衰变事件在探测器上的空间分布。对于注入-衰变关联，位置分布不能过于集中。若束流集中在少数几个 pixel 上，则同一位置在较短时间内可能有多个重离子注入，随机关联本底会明显增加。

按照 MWPC 信息对事件进行分类

根据上一节的定义：

• 重离子注入事件：me>0
• 衰变候选事件：me<0

需要注意，me<0 并不等价于纯净的衰变事件。它只表示该 DSSD 事件没有 MWPC prompt 符合，其中仍可能包含未被 veto 的轻粒子、本底或随机事件。后续还需要通过位置关联、时间关联和能量选择进一步筛选。

%%cpp -d
struct dssd
{
    Float_t energy;
    Int_t xstrip;
    Int_t ystrip;
};

dssd ds;

multimap<ULong64_t, dssd> mapimp; // implantation events
multimap<ULong64_t, dssd> mapdec; // decay-candidate events

Long64_t nentries = tree->GetEntriesFast();
for(Long64_t jentry=0; jentry<nentries; jentry++) {
    tree->GetEntry(jentry);

    ds.energy = de;
    ds.xstrip = xstrip;
    ds.ystrip = ystrip;
    if(me > 0) {
        mapimp.insert(pair<ULong64_t,dssd>(timestamp, ds));
    }
    else {
        mapdec.insert(pair<ULong64_t,dssd>(timestamp, ds));
    }
}

cout << "The number of implantation / decay-candidate events = "
     << mapimp.size() << " / " << mapdec.size() << endl;

The number of implantation / decay-candidate events = 77672 / 29461

这里仍然使用 multimap，是因为后续需要按照 timestamp 快速寻找某一注入事件前后一定时间范围内的衰变候选事件。

重离子与衰变事件的关联

位置关联

对于停止在 DSSD 中的重离子余核，其后续 $\alpha$ 衰变应发生在注入位置附近。DSSD 的一个 pixel 由一个正面条和一个背面条共同确定。

若注入事件位置为

$$ (hx,hy), $$

衰变候选事件位置为

$$ (bx,by), $$

最严格的位置关联条件为

$$ hx=bx,\quad hy=by. $$

在束流强度不高、DSSD pixel 足够小的情况下，可以近似认为：在一个合适的衰变时间窗内，同一 pixel 中的衰变事件主要来自该 pixel 中最近的重离子注入。

这也是衰变实验中利用 DSSD 做注入-衰变关联的基本物理依据。

时间关联

设重离子注入时间为 $ht$，衰变候选事件时间为 $bt$。定义衰变时间为

$$ \Delta t = bt - ht. $$

真实的注入-衰变关联应满足

$$ \Delta t > 0. $$

如果束流强度较低，在一个衰变时间窗内同一 pixel 中通常只有一个重离子注入，此时注入事件和衰变事件之间的对应关系比较清楚，$\Delta t$ 服从指数衰减分布。

当束流强度较高时，在同一 pixel 和同一时间窗内可能有多个重离子注入和多个衰变候选事件。此时，一部分注入-衰变组合是真实关联，另一部分只是随机组合。真实关联的 $\Delta t$ 服从指数衰减分布，随机组合在有限时间范围内近似形成平台本底。因此实验得到的衰变时间谱通常表现为：

$$ N(t)=B+N_0\exp\left(-\frac{t\ln2}{T_{1/2}}\right). $$

其中 $B$ 是随机关联本底，$T_{1/2}$ 是待提取的半衰期。

位置关联条件

对于 $\alpha$ 衰变，由于 $\alpha$ 粒子在硅中的射程较短，一般首先采用同一 pixel 的位置关联条件：

$$ \Delta x=0,\quad \Delta y=0. $$

在实际分析中，也可以扩大位置关联范围，例如

$$ |\Delta x|<2,\quad |\Delta y|<2. $$

这种条件包含同一 pixel 以及相邻 pixel。扩大位置范围可以增加效率，但也会增加随机关联本底。因此，位置关联范围必须通过数据检验确定，而不能任意选择。

对于本节数据，后面将比较同一 pixel 和相邻 pixel 的衰变时间谱，以确定合理的位置关联范围。

重新组织事件结构

为了后续分析，需要把每一个重离子注入事件与其附近的衰变候选事件重新组织成一个新的 tree。对每一个注入事件，寻找其前后一段时间范围内、位置满足条件的所有衰变候选事件。

这里保留负衰变时间区域。负衰变时间

$$ \Delta t < 0 $$

表示衰变候选事件发生在注入事件之前，不可能是真实的母核衰变。因此，负时间关联可以用来估计随机关联本底。

本节先采用较宽的位置条件

$$ |\Delta x|<2,\quad |\Delta y|<2, $$

把候选事件存入 decay.root。后续再通过 cut 选择同一 pixel 或相邻 pixel 进行比较。

下面实现的是 all-pairs correlation：每个注入与窗口内全部候选关联，不只选择最近一次注入。同一个衰变候选可能进入多个注入的列表。稳定束流和探测条件下，其随机关联部分可近似为平台；若改成“最近一次注入”或“第一个衰变”，选择本身就会改变时间分布，不能照搬同一平台模型。

// implantation event
ULong64_t hts;
Double_t he;
Int_t hx, hy;

// decay-candidate events within the decay-time window
const Int_t MAX_DEC_HIT = 1000;

Int_t bhit;
ULong64_t bts[MAX_DEC_HIT];
Double_t be[MAX_DEC_HIT];
Int_t bx[MAX_DEC_HIT], by[MAX_DEC_HIT];
Double_t decaytime[MAX_DEC_HIT]; // ms

// positive-time same-pixel decay candidates, used later for alpha-cascade analysis
Int_t bphit;
ULong64_t bpts[MAX_DEC_HIT];
Double_t bpe[MAX_DEC_HIT];

TFile *fout = new TFile("decay.root", "RECREATE");
TTree *tout = new TTree("tree", "decay");

tout->Branch("hts", &hts, "hts/l");
tout->Branch("he", &he, "he/D");
tout->Branch("hx", &hx, "hx/I");
tout->Branch("hy", &hy, "hy/I");

tout->Branch("bhit", &bhit, "bhit/I");
tout->Branch("bts", bts, "bts[bhit]/l");
tout->Branch("be", be, "be[bhit]/D");
tout->Branch("bx", bx, "bx[bhit]/I");
tout->Branch("by", by, "by[bhit]/I");
tout->Branch("decaytime", decaytime, "decaytime[bhit]/D");

tout->Branch("bphit", &bphit, "bphit/I");
tout->Branch("bpts", bpts, "bpts[bphit]/l");
tout->Branch("bpe", bpe, "bpe[bphit]/D");

衰变时间窗取为 20 s。该窗口足够覆盖 $^{211}\mathrm{Ac}$ 的多个半衰期，同时保留远离 prompt 衰变区的随机本底区域。

实际保存区间是相对注入时间的 −10 s 至 +20 s。run 起止附近的注入可观测时间较短；当运行时间与窗口长度可比时，应统计各时间 bin 的有效注入暴露量，或选择前后记录完整的注入。缺失的记录时间不能当成零衰变计数。

ULong64_t twindow = 20000000000ULL; // 20 s in ns

Int_t n = 0;
for(auto ia = mapimp.begin(); ia != mapimp.end(); ++ia) {
    hts = ia->first;
    he = ia->second.energy;
    hx = ia->second.xstrip;
    hy = ia->second.ystrip;

    // Include a negative-time region for random-background estimation.
    ULong64_t tmin = (ia->first > twindow/2) ? ia->first - twindow/2 : 0;
    ULong64_t tmax = ia->first + twindow;

    auto ib1 = mapdec.lower_bound(tmin);
    auto ib2 = mapdec.upper_bound(tmax);

    bhit = 0;
    bphit = 0;
    for( ; ib1 != ib2; ++ib1) {
        Int_t delx = TMath::Abs(ib1->second.xstrip - hx);
        Int_t dely = TMath::Abs(ib1->second.ystrip - hy);
        if (delx >= 2 || dely >= 2) continue; // 只保留同一及紧邻 pixel
        {
            if(bhit >= MAX_DEC_HIT) throw runtime_error("Increase MAX_DEC_HIT; candidate list exceeds capacity");

            bts[bhit] = ib1->first;
            be[bhit] = ib1->second.energy;
            bx[bhit] = ib1->second.xstrip;
            by[bhit] = ib1->second.ystrip;

            Long64_t dt_ns = Long64_t(ib1->first) - Long64_t(hts);
            decaytime[bhit] = dt_ns/1.e6; // ms

            // For alpha-cascade analysis:
            // keep only positive-time same-pixel decay candidates.
            if(decaytime[bhit] > 0 && delx == 0 && dely == 0) {
                // 正时间同 pixel 候选是 bhit 的子集，前面的容量检查也保护此数组。
                bpts[bphit] = ib1->first;
                bpe[bphit] = ib1->second.energy;
                bphit++;
            }

            bhit++;
        }
    }

    tout->Fill(); // 保留零候选的注入，避免改变注入计数分母

    n++;
    if(n%1000 == 0) cout << ".";
}

cout << endl;
cout << "Done!" << endl;

tout->Write();
fout->Close();

.............................................................................
Done!

这里要注意两点：

1. decaytime 必须用有符号整数计算。若直接用两个 ULong64_t 相减，再转换为 Long64_t，负时间会发生无符号下溢。
2. bhit 和 bphit 应设置最大数组长度，避免在高计数率或过宽时间窗下数组越界。

读取 decay.root

TFile *fdecay = new TFile("decay.root");
TTree *tree = (TTree*)fdecay->Get("tree");

tree->Print();

******************************************************************************
*Tree    :tree      : decay                                                  *
*Entries :    77672 : Total =         6837159 bytes  File  Size =    2731712 *
*        :          : Tree compression factor =   2.50                       *
******************************************************************************
*Br    0 :hts       : hts/l                                                  *
*Entries :    77672 : Total  Size=     623481 bytes  File Size  =     397618 *
*Baskets :       20 : Basket Size=      32000 bytes  Compression=   1.57     *
*............................................................................*
*Br    1 :he        : he/D                                                   *
*Entries :    77672 : Total  Size=     623457 bytes  File Size  =     351740 *
*Baskets :       20 : Basket Size=      32000 bytes  Compression=   1.77     *
*............................................................................*
*Br    2 :hx        : hx/I                                                   *
*Entries :    77672 : Total  Size=     311871 bytes  File Size  =      96589 *
*Baskets :       10 : Basket Size=      32000 bytes  Compression=   3.22     *
*............................................................................*
*Br    3 :hy        : hy/I                                                   *
*Entries :    77672 : Total  Size=     311871 bytes  File Size  =     108761 *
*Baskets :       10 : Basket Size=      32000 bytes  Compression=   2.86     *
*............................................................................*
*Br    4 :bhit      : bhit/I                                                 *
*Entries :    77672 : Total  Size=     311899 bytes  File Size  =      46357 *
*Baskets :       10 : Basket Size=      32000 bytes  Compression=   6.72     *
*............................................................................*
*Br    5 :bts       : bts[bhit]/l                                            *
*Entries :    77672 : Total  Size=     787638 bytes  File Size  =     318781 *
*Baskets :       35 : Basket Size=      32000 bytes  Compression=   2.47     *
*............................................................................*
*Br    6 :be        : be[bhit]/D                                             *
*Entries :    77672 : Total  Size=     787599 bytes  File Size  =     284322 *
*Baskets :       35 : Basket Size=      32000 bytes  Compression=   2.77     *
*............................................................................*
*Br    7 :bx        : bx[bhit]/I                                             *
*Entries :    77672 : Total  Size=     550244 bytes  File Size  =     159221 *
*Baskets :       27 : Basket Size=      32000 bytes  Compression=   3.45     *
*............................................................................*
*Br    8 :by        : by[bhit]/I                                             *
*Entries :    77672 : Total  Size=     550244 bytes  File Size  =     162920 *
*Baskets :       27 : Basket Size=      32000 bytes  Compression=   3.37     *
*............................................................................*
*Br    9 :decaytime : decaytime[bhit]/D                                      *
*Entries :    77672 : Total  Size=     787872 bytes  File Size  =     526490 *
*Baskets :       35 : Basket Size=      32000 bytes  Compression=   1.49     *
*............................................................................*
*Br   10 :bphit     : bphit/I                                                *
*Entries :    77672 : Total  Size=     311913 bytes  File Size  =      28349 *
*Baskets :       10 : Basket Size=      32000 bytes  Compression=  10.98     *
*............................................................................*
*Br   11 :bpts      : bpts[bphit]/l                                          *
*Entries :    77672 : Total  Size=     439600 bytes  File Size  =     132968 *
*Baskets :       24 : Basket Size=      32000 bytes  Compression=   3.30     *
*............................................................................*
*Br   12 :bpe       : bpe[bphit]/D                                           *
*Entries :    77672 : Total  Size=     439572 bytes  File Size  =     113694 *
*Baskets :       24 : Basket Size=      32000 bytes  Compression=   3.86     *
*............................................................................*

decay.root 中的一个 entry 对应一个重离子注入事件，以及在指定时间窗和位置范围内找到的所有衰变候选事件。

主要 branch 为：

• hts, he, hx, hy：注入事件的时间、能量和位置。
• bhit：与该注入事件关联的候选衰变事件数。
• bts, be, bx, by：候选衰变事件的时间、能量和位置。
• decaytime：候选衰变事件相对于注入事件的时间差，单位为 ms。
• bphit, bpts, bpe：正时间、同一 pixel 的衰变候选事件，用于后续两级 $\alpha$ 衰变关联。

Heavy ion and decay correlation in the same pixel $\Delta x=0$ and $\Delta y=0$

首先只选择同一 pixel 中的注入-衰变候选事件，画出衰变能量与衰变时间的二维图。

TCut cdxdy00 = "abs(bx-hx)==0 && abs(by-hy)==0";

tree->Draw("be:decaytime>>hEvsT00(200,0,2e4,300,5000,8000)",
           cdxdy00,
           "colz");

c1->SetLogy(0);
c1->SetLogz(0);
c1->Draw();

如果某个 $\alpha$ 能峰来自真实的注入后衰变，它应主要分布在正衰变时间区域，并随时间呈指数衰减。

$^{211}\mathrm{Ac}$ 的衰变时间分布

选择 $^{211}\mathrm{Ac}$ 的 $\alpha$ 能量区域：

$$ 7400~\mathrm{keV}<E_\alpha<7510~\mathrm{keV}. $$

TCut cAc211a = "be>7400 && be<7510 && decaytime>0";

tree->Draw("decaytime>>hdtAc211a(200,0,2e4)",
           cdxdy00 && cAc211a,
           "hist");

c1->SetLogy();
c1->Draw();

该时间谱由两部分组成：

• 真实的 $^{211}\mathrm{Ac}$ 衰变，服从指数衰减；
• 随机注入-衰变组合，形成近似平台本底。

因此，时间谱应使用“指数 + 常数平台”的形式拟合。

Heavy ion and decay correlation with neighboring pixels

$\Delta x=\pm 1$ and $\Delta y=\pm 1$

先检查对角相邻 pixel 中是否存在明显关联。

TCut cdxdy11 = "abs(bx-hx)==1 && abs(by-hy)==1 && decaytime>0";
TCut cAc211 = "be>7400 && be<7510";

tree->Draw("decaytime>>hdtAc211b(200,0,2e4)",
           cdxdy11 && cAc211,
           "hist");

c1->SetLogy();
c1->Draw();

$\Delta x=\pm 1$ or $\Delta y=\pm 1$

再检查上下左右相邻 pixel 中是否存在明显关联。

TCut cdxdy0110 = "abs(bx-hx) + abs(by-hy)==1 && decaytime>0";

tree->Draw("decaytime>>hdtAc211c(200,0,2e4)",
           cdxdy0110 && cAc211,
           "hist");

c1->SetLogy();
c1->Draw();

从同一 pixel 和相邻 pixel 的时间谱对比可以判断：本数据中 $^{211}\mathrm{Ac}$ 的 $\alpha$ 衰变关联主要出现在同一 pixel，相邻 pixel 中没有明显关联或贡献很小。

因此，下面的处理采用

$$ \Delta x=0,\quad \Delta y=0 $$

作为位置关联条件。

这个检查步骤很重要。不同粒子、不同能量和不同探测器厚度下，衰变粒子的射程和散射情况不同，合理的位置关联范围也可能不同。不能在所有实验中机械地使用同一个位置条件。

用远正时间区估计平台本底

对于 $^{211}\mathrm{Ac}$，同一 pixel、选定能量区间内的衰变时间谱可以写成

$$ N(t)=B+N_0\exp\left(-\frac{t\ln2}{T_{1/2}}\right). $$

当时间远大于半衰期时，指数衰变项已经很小，剩余计数主要来自随机关联本底。因此，可以先在远离衰变区的时间范围拟合常数平台。

这里以 $15$--$20$ s 区域估计平台本底。对于低计数的时间谱，拟合时应使用 likelihood 选项 "L"。

TF1 *fp01 = new TF1("fp01", "pol0", 1.5e4, 2.0e4);

hdtAc211a->Fit(fp01, "LR");

Double_t p01 = fp01->GetParameter(0);

hdtAc211a->Draw("hist");
fp01->Draw("same");
c1->SetLogy();
c1->Draw();

****************************************
Minimizer is Minuit2 / Migrad
MinFCN                    =      32.4205
Chi2                      =       64.841
NDf                       =           49
Edm                       =  2.27219e-07
NCalls                    =           51
p0                        =      1.77987   +/-   0.188666    

指数 + 平台拟合

先用远时间区估计平台 $B$，再固定该估计值拟合衰变部分。参数依次为本底高度、指数项初始高度和半衰期（ms）。本图 bin 宽为 100 ms，拟合从 100 ms 开始，跳过零点附近的整个首 bin；这不意味着仪器死时间为 100 ms。

L 用于未扣本底的计数谱；I 用 bin 内函数平均值比较计数，R 使用给定范围，S 返回拟合结果。固定 $B$ 后得到的半衰期误差只包含条件拟合误差。下面把 $B$ 改为 $B\pm\sigma_B$ 分别重拟合，用半衰期变化估计本底统计误差的贡献：

$$\sigma_T^2\simeq\sigma_{T\mid B}^2+\left[\frac{T(B+\sigma_B)-T(B-\sigma_B)}2\right]^2.$$

这里本底区与衰变拟合区不重叠，先忽略 all-pairs 共享事例带来的相关性，采用一阶误差传播。高注入率时需按完整注入—衰变记录评估这种相关性，不能只靠缩小拟合误差。更一般的做法是联合拟合信号区和本底区，让 $B$ 的不确定度及相关性一同进入拟合。

TF1 *fdecay1 = new TF1("fdecay1",
    "[0]+[1]*exp(-x*log(2.)/[2])",100,10000);
fdecay1->SetParNames("Background","Amplitude","Half-life (ms)");
fdecay1->FixParameter(0,p01);
fdecay1->SetParameter(1,600);
fdecay1->SetParameter(2,250); // 接近已知半衰期的初值，不是固定值
fdecay1->SetParLimits(1,0,1e6);
fdecay1->SetParLimits(2,1,20000);
TFitResultPtr result1 = hdtAc211a->Fit(fdecay1,"SLIR");
if (int(result1)!=0) throw runtime_error("Decay fit failed");
hdtAc211a->SetTitle("^{211}Ac;Correlation time (ms);Counts / 100 ms");
hdtAc211a->Draw("hist");
fdecay1->Draw("same");
c1->SetLogy();
c1->Draw();

****************************************
Minimizer is Minuit2 / Migrad
MinFCN                    =      52.2209
Chi2                      =      104.442
NDf                       =           97
Edm                       =  6.29568e-08
NCalls                    =           36
Background                =      1.77987                      	 (fixed)
Amplitude                 =      591.054   +/-   24.5674      	 (limited)
Half-life (ms)            =      256.363   +/-   7.00126      	 (limited)

再用本底估计值上下移动一个标准差的两次拟合，计算本底不确定度对半衰期的贡献。N 不保存试探拟合函数，因此不会替换上图的中心结果。

double half1=fdecay1->GetParameter(2), errFit1=fdecay1->GetParError(2);
double varied1[2];
for (int k=0;k<2;++k) {
    TF1 trial(*fdecay1);
    trial.FixParameter(0,p01+(2*k-1)*fp01->GetParError(0));
    if (int(hdtAc211a->Fit(&trial,"QLIRN"))!=0) throw runtime_error("Background variation fit failed");
    varied1[k]=trial.GetParameter(2);
}
double errB1=abs(varied1[1]-varied1[0])/2;
cout << "Half-life = " << half1 << " ms; conditional fit error = " << errFit1
     << ", background contribution = " << errB1
     << ", combined = " << hypot(errFit1,errB1) << " ms" << endl;

Half-life = 256.363 ms; conditional fit error = 7.00126, background contribution = 1.62651, combined = 7.18772 ms

通过该拟合可以得到 $^{211}\mathrm{Ac}$ 的半衰期。需要注意，拟合结果依赖于能量门、位置条件、拟合范围和本底估计范围，因此这些条件应在结果报告中明确给出。

$\alpha$ 能谱减本底

除了拟合时间谱，也可以利用时间窗对 $\alpha$ 能谱进行本底扣除。

基本思想是：

• 近时间窗：包含真实衰变信号和随机本底。
• 远时间窗：主要包含随机本底。
• 两者相减后得到衰变信号的近似净能谱。

对于半衰期约为 $250$ ms 的 $^{211}\mathrm{Ac}$，可取：

• 信号区：$0<t<5T_{1/2}$。
• 本底区：$20T_{1/2}<t<25T_{1/2}$。

两个时间窗宽度相同，因此可以直接相减。若时间窗宽度不同，则需要按时间窗宽度归一化。

对计数谱作本底扣除后，净 bin 可以为负。若两窗计数独立且缩放系数为 $a$，则 $\mathrm{Var}(N_{net})=N_{signal}+a^2N_{background}$；不能再对净谱直接使用 Poisson likelihood。这里保存原始计数谱并用 Sumw2 传播求差误差。

tree->Draw("be>>hbeall(300,5000,8000)",
           cdxdy00 && "decaytime>0 && decaytime<5*250", "hist");

tree->Draw("be>>hbebkg(300,5000,8000)",
           cdxdy00 && "decaytime>20*250 && decaytime<25*250", "hist");

TH1F *hbeall = (TH1F*)gROOT->FindObject("hbeall");
TH1F *hbebkg = (TH1F*)gROOT->FindObject("hbebkg");

TH1F *hbenet = (TH1F*)hbeall->Clone("hbenet");
hbenet->SetTitle("background-subtracted alpha spectrum;Energy (keV);Counts");

hbenet->Add(hbebkg, -1.0);

hbeall->SetLineColor(kGreen);
hbebkg->SetLineColor(kBlue);
hbenet->SetLineColor(kRed);

hbeall->Draw("hist");
hbebkg->Draw("hist same");
hbenet->Draw("hist same");

c1->SetLogy(0);
c1->Draw();

这种做法只适用于半衰期与所选时间窗匹配的衰变成分。对于半衰期更长的 $\alpha$ 峰，应取更长的信号时间窗，本底区也应相应后移。否则会把一部分真实衰变当成本底扣除掉。

减负衰变时间本底

负衰变时间表示候选衰变事件发生在注入事件之前：

$$ \Delta t = bt-ht < 0. $$

这在真实母核衰变中是不可能的。因此，负时间关联可以作为随机关联本底的实验估计。与远正时间本底相比，负时间本底不依赖待研究核素的半衰期，通常更适合用于检查随机关联。

TCut cdxdy00 = "abs(bx-hx)==0 && abs(by-hy)==0";

tree->Draw("be:decaytime>>hEvsTall(300,-1e4,2e4,300,5000,8000)",
           cdxdy00,
           "colz");

c1->SetLogy(0);
c1->SetLogz(0);
c1->Draw();

TCut cAc211a_alltime = "be>7400 && be<7510";

tree->Draw("decaytime>>hdtAc211d(300,-1e4,2e4)",
           cdxdy00 && cAc211a_alltime,
           "hist");

c1->SetLogy();
c1->Draw();

用负时间区估计平台本底

TF1 *fp01d = new TF1("fp01d", "pol0", -0.8e4, -0.3e4);

hdtAc211d->Fit(fp01d, "LR");

Double_t p01d = fp01d->GetParameter(0);

hdtAc211d->Draw("hist");
fp01d->Draw("same");
c1->SetLogy();
c1->Draw();

****************************************
Minimizer is Minuit2 / Migrad
MinFCN                    =      27.3895
Chi2                      =       54.779
NDf                       =           49
Edm                       =  6.99916e-08
NCalls                    =           51
p0                        =      1.57993   +/-   0.177756    

使用负时间本底进行指数 + 平台拟合

仍按上面的固定平台方法，并传播负时间区平台估计的统计误差。两次半衰期拟合使用同一批正时间事例，结果相关，不能把两者当作独立测量再取加权平均。

TF1 *fdecay2 = new TF1("fdecay2",
    "[0]+[1]*exp(-x*log(2.)/[2])",100,10000);
fdecay2->SetParNames("Background","Amplitude","Half-life (ms)");
fdecay2->FixParameter(0,p01d);
fdecay2->SetParameter(1,600);
fdecay2->SetParameter(2,250); // 接近已知半衰期的初值，不是固定值
fdecay2->SetParLimits(1,0,1e6);
fdecay2->SetParLimits(2,1,20000);
TFitResultPtr result2 = hdtAc211d->Fit(fdecay2,"SLIR");
if (int(result2)!=0) throw runtime_error("Decay fit failed");
hdtAc211d->SetTitle("^{211}Ac;Correlation time (ms);Counts / 100 ms");
hdtAc211d->Draw("hist");
fdecay2->Draw("same");
c1->SetLogy();
c1->Draw();

****************************************
Minimizer is Minuit2 / Migrad
MinFCN                    =      52.1304
Chi2                      =      104.261
NDf                       =           97
Edm                       =  3.06776e-07
NCalls                    =           36
Background                =      1.57993                      	 (fixed)
Amplitude                 =      587.596   +/-   24.3528      	 (limited)
Half-life (ms)            =      258.149   +/-   7.02636      	 (limited)

再用本底估计值上下移动一个标准差的两次拟合，计算本底不确定度对半衰期的贡献。N 不保存试探拟合函数，因此不会替换上图的中心结果。

double half2=fdecay2->GetParameter(2), errFit2=fdecay2->GetParError(2);
double varied2[2];
for (int k=0;k<2;++k) {
    TF1 trial(*fdecay2);
    trial.FixParameter(0,p01d+(2*k-1)*fp01d->GetParError(0));
    if (int(hdtAc211d->Fit(&trial,"QLIRN"))!=0) throw runtime_error("Background variation fit failed");
    varied2[k]=trial.GetParameter(2);
}
double errB2=abs(varied2[1]-varied2[0])/2;
cout << "Half-life = " << half2 << " ms; conditional fit error = " << errFit2
     << ", background contribution = " << errB2
     << ", combined = " << hypot(errFit2,errB2) << " ms" << endl;

Half-life = 258.149 ms; conditional fit error = 7.02636, background contribution = 1.65874, combined = 7.2195 ms

如果正时间远区本底和负时间本底给出的半衰期结果一致，说明随机本底估计比较稳定；若二者明显不同，则需要重新检查位置条件、能量门、时间窗、束流结构和可能的长寿命衰变贡献。

$\alpha$ 能谱减负时间本底

选择正时间信号区和等宽的负时间本底区。为了避开零点附近的死时间或事件重构效应，负时间本底区不直接取到 0，而是留出一定间隔。

TCut cdxdy00 = "abs(bx-hx)==0 && abs(by-hy)==0";

// signal: 0 < t < 5*T1/2
tree->Draw("be>>hbealln(300,5000,8000)",
           cdxdy00 && "decaytime>0 && decaytime<5*250", "hist");

TH1F *hbealln = (TH1F*)gROOT->FindObject("hbealln");

// background: an equal-width negative-time window, shifted away from zero
tree->Draw("be>>hbebkgn(300,5000,8000)",
           cdxdy00 && "decaytime>-5*250-300 && decaytime<-300", "hist");

TH1F *hbebkgn = (TH1F*)gROOT->FindObject("hbebkgn");

TH1F *hbenetn = (TH1F*)hbealln->Clone("hbenetn");
hbenetn->SetTitle("negative-time-background-subtracted alpha spectrum;Energy (keV);Counts");

hbenetn->Add(hbebkgn, -1.0);

hbealln->SetLineColor(kGreen);
hbebkgn->SetLineColor(kBlue);
hbenetn->SetLineColor(kRed);

hbealln->Draw("hist");
hbebkgn->Draw("hist same");
hbenetn->Draw("hist same");

c1->SetLogy(0);
c1->Draw();

负时间本底法的优点是：本底区不包含真实的正时间母核衰变，因此不容易受到半衰期选择的影响。它的前提是随机关联在正负时间两侧近似对称，并且束流强度和探测器状态在所选时间范围内稳定。

寻找下一级衰变事件

$^{211}\mathrm{Ac}$ 通过 $\alpha$ 衰变生成 $^{207}\mathrm{Fr}$，而 $^{207}\mathrm{Fr}$ 还可以继续通过 $\alpha$ 衰变生成 $^{203}\mathrm{At}$。

通过观察同一位置、连续发生的两级 $\alpha$ 衰变，可以进一步确认衰变链。若已知下一级衰变核的性质，也可以从衰变链下端向上推断未知母核的质量数和电荷数：

$$ A_N=A_p+4,\quad Z_N=Z_p+2. $$

由于重离子主要注入在 DSSD 表面附近，$\alpha$ 粒子可能向 DSSD 内部或外部发射。若只依赖 DSSD 测量，在浅注入、各向同性发射且只计算向硅内部发射的全能事件这一简化图像下，每一级的几何接受度约为 $50\%$，两级约为 $25\%$。实际全能效率还取决于注入深度、阈值、能量窗、分支比和时间窗。

在衰变时间窗内，将满足条件的所有衰变事件存入数组

在前面生成 decay.root 时，已经将正时间、同一 pixel 的候选衰变事件存入：

• bphit
• bpts[bphit]
• bpe[bphit]

其中：

• bphit：同一 pixel、正时间衰变候选事件数。
• bpts：这些衰变候选事件的 timestamp。
• bpe：这些衰变候选事件的能量。

这些数组可用于寻找同一注入事件之后的多级 $\alpha$ 衰变。

两级关联条件

对于 $^{211}\mathrm{Ac}\rightarrow ^{207}\mathrm{Fr}\rightarrow ^{203}\mathrm{At}$，可采用如下示范性条件：

第一级：

$$ t_{\alpha1}-t_\mathrm{imp}<700~\mathrm{ms}. $$

第二级：

$$ t_{\alpha2}-t_{\alpha1}<19000~\mathrm{ms}. $$

第二级时间范围原则上应覆盖 $^{207}\mathrm{Fr}$ 的多个半衰期。本例先展示两级关联。若用它提取半衰期，还需在拟合中计入有限时间窗和候选选择的影响。

tree->Draw("bphit", "bphit>0", "hist");

c1->SetLogy();
c1->Draw();

母核 $\alpha$-子核 $\alpha$ 能量关联图

下面使用时间顺序中的前两级正时间、同一 pixel 衰变候选事件，画出母核 $\alpha$ 与子核 $\alpha$ 的能量关联图。

TCut cs1 = "bphit>1 && (bpts[0]-hts)/1.e6<700";
TCut cs2 = "(bpts[1]-bpts[0])/1.e6<19000";

tree->Draw("bpe[1]:bpe[0]>>hbe12(200,6000,8000,120,6000,7200)",
           cs1 && cs2,
           "colz");

c1->SetLogy(0);
c1->SetLogz(0);
c1->Draw();

在二维能量关联图中，可以看到 $^{211}\mathrm{Ac}$ 的约 $7477$ keV $\alpha$ 与 $^{207}\mathrm{Fr}$ 的约 $6767$ keV $\alpha$ 之间存在关联。

需要注意，这里使用的是一个简化示范：bpe[0] 和 bpe[1] 表示时间顺序中的前两个正时间同 pixel 衰变候选事件。若同一注入事件后存在多个候选衰变，严格分析应遍历所有满足时间顺序的 $\alpha_i$-$\alpha_j$ 组合，而不是只使用前两个候选。

$^{207}\mathrm{Fr}$ 衰变时间

选择 $^{211}\mathrm{Ac}$ 和 $^{207}\mathrm{Fr}$ 的 $\alpha$ 能量门：

TCut cecut = "abs(bpe[1]-6767)<50 && abs(bpe[0]-7477)<50";

tree->Draw("(bpts[1]-bpts[0])/1.e6>>hdt12(100,0,20000)",
           cs1 && cs2 && cecut,
           "hist");

c1->SetLogy();
c1->Draw();

本例第二级候选统计较少，而且只选择时间顺序中的前两个候选，其时间分布受到候选选择和有限窗口影响。这里用于展示两级衰变关联；定量提取半衰期时，应建立与该选择一致的模型，而不是仅因窗口较短就认定原则上无法拟合。

作业

选择 $^{210}$Ra 的 α 能量区域，比较同一 pixel 与相邻 pixel 的关联时间谱，确定位置条件并提取半衰期。分别用远正时间区和负时间区估计本底，比较结果并给出半衰期及其误差。