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

6.1 衰变实验：基于 timestamp 的事件重构

学习目的：

1. 理解自触发数据中基于 timestamp 的事件重构方法。
2. 掌握 DSSD 正背面信号、相邻条信号和 MWPC 信号的基本关联方法。
3. 为后续注入事件与衰变事件的 position-time correlation 建立输入数据。

Experiment:

• 188.2 MeV $^{40}\mathrm{Ar}+^{175}\mathrm{Lu}$ fusion-evaporation @ SHANS, IMP
• 测量剩余核的 $\alpha$ 衰变：
  • 能量
  • 半衰期
  • 级联关系

实验装置：

Detectors

实验探测系统主要由 MWPC、DSSD、SSD 和 Si-BOX 组成。各探测器在实验中的作用不同：MWPC 主要标记重离子注入，DSSD 测量注入和衰变，SSD 与 Si-BOX 用于辅助排除或补充测量部分事件。

MWPC

• MWPC 放置在 DSSD 前方，用于探测经过 SHANS 分离后的重离子余核。
• 靶上产生的熔合蒸发反应余核经过 SHANS 分离后，首先穿过 MWPC，然后注入到 DSSD 表面。
• 因此，MWPC 信号主要用于标记重离子注入事件，并提供注入时间信息。

DSSD

• DSSD 是主要的 stop detector，尺寸为 $128 \times 48$ 条。
• 经过 SHANS 分离后的重离子余核最终停止在 DSSD 中。
• DSSD 同时测量两类信号：
  • 重离子注入信号；
  • 注入后发生的 $\alpha$ 衰变或其他带电粒子衰变信号。
• 通过 DSSD 正面和背面条号，可以确定注入和衰变发生的位置。

SSD

• SSD 放置在 DSSD 后方，作为 veto detector。
• 部分束流轻粒子可能穿过 DSSD，而不是停止在 DSSD 内。
• 这类穿透 DSSD 的轻粒子可以被后方 SSD 探测到，因此 SSD 可用于排除穿透事件。

Si-BOX

• Si-BOX 用于测量从 DSSD 表面逃逸出的带电粒子能量。
• 当衰变发生在 DSSD 表面附近时，部分衰变粒子可能不能把全部能量沉积在 DSSD 中，而是从 DSSD 逃逸出来。
• Si-BOX 可以补充测量这部分逃逸能量。

需要注意的是，MWPC 主要用于探测重离子余核，而不是束流中的轻粒子。轻粒子质量小、电荷低，在 MWPC 气体中的能量损失较小，产生的电离信号弱，通常低于有效触发或鉴别阈值。因此，这部分轻粒子可能不能被 MWPC 有效探测；若其能量足够高，则可能穿过 DSSD，并被后方 SSD 探测到。

本课件使用的数据只包含 MWPC 和 DSSD 信息，未包含 SSD 和 Si-BOX 数据。

参考文献：

• Nuclear Instruments and Methods in Physics Research B 317 (2013) 315–318

DAQ:

• CAEN V1724 100 MHz, 14 bit digitizer

Trigger: 自触发模式

每个通道独立判断自身信号是否超过阈值；一旦触发，就记录该通道信号的能量和 timestamp。

这种采集方式不会在数据记录阶段直接形成完整物理事件。DSSD 正面、背面和 MWPC 的信号最初都是独立 hit，需要在离线分析中根据 timestamp 重新组合。

Data:

本节使用两个 ROOT 文件：

alpha.root

• Branch：
  • energy / keV
  • timestamp / ns
  • side: 0-front, 1-back
  • strip

mwpc.root

• Branch：
  • energy / ch
  • timestamp / ns

衰变实验的目标是研究剩余核注入到 DSSD 后产生的 $\alpha$ 衰变。对于一个完整的物理事件，通常需要同时知道粒子在 DSSD 正面和背面的条号、能量和时间，并判断该事件是否与 MWPC 有 prompt 符合。

但是在自触发采集方式中，每一个通道独立触发并独立记录数据。因此，原始 ROOT tree 中的一条 entry 只是某一个探测器通道的一个 hit，而不是一个完整的物理事件。本节的任务就是根据 timestamp，将这些分散的 hit 重新组织成后续衰变分析可以使用的事件结构。

alpha.root 中 DSSD 能量单位为 keV，mwpc.root 的幅度仍为 channel；两者 timestamp 的单位均为 ns。这里直接使用已刻度的 DSSD 能量，MWPC 幅度仅用于符合标记，不另作连续能量刻度。整数道值的 dithering 应在相应通道的连续刻度之前做一次，不在这些下游关联步骤中重复。

%jsroot on

TFile *fds = new TFile("alpha.root");
TTree *tree = (TTree*)fds->Get("tree");

Double_t energy;
ULong64_t timestamp;
Int_t side, strip;

