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

3.5 补充：多路径刻度与误差传播

相对刻度的多路径表述

与 3.5 相同，选定一条正面条 $X_r$ 为参考，用其原始幅度定义探测器内部的相对尺度：

$$ E^{(\mathrm{rel})} \equiv A_{x,r}, \qquad k_{x,r}=1, \qquad b_{x,r}=0, \qquad \sigma(k_{x,r})=0, \qquad \sigma(b_{x,r})=0. $$

目标仍是确定每条的共同尺度参数：

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

使正反面所有条的信号都能在同一相对尺度上比较。

这里改变的是求解途径，而不是刻度的物理目标。不再只沿固定的三步路线传播，而是把条与交叉 pixel 的关联写成一个网络：每条路径产生一组候选参数，再传播误差、合并多条路径的信息。

---

以交叉 pixel 连接两面的刻度

3.5 从统计量较好的条组逐步传递参考尺度。这里把每一对交叉条的局部拟合视为连接两面的基本关系。

对 $(i,j)$ 处选出的 single-pixel 候选，局部直线拟合给出：

$$ A_{x,i}=s_{ij}A_{y,j}+t_{ij}. $$

这条关系把 $X_i$ 与 $Y_j$ 联系起来。若一侧条在共同尺度中的增益和截距已知，就可借此求另一侧条的一组候选参数。

为区分计算中的两种信息，分别记录：

• 已刻度参考条当前采用的参数。
• 目标条通过某个交叉 pixel 得到的候选参数。同一目标条可以有多组候选。

正向传播：$X\rightarrow Y$

将 $X_i$ 的已知刻度代入局部关联式，得到 $Y_j$ 的候选参数： $$ k^{[i]}_{y,j}=k_{x,i}s_{ij}, \qquad b^{[i]}_{y,j}=k_{x,i}t_{ij}+b_{x,i}. $$

反向传播：$Y\rightarrow X$

将局部关联反解，并代入 $Y_j$ 的已知刻度，得到 $X_i$ 的候选参数： $$ k^{[j]}_{x,i}=\frac{k_{y,j}}{s_{ij}}, \qquad b^{[j]}_{x,i}=b_{y,j}-\frac{k_{y,j}t_{ij}}{s_{ij}}. $$

这样，每个有足够数据的交叉 pixel 都可传递刻度；多个交叉区域构成覆盖探测器的传播路径。

沿一条路径逐步传播

与三步法一样，先固定可靠的参考条，再依次覆盖尚未刻度的条。

例如，以 $X_r$ 为参考，通过 pixel $(r,j_1)$ 的正向关系刻度 $Y_{j_1}$；再通过 pixel $(i_1,j_1)$ 的反向关系刻度 $X_{i_1}$。新得到的 $X_{i_1}$ 又可以为其他 Y 条提供参考。每一步使用的是源条已有的全局参数与当前 pixel 的局部拟合参数。

因此，已刻度条提供参考，局部拟合提供转换关系，新刻度条继续传递尺度。它与主讲义完成的是同一项工作，但把每一步的来源写得更明确，便于比较不同路径。

候选参数的误差传播

下面四个逐项平方和是忽略相关项的近似，用于读懂参数如何传递。一般形式是 $C_v=J C_uJ^T$；同一次直线拟合的 slope 与 intercept 通常相关，完整计算需要它们的 covariance。参考条的 σ=0 仅表示相对单位被固定，不表示探测器没有测量误差。

局部拟合本身有 $\sigma(s_{ij})$ 和 $\sigma(t_{ij})$，参考条的参数也有不确定度。候选参数的方差由这两部分共同决定：

从 $X_i$ 传播到 $Y_j$： $$\sigma^2\!\left(k_{y,j}^{[i]}\right) \approx s_{ij}^2\,\sigma^2(k_{x,i}) + k_{x,i}^2\,\sigma^2(s_{ij})$$ $$\sigma^2\!\left(b_{y,j}^{[i]}\right) \approx t_{ij}^2\,\sigma^2(k_{x,i}) + k_{x,i}^2\,\sigma^2(t_{ij}) + \sigma^2(b_{x,i})$$

