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

3.5 DSSD 正反面关联（II）：逐条归一

3.4 的 pixel-pixel 方法从一对交叉条的关联求刻度，物理含义直接。但不同 pixel 的事例数往往相差很大：少数中心 pixel 可以稳定拟合，边缘 pixel 的统计量和幅度覆盖却不足。因此，不能只选一条固定参考条，就假定它与所有对面条的交叉区域都能给出可靠参数。

从单个 pixel 到 strip normalization

逐条归一先刻度一组统计量较好的条，再把它们转换到共同尺度，合并作为另一面的参考。同一目标条因此可利用多个交叉 pixel 的事例。顺序是：选择参考条；刻度对面的高统计条组；以该组刻度全部本面条；最后返回对面，补齐其余条。

本例采用 $X[16]\rightarrow Y[8\text{–}16]\rightarrow\text{全部 X}\rightarrow\text{全部 Y}$。固定参考条 $b=0,k=1$，其他系数均相对于它确定。得到的是探测器内部统一的相对幅度，而不是 MeV。只有与参考区域通过有效数据相连的条才能得到刻度；没有足够数据的条保留未刻度状态。

%jsroot on

#include <iostream>
#include <fstream>
#include <cmath>

#include "TFile.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TGraph.h"
#include "TF1.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TRandom.h"
#include "TMath.h"
#include "TString.h"

using namespace std;

// ---------- notebook-global objects ----------
TFile   *fin  = nullptr;
TTree   *tree = nullptr;
TCanvas *c1   = nullptr;

Int_t ix, iy;
Int_t xe, ye;

double parx[32][2], pary[32][2];

// Residual objects kept in memory for direct display
TH2F   *gResiduals[32] = {nullptr};
TCanvas *gResidualCanvas = nullptr;

%%cpp -d
    
void DeleteIfExists(const char* name) {
    TObject *obj = gROOT->FindObject(name);
    if (obj) delete obj;
}

void ClearResiduals() {
    for (int i = 0; i < 32; ++i) {
        if (gResiduals[i]) {
            delete gResiduals[i];
            gResiduals[i] = nullptr;
        }
    }
}

void InitPars() {
    for (int i = 0; i < 32; ++i) {
        parx[i][0] = 0.0; parx[i][1] = 0.0;   // k=0 表示尚未刻度
        pary[i][0] = 0.0; pary[i][1] = 0.0;   // b_y, k_y
    }
    parx[16][1] = 1.0; // 唯一参考尺度
}

inline double Ex(int ix, double xe) {
    return parx[ix][1] * xe + parx[ix][0];
}

inline double Ey(int iy, double ye) {
    return pary[iy][1] * ye + pary[iy][0];
}

void DrawResiduals(int id1, int id2, const char* cname = "c_residuals")
{
    static int stage = 0;
    // 各阶段保留自己的 canvas，后续拟合不更新前一单元的图。
    gResidualCanvas = new TCanvas(Form("%s_%d",cname,++stage),cname,1000,1400);

    gResidualCanvas->Divide(4, 8);

    for (int id = id1; id <= id2; ++id) {
        if (!gResiduals[id]) continue;
        gResidualCanvas->cd(id + 1);
        gResiduals[id]->Draw("colz");
    }

    gResidualCanvas->Update();
}

输入与参数表

data/d1xy.root 保存前一节方法得到的候选：ix,iy 是条号，xe,ye 是两条的 raw amplitude。选择降低了可见 sharing 和远离条带的组合，但样本仍需 residual 检查。

每行 (ix,iy,xe,ye) 描述一个候选交叉 pixel，其中 ix、iy 是条号，xe、ye 是这两条的原始幅度。程序分别保存两张参数表 parx、pary，每条的第 0 项为截距 b，第 1 项为增益 k。这样能从一个参考条开始逐步更新，而不修改输入的原始幅度。

DeleteIfExists("c1");
c1 = new TCanvas("c1", "c1", 900, 700);