tree->SetBranchAddress("energy", &energy);
tree->SetBranchAddress("timestamp", &timestamp);
tree->SetBranchAddress("side", &side);
tree->SetBranchAddress("strip", &strip);
TH1::AddDirectory(kFALSE); // 关闭文件后仍保留本节创建的检查谱
tree->Print();

******************************************************************************
*Tree    :tree      : dssd                                                   *
*Entries :   269955 : Total =         6499612 bytes  File  Size =    4793737 *
*        :          : Tree compression factor =   1.36                       *
******************************************************************************
*Br    0 :energy    : energy/D                                               *
*Entries :   269955 : Total  Size=    2166281 bytes  File Size  =    2163446 *
*Baskets :       68 : Basket Size=      32000 bytes  Compression=   1.00     *
*............................................................................*
*Br    1 :timestamp : timestamp/l                                            *
*Entries :   269955 : Total  Size=    2166497 bytes  File Size  =    1904904 *
*Baskets :       68 : Basket Size=      32000 bytes  Compression=   1.14     *
*............................................................................*
*Br    2 :strip     : strip/I                                                *
*Entries :   269955 : Total  Size=    1083253 bytes  File Size  =     518726 *
*Baskets :       34 : Basket Size=      32000 bytes  Compression=   2.09     *
*............................................................................*
*Br    3 :side      : side/I                                                 *
*Entries :   269955 : Total  Size=    1083215 bytes  File Size  =     203601 *
*Baskets :       34 : Basket Size=      32000 bytes  Compression=   5.32     *
*............................................................................*

正/背面的能量分布

首先分别观察 DSSD 正面和背面的能量分布。

TCanvas *c2 = new TCanvas("c2","c2",1000,400);
c2->Divide(2,1);
c2->cd(1);
tree->Draw("energy>>h1(3000,0,30000)","side==0","hist");
c2->cd(2);
tree->Draw("energy>>h2(3000,0,30000)","side==1","hist");
c2->Draw();

DSSD 事件结构

原始数据中的一个 entry 只对应 DSSD 某一面某一条的一个 hit，而不是一个已经完成正背面匹配的 DSSD 事件。后续分析需要根据 timestamp 将正面和背面的 hit 关联起来，从而得到带电粒子在 DSSD 中的能量、位置和时间信息。

在自触发采集方式下，单通道 hit 先写入前端电子学的 buffer，再传输到后端数据流中。由于不同通道、不同插件的 buffer 独立工作，数据写入 ROOT tree 的顺序不一定等同于粒子实际到达 DSSD 的时间顺序。因此，entry 编号不能直接作为事件顺序使用。

直接检查原始数据可以发现，timestamp 并不随 entry 编号增加而单调增加：

tree->Scan("side:strip:timestamp:int(energy)","","colsize=15",10,1);

************************************************************************************
*    Row   *            side *           strip *       timestamp *     int(energy) *
************************************************************************************
*        1 *               0 *              14 *   2484394953520 *            4293 *
*        2 *               1 *             110 *   3553975013340 *            6755 *
*        3 *               1 *              41 *   1149167130120 *            2165 *
*        4 *               0 *              15 *   3469034044010 *            6763 *
*        5 *               1 *              14 *    217811466190 *             994 *
*        6 *               1 *             102 *   2051660598330 *            5173 *
*        7 *               0 *              14 *   1928943248620 *            3068 *
*        8 *               1 *             106 *   2414539135320 *            5698 *
*        9 *               1 *             102 *   2897908302230 *            7489 *
*       10 *               1 *              84 *   3064794615050 *           10006 *
************************************************************************************

正面和背面事件分别按照时间增加的顺序排序

分别定义正面和背面信号的 multimap 结构。将 timestamp 作为 key，将条号和能量作为 value。当数据插入 multimap 后，数据会自动按照 key 的大小排序。 使用 multimap，是因为不同通道可能具有相同或非常接近的 timestamp，不能假定 timestamp 是唯一的。

%%cpp -d
struct dssd
{
    Double_t energy;
    Int_t strip;
};

multimap<ULong64_t, dssd> mapf1, mapb; // map for front and back side

dssd ds;

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

    ds.energy = energy;
    ds.strip = strip;
    if(side==0) mapf1.insert(make_pair(timestamp, ds));
    if(side==1) mapb.insert(make_pair(timestamp, ds));
    if(jentry%100000==0) cout << jentry << endl;
}

cout << "Total number of frontside/backside events = "
     << mapf1.size() << " / " << mapb.size() << endl;

fds->Close();

0
100000
200000
Total number of frontside/backside events = 111442 / 158513

查看事件的时间排序结果

排序后，可以比较正面和背面的前若干个 hit。需要注意，此时只是分别按时间排序，并没有建立正背面对应关系。 如果简单地把正面第 $i$ 个 hit 和背面第 $i$ 个 hit 配对，通常不会得到正确的物理事件。因为正面和背面 hit 的数目不同，触发阈值不同，噪声和随机事件也不同。

TString sout;

auto ifts = mapf1.begin();
auto ibts = mapb.begin();