从 $Y_j$ 传播到 $X_i$： $$\sigma^2\!\left(k_{x,i}^{[j]}\right) \approx \frac{\sigma^2(k_{y,j})}{s_{ij}^2} + \frac{k_{y,j}^2}{s_{ij}^4}\sigma^2(s_{ij})$$ $$\sigma^2\!\left(b_{x,i}^{[j]}\right) \approx \sigma^2(b_{y,j}) + \left(\frac{t_{ij}}{s_{ij}}\right)^2\!\sigma^2(k_{y,j}) + \left(\frac{k_{y,j}}{s_{ij}}\right)^2\!\sigma^2(t_{ij}) + \left(\frac{k_{y,j}t_{ij}}{s_{ij}^2}\right)^2\!\sigma^2(s_{ij})$$

右端的 $k_{x,i}$、$\sigma^2(k_{x,i})$ 等均指传播前已经得到的参考条参数及其方差，不是本次局部拟合重新测得的量。

---

用不同路径合并候选参数

固定三步法只采用选定的传播路线，而完整的交叉网络通常允许一条目标条沿多条路线获得刻度。

一个目标条可同时与几条已刻度条相交。各交叉区域的事例数、幅度范围和参考误差不同，所以各路径得到的参数及精度也不同。合并之前应先检查候选是否相容。

若各路径候选确实独立，可以采用 inverse-variance 标量组合。路径复用了事例或共享有误差的参考时，应纳入路径间的 covariance；“来自不同路径”本身不等于独立。

以 Y 条为例，用 $\mathcal C_j$ 表示能够提供候选的已刻度 X 条集合。在独立路径近似下，分别定义增益和截距的 inverse-variance 权重： $$ w^{(k)}_{ij}=\frac{1}{\sigma^2\!\left(k^{[i]}_{y,j}\right)}, \qquad w^{(b)}_{ij}=\frac{1}{\sigma^2\!\left(b^{[i]}_{y,j}\right)}. $$

分别加权得到的参数及其方差为： $$ k_{y,j} = \frac{\sum_{i\in\mathcal{C}_j} w^{(k)}_{ij}k^{[i]}_{y,j}}{\sum_{i\in\mathcal{C}_j} w^{(k)}_{ij}}, \qquad \sigma^2(k_{y,j}) = \frac{1}{\sum_{i\in\mathcal{C}_j} w^{(k)}_{ij}}, $$ $$ b_{y,j} = \frac{\sum_{i\in\mathcal{C}_j} w^{(b)}_{ij}b^{[i]}_{y,j}}{\sum_{i\in\mathcal{C}_j} w^{(b)}_{ij}}, \qquad \sigma^2(b_{y,j}) = \frac{1}{\sum_{i\in\mathcal{C}_j} w^{(b)}_{ij}}. $$

对 X 条也可作相同组合。计算按“已知参数 → 各路径候选 → 加权组合 → 更新参数”的顺序向未刻度区域扩展。但同一批数据不能在循环中被反复当作新的独立测量，否则误差会被人为缩小。

Covariance 与独立性

上面的逐项平方和仅适用于忽略相关项的近似。令输入参数向量为 u，传播后的参数为 v(u)，一般公式为 \(C_v=J C_u J^T\)，其中 \(J_{ij}=\partial v_i/\partial u_j\)。同一次直线拟合的 slope 与 intercept 通常相关，不能遗漏它们的 covariance。

下面的小例子固定 X15 为相对单位的定义，三条路径使用互不重叠的 pixel 事例，并近似忽略探测器共同系统误差。代码保留每条路径内部 k、b 的 covariance。分别对 k、b 作 inverse-variance 平均只是便于演示的组合，不是一般情况下的最优联合估计。若循环使用同一 pixel 或共享一个有误差的参考，不能反复当作独立新信息迭代，否则会低估误差。

%jsroot on

TCanvas *c1 = new TCanvas("c1","c1");
TFile *fin = new TFile("./data/d1xy.root");
if (!fin || fin->IsZombie()) throw std::runtime_error("无法打开输入 ROOT 文件");
TTree *tree =(TTree *)fin->Get("tree");
if (!tree) throw std::runtime_error("输入文件中缺少 tree");
tree->Draw("ye:xe>>(1000,0,8000,1000,0,8000)","","colz");
c1->SetLogz();
c1->Draw();

多路径传播的具体实例

仍使用 data/d1xy.root 中的候选事例。下面不是替换 3.5 的三步法，而是用一个较小的网络演示局部拟合、候选传播和组合：

