望远镜法粒子鉴别：已刻度实验数据¶

本页使用已刻度实验数据说明三层 Si 望远镜中的停止、穿透和粒子鉴别。

实验设置¶

  • 束流：25 MeV/u 的 $^{20}\mathrm{O}$ 束流轰击 $^9\mathrm{Be}$ 靶，产生碎裂反应；
  • 探测器：厚度依次为 $1000\,\mu\mathrm{m}$（$\Delta E_1$）、$500\,\mu\mathrm{m}$（$\Delta E_2$）和 $1000\,\mu\mathrm{m}$（$\Delta E_3$）的三层位置灵敏 Si 探测器。

在本实验设置中，所关注的粒子最终停在 D3 或前面的探测器中，不穿出 D3。三层能量均已刻度，单位为 MeV。

三层 Si 望远镜示意图¶

示意图画出束流、靶和三层 Si 探测器的先后关系。D1、D2、D3 的画面宽度只用于提示相对厚度，探测器间距不按比例。蓝色路径表示粒子停止在 D2；红色路径表示粒子穿透 D2 后停止在 D3。

三层 Si 望远镜立体示意图

读入数据¶

分析采用事件选择 d1.hit==2 && ssd1<1。它只确定本页使用的事件样本，不区分粒子种类。dE1、dE2、dE3 分别表示三层探测器的已刻度能量。

In [1]:
TFile *input = TFile::Open("telenew.root", "READ");
TTree *tree = (TTree *)input->Get("tree");

tree->SetAlias("dE1", "d1.e[0]");
tree->SetAlias("dE2", "d2.e[0]");
tree->SetAlias("dE3", "d3.e[0]");

TCut quality = "d1.hit==2 && ssd1<1";
TCut stopInD2 = quality && "dE3<1.0";
TCut reachD3 = quality && "dE3>2.0";
gStyle->SetOptStat(0);

std::cout << "all entries: " << tree->GetEntries() << std::endl;

all entries: 2736751

Warning in <TClass::Init>: no dictionary for class ppac is available
Warning in <TClass::Init>: no dictionary for class dssd is available

D1 上的入射位置¶

D1 是位置灵敏探测器。条带编号的二维分布显示所选事件在探测器表面的覆盖区域。

In [2]:
TH2D *hXY = new TH2D("hXY", "D1 hit position;D1 y strip;D1 x strip",
                      32, 0, 32, 32, 0, 32);
tree->Draw("d1.x:d1.y>>hXY", quality, "goff");

TCanvas *position = new TCanvas("position", "D1 position", 700, 620);
hXY->Draw("colz");
position->Draw();

$\Delta E_1$–$\Delta E_2$ 能量关联¶

不同粒子在前两层中的能量沉积形成多条弯曲的带，并出现分支、折点和截止位置。

In [3]:
TH2D *hD1D2 = new TH2D("hD1D2", "D1-D2 energy correlation;#Delta E_{2} (MeV);#Delta E_{1} (MeV)",
                        600, 0, 300, 840, 0, 420);
tree->Draw("dE1:dE2>>hD1D2", quality, "goff");
std::cout << "selected entries: " << (Long64_t)hD1D2->GetEntries() << std::endl;

TCanvas *d1d2 = new TCanvas("d1d2", "D1-D2 correlation", 760, 620);
gPad->SetLogz();
hD1D2->Draw("colz");
d1d2->Draw();

selected entries: 21807

D2 的停止、穿透与穿透点¶

“停止在 D2”表示粒子在 D2 中耗尽剩余动能，没有到达 D3；“穿透 D2”表示粒子离开 D2 时仍有剩余动能，并进入 D3。实验中不能直接读取粒子的剩余动能，因此这里用 D3 信号作操作性判断：

  • 蓝色，dE3 < 1 MeV：D3 没有有效信号，与粒子停止在 D2 相符；
  • 红色，dE3 > 2 MeV：D3 有信号，表示粒子到达 D3，即穿透 D2。

当前选择下没有 $1\leq dE3\leq2$ MeV 的事件。对同一粒子带，停止分支与穿透分支相接的位置称为穿透点。

沿蓝色的停止分支，$\Delta E_2$ 增加时 $\Delta E_1$ 总体减小；沿红色的穿透分支，$\Delta E_2$ 增加时 $\Delta E_1$ 总体增加。两种趋势在穿透点附近相接，反映粒子从“在 D2 中耗尽能量”过渡到“带着剩余能量进入 D3”。

In [4]:
TH2D *hStopD2 = new TH2D("hStopD2", "Stop in D2;#Delta E_{2} (MeV);#Delta E_{1} (MeV)",
                          300, 0, 300, 420, 0, 420);