TH2F *hxyraw = new TH2F("hxyraw",
                        "front-back correlation for unmatched data",
                        3000,0,30000,3000,0,30000);

sout.Form("%4s%15s%15s%5s%5s%5s%5s",
          "i","front-ts","back-ts",
          "     front-strip", "    back-strip",
          "    front-energy", "    back-energy");
cout << sout << endl;
for(int i=0; i<TMath::Min(mapf1.size(),mapb.size()); i++) {
    if(i<20) {
        sout.Form("%2d  %15llu  %15llu %10d    %10d   %10.1f    %10.1f",
                  i,
                  ifts->first,
                  ibts->first,
                  ifts->second.strip,
                  ibts->second.strip,
                  ifts->second.energy,
                  ibts->second.energy);
        cout << sout << endl;
    }

    hxyraw->Fill(ifts->second.energy, ibts->second.energy);

    ifts++;
    ibts++;
}

   i       front-ts        back-ts     front-strip    back-strip    front-energy    back-energy
 0       1339474350       1243591410         11            97       6141.5         609.7
 1       9022757640       1339474900          3            94      10898.4        6140.3
 2       9122771170       9022758220         24            48      18331.7       12394.7
 3       9202986900       9122771720         31           121      13594.9       18350.6
 4       9305276750       9202987450         27            76      10988.4        3868.6
 5       9418690000       9202987450         29            77      16461.8        9703.6
 6       9423589230       9305277310         37            91      18753.1       10993.2
 7       9441338680       9418690560         43            79      23666.8       16449.1
 8       9470639040       9423589790          8           100      15354.3       18779.6
 9       9531671070       9441339240         22            94      15308.4       23665.6
10       9655957710       9445831200         26            42       7575.8        1169.2
11       9664140820       9449732290         46            74       4508.0         813.0
12       9703733650       9470639600         40           124      13927.8       15331.8
13       9765656470       9492181880         40             4       7464.1        2167.5
14       9974541640       9531671620         34           108      15086.1       15325.0
15       9977849690       9534153370         28           111       5240.2        1001.0
16       9994473120       9655958280         33            67      10240.3        7563.0
17      10017904000       9664141380         11            94      10032.1        4496.9
18      10024458920       9703734210          7           127       2146.2       13927.4
19      10043706070       9765657030         37           127       6581.2        7476.1

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

从图中可以看到，直接按照排序后的序号建立正背面能量关联，并不能得到清晰的对角线结构。这说明这种配对方式不是正确的物理事件重构方法。 对于同一个带电粒子，它在 DSSD 正面和背面应几乎同时产生信号。因此，正确的正背面匹配应当基于 timestamp 的时间邻近关系，而不是 entry 编号或排序后的序号。

寻找正面和背面的符合时间窗

当一个粒子同时在 DSSD 正面和背面产生信号时，两个信号之间应具有稳定的时间关系。由于电子学延迟、信号形成时间和时间提取方法不同，正面和背面的 timestamp 不一定完全相同，但时间差应集中在一个窄峰内。对于不相关的 hit，正面和背面的时间差是随机的，因此会形成较平缓的随机本底。

为了确定正背面的符合时间窗，先在一个较宽的时间范围内搜索背面 hit。对每一个正面 hit，在背面 timestamp 表中寻找

$$ [t_\mathrm{front}-T,\ t_\mathrm{front}+T] $$

范围内的所有背面 hit，并把每一个组合的时间差填入直方图：

$$ \Delta t = t_\mathrm{front}-t_\mathrm{back}. $$

这里使用 lower_bound() 和 upper_bound() 在已排序的 multimap 中快速寻找时间范围内的 hit：

• map::lower_bound(key)：返回 map 中第一个大于或等于 key 的迭代器。
• map::upper_bound(key)：返回 map 中第一个大于 key 的迭代器。

窗口下边界先限制到 0，避免 unsigned timestamp 减去较大窗口后下溢。时间差则分别转为 Long64_t 后相减。

ULong64_t twindow = 100000; // ns

TH1I *hdt = new TH1I("hdt",
                     "fts-bts distribution(ns)",
                     20000,-100000,100000);

// time window: [ts-twindow, ts+twindow]
for(auto ifts=mapf1.begin(); ifts!=mapf1.end(); ifts++) {
    auto ibts1 = mapb.lower_bound((ifts->first > twindow ? ifts->first - twindow : 0));
    auto ibts2 = mapb.upper_bound(ifts->first + twindow);
    for( ; ibts1 != ibts2; ++ibts1) {
        Long64_t dt = Long64_t(ifts->first) - Long64_t(ibts1->first);
        hdt->Fill(dt);
    }
}

hdt->Draw("hist");
c1->SetLogy();
c1->Draw();
hdt->GetXaxis()->SetRangeUser(-1000,500);
hdt->Draw("hist");
c1->Draw();