fin = new TFile("./data/d1xy.root");
tree = (TTree*)fin->Get("tree");

tree->SetBranchAddress("ix", &ix);
tree->SetBranchAddress("iy", &iy);
tree->SetBranchAddress("xe", &xe);
tree->SetBranchAddress("ye", &ye);

InitPars();
gRandom->SetSeed(20260418); // 固定种子，使整本 notebook 的抽样可重复

cout << "Entries = " << tree->GetEntries() << endl;

Entries = 210820

先检查整体关联

查看所有候选的 x-y 主条带和离群区域。这一步判断是否有可用的线性关联，不直接给出每条的刻度，也不能证明每个候选都正确。

多条原始响应叠加后，不必形成一条很窄的直线：各条本来就有不同的增益和 offset。先确认主要相关区域及离群点的位置；逐条校正后，再检查这些条带能否集中到共同对角线。

DeleteIfExists("hxy_global");
tree->Draw("ye:xe>>hxy_global(1000,0,8000,1000,0,8000)", "", "colz");
c1->SetLogz();
c1->SetGrid(1,1);
c1->Draw();

选择起始参考条

统计各条的候选数。本例 X[16] 的统计量较好，并与 Y[8–16] 有较多交叉事例，因而采用这一区域开始传播。除了数量，也应查看幅度覆盖范围和 residual。

DeleteIfExists("hxy_stripmap");
tree->Draw("iy:ix>>hxy_stripmap(32,0,32,32,0,32)", "", "colz");
c1->SetLogz(0);
c1->Draw();

DeleteIfExists("hx");
DeleteIfExists("hy");

tree->Draw("ix>>hx(32,0,32)", "", "goff");
tree->Draw("iy>>hy(32,0,32)", "", "goff");

auto hx = (TH1F*)gROOT->FindObject("hx");
auto hy = (TH1F*)gROOT->FindObject("hy");

DeleteIfExists("c_strip_counts");
auto c_strip_counts = new TCanvas("c_strip_counts", "Strip counts", 900, 700);
hx->SetLineColor(kBlue+1);
hy->SetLineColor(kRed+1);
hx->SetStats(0);
hy->SetStats(0);
hx->Draw("hist");
hy->Draw("hist same");
c_strip_counts->BuildLegend();
c_strip_counts->SetGrid(1,1);
c_strip_counts->Draw();

cout << "Highest x-strip bin = " << hx->GetMaximumBin() - 1 << endl;
cout << "Highest y-strip bin = " << hy->GetMaximumBin() - 1 << endl;

Highest x-strip bin = 16
Highest y-strip bin = 12

DeleteIfExists("hy_ix16");
tree->Draw("iy>>hy_ix16(32,0,32)", "ix==16", "hist");

DeleteIfExists("c_ix16");
auto c_ix16 = new TCanvas("c_ix16", "iy distribution under ix==16", 900, 700);
auto hy_ix16 = (TH1F*)gROOT->FindObject("hy_ix16");
hy_ix16->SetStats(0);
hy_ix16->SetLineColor(kRed+1);
hy_ix16->Draw("hist");
c_ix16->SetGrid(1,1);
c_ix16->Draw();

以下固定 X[16] 为参考，第一步刻度 Y[8–16]。

拟合区间内的事例分布

拟合前先查看所选事例的幅度谱。事例主要集中在约 400–700 channel；逐事例拟合时，密集区的大量 residual 会累加到目标函数中，拟合容易更侧重这一段的符合程度，而稀疏的高幅度区约束不足。实际刻度若存在轻微非线性或不同幅度区的响应差异，得到的系数便可能偏重低幅度区，这种偏移还会随逐条归一向后传播。

本例采用以下抽样处理：对密集区中的 400 < target_raw < 600 channel 随机保留 20% 的事例，其余幅度全部保留。拟合 Y 条时取 target_raw = ye，拟合 X 条时取 target_raw = xe，即始终按待刻度条的原始幅度抽样。这样降低这段密集区对拟合的主导作用，同时保留高幅度区的统计量。