TH2D *hReachD3 = new TH2D("hReachD3", "Punch through D2;#Delta E_{2} (MeV);#Delta E_{1} (MeV)",
                           300, 0, 300, 420, 0, 420);
tree->Draw("dE1:dE2>>hStopD2", stopInD2, "goff");
tree->Draw("dE1:dE2>>hReachD3", reachD3, "goff");
std::cout << "D3 < 1 MeV: " << (Long64_t)hStopD2->GetEntries() << std::endl;
std::cout << "D3 > 2 MeV: " << (Long64_t)hReachD3->GetEntries() << std::endl;

hStopD2->SetLineColor(kBlue + 1);
hStopD2->SetFillColorAlpha(kBlue + 1, 0.55);
hReachD3->SetLineColor(kRed + 1);
hReachD3->SetFillColorAlpha(kRed + 1, 0.55);

TCanvas *d2Overlay = new TCanvas("d2Overlay", "D2 stopping and punch-through", 760, 620);
hStopD2->SetTitle("D2 stopping (blue) and punch-through (red);#Delta E_{2} (MeV);#Delta E_{1} (MeV)");
hStopD2->Draw("box");
hReachD3->Draw("box same");
d2Overlay->Draw();

D3 < 1 MeV: 14946
D3 > 2 MeV: 6861

分开观察 D2 停止与穿透事件¶

分开显示两类事件后，可以更清楚地看到：停止在 D2 的粒子能够使用 $\Delta E_1$–$\Delta E_2$ 关系进行同位素鉴别；穿透 D2 的粒子还带有未计入的剩余能量，仅用 D1–D2 不能完成同位素鉴别，需要结合 D3。

In [5]:
TCanvas *d2Classes = new TCanvas("d2Classes", "D2 classes", 1200, 560);
d2Classes->Divide(2, 1);
d2Classes->cd(1);
gPad->SetLogz();
hStopD2->SetTitle("Stop in D2: D3 < 1 MeV;#Delta E_{2} (MeV);#Delta E_{1} (MeV)");
hStopD2->Draw("colz");
d2Classes->cd(2);
gPad->SetLogz();
hReachD3->SetTitle("Punch through D2: D3 > 2 MeV;#Delta E_{2} (MeV);#Delta E_{1} (MeV)");
hReachD3->Draw("colz");
d2Classes->Draw();

$\Delta E_2$–$\Delta E_3$ 能量关联¶

对穿透 D2 的事件，D2 是前方的 $\Delta E$ 探测器，D3 记录后续能量沉积，因此可以用 $\Delta E_2$–$\Delta E_3$ 关系继续鉴别粒子。

In [6]:
TH2D *hD2D3 = new TH2D("hD2D3", "D2-D3 energy correlation;#Delta E_{3} (MeV);#Delta E_{2} (MeV)",
                        560, 0, 280, 600, 0, 300);
tree->Draw("dE2:dE3>>hD2D3", quality, "goff");

TCanvas *d2d3 = new TCanvas("d2d3", "D2-D3 correlation", 760, 620);
gPad->SetLogz();
hD2D3->Draw("colz");
d2d3->Draw();

实验条带宽度与有效 $\sigma(E)$¶

实验中的粒子带具有有限宽度。这里用常数项和能量相关项组成的有效分辨率描述其大致变化：

$$ \sigma(E)=\sqrt{\sigma_0^2+(rE)^2} $$

其中 $\sigma_0$ 近似表示与能量无关的基线和电子学噪声，$rE$ 表示相对响应变化。本数据的条带宽度可从 $\sigma_0\approx0.5$ MeV、$r\approx1\%$ 的量级开始描述。这里的参数包含探测器响应和实验条件的共同影响，不是单片 Si 探测器分辨率的测量值。

下一章将进一步区分 $\sqrt{E}$ 项和 $E$ 项的物理与统计来源，并学习为什么独立贡献按方差相加，以及 FWHM 与 $\sigma$ 的关系。

从薄探测器近似出发¶

在非相对论能区的简化近似中，同一种粒子在很薄的探测器内可写成 $\Delta E\propto1/E$，因此

$$ E_f^{(0)}=\sqrt{\Delta E\,E} $$

对同一种粒子的不同能量近似保持不变。把这个量用于实验数据后，条带仍有明显弯曲，说明 D1 和 D2 的厚度不能忽略，简单的薄层近似还不足以得到良好的粒子鉴别参数。