从图中可以看到，正背面符合峰的中心约为 $-550$ ns。这个 offset 不是物理常数，而是由当前电子学延迟、通道设置和 timestamp 定义共同决定的实验参数。为了后续用以 0 为中心的时间窗做正背面关联，需要对正面 timestamp 进行平移。

本节选取正背面符合时间窗为峰中心附近约 $\pm 100$ ns。

正面 timestamp 减去 offset，将正背面时间位置对齐

原始时间差峰中心在

$$ t_\mathrm{front}-t_\mathrm{back}\approx -550~\mathrm{ns}. $$

因此，将正面 timestamp 修正到与背面 timestamp 对齐：

$$ t'_\mathrm{front}=t_\mathrm{front}-(-550~\mathrm{ns}). $$

代码中涉及负 offset 时，应使用有符号整数，避免把负数赋给 ULong64_t。

Long64_t toffset = -550;

multimap<ULong64_t,dssd> mapf;
for(auto ifts=mapf1.begin(); ifts!=mapf1.end(); ifts++) {
    ULong64_t tcorr = ULong64_t(Long64_t(ifts->first) - toffset);
    mapf.insert(pair<ULong64_t,dssd>(tcorr, ifts->second));
}
cout << mapf.size() << endl;

111442

在符合时间窗内观察正背面关联

对齐正背面时间后，在 $\pm 100$ ns 时间窗内寻找正背面符合事件，并重新画正背面能量关联图。

TH2F *hxyc = new TH2F("hxyc",
                      "front-back correlation within coincidence window",
                      3000,0,30000,3000,0,30000);
twindow = 100; // ns
for(auto ifts=mapf.begin(); ifts!=mapf.end(); ifts++) {
    auto ibts1 = mapb.lower_bound((ifts->first > twindow ? ifts->first - twindow : 0));
    auto ibts2 = mapb.upper_bound(ifts->first + twindow);
    for( ; ibts1 != ibts2; ++ibts1) {
        Long64_t dt = Long64_t(ifts->first) - Long64_t(ibts1->first);
        hxyc->Fill(ifts->second.energy, ibts1->second.energy);
    }
}

hxyc->Draw("colz");
c1->SetLogy(0);
c1->SetLogz();
c1->Draw();

此时可以看到明显的正背面能量关联。对于同一个带电粒子，正面和背面测量的是同一粒子在同一片硅中的能量沉积，因此正确匹配事件应接近

$$ E_\mathrm{front}=E_\mathrm{back}. $$

图中偏离对角线的事件，一部分来自错误匹配和随机符合，另一部分来自相邻条之间的能量共享。

背面相邻条之间事件关联

带电粒子入射到 DSSD 条之间的边界附近时，产生的电荷可能被相邻两条共同收集。此时单独使用其中一条的能量会低估粒子能量，并在正背面能量关联图中形成偏离对角线的事件。 因此，在建立最终 DSSD 事件前，需要检查同一面相邻条之间是否存在能量共享。先从当前 hit 开始，在 200 ns 范围内寻找同一面其他 hit。由于 map 中事件已经按 timestamp 排序，可以顺序扫描并统计局部时间窗内的多重性。

Int_t nhit;

TH1I *hnhit = new TH1I("hnhit",
                       "back-side hit multiplicity within 200 ns;Multiplicity;Counts",
                       6,-0.5,5.5);
twindow = 200; // ns

ifts = mapb.begin();
while(ifts != mapb.end()) {
    nhit = 0;
    ULong64_t t0 = ifts->first;
    while(ifts != mapb.end()) {
        if(ifts->first > t0 + twindow) break;
        nhit++;
        ifts++;
    }

    hnhit->Fill(nhit);
}

hnhit->Draw("hist");
gPad->SetLogy();
c1->Draw();

观察背面相邻条关联

根据多重性分布，选择两重事件作为相邻条能量共享的主要候选。正面的多重性分布和时间差分布与背面类似，因此这里只以背面为例说明。

Int_t xhit, yhit;

Int_t xstrip[2], ystrip[2];
ULong64_t xtimestamp[2], ytimestamp[2];
Float_t xenergy[2], yenergy[2];

TH2F *hx2 = new TH2F("hx2",
                     " x0-x1 energy correlation",
                     500,0,10000,500,0,10000);

TH1I *hdtx = new TH1I("hdtx"," tx0-tx1 ",20,-150,50);
ifts = mapb.begin();
while(ifts != mapb.end()) {
    xhit = 0;
    ULong64_t t0 = ifts->first;
    while(ifts != mapb.end()) {
        if(ifts->first > t0 + twindow) break;
        if (xhit < 2) { // 只存前两条，但遍历完整时间窗
        xtimestamp[xhit] = ifts->first;
        xstrip[xhit] = ifts->second.strip;
        xenergy[xhit] = ifts->second.energy;
        }

        xhit++;
        ifts++;
    }
    if(xhit == 2) {
        hx2->Fill(xenergy[0], xenergy[1]);

        Long64_t dt = Long64_t(xtimestamp[0]) - Long64_t(xtimestamp[1]);
        hdtx->Fill(dt);
    }
}
c1->Clear();
hx2->Draw("colz");
c1->SetLogz(0);
c1->SetLogy(0);
c1->Draw();
hdtx->Draw("hist");
c1->SetLogy();
c1->Draw();