这一规则只用于构建刻度拟合样本，不改动原始数据，也不用于后面的能谱统计。400–600 channel 和 20% 是本例的设置，其他数据应根据幅度分布选择。拟合后仍要用完整候选样本检查整个幅度范围的 residual；ROB 抑制离群点，不能代替这里的幅度抽样。

DeleteIfExists("hxe");
DeleteIfExists("hye");

tree->Draw("xe>>hxe(100,0,1000)", "", "goff");
tree->Draw("ye>>hye(100,0,1000)", "", "goff");

auto hxe = (TH1F*)gROOT->FindObject("hxe");
auto hye = (TH1F*)gROOT->FindObject("hye");

DeleteIfExists("c_amp");
auto c_amp = new TCanvas("c_amp", "Amplitude spectra", 900, 700);
hxe->SetStats(0);
hye->SetStats(0);
hxe->SetLineColor(kBlue+1);
hye->SetLineColor(kRed+1);
hxe->Draw("hist");
hye->Draw("hist same");
c_amp->SetLogy();
c_amp->SetGrid(1,1);
c_amp->BuildLegend();
c_amp->Draw();

%%cpp -d
bool KeepPoint(double target_raw,
               double center = 500.0,
               double halfWidth = 100.0,
               double keepProb = 0.20)
{
    // 区间外全部保留；区间内按给定概率保留。
    if (TMath::Abs(target_raw - center) >= halfWidth) return true;
    return gRandom->Rndm() < keepProb;
}

用 Y[12] 演示一条的刻度

选择 ix==16 && iy==12，拟合横轴为 Y[12] 原始幅度、纵轴为参考 X[16] 幅度的关联。本例把 X[16] 定义为参考，所以纵轴直接用 xe。

拟合关系是 $xe=b_y+k_y\,ye$。因此 pol1 的参数 0、1 就是 Y[12] 的 b、k，后续用 $E_y^{(\rm rel)}=b_y+k_yA_y$ 转换幅度。调换横纵轴后不能继续直接套用同一组参数。

按上面的规则抽样后建立 TGraph；输出的 N 是实际进入拟合的点数。

tree->SetEstimate(tree->GetEntries()+1);
tree->Draw("xe:ye","ix==16 && iy==12","goff");
TGraph *gExample=new TGraph();
for (Long64_t i=0; i<tree->GetSelectedRows(); ++i) {
    double target_raw=tree->GetV2()[i]; // 横轴：待刻度 Y 条的幅度
    if (!KeepPoint(target_raw)) continue;
    gExample->SetPoint(gExample->GetN(),target_raw,tree->GetV1()[i]);
}
TCanvas *cStripExample=new TCanvas("cStripExample","One-strip normalization",700,430);
gExample->SetTitle("X[16] reference vs Y[12];Raw Y[12];Reference amplitude X[16]");
gExample->SetMarkerStyle(7);
gExample->Draw("AP");
gExample->Fit("pol1","Q ROB");
TF1 *fExample=gExample->GetFunction("pol1");
cout << "Y[12]: b=" << fExample->GetParameter(0) << ", k=" << fExample->GetParameter(1)
     << ", N=" << gExample->GetN() << endl;
cStripExample->Draw();

Y[12]: b=14.1673, k=0.999621, N=587

对多条重复相同操作

parx[i][0/1]、pary[j][0/1] 分别保存 b、k。Ex、Ey 按这两张表转换幅度。下面把例行操作分为收集关联点、拟合、画 residual 三部分。

fity 使用已刻度 X 面作纵轴参考，求 Y 系数；fitx 则相反。无逐点误差的 TGraph 拟合输出 RSS/NDF 带有幅度平方单位，不能要求它接近 1 来评价拟合。

两个拟合方向