$$X_{15}\rightarrow\{Y_7,Y_{12},Y_{20}\}\rightarrow X_5.$$

固定 X[15] 为参考，先分别刻度三条 Y，再沿三条 Y 各自回到 X[5]，得到三组 $(k_{x,5},b_{x,5})$。比较各组参数与传播误差后进行组合，并用刻度前后的关联图检查结果。

代码的组织

下面依次定义参数结构、正反向传播、候选组合、局部拟合及结果显示。各部分都对应上面的公式，逐段执行即可看到中间参数。

局部拟合先用 ROB 找初始直线，再依据 residual 的稳健尺度选择主体条带。最终 least squares 的误差是在所选样本和模型下的条件误差，不包含改变 cut 或误认组合的系统影响。三条路径使用互不重叠的 pixel 事例，并近似忽略共同系统误差；代码保留每条路径内 k、b 的 covariance。

c1->Close(); // 初始图已保存在上一个输出单元
#include <iostream>
#include <vector>
#include <cmath>
#include "TFile.h"
#include "TTree.h"
#include "TGraph.h"
#include "TF1.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TStyle.h"

using namespace std;

%%cpp -d
// Block 1: Data Structures

// Current Global Parameters for a strip
struct StripParam {
    double k = 1.0, b = 0.0;
    double vk = 1e9, vb = 1e9;  // Variance (large = uncalibrated)
    bool isCalibrated = false;
    double ckb = 0; // Cov(k,b)
};

// Local Pixel Bridge from fit
struct PixelFit {
    double s = 0, t = 0;
    double vs = 0, vt = 0;
    bool isValid = false;
    double cst = 0; // Cov(s,t)
};

// Global arrays
StripParam calX[32], calY[32];

%%cpp -d
// Block 2: Propagation Engine

// Forward: X -> Y (Eq.4 + Error Propagation)
StripParam ForwardPropagate(const StripParam& X, const PixelFit& p) {
    StripParam Y;
    Y.k  = X.k * p.s;
    Y.b  = X.k * p.t + X.b;
    Y.vk = p.s*p.s * X.vk + X.k*X.k * p.vs;
    Y.vb = p.t*p.t * X.vk + X.k*X.k * p.vt + X.vb + 2*p.t*X.ckb;
    Y.ckb = p.s*p.t*X.vk + p.s*X.ckb + X.k*X.k*p.cst;
    Y.isCalibrated = true;
    return Y;
}

%%cpp -d
// Reverse: Y -> X (Eq.5 + Error Propagation)
StripParam ReversePropagate(const StripParam& Y, const PixelFit& p) {
    StripParam X;
    if (p.s == 0) return X;
    
    double s2 = p.s * p.s;
    double s4 = s2 * s2;
    
    X.k  = Y.k / p.s;
    X.b  = Y.b - Y.k * p.t / p.s;
    X.vk = Y.vk / s2 + Y.k*Y.k / s4 * p.vs;
    X.vb = Y.vb + pow(p.t/p.s, 2) * Y.vk 
               + pow(Y.k/p.s, 2) * p.vt 
               + pow(Y.k*p.t/s2, 2) * p.vs
               - 2*p.t/p.s*Y.ckb - 2*Y.k*Y.k*p.t/(s2*p.s)*p.cst;
    X.ckb = Y.ckb/p.s - p.t/s2*Y.vk
          + Y.k*Y.k/(s2*p.s)*p.cst - Y.k*Y.k*p.t/s4*p.vs;
    X.isCalibrated = true;
    return X;
}

%%cpp -d
// Block 3: Fusion Engine (Inverse-Variance Weighting)

StripParam FuseCandidates(const vector<StripParam>& cands) {
    double swk = 0, swkk = 0;
    double swb = 0, swbb = 0;

    for (const auto& c : cands) {
        if (!c.isCalibrated || c.vk<=0 || c.vb<=0) continue;
        double wk = 1.0 / c.vk;
        double wb = 1.0 / c.vb;
        swkk += wk * c.k;  swk += wk;
        swbb += wb * c.b;  swb += wb;
    }

    StripParam result;
    if (swk<=0 || swb<=0) return result;
    result.k  = swkk / swk;
    result.b  = swbb / swb;
    result.vk = 1.0 / swk;
    result.vb = 1.0 / swb;
    for (const auto& c : cands) {
        if (c.isCalibrated && c.vk>0 && c.vb>0)
            result.ckb += c.ckb / (c.vk*c.vb*swk*swb);
    }
    result.isCalibrated = true;
    return result;
}