同一面不同条之间的时间差集中在约 $\pm 100$ ns 范围内，说明这些两重事件中存在真实的相邻条能量共享成分。

合并正背面相邻条的两重事件，忽略三重事件

本节采用简单处理：

• 对同一面 200 ns 内的单重事件，直接保留。
• 对同一面 200 ns 内的两重事件：
  • 若两个 strip 相邻，则认为可能来自同一粒子的相邻条能量共享；
  • 将两条能量相加；
  • 使用能量较大的一条作为位置代表；
  • timestamp 也采用能量较大的一条的 timestamp。
• 若两重事件的 strip 不相邻（或重复同一条），本节不将其作为单粒子候选写入下一阶段，数量单独统计。
• 三重及更高多重事件在本节中忽略。

这种处理保持了算法的简单性，适合当前数据的基本重构。它的近似条件是：真实能量共享主要发生在相邻条之间，并且共享信号均超过阈值。若其中一条低于阈值，则能量无法完全恢复。

内层循环遍历完整的 200 ns 时间窗，再按总 hit 数分类。不能读到两个 hit 就结束，否则一个三重组会被错误地拆成“双重 + 单重”。

multimap<ULong64_t,dssd> mapfn, mapbn;

ULong64_t t0;

int mby[4] = {0};
int mfx[4] = {0};

先处理背面：

twindow = 200; // ns
ibts = mapb.begin();
while (ibts != mapb.end()) {
    yhit = 0;
    t0 = ibts->first;
    while (ibts != mapb.end() && ibts->first <= t0+twindow) {
        if (yhit < 2) { // 只存前两条，但完整数完本时间窗中的 hit
            ytimestamp[yhit] = ibts->first;
            ystrip[yhit] = ibts->second.strip;
            yenergy[yhit] = ibts->second.energy;
        }
        ++yhit;
        ++ibts;
    }
    if (yhit > 2) { ++mby[3]; continue; }
    if (yhit == 2 && abs(ystrip[0]-ystrip[1]) != 1) {
        ++mby[2]; // 两条不相邻，不能按相邻条能量共享处理
        continue;
    }
    int peak = 0;
    ds.energy = yenergy[0];
    if (yhit == 2) {
        ds.energy = yenergy[0] + yenergy[1];
        peak = (yenergy[0] > yenergy[1]) ? 0 : 1;
    }
    ds.strip = ystrip[peak]; // 较大能量条给出位置和时间
    mapbn.emplace(ytimestamp[peak], ds);
    ++mby[yhit-1]; // [0] 单条，[1] 相邻双条
}
cout << "single / adjacent double / other double / >=3 groups: "
     << mby[0] << " / " << mby[1] << " / "
     << mby[2] << " / " << mby[3] << endl;

single / adjacent double / other double / >=3 groups: 126371 / 15935 / 106 / 17

再处理正面：

twindow = 200; // ns
ifts = mapf.begin();
while (ifts != mapf.end()) {
    xhit = 0;
    t0 = ifts->first;
    while (ifts != mapf.end() && ifts->first <= t0+twindow) {
        if (xhit < 2) { // 只存前两条，但完整数完本时间窗中的 hit
            xtimestamp[xhit] = ifts->first;
            xstrip[xhit] = ifts->second.strip;
            xenergy[xhit] = ifts->second.energy;
        }
        ++xhit;
        ++ifts;
    }
    if (xhit > 2) { ++mfx[3]; continue; }
    if (xhit == 2 && abs(xstrip[0]-xstrip[1]) != 1) {
        ++mfx[2]; // 两条不相邻，不能按相邻条能量共享处理
        continue;
    }
    int peak = 0;
    ds.energy = xenergy[0];
    if (xhit == 2) {
        ds.energy = xenergy[0] + xenergy[1];
        peak = (xenergy[0] > xenergy[1]) ? 0 : 1;
    }
    ds.strip = xstrip[peak]; // 较大能量条给出位置和时间
    mapfn.emplace(xtimestamp[peak], ds);
    ++mfx[xhit-1]; // [0] 单条，[1] 相邻双条
}
cout << "single / adjacent double / other double / >=3 groups: "
     << mfx[0] << " / " << mfx[1] << " / "
     << mfx[2] << " / " << mfx[3] << endl;

single / adjacent double / other double / >=3 groups: 109581 / 902 / 21 / 5

相邻条修正完成后的正背面关联

相邻条修正后，再次进行正背面时间符合，并画出正背面能量关联图。

先画时间窗内全部候选的能量关联，再应用 |Efront−Eback|<500 keV。写入事件时，还要求该 front 和 back 在双方的候选列表中唯一对应，避免同一个 hit 被重复用于多个 DSSD 事件。歧义候选单独计数；高计数率数据需进一步用第 3 章的多粒子重建方法处理。

