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

3.1 DSSD 能量刻度

DSSD 探测器

• 探测器： S4

• 128(Pie) x 128(Ring)

• Pie: $1.4^{\circ}$/条, HV=-110V

• Ring: 0.43 mm 条宽，~0.03 mm 条间距, HV=0V

• 厚度1014 $\mu$m

放射源

• $\alpha$ 放射源： $^{232}U$

• $\alpha$粒子从Pie面入射，考虑放射源表面的金膜和Si前表面的死层后，最终沉积在Si里的能量分别为 $5.5658$ MeV, $6.1748$ MeV, $6.6708$ MeV, $8.6931$ MeV.

• 使用下图中标红色的4个峰进行刻度.

• 数据: s4.root
• 数据由数字化模块在自触发模式下记录， 数据截取了pie和ring的前48条信息。

 Int_t pe[48];//Pie 0-47
 Int_t re[48];//Ring 0-47

%jsroot on

TFile *f = new TFile("s4.root");
if (!f || f->IsZombie()) throw std::runtime_error("无法打开输入 ROOT 文件");
TTree *tree = (TTree*)f->Get("tree");
if (!tree) throw std::runtime_error("输入文件中缺少 tree");
TCanvas *c1 = new TCanvas;

tree->Draw("pe[0]>>hx0(1000, 0, 2000)", "", "goff");//单条
tree->Draw("pe[1]>>hx1(1000, 0, 2000)", "", "hist");//单条
TH1 *hx0 = (TH1*)gROOT->FindObject("hx0");
TH1 *hx1 = (TH1*)gROOT->FindObject("hx1");
hx0->SetLineColor(kGreen);
hx1->Draw("hist");
hx0->Draw("hist same");
c1->Draw();

tree->Draw("pe>>(1000,0,2000)", "", "hist");//所有条pe[0-47]累积
c1->Draw();

tree->Draw("re[0]>>hy0(1000, 0, 2000)", "", "goff");//单条
tree->Draw("re[1]>>hy1(1000, 0, 2000)", "", "goff");//单条
TH1 *hy0 = (TH1*)gROOT->FindObject("hy0");
TH1 *hy1 = (TH1*)gROOT->FindObject("hy1");
gPad->SetLogy();
hy0->SetLineColor(kGreen);
hy1->Draw("hist");
hy0->Draw("hist same");
c1->Draw();

TSpectrum

TSpectrum 用于一维谱分析中的几个核心问题：本底估计、平滑、反卷积和寻峰。在刻度问题中，最直接相关的是本底估计和寻峰。

创建对象时需要先给出允许搜索的最大峰数，例如：

TSpectrum *s = new TSpectrum(500);

这个数字不是实际峰数，而是允许搜索的上限。

Background

用于本底估计的接口可以写成：

TH1 *TSpectrum::Background(const TH1 *h,
                           Int_t niter = 20,
                           Option_t *option = "");

niter 是迭代次数，控制 clipping 窗口覆盖的尺度，不是简单的平滑强度。窗口需要跨过待分离的峰。将估计本底叠加在原谱上检查；本底扣除后的谱用于寻找候选峰，峰位拟合仍使用原始计数谱。

先由 Background 构造本底，再从原谱中减去本底，得到更适合后续寻峰的谱。

接口与参数说明：ROOT TSpectrum reference。

TFile *fgamma = new TFile("gamma.root");
if (!fgamma || fgamma->IsZombie()) throw std::runtime_error("无法打开输入 ROOT 文件");
TH1F *h0 = (TH1F*)fgamma->Get("h0");

TSpectrum *s = new TSpectrum(500);
TH1 *hbg = s->Background(h0, 30, "same");
TH1F *hpeak  = (TH1F*)h0->Clone("hpeak");
hpeak->Add(hbg, -1);
c1->Clear();
h0->SetLineColor(kBlack);
hbg->SetLineColor(kRed);
hpeak->SetLineColor(kBlue);
h0->Draw("hist");
hpeak->Draw("hist same");
hbg->Draw("hist same");
gPad->SetLogy();
c1->Draw();

Search

用于自动找峰的接口可以写成：

Int_t TSpectrum::Search(const TH1 *hin,
                        Double_t sigma = 2,
                        Option_t *option = "",
                        Double_t threshold = 0.05);

sigma 以 bin 为单位描述预期峰宽；threshold 是相对于搜索结果中最高峰的阈值。它们控制候选峰集，不能独立判断某个结构是否为真实谱线。

Search 默认包含本底处理。上一步已显式扣除本底，所以这里用 nobackground，避免重复扣除；nodraw 关闭自动绘图，由 hist 画谱，再显式叠加候选峰标记。hist 只改变显示，不删除直方图误差。参见 TSpectrum::Search。

