作业 5.1：液体闪烁体的脉冲形状甄别¶
本作业使用 $2\times2$ 英寸 BC501A 液体闪烁体探测器测得的钚碳中子源数据。信号由 XIA 500 MS/s、14 bit 数字化采集卡记录，因此相邻采样点间隔为 2 ns。
数据文件 liquidpsd_tree.root 中的 TTree wave 包含 10 000 个 entries，每个 entry 对应一次触发，branch adc[250] 保存该次触发的 250 个 ADC samples。波形保留了基线偏置和触发位置的变化，需要在积分前完成基线修正和时间对齐。
方法¶
液体有机闪烁体中，中子主要通过反冲质子沉积能量，gamma 主要通过电子沉积能量。中子脉冲的慢成分占比通常更大，因此可用脉冲尾部积分与总积分的比值进行脉冲形状甄别（PSD）。
处理时先修正每条波形的基线，再统一脉冲的时间位置，最后计算积分和 PSD。
import ROOT
import math
from array import array
%jsroot on
ROOT.gStyle.SetOptStat(0)
input_file = ROOT.TFile.Open("liquidpsd_tree.root")
tree = input_file.Get("wave")
tree.Print()
pulse_count, sample_count = tree.GetEntries(), 250
adc = array("d", [0.0] * sample_count)
tree.SetBranchAddress("adc", adc)
tree.GetEntry(20) # entry 从 0 编号；读取 event 20
//%jsroot on
gStyle->SetOptStat(0);
TFile* inputFile = TFile::Open("liquidpsd_tree.root");
TTree* tree = (TTree*)inputFile->Get("wave");
tree->Print();
int pulseCount = tree->GetEntries();
const int sampleCount = 250;
Double_t adc[250];
tree->SetBranchAddress("adc", adc);
tree->GetEntry(20); // entry 从 0 编号；读取 event 20
****************************************************************************** *Tree :wave : BC501A digitized waveforms * *Entries : 10000 : Total = 20060726 bytes File Size = 4626311 * * : : Tree compression factor = 4.34 * ****************************************************************************** *Br 0 :adc : adc[250]/D * *Entries : 10000 : Total Size= 20060345 bytes File Size = 4620663 * *Baskets : 667 : Basket Size= 32000 bytes Compression= 4.34 * *............................................................................*
基线修正¶
原始波形叠加了约 1665 ADC 的基线偏置，各次记录的偏置略有不同。积分前要减去每条波形自己的基线，否则偏置会随积分门长度累加到电荷中。
用脉冲到来前的 40 个采样点估计基线，再从整条波形中减去它：
$$b=\frac{1}{40}\sum_{i=0}^{39}y_i,\qquad y'_i=y_i-b.$$
tree.GetEntry(20)
baseline = sum(adc[:40]) / 40
corrected = array("d", [0.0] * sample_count)
for i in range(sample_count):
corrected[i] = adc[i] - baseline
print(f"baseline: event 20 = {baseline:.2f} ADC")
tree->GetEntry(20);
double baseline = 0;
for (int i = 0; i < 40; ++i) baseline += adc[i];
baseline /= 40;
double corrected[250];
for (int i = 0; i < sampleCount; ++i) corrected[i] = adc[i] - baseline;
std::cout << std::fixed << std::setprecision(2)
<< "baseline: event 20 = " << baseline << " ADC" << std::endl;
baseline: event 20 = 1666.78 ADC
图中放大了触发前的基线区间。修正后，两条波形都围绕零涨落；减去基线的操作应用于全部 250 个采样点。
时间对齐¶
触发时刻相对于采样时钟会有变化，因此脉冲在各条记录中的位置并不完全相同。若直接使用固定积分门，同一个门的边界会落在不同的脉冲部位，尾部积分随之改变，给 PSD 带来额外的变化。时间对齐使积分门相对于脉冲的位置一致。
通常应以脉冲前沿的 CFD（constant-fraction discriminator）定时位置作为时间起点。对于前沿形状相近、幅度不同的脉冲，CFD 可以减小幅度变化引起的定时偏移；数字处理中可用前沿达到固定幅度比例的位置确定时间参考。CFD 定时简介
这组液体闪烁体（LS）脉冲的上升沿很快、上升时间近似不变，主要形状差异在下降尾部。因此，前沿参考点到峰顶的时间间隔近似固定；峰顶又比较尖锐，找最大采样点就能方便地完成本例的积分门对齐。这是针对这组波形的简化，不替代精密时间测量中的 CFD 定时。
若峰顶较宽、噪声较大，最大采样点容易跳动；若上升时间随事例变化，即使峰顶重合，脉冲前沿也可能没有对齐。一般情况下应使用适合该探测器脉冲的前沿定时方法，而不能直接照用峰顶对齐。
取第一条波形的峰位 $k_0$ 为参考，找出当前波形的峰位 $k_{\max}$。减去常数基线不改变最大点的位置。
tree.GetEntry(0)
alignment_target = adc.index(max(adc))
peak = corrected.index(max(corrected)) # event 20 的基线修正波形
shift = alignment_target - peak
print(f"reference maximum = {alignment_target}; event 20 maximum = {peak}")
print(f"shift = {shift} samples")
tree->GetEntry(0);
int alignmentTarget = 0, peak = 0;
for (int i = 1; i < sampleCount; ++i) {
if (adc[i] > adc[alignmentTarget]) alignmentTarget = i;
if (corrected[i] > corrected[peak]) peak = i; // event 20
}
int shift = alignmentTarget - peak;
std::cout << "reference maximum = " << alignmentTarget
<< "; event 20 maximum = " << peak << "\n"
<< "shift = " << shift << " samples" << std::endl;
reference maximum = 61; event 20 maximum = 58 shift = 3 samples
按 $\Delta=k_0-k_{\max}$ 平移：$s_j=y'_{j-\Delta}$。$\Delta>0$ 时向右移，$\Delta<0$ 时向左移；空出的采样点补零。
aligned = array("d", [0.0] * sample_count)
for i in range(sample_count):
source = i - shift
if 0 <= source < sample_count:
aligned[i] = corrected[source]
double aligned[250] = {0};
for (int i = 0; i < sampleCount; ++i) {
int source = i - shift;
if (source >= 0 && source < sampleCount) aligned[i] = corrected[source];
}
第一条波形的峰位为 61。event 20 的峰位由 58 向右移到 61，event 30 已在 61，不需平移。图中两侧都已减去基线，仅比较时间位置。
积分门与脉冲的相对位置¶
对基线修正并完成时间对齐的波形 $s_i$，分别对整段脉冲和尾部求和：
$$Q_{\rm total}=\sum_{i=t_0}^{t_2-1}s_i,\qquad Q_{\rm tail}=\sum_{i=t_1}^{t_2-1}s_i,\qquad PSD=\frac{Q_{\rm tail}}{Q_{\rm total}}.$$
total gate 从上升沿之前开始，覆盖主要脉冲及其尾部；tail gate 从峰后的下降段开始，覆盖其中的慢成分。这里 $t_0,t_1,t_2$ 表示采样点位置，每个采样间隔为 2 ns。
积分门位置示意图。tail gate 是 total gate 内的尾部区间。
电荷积分与 PSD 二维图¶
对每条处理后的波形，在两个积分门内分别求和，再用 $Q_{\rm tail}/Q_{\rm total}$ 计算 PSD。每条波形在二维图中贡献一个事例：左图比较尾部积分与总积分，右图比较尾部占比与总积分。
以下示例取 total gate 为 $[50,240)$、tail gate 为 $[72,240)$，长度分别为 380 ns 和 336 ns。在相近的 $Q_{\rm total}$ 下，中子脉冲的尾部占比更大，因此形成上方条带；gamma 位于下方。
processed 10000 waveforms total gate: [50, 240), 380 ns tail gate: [72, 240), 336 ns
固定 $Q_{\rm total}$ 区间计算 FoM¶
PSD 条带的位置和宽度随脉冲总积分变化，需要在相近的脉冲大小下比较分离效果。选择 $40000\le Q_{\rm total}<60000$ 的事例，将其 PSD 填入一维直方图，再对 gamma 和 neutron 两个峰分别作 Gaussian fit。
用拟合得到的峰位 $\mu$ 和 $FWHM=2.355\sigma$ 计算
$$FoM=\frac{|\mu_n-\mu_\gamma|}{FWHM_n+FWHM_\gamma}.$$
两峰相距越远、各自越窄，FoM 越大。左图标出所选总积分区间，右图给出该区间的 PSD 分布及拟合。
gamma: mu=0.0760, FWHM=0.0234 neutron: mu=0.2862, FWHM=0.0472 FoM = 2.977; fit status = 0/0
该示例得到 $FoM\approx2.98$。这个数值只对应所选 $Q_{\rm total}$ 区间和积分门；不能直接当作整个二维分布的单一性能指标。
平均脉冲与积分门¶
单条波形既有幅度差异，也有统计涨落。为了看清两类脉冲的形状差异，在上面同一个 $Q_{\rm total}$ 区间内，用两个 PSD 峰位的中点将事例分为两组。先把每条已对齐的波形除以自己的 $Q_{\rm total}$，再在每个采样点分别求平均：
$$\bar u_a(i)=\frac{1}{N_a}\sum_{j\in a}\frac{s_{j,i}}{Q_{{\rm total},j}}, \qquad a=\gamma,\ n.$$
归一化使每条波形在 total gate 内的面积都为 1，避免大脉冲主导平均值；逐点平均减小随机涨落，留下较清楚的脉冲形状。上图比较两类平均波形，下图给出差值 $D_i=\bar u_n(i)-\bar u_\gamma(i)$。
差值为正的尾部表示中子脉冲在这里占有更大的面积比例。将差值在 tail gate 内求和，恰好得到两组事例平均 PSD 的差：
$$\sum_{i=t_1}^{t_2-1}D_i=\overline{PSD}_n-\overline{PSD}_\gamma.$$
这说明了尾部积分为什么能够甄别两类脉冲。平均波形帮助选择尾部区间，但门也不能无限延长：脉冲回到基线后的采样点主要增加噪声，最终门宽仍通过 FoM 比较。
averaged 267 gamma-like and 588 neutron-like pulses
作业要求¶
- 参照上述方法处理全部波形，完成基线修正和时间对齐，绘制 $Q_{\rm tail}$–$Q_{\rm total}$ 与 $PSD$–$Q_{\rm total}$ 关联图。
- 选择一个 $Q_{\rm total}$ 区间，拟合其中的 PSD 分布，计算 FoM。
- 比较两类平均脉冲，选择积分门；在同一 $Q_{\rm total}$ 区间内比较调整前后的 FoM，给出所选门的范围。
进度对应第 5 章。