本文件约定 side=0 为 front（48 条），side=1 为 back（128 条）；输出 xstrip 对应 front，ystrip 对应 back。位置图沿用这一顺序，每条占一个 bin。

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

TH2F *hxy = new TH2F("hxy",
                     "front-back correlation for corrected data",
                     3000,0,30000,3000,0,30000);

TH2F *hxyfbc = new TH2F("hxyfbc",
                        "front-back correlation after energy-consistency cut",
                        3000,0,30000,3000,0,30000);

TH2I *hxystrip = new TH2I("hxystrip",
                          "front-back position;front strip;back strip",
                          48,-0.5,47.5,128,-0.5,127.5);

dssdxy xy;

multimap<ULong64_t, dssdxy> mapdssd;
Long64_t ambiguousPairs = 0;

twindow = 200; // ns

先统计双方的候选数，再保留唯一匹配。全部候选图仍在施加唯一性条件之前填充。

for(ifts=mapfn.begin(); ifts!=mapfn.end(); ++ifts) {
    ULong64_t tmin = (ifts->first > twindow) ? ifts->first - twindow : 0;
    ULong64_t tmax = ifts->first + twindow;

    auto ibts1 = mapbn.lower_bound(tmin);
    auto ibts2 = mapbn.upper_bound(tmax);

    // 同一 front 的 back 候选数不随当前 pair 改变，只统计一次。
    int nback = 0;
    for (auto ib=ibts1; ib!=ibts2; ++ib)
        if (abs(ifts->second.energy-ib->second.energy)<500) ++nback;
    for( ; ibts1 != ibts2; ++ibts1) {
        hxy->Fill(ifts->second.energy, ibts1->second.energy);
        if(TMath::Abs(ifts->second.energy - ibts1->second.energy) < 500) {
            // 两个方向都只有一个能量相容的候选，才写入一个 DSSD 事件。
            int nfront = 0;
            ULong64_t lo = ibts1->first>twindow ? ibts1->first-twindow : 0;
            auto fend = mapfn.upper_bound(ibts1->first+twindow);
            for (auto jf=mapfn.lower_bound(lo); jf!=fend; ++jf)
                if (abs(jf->second.energy-ibts1->second.energy)<500) ++nfront;
            if (nback!=1 || nfront!=1) { ++ambiguousPairs; continue; }
            hxyfbc->Fill(ifts->second.energy, ibts1->second.energy);

            xy.xstrip = ifts->second.strip;
            xy.ystrip = ibts1->second.strip;

            // The back-side energy is used as the DSSD event energy here.
            // After the front-back energy-consistency cut, the two energies
            // should be close to each other.
            xy.energy = ibts1->second.energy;

            mapdssd.insert(pair<ULong64_t, dssdxy>(ibts1->first, xy));
            hxystrip->Fill(xy.xstrip, xy.ystrip);
        }
    }
}
cout << "Unique DSSD events = " << mapdssd.size()
     << ", ambiguous candidate pairs = " << ambiguousPairs << endl;

Unique DSSD events = 107133, ambiguous candidate pairs = 0

c2->Clear();
c2->Divide(2,1);
c2->cd(1);
hxy->Draw("colz");
c2->cd(2);
hxyfbc->Draw("colz");
c2->Draw();

c1->Clear();
c1->SetLogx(0);
c1->SetLogy(0);
hxystrip->SetStats(0);
hxystrip->Draw("colz");
c1->Draw();
cout << "Position plot: entries / regular-bin sum = "
     << hxystrip->GetEntries() << " / " << hxystrip->Integral() << endl;

Position plot: entries / regular-bin sum = 107133 / 107133

对比相邻条修正前后的正背面能量关联图，可以看到杂散点明显减少。说明部分偏离对角线的事件确实来自相邻条能量共享。 但是图中仍可能保留少量相邻条关联的迹象。这可能来自探测器阈值：如果相邻条中的小幅度信号低于阈值，则该部分能量无法被恢复。

正背面能量关联条件为

$$ |E_\mathrm{front}-E_\mathrm{back}|<500~\mathrm{keV}. $$

熔合蒸发反应余核的 $\alpha$ 衰变产物

应用正背面能量关联条件后，可以得到 DSSD 中较干净的带电粒子能谱。

c1->Clear();
TH1F *hfe = new TH1F("hfe","Reconstructed DSSD;Energy (keV);Counts / 10 keV",3000,0,30000);
for (const auto &hit : mapdssd) hfe->Fill(hit.second.energy);
hfe->Draw("hist");
c1->SetLogy(0);
c1->Draw();

到这里，DSSD 正面和背面的信号已经被重构成带有位置、能量和 timestamp 的 DSSD 事件。下一步需要判断这些 DSSD 事件是否与 MWPC 有 prompt 符合。

MWPC