自动寻峰给的是候选位置，不是精确 centroid。后面仍沿用“识别谱线 → 局部拟合 centroid → 能量刻度 → residual 检查”的顺序。

Double_t *xpeaks = nullptr;
Double_t *ypeaks = nullptr;

Int_t nfound = s->Search(hpeak, 2, "nobackground nodraw", 0.01);

xpeaks = s->GetPositionX();
ypeaks = s->GetPositionY();

cout << nfound << endl;
for(int i=0; i<nfound; i++)
    cout << i << ": " << int(xpeaks[i]) 
         << ", " << ypeaks[i] << endl;
hpeak->Draw("hist");
TPolyMarker *peakMarkers=(TPolyMarker*)hpeak->GetListOfFunctions()->FindObject("TPolyMarker");
if(peakMarkers) peakMarkers->Draw();
c1->Draw();

22
0: 332, 984052
1: 135, 786091
2: 101, 664539
3: 322, 538447
4: 67, 355974
5: 71, 117777
6: 63, 71649.2
7: 96, 61243.9
8: 59, 57259.4
9: 94, 43959.7
10: 98, 39411.4
11: 87, 29171.4
12: 132, 17675.6
13: 168, 23274.4
14: 104, 18971.5
15: 78, 18438.1
16: 160, 17473.4
17: 129, 19779
18: 165, 11409.6
19: 161, 12014.4
20: 106, 12859.6
21: 121, 10732.4

保存与排序候选峰

vector<Double_t> 存候选位置，push_back 添加一个值。下面的 helper 用 multimap<Double_t,Double_t> 保存“峰高、峰位”，反向遍历得到峰高从大到小的序列；允许相同峰高，避免漏掉等高峰。实际对应哪条 α 谱线，还要结合能谱判断。

vector<Double_t> pe;
pe.push_back(993.0);
pe.push_back(1079.0);
// 若要按位置从小到大排列：
sort(pe.begin(),pe.end());

vector 用下标访问元素，size() 给出已保存的元素个数，clear() 清空上一次调用的结果。下面的函数以引用 vector<Double_t>& 接收容器，因而调用者能得到函数写入的候选峰位。

map 按 key 排序且每个 key 只能保留一个值；multimap 允许重复 key。这里把峰高作为 key、峰位作为 value，用 rbegin() 到 rend() 的反向迭代器读取，得到从高峰到低峰的顺序。随后按峰位排序是另一项操作，用于与已知谱线的能量顺序对应。

可重复调用的寻峰函数

函数只在直方图副本上扣本底，保留原始计数用于后面的峰拟合。backsub=0 表示输入已经扣过本底。

#include <TSpectrum.h>
#include <TH1.h>
#include <TROOT.h>
#include <TPolyMarker.h>
#include <TLatex.h>
#include <TString.h>
#include <map>
#include <vector>
#include <iostream>
using namespace std;
TH1 *h = nullptr, *hb = nullptr;

void peaks(TString hname, vector<Double_t> &pe, Double_t thres=0.05, int backsub=1)
{
Int_t nfound;
Double_t *xpeaks = NULL, *ypeaks = NULL;
pe.clear();
multimap<Double_t, Double_t> me;   // key: ypeaks, value: xpeaks
TH1 source = dynamic_cast<TH1>(gROOT->FindObject(hname));
h = source ? static_cast<TH1*>(source->Clone(hname+"_search")) : nullptr;
if(!h) {
cout << "Histogram " << hname << " not found." << endl;
return;
}
TSpectrum *sp = new TSpectrum(500);
if(backsub) {
hb = (TH1F*)sp->Background(h, 80);
h->Add(hb, -1);
}
// 本底已在上面处理，不再使用 Search 的默认本底扣除。
nfound = sp->Search(h, 2, "nobackground nodraw", thres);
h->Draw("hist");
TPolyMarker pm = (TPolyMarker)h->GetListOfFunctions()->FindObject("TPolyMarker");
if(pm) {
pm->SetMarkerStyle(32);
pm->SetMarkerColor(kGreen);
pm->SetMarkerSize(0.4);
pm->Draw(); // HIST 不画附加对象，显式叠加候选峰标记
}
xpeaks = sp->GetPositionX();
ypeaks = sp->GetPositionY();
for(int j=0; j<nfound; j++) {
me.emplace(ypeaks[j], xpeaks[j]);
TLatex *tex = new TLatex(xpeaks[j], ypeaks[j], Form("%.0f", xpeaks[j]));
tex-&gt;SetTextFont(42);
tex-&gt;SetTextSize(0.025);
tex-&gt;SetTextAlign(12);
tex-&gt;SetTextAngle(90);
tex-&gt;SetTextColor(kRed);
tex-&gt;Draw();
}
for(auto ie = me.rbegin(); ie != me.rend(); ++ie) {
cout << ie->second << " " << ie->first << endl;
pe.push_back(ie->second);
}
delete sp;
}