fity：用已刻度 X 面作参考，横轴仍为待刻度 Y 条的原始幅度，拟合

$$E_{x,i}^{(\mathrm{rel})}=k_{y,j}A_{y,j}+b_{y,j}.$$

fitx：把方向调换，用已刻度 Y 面作参考，拟合

$$E_{y,j}^{(\mathrm{rel})}=k_{x,i}A_{x,i}+b_{x,i}.$$

每一条直线的横轴都是目标条的原始幅度，纵轴都是参考面转换后的幅度。这一约定使提取参数的代码在两个方向上保持一致。

%%cpp -d
TGraph *g[32] = {nullptr};
Int_t npt[32];
void collectCorrelations(Int_t ix1, Int_t ix2, Int_t iy1, Int_t iy2, TString fitmethod)
{
    if (fitmethod!="fity" && fitmethod!="fitx")
        throw std::runtime_error("fitmethod 应为 fity 或 fitx");

    for (Int_t i = 0; i < 32; ++i) {
        delete g[i]; // 重收集之前清理旧点；显示 residual 不删除拟合点
        g[i] = new TGraph();
        npt[i] = 0;
    }

    ClearResiduals();

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

        if (ix < ix1 || ix > ix2 || iy < iy1 || iy > iy2) continue;

        if (fitmethod == "fity") {
            if (parx[ix][1] <= 0) continue;
            double Eref = Ex(ix, xe);
            double target_raw = ye;
            if (!KeepPoint(target_raw)) continue;
            g[iy]->SetPoint(npt[iy]++, target_raw, Eref);
        } else if (fitmethod == "fitx") {
            if (pary[iy][1] <= 0) continue;
            double Eref = Ey(iy, ye);
            double target_raw = xe;
            if (!KeepPoint(target_raw)) continue;
            g[ix]->SetPoint(npt[ix]++, target_raw, Eref);
        }
    }

}

拟合各条并更新系数

只拟合至少有 20 个候选点的条，这是本例的最低数量检查，不保证幅度范围或拟合质量。参考 X[16] 始终保持 b=0、k=1，避免改变尺度定义。

%%cpp -d
void fitCorrelations(Int_t id1, Int_t id2, TString fitmethod) {
    TString sid = fitmethod=="fity" ? "iy":"ix";
    cout << Form("%4s%12s%12s%12s%10s",
                 sid.Data(), "b", "k", "RSS/NDF", "Ncounts") << endl;

    for (Int_t id = id1; id <= id2; ++id) {
        if (npt[id] < 20) {
            printf("Warning: too few points for %s-strip %d\n", sid.Data(), id);
            continue;
        }

        int status = g[id]->Fit("pol1", "Q ROB");
        if (status != 0) throw std::runtime_error("strip fit failed");
        TF1 *f = g[id]->GetFunction("pol1");
        if (!f) continue;

        double b = f->GetParameter(0);
        double k = f->GetParameter(1);
        double c2n = (f->GetNDF() > 0) ? f->GetChisquare() / f->GetNDF() : 0.0;

        if (fitmethod == "fity") {
            pary[id][0] = b;
            pary[id][1] = k;
        } else {
            if (id == 16) { b = 0; k = 1; } // 固定参考条定义
            parx[id][0] = b;
            parx[id][1] = k;
        }

        cout << Form("%4d%12.3f%12.6f%12.2f%10d",
                     id, b, k, c2n, npt[id]) << endl;
    }

}

查看 residual

每个点计算“参考幅度 − 新刻度幅度”，在各条的图中检查是否有随幅度变化的弯曲或偏移。绘图使用拟合过的同一批点，属于一致性检查；独立 run 可用于进一步验证。

%%cpp -d