MWPC 用于标记穿过分离器并进入 DSSD 的重离子注入事件。重离子注入时，MWPC 和 DSSD 应在 prompt 时间范围内同时产生信号。 而 $\alpha$ 衰变发生在余核注入之后，通常不应与 MWPC 同时符合。因此，MWPC-DSSD prompt 符合主要对应注入事件；无 MWPC prompt 符合的 DSSD 事件则主要包含后续衰变候选事件。MWPC 数据也需要先按 timestamp 排序，再寻找与 DSSD 的符合时间窗。

TFile *fmw = new TFile("mwpc.root");
TTree *tmw = (TTree*)fmw->Get("tree");

tmw->SetBranchAddress("energy",&energy);
tmw->SetBranchAddress("timestamp",&timestamp);

tmw->Print();
tmw->Scan("timestamp:energy","","colsize=15",10,1);

******************************************************************************
*Tree    :tree      : mwpc                                                   *
*Entries :   353117 : Total =         5667688 bytes  File  Size =    3824341 *
*        :          : Tree compression factor =   1.48                       *
******************************************************************************
*Br    0 :energy    : energy/D                                               *
*Entries :   353117 : Total  Size=    2833530 bytes  File Size  =    1324555 *
*Baskets :       89 : Basket Size=      32000 bytes  Compression=   2.14     *
*............................................................................*
*Br    1 :timestamp : timestamp/l                                            *
*Entries :   353117 : Total  Size=    2833809 bytes  File Size  =    2497325 *
*Baskets :       89 : Basket Size=      32000 bytes  Compression=   1.13     *
*............................................................................*
************************************************
*    Row   *       timestamp *          energy *
************************************************
*        1 *   1967773996860 *           15365 *
*        2 *   1018727462730 *           15699 *
*        3 *   2771711086400 *           15097 *
*        4 *    934774734690 *              24 *
*        5 *   2435449424700 *           16088 *
*        6 *    422204325350 *            2796 *
*        7 *   3178976019570 *            3133 *
*        8 *   2887172375060 *            7291 *
*        9 *   1924594027780 *           12244 *
*       10 *   2430242405260 *            3944 *
************************************************

将 MWPC 记录按 timestamp 排序，然后查看幅度谱。

multimap<ULong64_t, Float_t> mapmw1; // map for mwpc

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

    mapmw1.insert(pair<ULong64_t,Float_t>(timestamp, energy));
    if(jentry%100000==0) cout << jentry << endl;
}

cout << "Total number of mwpc events = " << mapmw1.size() << endl;

fmw->Close();
auto imw = mapmw1.begin();

TH1F *hmwe = new TH1F("hmwe","energy of mwpc",3500,0,35000);
for(int i=0; i<mapmw1.size(); i++) {
    if(i<10) {
        sout.Form("%2d  %15llu   %5.1f",
                  i, imw->first, imw->second);
        cout << sout << endl;
    }

    hmwe->Fill(imw->second);
    imw++;
}

hmwe->Draw("hist");
c1->SetLogy();
c1->Draw(); // threshold: 1000

0
100000
200000
300000
Total number of mwpc events = 353117
 0         78196010     3.0
 1        259521960   32763.0
 2        574208200   32726.0
 3        750230460    24.0
 4       1086207510    15.0
 5       1363535160    24.0
 6       1364841100   32749.0
 7       1426179500     5.0
 8       1495511510    34.0
 9       1731485990     8.0

MWPC 与 DSSD 符合时间窗

与正背面关联类似，先用较宽时间窗观察 MWPC 与 DSSD 的时间差分布：

$$ \Delta t = t_\mathrm{MWPC}-t_\mathrm{DSSD}. $$

twindow = 5000; // ns

TH1I *hmdt = new TH1I("hmdt",
                      "mts-fts distribution(ns)",
                      1000,-5000,5000);
for(auto idts=mapdssd.begin(); idts!=mapdssd.end(); idts++) {
    auto imts1 = mapmw1.lower_bound((idts->first > twindow ? idts->first - twindow : 0));
    auto imts2 = mapmw1.upper_bound(idts->first + twindow);
    for( ; imts1 != imts2; ++imts1) {
        Long64_t dt = Long64_t(imts1->first) - Long64_t(idts->first);
        hmdt->Fill(dt);
    }
}

hmdt->Draw("hist");
c1->SetLogy();
c1->Draw();
hmdt->GetXaxis()->SetRangeUser(400,1000);
hmdt->Draw("hist");
c1->Draw();

从图中可以看到，MWPC-DSSD 符合峰中心约为 $650$ ns。下面以峰中心附近 $\pm 200$ ns 作为 MWPC-DSSD 符合窗口。

MWPC 的 timestamp 减去 offset，与 DSSD 时间对齐

toffset = 650;

multimap<ULong64_t, Float_t> mapmw;
for(auto imts=mapmw1.begin(); imts!=mapmw1.end(); imts++) {
    ULong64_t tcorr = ULong64_t(Long64_t(imts->first) - toffset);
    mapmw.insert(pair<ULong64_t,Float_t>(tcorr, imts->second));
}