cout<<hpeak->GetName()<<endl;

hpeak

gROOT->ProcessLine(".L peaks.C");
vector<Double_t> ge;
peaks("hpeak", ge, 0.05, 0); // 已扣除本底，不重复扣除
gPad->SetLogy(0);
c1->Draw();

332.7 984052
135.9 786091
101.7 664539
322.9 538447
67.1 355974
71.9 117777
63.3 71649.2
96.7 61243.9
59.5 57259.4

DSSD 能量刻度

对某一条 Pie 谱执行：

TFile *fs4 = new TFile("s4.root");
if (!fs4 || fs4->IsZombie()) throw std::runtime_error("无法打开输入 ROOT 文件");
TTree *tree = (TTree*)fs4->Get("tree");
if (!tree) throw std::runtime_error("输入文件中缺少 tree");
tree->Draw("pe[1]>>hx1(1000, 0, 2000)", "", "hist");//单条
vector<Double_t> pe;
Double_t e[4] = {5.5658, 6.1748, 6.6708, 8.6931};//alpha energy

peaks("hx1", pe);
gPad->SetLogy(0);
c1->Draw();

941 1750.87
1079 1558.92
993 1549.89
1147 1493.96
957 1358.88
1429 604
1047 479.911
277 238.968
1405 133
915 119.867

寻峰给出候选位置，不会自动识别核素。本例中先按峰高取前六个候选，再按 channel 排序。下面的 pe[2+i] 取其中后四个，与四条已知 α 能量依次对应；这一选择来自当前谱形，不能原样套用到全部条。

候选峰与已知能量的匹配

推广到其他条时，可以按以下方法寻找对应关系：

1. 将候选峰按 channel 排序，枚举四个峰的不同组合，按相同顺序与四个已知能量配对。
2. 对每组拟合线性刻度 E=b+kA，比较四个刻度点的 residual。初步峰位尚无可靠误差时，比较残差平方和；得到峰心及其误差后，再作相应的误差加权拟合。
3. 将较好的候选组合带回原谱，检查谱线身份、能量覆盖及其他可识别峰的位置。若不同组合都能符合，就保留待判定状态，不能仅因某组残差最小就认定匹配正确。

确定对应关系后，再从原始计数谱的局部拟合中提取精确峰心，更新刻度并检查 residual。

if (pe.size() < 6) throw std::runtime_error("候选峰不足六个，请检查寻峰结果");
sort(pe.begin(), pe.begin()+6); 
Double_t mpe[4], cpe[4];//peak,count
for(int i=0; i<4; i++) {
    mpe[i] = pe[2+i];
    cpe[i] = hx1->GetBinContent(hx1->FindBin(mpe[i]));
    cout << mpe[i] << " " << cpe[i] << endl;
}

993 1550
1079 1559
1147 1494
1429 604

粗刻度

选定 4 个候选峰位以后，可以先做一轮峰位—能量线性拟合，得到粗刻度关系：

TGraph *gr = new TGraph(4, mpe, e);
gr->Draw("A*");
gr->Fit("pol1");
TF1 *f1 = gr->GetFunction("pol1");
c1->Draw();

****************************************
Minimizer is Linear / Migrad
Chi2                      =  4.40562e-05
NDf                       =            2
p0                        =     -1.56557   +/-   0.0168397   
p1                        =   0.00717874   +/-   1.43506e-05 

通过残差图评估拟合结果

得到线性拟合后，还要快速检查拟合是否自然。这里引入残差

$$ e_i = y_i - f(x_i|\alpha). $$

如果模型和匹配都合理，残差应当围绕 0 随机分布，不应表现出明显的系统趋势。

Double_t p0 = f1->GetParameter(0);
Double_t p1 = f1->GetParameter(1);
Double_t dpe[4];
for(int i=0; i<4; i++) 
    dpe[i] = e[i] - (p0 + p1 * mpe[i]);
TGraph *rgr = new TGraph(4,mpe,dpe);
rgr->Draw("A*");
c1->Draw();

高斯峰拟合

粗刻度之后，不能直接把 TSpectrum 返回的峰位当作最终峰心。更合理的做法是对每个峰分别进行高斯拟合，用拟合得到的均值作为最终峰位。

在每个峰附近设置初值：