%%cpp -d
PixelFit GetPixelRelation(TTree* tree, int xId, int yId, double nSigma = 3.0) {
    PixelFit fit;
    TString cut = Form("ix==%d && iy==%d", xId, yId);

    tree->SetEstimate(tree->GetEntries()+1);
    tree->Draw("xe:ye", cut, "goff");
    int n = tree->GetSelectedRows();
    if (n < 50) return fit; // Draw 已完成选择，不再扫描整棵 Tree 计数
    vector<double> x(tree->GetV2(), tree->GetV2() + n);
    vector<double> y(tree->GetV1(), tree->GetV1() + n);

    // Stage 1: ROB fit for initial estimate
    TGraph gAll(n, &x[0], &y[0]);
    int status = gAll.Fit("pol1", "ROB Q 0");
    if (status!=0 || !gAll.GetFunction("pol1")) return fit;
    double s0 = gAll.GetFunction("pol1")->GetParameter(1);
    double t0 = gAll.GetFunction("pol1")->GetParameter(0);

    // Stage 2: MAD-based outlier rejection
    vector<double> absRes;
    for (int i = 0; i < n; i++) 
        absRes.push_back(fabs(y[i] - s0 * x[i] - t0));
    vector<double> sortedRes = absRes; // 只排序副本，保留事例对应
    sort(sortedRes.begin(), sortedRes.end());
    double sigma = 1.4826 * sortedRes[n / 2];
    if (sigma <= 0) return fit;

    // Stage 3: Refit cleaned data with standard LS
    TGraph gClean;
    for (int i = 0; i < n; i++)
        if (absRes[i] < nSigma * sigma) 
            gClean.SetPoint(gClean.GetN(), x[i], y[i]);

    if (gClean.GetN() < 3) return fit;
    TFitResultPtr r = gClean.Fit("pol1", "Q S 0");
    if (r.Get() && r->IsValid()) {
        fit = {r->Parameter(1), r->Parameter(0), 
               pow(r->ParError(1), 2), pow(r->ParError(0), 2), true};
        fit.cst = r->CovMatrix(1,0);
        printf("  Pixel(%2d,%2d): s=%.4f±%.4f, t=%5.1f±%.1f [N=%d, cut=%d]\n",
               xId, yId, fit.s, sqrt(fit.vs), fit.t, sqrt(fit.vt), 
               gClean.GetN(), n - gClean.GetN());
    }
    return fit;
}

%%cpp -d
// Block 5: Visualization

void DrawValidation(TTree* tree, 
                    const vector<pair<int,int>>& pixels,
                    const char* title = "Validation") {
    TH2F* hRaw = new TH2F("hRaw", "Raw ADC;X (ADC);Y (ADC)", 
                          1000, 0, 8000, 1000, 0, 8000);
    TH2F* hCal = new TH2F("hCal", "Calibrated;X (Rel.Energy);Y (Rel.Energy)", 
                          1000, 0, 8000, 1000, 0, 8000);
    // 保留前面阶段的 canvas。

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

    for (Long64_t i = 0; i < tree->GetEntries(); i++) {
        tree->GetEntry(i);
        for (const auto& pix : pixels) {
            if (ix == pix.first && iy == pix.second) {
                hRaw->Fill(xe, ye);
                if (calX[ix].isCalibrated && calY[iy].isCalibrated) {
                    double eX = calX[ix].k * xe + calX[ix].b;
                    double eY = calY[iy].k * ye + calY[iy].b;
                    hCal->Fill(eX, eY);
                }
            }
        }
    }

    gStyle->SetOptStat(0);
    TCanvas* c = new TCanvas("cVal", title, 1000, 450);
    c->Divide(2, 1);
    c->cd(1); hRaw->Draw("colz");
    c->cd(2); hCal->Draw("colz");
    c->Draw();
    tree->ResetBranchAddresses(); // 局部变量即将离开作用域
}

%%cpp -d
// Main Demo: Multi-path Propagation Example

