衰变实验的数据分析 II - 重离子和衰变事件的关联¶
学习目的：¶
- 理解重离子注入事件与后续衰变事件之间的 position-time correlation 方法。
- 掌握利用 DSSD 位置和 timestamp 建立注入-衰变关联的基本流程。
- 理解随机关联本底的来源，并掌握正时间本底和负时间本底的基本扣除方法。
- 通过 $^{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¶
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", ×tamp);
tree->SetBranchAddress("xstrip", &xstrip);
tree->SetBranchAddress("ystrip", &ystrip);
tree->SetBranchAddress("de", &de);
tree->SetBranchAddress("me", &me);
****************************************************************************** *Tree :tree : sorted events * *Entries : 107136 : Total = 2580874 bytes File Size = 1454827 * * : : Tree compression factor = 1.77 * ****************************************************************************** *Br 0 :timestamp : timestamp/l * *Entries : 107136 : Total Size= 860009 bytes File Size = 544219 * *Baskets : 27 : Basket Size= 32000 bytes Compression= 1.58 * *............................................................................* *Br 1 :me : me/F * *Entries : 107136 : Total Size= 430083 bytes File Size = 246434 * *Baskets : 14 : Basket Size= 32000 bytes Compression= 1.74 * *............................................................................* *Br 2 :de : de/F * *Entries : 107136 : Total Size= 430083 bytes File Size = 379997 * *Baskets : 14 : Basket Size= 32000 bytes Compression= 1.13 * *............................................................................* *Br 3 :xstrip : xstrip/I * *Entries : 107136 : Total Size= 430155 bytes File Size = 133106 * *Baskets : 14 : Basket Size= 32000 bytes Compression= 3.23 * *............................................................................* *Br 4 :ystrip : ystrip/I * *Entries : 107136 : Total Size= 430155 bytes File Size = 149610 * *Baskets : 14 : Basket Size= 32000 bytes Compression= 2.87 * *............................................................................*
tree->Draw("de>>hde(3000,0,30000)");
tree->Draw("de>>hdem(3000,0,30000)", "me>0");
tree->Draw("de>>hdenm(3000,0,30000)", "me<0");
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();
hdem->SetLineColor(kGreen);
hdenm->SetLineColor(kRed);
hdem->Draw("same");
hdenm->Draw("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();
c1->SetLogy(0);
c1->Draw();
DSSD 位置分布¶
在建立 position-time correlation 之前，先检查 DSSD 上所有事件的位置分布。
tree->Draw("xstrip:ystrip>>hxypos(128,0,128,48,0,48)", "", "colz");
c1->SetLogy(0);
c1->SetLogz(0);
c1->Draw();
DSSD 的位置分布反映了束流注入和衰变事件在探测器上的空间分布。对于注入-衰变关联，位置分布不能过于集中。若束流集中在少数几个 pixel 上，则同一位置在较短时间内可能有多个重离子注入，随机关联本底会明显增加。
按照 MWPC 信息对事件进行分类¶
根据上一节的定义：
- 重离子注入事件：
me>0 - 衰变候选事件：
me<0
需要注意，me<0 并不等价于纯净的衰变事件。它只表示该 DSSD 事件没有 MWPC prompt 符合，其中仍可能包含未被 veto 的轻粒子、本底或随机事件。后续还需要通过位置关联、时间关联和能量选择进一步筛选。
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 = 77675 / 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 进行比较。
// 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 衰变区的随机本底区域。
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) {
if(bhit >= MAX_DEC_HIT) break;
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) {
if(bphit < MAX_DEC_HIT) {
bpts[bphit] = ib1->first;
bpe[bphit] = ib1->second.energy;
bphit++;
}
}
bhit++;
}
}
if(bhit > 0) tout->Fill();
n++;
if(n%1000 == 0) cout << ".";
}
cout << endl;
cout << "Done!" << endl;
tout->Write();
fout->Close();
............................................................................. Done!
这里要注意两点：
decaytime必须用有符号整数计算。若直接用两个ULong64_t相减，再转换为Long64_t，负时间会发生无符号下溢。bhit和bphit应设置最大数组长度，避免在高计数率或过宽时间窗下数组越界。
读取 decay.root¶
TFile *fdecay = new TFile("decay.root");
TTree *tree = (TTree*)fdecay->Get("tree");
tree->Print();
****************************************************************************** *Tree :tree : decay * *Entries : 39015 : Total = 4507918 bytes File Size = 2059645 * * : : Tree compression factor = 2.19 * ****************************************************************************** *Br 0 :hts : hts/l * *Entries : 39015 : Total Size= 313325 bytes File Size = 202856 * *Baskets : 10 : Basket Size= 32000 bytes Compression= 1.54 * *............................................................................* *Br 1 :he : he/D * *Entries : 39015 : Total Size= 313311 bytes File Size = 175806 * *Baskets : 10 : Basket Size= 32000 bytes Compression= 1.78 * *............................................................................* *Br 2 :hx : hx/I * *Entries : 39015 : Total Size= 156878 bytes File Size = 48252 * *Baskets : 5 : Basket Size= 32000 bytes Compression= 3.24 * *............................................................................* *Br 3 :hy : hy/I * *Entries : 39015 : Total Size= 156878 bytes File Size = 52323 * *Baskets : 5 : Basket Size= 32000 bytes Compression= 2.99 * *............................................................................* *Br 4 :bhit : bhit/I * *Entries : 39015 : Total Size= 156896 bytes File Size = 19722 * *Baskets : 5 : Basket Size= 32000 bytes Compression= 7.93 * *............................................................................* *Br 5 :bts : bts[bhit]/l * *Entries : 39015 : Total Size= 632070 bytes File Size = 295970 * *Baskets : 25 : Basket Size= 32000 bytes Compression= 2.13 * *............................................................................* *Br 6 :be : be[bhit]/D * *Entries : 39015 : Total Size= 632041 bytes File Size = 262662 * *Baskets : 25 : Basket Size= 32000 bytes Compression= 2.40 * *............................................................................* *Br 7 :bx : bx[bhit]/I * *Entries : 39015 : Total Size= 394763 bytes File Size = 137012 * *Baskets : 18 : Basket Size= 32000 bytes Compression= 2.88 * *............................................................................* *Br 8 :by : by[bhit]/I * *Entries : 39015 : Total Size= 394763 bytes File Size = 140957 * *Baskets : 18 : Basket Size= 32000 bytes Compression= 2.80 * *............................................................................* *Br 9 :decaytime : decaytime[bhit]/D * *Entries : 39015 : Total Size= 632244 bytes File Size = 476562 * *Baskets : 25 : Basket Size= 32000 bytes Compression= 1.32 * *............................................................................* *Br 10 :bphit : bphit/I * *Entries : 39015 : Total Size= 156905 bytes File Size = 18385 * *Baskets : 5 : Basket Size= 32000 bytes Compression= 8.51 * *............................................................................* *Br 11 :bpts : bpts[bphit]/l * *Entries : 39015 : Total Size= 283982 bytes File Size = 123092 * *Baskets : 14 : Basket Size= 32000 bytes Compression= 2.30 * *............................................................................* *Br 12 :bpe : bpe[bphit]/D * *Entries : 39015 : Total Size= 283964 bytes File Size = 103080 * *Baskets : 14 : Basket Size= 32000 bytes Compression= 2.75 * *............................................................................*
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,
"colz");
c1->SetLogy();
c1->Draw();
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,
"colz");
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,
"colz");
c1->SetLogy();
c1->Draw();
从同一 pixel 和相邻 pixel 的时间谱对比可以判断：本数据中 $^{211}\mathrm{Ac}$ 的 $\alpha$ 衰变关联主要出现在同一 pixel，相邻 pixel 中没有明显关联或贡献很小。
因此，下面的处理采用
$$ \Delta x=0,\quad \Delta y=0 $$
作为位置关联条件。
这个检查步骤很重要。不同粒子、不同能量和不同探测器厚度下，衰变粒子的射程和散射情况不同，合理的位置关联范围也可能不同。不能在所有实验中机械地使用同一个位置条件。
减去 20 个半衰期以后的平台本底¶
对于 $^{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);
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
指数 + 平台拟合¶
拟合 $\alpha$ 或 $\beta$ 衰变时间谱时，通常避开非常靠近零时间的区域，以减小电子学死时间、pile-up 或事件重构不完全带来的影响。这里从 decaytime>1 ms 开始拟合。
TF1 *fdecay1 = new TF1("fdecay1",
"[0]+[1]*TMath::Exp(-x*TMath::Log(2.)/[2])",
1, 1.0e4);
fdecay1->FixParameter(0, p01);
fdecay1->SetParameter(1, 600);
fdecay1->SetParameter(2, 250); // ms
hdtAc211a->Fit(fdecay1, "RL");
cout << endl;
cout << "half-life of 211Ac is "
<< fdecay1->GetParameter(2)
<< " +/- "
<< fdecay1->GetParError(2)
<< " ms" << endl;
c1->SetLogy();
c1->Draw();
**************************************** Minimizer is Minuit2 / Migrad MinFCN = 52.6435 Chi2 = 105.287 NDf = 98 Edm = 1.4813e-09 NCalls = 36 p0 = 1.77987 (fixed) p1 = 608.163 +/- 18.9374 p2 = 253.085 +/- 5.91764 half-life of 211Ac is 253.085 +/- 5.91764 ms
通过该拟合可以得到 $^{211}\mathrm{Ac}$ 的半衰期。需要注意，拟合结果依赖于能量门、位置条件、拟合范围和本底估计范围，因此这些条件应在结果报告中明确给出。
$\alpha$ 能谱减本底¶
除了拟合时间谱，也可以利用时间窗对 $\alpha$ 能谱进行本底扣除。
基本思想是：
- 近时间窗：包含真实衰变信号和随机本底。
- 远时间窗：主要包含随机本底。
- 两者相减后得到衰变信号的近似净能谱。
对于半衰期约为 $250$ ms 的 $^{211}\mathrm{Ac}$，可取：
- 信号区：$0<t<5T_{1/2}$。
- 本底区：$20T_{1/2}<t<25T_{1/2}$。
两个时间窗宽度相同，因此可以直接相减。若时间窗宽度不同，则需要按时间窗宽度归一化。
tree->Draw("be>>hbeall(300,5000,8000)",
cdxdy00 && "decaytime>0 && decaytime<5*250");
tree->Draw("be>>hbebkg(300,5000,8000)",
cdxdy00 && "decaytime>20*250 && decaytime<25*250");
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();
hbebkg->Draw("same");
hbenet->Draw("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,
"colz");
c1->SetLogy();
c1->Draw();
用负时间区估计平台本底
TF1 *fp01d = new TF1("fp01d", "pol0", -0.8e4, -0.3e4);
hdtAc211d->Fit(fp01d, "LR");
Double_t p01d = fp01d->GetParameter(0);
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]*TMath::Exp(-x*TMath::Log(2.)/[2])",
1, 1.0e4);
fdecay2->FixParameter(0, p01d);
fdecay2->SetParameter(1, 600);
fdecay2->SetParameter(2, 250); // ms
hdtAc211d->Fit(fdecay2, "RL");
hdtAc211d->SetMaximum(1000);
cout << endl;
cout << "half-life of 211Ac is "
<< fdecay2->GetParameter(2)
<< " +/- "
<< fdecay2->GetParError(2)
<< " ms" << endl;
c1->SetLogy();
c1->Draw();
**************************************** Minimizer is Minuit2 / Migrad MinFCN = 52.6433 Chi2 = 105.287 NDf = 98 Edm = 4.05138e-09 NCalls = 36 p0 = 1.57993 (fixed) p1 = 606.06 +/- 18.8425 p2 = 254.538 +/- 5.93313 half-life of 211Ac is 254.538 +/- 5.93313 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");
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");
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();
hbebkgn->Draw("same");
hbenetn->Draw("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 测量，每一级 $\alpha$ 的几何探测效率约为 $50\%$，两级 $\alpha$ 衰变同时被探测到的效率约为 $25\%$。
在衰变时间窗内，将满足条件的所有衰变事件存入数组¶
在前面生成 decay.root 时，已经将正时间、同一 pixel 的候选衰变事件存入：
bphitbpts[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}$ 的多个半衰期。当前示范中时间范围仍偏短，因此只能用于观察关联，不能可靠拟合 $^{207}\mathrm{Fr}$ 半衰期。
tree->Draw("bphit", "bphit>0");
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,
"colz");
c1->SetLogy();
c1->Draw();
由于当前时间范围小于或接近 $^{207}\mathrm{Fr}$ 半衰期的有效拟合要求，不能通过该谱得到可靠半衰期。这里的目的主要是展示两级 $\alpha$ 衰变之间的时间和能量关联。
作业：¶
- 求210Ra的半衰期。