本例选取彼此分开的 α 峰，窄区间内先用 Gaussian 描述峰心。若局部本底或低能尾明显，应检查残差后调整模型和区间。下面的 S 返回拟合结果，Q 减少日志，+ 保留前面各峰的曲线。

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

double centroid[4], centroidError[4];
TF1 *fg[4];
TFitResultPtr fr;
for(int i=0; i<4; i++) {
    fg[i] = new TF1(Form("fg%d", i), "gaus");
    fg[i]->SetParameters(cpe[i], mpe[i], 4);
    fg[i]->SetRange(mpe[i]-5,mpe[i]+8);
    fg[i]->SetLineColor(kRed+1);
    fg[i]->SetLineWidth(2);
    fg[i]->SetParLimits(2, 0.2, 20);// 参数 2 是 sigma；限制为正值
    fr = hx1->Fit(fg[i], "SQ+", "", mpe[i]-5, mpe[i]+8);// S：返回拟合结果；Q：简短输出；+：保留前面的峰函数
    if (int(fr) != 0 || !fr.Get() || !fr->IsValid()) throw std::runtime_error("峰拟合失败");
    centroid[i] = fg[i]->GetParameter(1);
    double sigma = fg[i]->GetParameter(2);
    double chi2Ndf = fr->Chi2()/fr->Ndf();
    centroidError[i] = fg[i]->GetParError(1);
    cout << Form("peaks = %4.1f, sigma = %.2f, chi2/ndf = %.2f", 
                 centroid[i], sigma, chi2Ndf) << endl;
}
hx1->Draw("hist");
for(int i=0;i<4;++i) fg[i]->Draw("same");
c1->Draw();

peaks = 992.9, sigma = 3.98, chi2/ndf = 12.91
peaks = 1078.3, sigma = 4.00, chi2/ndf = 11.31
peaks = 1147.3, sigma = 3.98, chi2/ndf = 13.35
peaks = 1429.0, sigma = 4.83, chi2/ndf = 1.27

TGraphErrors *gr1 = new TGraphErrors(4, centroid, e, centroidError, nullptr);
gr1->SetTitle("Energy calibration;Channel;Energy (MeV)");
gr1->Draw("AP");
TFitResultPtr calibrationFit = gr1->Fit("pol1", "S");
if (int(calibrationFit)!=0 || !calibrationFit.Get() || !calibrationFit->IsValid())
    throw std::runtime_error("能量刻度拟合失败");
TF1 *f2 = gr1->GetFunction("pol1");
c1->Draw();

****************************************
Minimizer is Minuit2 / Migrad
Chi2                      =      16.6311
NDf                       =            2
Edm                       =  1.88825e-08
NCalls                    =           42
p0                        =     -1.55384   +/-   0.00292773  
p1                        =   0.00716924   +/-   2.65536e-06 

p0 = f2->GetParameter(0);
p1 = f2->GetParameter(1);
for(int i=0; i<4; i++) 
    dpe[i] = e[i] - (p0 + p1 * centroid[i]);
TGraph *rgr1 = new TGraph(4, centroid, dpe);
rgr1->Draw("A*");
c1->Draw();

逐事件应用能量刻度

刻度函数为 $E=p_0+p_1A$。参照 1.3，先对有效的整数道值进行 dithering，再计算能量。前面的 channel 谱使用整数边界，按 $N\leftrightarrow[N,N+1)$ 的约定处理；在这些整数边界内随机化不改变原谱的 bin 计数。转换到 MeV 后，新的 bin 边界通常不再与整数道对齐，dithering 可以避免直接映射离散码造成的锯齿。

下面读取与峰拟合相同的 Pie 第 1 条；本数据以 −10 标记未记录到有效信号的条，按 pe[1]>0 选择有效读数。adc 每个事例只生成一次，再由它求能量。rawPie 保留原始整数读数，p0,p1 使用上面精确峰心拟合后的结果。

Int_t rawPie[48];
tree->SetBranchAddress("pe", rawPie);
TRandom3 dither(3101); // 固定种子，重复运行得到相同结果
TH1D *hx1e = new TH1D("hx1e", "Pie strip 1;Energy (MeV);Counts", 1000, -1, 10);
for (Long64_t i=0; i<tree->GetEntries(); ++i) {
    tree->GetEntry(i);
    if (rawPie[1]<=0) continue; // 先排除无效道值
    double adc = rawPie[1] + dither.Uniform();
    double energy = p0 + p1*adc;
    hx1e->Fill(energy);
}
hx1e->Draw("hist");
c1->Draw();
cout << "Calibrated hits = " << hx1e->GetEntries() << endl;

Calibrated hits = 68119