将 DSSD 和 MWPC 事件存入 ROOT 文件

最终生成 sort1.root。其中每个 entry 对应一个已经完成正背面匹配和相邻条修正的 DSSD 事件。

输出变量包括：

• timestamp：DSSD 事件 timestamp
• de：DSSD 能量
• xstrip：DSSD 正面条号
• ystrip：DSSD 背面条号
• me：MWPC 能量
  • 若 DSSD 事件与 MWPC 在 prompt 时间窗内符合，则 me>0
  • 若无 MWPC prompt 符合，则记为 me=-1

这里 me>0 的事件主要对应重离子注入事件，me<0 的事件主要作为后续衰变候选事件。真正的注入-衰变关联将在下一节通过位置和时间关联进一步建立。

若窗口内有多个 MWPC hit，me 取时间差最小的一路；一个 DSSD 事件只填一次符合能谱。这里的分类依据是有没有 prompt MWPC 信号，不把最邻近关系当作唯一的物理身份判定。

TH1F *hmwec = new TH1F("hmwec",
                       "energy of mwpc coincided with DSSD ",
                       3500,0,35000);

TH1F *hdsec = new TH1F("hdsec",
                       "energy of DSSD coincided with MWPC ",
                       3000,0,30000);

TH1F *hdsenc = new TH1F("hdsenc",
                        "energy of DSSD no coincided with MWPC ",
                        3000,0,30000);

TFile *ff = new TFile("sort1.root","recreate");
TTree *t1 = new TTree("tree","sorted events");

Float_t de, me;
Int_t x, y;
ULong64_t tt;

t1->Branch("timestamp",&tt,"timestamp/l");
t1->Branch("me",&me,"me/F");
t1->Branch("de",&de,"de/F");
t1->Branch("xstrip",&x,"xstrip/I");
t1->Branch("ystrip",&y,"ystrip/I");
twindow = 200; // ns

逐个 DSSD 事件寻找时间窗内最近的 MWPC 信号；每个 DSSD 事件只写一次。

for(auto idts=mapdssd.begin(); idts!=mapdssd.end(); idts++) {
    tt = idts->first;
    x = idts->second.xstrip;
    y = idts->second.ystrip;
    de = idts->second.energy;
    me = -1;

    auto imts1 = mapmw.lower_bound((idts->first > twindow ? idts->first - twindow : 0));
    auto imts2 = mapmw.upper_bound(idts->first + twindow);

    Long64_t nearest = twindow+1;
    for( ; imts1 != imts2; imts1++) {
        Long64_t dt = llabs(Long64_t(imts1->first)-Long64_t(tt));
        if (dt<nearest) { nearest=dt; me=imts1->second; }
    }
    if (me>0) { hmwec->Fill(me); hdsec->Fill(de); }
    if(me < 0) hdsenc->Fill(idts->second.energy);

    t1->Fill();
}

hmwec->Draw("hist");
c1->SetLogy();
c1->Draw();

t1->Write();
ff->Close();
cout << "DSSD total / with MWPC / without MWPC = "
     << hfe->GetEntries() << " / " << hdsec->GetEntries()
     << " / " << hdsenc->GetEntries() << endl;

DSSD total / with MWPC / without MWPC = 107133 / 77672 / 29461

与 DSSD 符合的 MWPC 能谱

hmwe->Draw("hist");
hmwec->SetLineColor(kRed);
hmwec->Draw("hist same");

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

与 DSSD 符合后，MWPC 谱中的部分高道值计数受到明显抑制。可进一步结合模块 overflow 状态判断其来源；仅凭符合谱中的缺失不能确定它们都是无效信号。

与 MWPC 符合的 DSSD 能谱

hfe->Draw("hist"); // dssd only

hdsec->SetLineColor(kGreen);
hdsec->Draw("hist same"); // coincide with MWPC

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

从图中可以看到，明显的 $\alpha$ 峰主要出现在无 MWPC prompt 符合的 DSSD 能谱中。其物理原因是：$\alpha$ 衰变发生在重离子注入之后，不与注入时刻同步。因此，$\alpha$ 衰变事件通常没有 MWPC prompt 符合。这说明 MWPC 可以作为一个有效的注入标记。通过区分有无 MWPC prompt 符合，可以将 DSSD 事件初步分为注入事件和衰变候选事件。

c1->Clear();

hdsenc->Draw("hist");
c1->Draw();

0--2000 keV 范围内的 DSSD 能谱还可能包含低能衰变产物和入射轻粒子，例如 $p$、$d$ 等。这些轻粒子在 MWPC 中能量沉积较小，可能不产生有效 MWPC 信号，但通常有足够能量穿透 DSSD。实际分析中，这部分背景可由 DSSD 后方的其他 Si 探测器进行 veto。本节使用的数据没有包含 SSD 和 Si-BOX 信息，因此这里只能说明这一局限。