void MultiPathDemo() {
    TFile* fin = TFile::Open("./data/d1xy.root");
    if (!fin || fin->IsZombie()) { 
        cout << "Error: File missing!" << endl; 
        return; 
    }
    TTree* tree = (TTree*)fin->Get("tree");
if (!tree) throw std::runtime_error("输入文件中缺少 tree");

    // X[15] 定义本例的相对幅度尺度，不是绝对能量刻度
    for (int i=0; i<32; ++i) { calX[i]=StripParam{}; calY[i]=StripParam{}; }
    int refX = 15;
    calX[refX] = {1.0, 0.0, 0.0, 0.0, true};

    // Step 1: Forward Propagation (X15 -> Y7, Y12, Y20)
    cout << "=== Step 1: Forward X[15] -> Y[7,12,20] ===" << endl;
    vector<int> targetY = {7, 12, 20};
    
    for (int y : targetY) {
        PixelFit pf = GetPixelRelation(tree, refX, y);
        if (pf.isValid) {
            calY[y] = ForwardPropagate(calX[refX], pf);
            printf("  -> Y[%2d]: k = %.4f ± %.4f, b = %6.1f ± %.1f\n",
                   y, calY[y].k, sqrt(calY[y].vk), calY[y].b, sqrt(calY[y].vb));
        }
    }

    // Step 2: Reverse Propagation (Y7, Y12, Y20 -> X5)
    cout << "\n=== Step 2: Reverse Y[7,12,20] -> X[5] ===" << endl;
    int targetX = 5;
    vector<StripParam> candidates;

    for (int y : targetY) {
        PixelFit pf = GetPixelRelation(tree, targetX, y);
        if (pf.isValid && calY[y].isCalibrated) {
            StripParam cand = ReversePropagate(calY[y], pf);
            candidates.push_back(cand);
            printf("  -> Candidate from Y[%2d]: k = %.4f ± %.4f, b = %6.1f ± %.1f\n",
                   y, cand.k, sqrt(cand.vk), cand.b, sqrt(cand.vb));
        }
    }

    // Step 3: Inverse-Variance Weighted Fusion
    cout << "\n=== Step 3: Fusion ===" << endl;
    calX[targetX] = FuseCandidates(candidates);
    printf(">>> X[5] Final: k = %.4f ± %.4f, b = %6.1f ± %.1f <<<\n",
           calX[targetX].k, sqrt(calX[targetX].vk), 
           calX[targetX].b, sqrt(calX[targetX].vb));

    // Step 4: Visualization
    cout << "\n=== Step 4: Drawing Validation ===" << endl;
    vector<pair<int,int>> pixels = {{15,7}, {15,12}, {15,20}, {5,7}, {5,12}, {5,20}};
    DrawValidation(tree, pixels, "Multi-path Propagation Validation");
}

MultiPathDemo();

=== Step 1: Forward X[15] -> Y[7,12,20] ===
  Pixel(15, 7): s=0.9892±0.0001, t= 11.0±0.3 [N=284, cut=121]
  -> Y[ 7]: k = 0.9892 ± 0.0001, b =   11.0 ± 0.3
  Pixel(15,12): s=0.9971±0.0001, t=  4.9±0.4 [N=333, cut=162]
  -> Y[12]: k = 0.9971 ± 0.0001, b =    4.9 ± 0.4
  Pixel(15,20): s=0.9797±0.0001, t=-17.9±0.3 [N=224, cut=111]
  -> Y[20]: k = 0.9797 ± 0.0001, b =  -17.9 ± 0.3

=== Step 2: Reverse Y[7,12,20] -> X[5] ===
  Pixel( 5, 7): s=0.9736±0.0003, t= 13.4±0.3 [N=145, cut=71]
  -> Candidate from Y[ 7]: k = 1.0160 ± 0.0003, b =   -2.6 ± 0.4
  Pixel( 5,12): s=0.9814±0.0003, t=  7.8±0.3 [N=138, cut=62]
  -> Candidate from Y[12]: k = 1.0160 ± 0.0003, b =   -3.1 ± 0.5
  Pixel( 5,20): s=0.9632±0.0002, t=-13.9±0.3 [N=126, cut=62]
  -> Candidate from Y[20]: k = 1.0172 ± 0.0003, b =   -3.7 ± 0.4

=== Step 3: Fusion ===
>>> X[5] Final: k = 1.0164 ± 0.0002, b =   -3.2 ± 0.2 <<<

=== Step 4: Drawing Validation ===