void showResiduals(Int_t id1, Int_t id2, TString fitmethod) {
    ClearResiduals();
    TString sid=fitmethod=="fity" ? "iy":"ix";
    for(int id=id1;id<=id2;++id) {
        double b=fitmethod=="fity"?pary[id][0]:parx[id][0];
        double k=fitmethod=="fity"?pary[id][1]:parx[id][1];
        if(k<=0 || g[id]->GetN()<20) continue;
        gResiduals[id] = new TH2F(
            Form("h2res_%s_%02d", sid.Data(), id),
            Form("Residuals %s-strip %d;Residual;Raw amplitude", sid.Data(), id),
            100, -50, 50,
            800, 0, 8000
        );

        double *x = g[id]->GetX();
        double *y = g[id]->GetY();
        for (Int_t i = 0; i < g[id]->GetN(); ++i) {
            double r = y[i] - (b + k * x[i]);
            gResiduals[id]->Fill(r, x[i]);
        }

    }
    DrawResiduals(id1,id2,Form("c_res_%s",fitmethod.Data()));
}
void fit(Int_t ix1,Int_t ix2,Int_t iy1,Int_t iy2,TString method) {
    int first=method=="fity"?iy1:ix1, last=method=="fity"?iy2:ix2;
    collectCorrelations(ix1,ix2,iy1,iy2,method);
    fitCorrelations(first,last,method);
    showResiduals(first,last,method);
}

现在三个调用按同一顺序收集、拟合并显示结果，物理选择与前面的单条例子相同。

三步刻度

1. X[16] → Y[8–16]

用初始参考条建立中央 Y 区域的共同尺度。

此时只允许 X[16] 提供参考，避免把尚未刻度的 X 条混进纵轴。每条 Y 的拟合图包含它与 X[16] 的交叉事例，拟合后将结果写入 pary。

fit(16,16,8,16,"fity");
gResidualCanvas->Draw();

  iy           b           k     RSS/NDF   Ncounts
   8       9.212    0.992384     3093.02       469
   9      11.750    0.963387     3061.26       508
  10      12.719    0.987965     3144.65       511
  11      15.792    0.981848     4327.23       525
  12      14.240    0.999619     5088.55       590
  13      15.911    0.979801     3701.72       519
  14       8.777    0.980473     3075.19       503
  15      18.704    0.988394     3480.36       487
  16      -2.730    0.980747     2656.77       477

2. Y[8–16] → 全部 X

合并已经刻度的 Y 条，增加每条 X 的统计量。X[16] 不改变。

现在参考不再是某一条 Y 的原始幅度，而是 Y[8–16] 各自刻度后的幅度。对目标 X[i]，这些不同 Y 条的事例合并到同一个 TGraph 中，增加可用统计量和幅度覆盖。

fit(0,31,8,16,"fitx");
gResidualCanvas->Draw();

  ix           b           k     RSS/NDF   Ncounts
   0      10.333    1.001017     7239.54       486
   1       4.859    1.008147     4928.05       586
   2       0.175    1.009344     4816.26       754
   3       5.585    0.991565     4709.86       971
   4       0.500    1.001091     4161.21      1169
   5       5.897    1.018953     3754.75      1176
   6       1.688    1.013953     3836.93      1315
   7       5.713    1.011026     3799.34      1419
   8       0.788    1.008791     3575.19      1546
   9       6.709    1.005999     3991.97      1715
  10       2.775    1.022553     3899.48      1864
  11       6.389    1.016543     3754.57      1887
  12       2.091    0.995608     4539.31      2205
  13       6.657    1.020993     4335.05      2280
  14      -0.180    1.027824     3755.34      2542
  15       9.087    1.002689     4054.74      3045
  16       0.000    1.000000     3186.78      4553
  17       6.837    1.006057     3476.73      3709
  18      -0.429    1.014302     3373.57      3855
  19       3.248    0.982068     3661.03      3673
  20      -1.053    1.011768     3213.99      3411
  21       3.400    1.009598     3410.58      2650
  22      -0.878    1.001278     3148.34      2514
  23       6.202    1.006310     3247.06      2208
  24       2.756    0.990219     3193.75      1993
  25       2.380    0.978621     3370.83      1650
  26       1.207    0.985333     4115.24      1477
  27       3.596    0.972839     3890.81      1411
  28       1.391    0.965847     3946.58      1113
  29       1.287    0.999668     4542.31       932
  30       2.459    1.004284     4686.67       653
  31       3.165    0.969239     7631.50       618