In [7]:
tree->SetAlias("P0_12", "sqrt(dE1*dE2)");
tree->SetAlias("P0_23", "sqrt(dE2*dE3)");

TH2D *hP0_12 = new TH2D("hP0_12", "Thin-detector approximation;#Delta E_{2} (MeV);#sqrt{#Delta E_{1}#Delta E_{2}} (MeV)",
                         600, 0, 300, 800, 0, 400);
TH2D *hP0_23 = new TH2D("hP0_23", "Thin-detector approximation;#Delta E_{3} (MeV);#sqrt{#Delta E_{2}#Delta E_{3}} (MeV)",
                         560, 0, 280, 640, 0, 320);
tree->Draw("P0_12:dE2>>hP0_12", stopInD2, "goff");
tree->Draw("P0_23:dE3>>hP0_23", quality, "goff");

TCanvas *thinPid = new TCanvas("thinPid", "Thin-detector PID", 1200, 560);
thinPid->Divide(2, 1);
thinPid->cd(1);
gPad->SetLogz();
hP0_12->Draw("colz");
thinPid->cd(2);
gPad->SetLogz();
hP0_23->Draw("colz");
thinPid->Draw();

进阶：有限厚度的经验 PID 参数¶

把能量从 $E$ 到 $E+\Delta E$ 的变化考虑进去，可以得到

$$ E_f^2=\int_E^{E+\Delta E}E\,dE =\Delta E\,E+\frac{1}{2}\Delta E^2. $$

积分给出的二次项系数是 $0.5$。实际数据中，$0.7$ 比 $0.5$ 更能补偿 D1 有限厚度带来的弯曲，同时还存在对后层能量 $E$ 的轻微线性依赖，因此再加入 $\beta E$ 项：

$$ P_{12}=\sqrt{\Delta E_1\Delta E_2+0.70\Delta E_1^2}+0.05\Delta E_2, \qquad P_{23}=\sqrt{\Delta E_2\Delta E_3+0.62\Delta E_2^2}-0.01\Delta E_3. $$

对 D1–D2，$\alpha=0.7$、$\beta=0.05$ 使同一种粒子在不同能量下的 $P_{12}$ 更接近恒定值；D2–D3 的系数按同样方法得到。式中的 $E$、$\Delta E$ 和 PID 参数单位均为 MeV。

这些系数来自当前数据的经验调整，不是普适常数；$P_{12}$ 和 $P_{23}$ 也不表示真实总能量或质量。

In [8]:
tree->SetAlias("P12", "sqrt(dE1*dE2+0.70*dE1*dE1)+0.05*dE2");
tree->SetAlias("P23", "sqrt(dE2*dE3+0.62*dE2*dE2)-0.01*dE3");

TH2D *hP12 = new TH2D("hP12", "Empirical P_{12};#Delta E_{2} (MeV);P_{12} (MeV)",
                       600, 0, 300, 900, 0, 450);
TH2D *hP23 = new TH2D("hP23", "Empirical P_{23};#Delta E_{3} (MeV);P_{23} (MeV)",
                       560, 0, 280, 640, 0, 320);
tree->Draw("P12:dE2>>hP12", stopInD2, "goff");
tree->Draw("P23:dE3>>hP23", reachD3, "goff");

TCanvas *empiricalPid = new TCanvas("empiricalPid", "Empirical PID", 1200, 560);
empiricalPid->Divide(2, 1);
empiricalPid->cd(1);
gPad->SetLogz();
hP12->Draw("colz");
empiricalPid->cd(2);
gPad->SetLogz();
hP23->Draw("colz");
empiricalPid->Draw();

进阶：PID 参数的一维投影¶

经验变换把同一粒子带压缩到近似相同的 PID 值，因此一维投影中的峰对应不同的粒子群。具体核素仍需通过已知参考确定。

In [9]:
TH1D *hP12Projection = new TH1D("hP12Projection", "P_{12} projection;P_{12} (MeV);Counts / 0.5 MeV",
                                900, 0, 450);
TH1D *hP23Projection = new TH1D("hP23Projection", "P_{23} projection;P_{23} (MeV);Counts / 0.5 MeV",
                                640, 0, 320);
tree->Draw("P12>>hP12Projection", stopInD2, "goff");
tree->Draw("P23>>hP23Projection", reachD3, "goff");

TCanvas *pidProjection = new TCanvas("pidProjection", "PID projections", 1200, 560);
pidProjection->Divide(2, 1);
pidProjection->cd(1);
gPad->SetLogy();
hP12Projection->Draw("hist");
pidProjection->cd(2);
gPad->SetLogy();
hP23Projection->Draw("hist");
pidProjection->Draw();