3. 全部已刻度 X → 全部 Y

再用有可用系数的 X 条完成其余 Y 条。

此时一个 Y 条可以利用多条已刻度 X 的事例。每步图中的空白面板对应没有得到可用拟合的条，应结合计数输出检查；不要把空白误当成零 residual。

collectCorrelations(0,31,0,31,"fity");

fitCorrelations(0,31,"fity");

  iy           b           k     RSS/NDF   Ncounts
   0      23.107    0.955886     5094.34      2125
   1      16.737    0.975259     4529.77      2770
   2      16.124    0.971054     3948.14      3500
   3      21.220    0.975471     4191.38      4107
   4       8.192    0.978029     3670.25      4735
   5      15.850    0.986296     4166.03      5116
   6      10.222    0.982382     3662.05      5537
   7      19.977    0.992016     4224.47      5996
   8       9.272    0.992379     3831.63      6568
   9      12.121    0.963360     3799.20      6842
  10      12.542    0.988041     4133.43      6980
  11      16.011    0.981685     4026.06      6952
  12      13.800    0.999753     4414.65      7262
  13      15.934    0.979866     3680.22      7029
  14       8.826    0.980464     3809.90      7093
  15      18.779    0.988287     4213.65      6543
  16      -2.835    0.980683     3581.55      6546
  17      11.591    0.969758     3726.28      6171
  18       1.033    1.000564     4119.22      5771
  19       8.278    0.996615     3948.44      5222
  20      -9.409    0.982335     3460.19      4904
  21       2.832    0.967582     4736.64      5320
Warning: too few points for iy-strip 22
  23       6.545    0.956047     5131.98      4957
  24       4.900    0.992193     3744.83      3499
  25      -0.265    0.974654     3620.53      3510
  26      -0.221    0.993838     3905.24      3086
  27       8.756    0.981591     3801.50      2624
  28      -4.470    0.984534     4082.06      2152
  29      -3.107    0.969547     4225.63      1746
  30       0.757    0.976658     4482.74      1218
  31       1.938    0.969067     5252.99      1146

showResiduals(0,31,"fity");
gResidualCanvas->Draw();

保存与读取参数

文件前 32 行是 X，后 32 行是 Y，每行为 strip b k。无可用数据的条保持 k=0，后续跳过并报告；不能把恒等变换当成它的刻度。

%%cpp -d
    
void saveParameters(const char *fname="cal_dssd1.txt")
{
    ofstream fout(fname);
    if (!fout) throw std::runtime_error("无法写入刻度参数文件");
    for (int i = 0; i < 32; ++i)
        fout << i << " " << parx[i][0] << " " << parx[i][1] << endl;
    for (int i = 0; i < 32; ++i)
        fout << i << " " << pary[i][0] << " " << pary[i][1] << endl;
    fout.close();
}

void readParameters(const char *fname, double px[32][2], double py[32][2])
{
    ifstream fin(fname);
    int id;
    for (int i = 0; i < 32; ++i) {
        if (!(fin >> id >> px[i][0] >> px[i][1]) || id!=i)
            throw std::runtime_error("X 刻度参数缺失或条号顺序错误");
    }
    for (int i = 0; i < 32; ++i) {
        if (!(fin >> id >> py[i][0] >> py[i][1]) || id!=i)
            throw std::runtime_error("Y 刻度参数缺失或条号顺序错误");
    }
    fin.close();
}

saveParameters("cal_dssd1.txt");

检查整体归一结果

这里不再调用 KeepPoint，包括抽样时未进入拟合的事例。重新计算全部有可用刻度参数的候选的 $E_x^{(\rm rel)}$ 和 $E_y^{(\rm rel)}$，查看二维主条带及 $E_y^{(\rm rel)}-E_x^{(\rm rel)}$ 是否在零附近、是否随幅度发生漂移。整体主条带正确不替代逐条检查；少数错误系数可能被高统计条掩盖。

%%cpp -d
void calibrateAndPlot(const char* fname="cal_dssd1.txt")
{
    double px[32][2], py[32][2];
    readParameters(fname, px, py);

    DeleteIfExists("hxy_check");
    DeleteIfExists("hdiff_check");

    auto hxy = new TH2F("hxy_check", "Normalized front-back;X relative amplitude;Y relative amplitude", 400, 0, 8000, 400, 0, 8000);
    auto hdiff = new TH2F("hdiff_check", "Normalized residual;X relative amplitude;Y-X relative amplitude", 400, 0, 8000, 200, -200, 200);

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

        if (px[ix][1]<=0 || py[iy][1]<=0) { ++skipped; continue; }
        double cxe = px[ix][1] * xe + px[ix][0];
        double cye = py[iy][1] * ye + py[iy][0];

        hxy->Fill(cxe, cye);
        hdiff->Fill(cxe, cye - cxe);
    }

    DeleteIfExists("c_check");
    auto c_check = new TCanvas("c_check", "Global check", 1000, 430);
    c_check->Divide(2,1);
    c_check->cd(1); gPad->SetLeftMargin(0.14); gPad->SetRightMargin(0.16); hxy->SetStats(0); hxy->Draw("colz");
    c_check->cd(2); gPad->SetLeftMargin(0.14); gPad->SetRightMargin(0.16); hdiff->SetStats(0); hdiff->Draw("colz");
    c_check->Update();
    cout << "Input=" << nentries << ", calibrated=" << nentries-skipped
         << ", unavailable calibration=" << skipped << endl;
}

calibrateAndPlot("cal_dssd1.txt");
TCanvas *cCheck=(TCanvas*)gROOT->GetListOfCanvases()->FindObject("c_check");
cCheck->Draw();

Input=210820, calibrated=210820, unavailable calibration=0

讨论：传播路径与参数误差

三步法把复杂的条间关联拆成几轮直线拟合，容易实施，也便于逐条检查。它的限制是校准路线由所选参考条和参考区域决定：起始参数的偏移会传给后续条，不同条也可能共享同一参考误差。传播只有一两步，并不意味着这些相关误差可以忽略。

同一目标条通常有多个可用交叉 pixel。沿不同路径求得的参数可以用来检查刻度稳定性；若要合并路径，应同时考虑局部拟合质量、误差传播和路径之间共享的信息。补充：多路径刻度与误差传播保留了这一推广的公式和实例。

作业：三层 DSSD 的相对归一与 hit 列表

参照三步方法求 DSSD1–3 的参数，检查各条 residual，保存参数。把 raw amplitude 转为相对幅度，选择幅度≥100 的条，分别按两面幅度从高到低排列，条号与幅度一起移动，写入新的 ROOT 文件。

Int_t x1hit, y1hit;
Int_t x1[32], y1[32];
Double_t x1e[32], y1e[32];
Double_t x1es, y1es;
Int_t x1m, y1m;

branch 的 leaflist 使用 x1[x1hit]/I、x1e[x1hit]/D 等形式。x1es,y1es 是通过阈值的相对幅度之和；x1m,y1m 只统计相邻有效条的连通组数，不代表已确定的粒子数。DSSD2、3 使用同样结构，保留源事例编号。3.6 再根据两面能量关系判断如何合并或配对。

x1hit,y1hit 是通过阈值的条数；x1,y1 保存这些条的条号；x1e,y1e 保存对应相对幅度。排序移动的是完整 hit，而不是独立的能量数组。相邻条的连通组计数 x1m,y1m 用于描述空间拓扑，不能代替粒子数判定。