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

不变质量谱重建

在逆运动学中，α 和 $^{10}\mathrm{Be}$ 碎片集中在前角，反冲核的角度分布更宽。这里研究只使用两个前向碎片的分析方法：先识别靶核成分，再用 Q 门区分衰变分支，最后重建不变质量。能否同时测得反冲核还取决于探测器覆盖和阈值；本节不据无几何的 toy MC 推断实际的两体、三体探测效率。

读取 4.4 生成的同一混合样本。本节不使用反冲核的测量动能和角度，仍保留样本已有的生成选择。后面的真值谱仅用于检验，不参与门选。

靶核类型判断

对于 $(CH)_n$ 复合靶中的反应事例，若只测得碎裂末态 $\alpha$ 和 $^{10}\mathrm{Be}$，而未测得反冲靶核，则可利用 $x$-$y$ 关联图对靶核类型进行区分，并据此选择相应反应道。

反冲核 $A'$ 的动量可写为

$$ p_{A'}^2=(\vec p_a-\vec p_c-\vec p_d)^2 . $$

在低能条件下，反冲核动能近似为

$$ E_{K,A'} \approx \frac{p_{A'}^2}{2m_u A}, $$

其中 $m_u$ 为原子质量单位，$A$ 为反冲核质量数。另一方面，由能量守恒有

$$ E_{K,A'}=E_{K,a}-E_{K,c}-E_{K,d}+Q . $$

定义

$$ x=\frac{p_{A'}^2}{2m_u}, \qquad y=E_{K,a}-E_{K,c}-E_{K,d}, $$

则有近似线性关系

$$ y=\frac{x}{A}-Q . $$

因此 H、C 的斜率分别接近 1 和 $1/12$。这里 $x,y$ 均用 MeV，$x$ 使用平均束流动量计算；对于本例的负 Q，截距 $-Q$ 为正。带的宽度还包含束流能散、有限测量分辨和未探测 γ 的影响。

后续分析中，取两条经验直线之间的区域作为 H 靶候选事例： $$ y=9.5+\frac{58}{60}x, \qquad y=18.5+\frac{58}{60}x. $$

斜率 $58/60$ 及两条边界是本例的经验门，不是 H 核质量对应的精确运动学系数。程序中的 $Q_{\mathrm{rec}}=(58/60)x-y$ 是沿此带定义的分支识别量，与 4.5 的相对论 Q2 不严格相同。这些门只在本例单位、束流和分辨设置下使用；门内仍可能有 C 本底，不能把候选事件全部视为纯 H。

先按代码构造两碎片动量，再用 pvBeam-pvAlpha-pvBe 得到缺失动量。continue 跳过门外事件；只有门内事件进入右侧 Q 谱。运行完整 macro 前，先在相同目录运行 4.4 的模拟。

%jsroot on

动能、质量和测量角度到三动量的转换会重复使用，先写成一个短函数。

%%cpp -d
TVector3 MeasuredMomentum(double kineticEnergy, double mass,
                          double thetaDeg, double phiDeg)
{
    double p = std::sqrt(kineticEnergy*kineticEnergy + 2.0*kineticEnergy*mass);
    TVector3 momentum;
    momentum.SetMagThetaPhi(p, thetaDeg*TMath::DegToRad(), phiDeg*TMath::DegToRad());
    return momentum; // MeV/c；动能与质量均用 MeV
}

%%cpp -d
void TargetID_CHn()
{
    TFile *f = TFile::Open("C14_CHn_He4Be10x.root");
    if (!f || f->IsZombie()) throw std::runtime_error("Run Section 4.4 first");
    TTree *tree = f->Get<TTree>("tree");
    if (!tree) throw std::runtime_error("Missing tree");

    gStyle->SetPadLeftMargin(0.15);
    gStyle->SetPadRightMargin(0.16);

    const double u  = 1000.0;          // mass: GeV -> MeV
    const double mu = 931.494;
    const double k  = 58.0 / 60.0;

    // masses and beam mean from TParameter
    double mHe4       = ((TParameter<Double_t>*)f->Get("massHe4"))->GetVal() * u;
    double mBe10      = ((TParameter<Double_t>*)f->Get("massBe10"))->GetVal() * u;
    double mC14       = ((TParameter<Double_t>*)f->Get("massC14"))->GetVal() * u;
    double ekBeamMean = ((TParameter<Double_t>*)f->Get("ekBeamMean"))->GetVal(); // MeV

    // branches
    Double_t ek[4], theta[4], phi[4];
    Double_t weight_total;

    tree->SetBranchAddress("ek",           ek);
    tree->SetBranchAddress("theta",        theta);
    tree->SetBranchAddress("phi",          phi);
    tree->SetBranchAddress("weight_total", &weight_total);

    TH2F *hEP = new TH2F("hEP", "Target identification;x = p_{miss}^{2}/(2m_{u}) [MeV];y = T_{beam}-T_{#alpha}-T_{Be} [MeV]", 400, 0, 80, 400, 0, 80);
    TH1F *hQ  = new TH1F("hQ",  "H-candidate gate;Q_{rec} [MeV];Sum of weights", 400, -20, -5);

    // mean beam momentum
    double pBeamMean = std::sqrt(ekBeamMean*ekBeamMean + 2.0*ekBeamMean*mC14);
    TVector3 pvBeam;
    pvBeam.SetMagThetaPhi(pBeamMean, 0.0, 0.0);

    Long64_t nentries = tree->GetEntriesFast();
    for (Long64_t i = 0; i < nentries; ++i) {
        tree->GetEntry(i);

        // alpha
        TVector3 pvAlpha = MeasuredMomentum(ek[0], mHe4, theta[0], phi[0]);

        // final 10Be
        TVector3 pvBe = MeasuredMomentum(ek[1], mBe10, theta[1], phi[1]);

        // missing recoil momentum
        double pRec = (pvBeam - pvAlpha - pvBe).Mag();

        double x = pRec * pRec / (2.0 * mu);
        double y = ekBeamMean - ek[0] - ek[1];

        hEP->Fill(x, y, weight_total);

        // H-target gate
        if (y <  9.5 + k*x) continue;
        if (y > 18.5 + k*x) continue;

        double Q = k*x - y;
        hQ->Fill(Q, weight_total);
    }

    TCanvas *c1 = new TCanvas("c1", "Target identification", 800, 400);
    c1->Divide(2, 1);

    c1->cd(1);
    hEP->Draw("colz");
    TLine *l1 = new TLine(0,  9.5, 60, 67.5);
    TLine *l2 = new TLine(0, 18.5, 60, 76.5);
    l1->Draw();
    l2->Draw();
    gPad->SetLogz();

    c1->cd(2);
    hQ->Draw("hist");

    c1->Draw();
}

TargetID_CHn();

---

不变质量谱

在 H 候选事件中构造两碎片不变质量： $$M_{\mathrm{inv}}^2=(E_\alpha+E_{\mathrm{Be}})^2-|\vec p_\alpha+\vec p_{\mathrm{Be}}|^2.$$ 若包含母核的全部衰变产物，则 $$E_x(^{14}\mathrm C)=M_{\mathrm{inv}}-m(^{14}\mathrm C_{\mathrm{gs}}).$$ 探测器测到的 $^{10}\mathrm{Be}$ 已是 γ 退激后的基态核，但它可能来自 3.368 MeV 激发态分支。若 γ 未测得，两碎片四动量之和并不包含完整母核能量。此时基态质量是对测到的粒子的正确描述；缺失的是退激前母核的一部分四动量，而不是测量粒子的质量写错。因此先用 Q 门识别分支，再讨论如何近似恢复母核激发能。

在本节中，为保持与前述靶核甄别的处理一致，仍采用 H 带对应的线性关系定义 $$ Q_{\mathrm{rec}}=\frac{58}{60}x-y. $$ 然后按经验门对 $^{10}\mathrm{Be}$ 分支进行区分：

• 当 $Q_{\mathrm{rec}}<-14.5\ \mathrm{MeV}$ 时，视为 $^{10}\mathrm{Be}^{*}(3.368)$ 分支；
• 当 $Q_{\mathrm{rec}}>-13.5\ \mathrm{MeV}$ 时，视为 $^{10}\mathrm{Be}$ 基态分支；
• 中间过渡区事例暂不参与不变质量重建，以减少分支混叠。

保留本例的近似：对于激发分支，把测得的 $^{10}\mathrm{Be}$ 动能和方向近似作为退激前的量，同时采用 $m_{\mathrm{Be}}+3.368$ MeV 构造四动量。这忽略了 γ 发射造成的动能变化和反冲，不能精确恢复母核四动量。随后与相同选中事件的 exC14 比较，检查峰位恢复和剩余展宽。

因此，本节的分析流程可以概括为：

1. 从统一模拟文件中读取前向可测粒子信息；
2. 利用 $x$-$y$ 关联图进行靶核甄别，选出 H 靶候选事例；
3. 在 H 门内利用重建 $Q$ 值区分 $^{10}\mathrm{Be}$ 基态与激发态分支；
4. 对不同分支分别采用相应质量，构建 $^{14}\mathrm{C}$ 的两体不变质量谱。

%%cpp -d
void MassSpectrum_CHn()
{
    TFile *f = TFile::Open("C14_CHn_He4Be10x.root");
    if (!f || f->IsZombie()) throw std::runtime_error("Run Section 4.4 first");
    TTree *tree = f->Get<TTree>("tree");
    if (!tree) throw std::runtime_error("Missing tree");

    gStyle->SetPadLeftMargin(0.15);
    gStyle->SetPadRightMargin(0.16);

    const double u  = 1000.0;          // mass: GeV -> MeV
    const double mu = 931.494;
    const double k  = 58.0 / 60.0;

    // masses and beam mean from TParameter
    double mHe4       = ((TParameter<Double_t>*)f->Get("massHe4"))->GetVal() * u;
    double mBe10      = ((TParameter<Double_t>*)f->Get("massBe10"))->GetVal() * u;
    double mC14       = ((TParameter<Double_t>*)f->Get("massC14"))->GetVal() * u;
    double ekBeamMean = ((TParameter<Double_t>*)f->Get("ekBeamMean"))->GetVal(); // MeV

    // branches
    Double_t ek[4], theta[4], phi[4];
    Double_t exC14, weight_total;

    tree->SetBranchAddress("ek",           ek);
    tree->SetBranchAddress("theta",        theta);
    tree->SetBranchAddress("phi",          phi);
    tree->SetBranchAddress("exC14",        &exC14);
    tree->SetBranchAddress("weight_total", &weight_total);

    // 真值只用于评价下述选择，不参与门选。
    Int_t processID;
    Double_t exBe10;
    tree->SetBranchAddress("processID", &processID);
    tree->SetBranchAddress("exBe10", &exBe10);
    double selectedWeight=0, hWeight=0, correctBranchWeight=0;

    TH1F *hExcal  = new TH1F("hExcal",  "Reconstructed;E_{x}(^{14}C) [MeV];Sum of weights",   400, 12, 22);
    TH2F *hExcalQ = new TH2F("hExcalQ", "Reconstructed;E_{x}(^{14}C) [MeV];Q_{rec} [MeV]", 400, 12, 22, 400, -20, -5);
    TH1F *hEx     = new TH1F("hEx",     "Truth, same selected events;E_{x}(^{14}C) [MeV];Sum of weights", 400, 12, 22);
    TH2F *hExQ    = new TH2F("hExQ",    "Truth, same selected events;E_{x,true}(^{14}C) [MeV];Q_{rec} [MeV]",    400, 12, 22, 400, -20, -5);

    // mean beam momentum
    double pBeamMean = std::sqrt(ekBeamMean*ekBeamMean + 2.0*ekBeamMean*mC14);
    TVector3 pvBeam;
    pvBeam.SetMagThetaPhi(pBeamMean, 0.0, 0.0);

    Long64_t nentries = tree->GetEntriesFast();
    for (Long64_t i = 0; i < nentries; ++i) {
        tree->GetEntry(i);

        // step 1: use final alpha + final 10Be to do H gate
        TVector3 pvAlpha0 = MeasuredMomentum(ek[0], mHe4, theta[0], phi[0]);

        TVector3 pvBe0 = MeasuredMomentum(ek[1], mBe10, theta[1], phi[1]);

        double pRec = (pvBeam - pvAlpha0 - pvBe0).Mag();

        double x = pRec * pRec / (2.0 * mu);
        double y = ekBeamMean - ek[0] - ek[1];

        // H-target gate
        if (y <  9.5 + k*x) continue;
        if (y > 18.5 + k*x) continue;

        // reconstructed Q from x-y band
        double Q = k*x - y;

        // step 2: choose 10Be mass according to Q branch
        double mBeUse = -1.0;
        if (Q < -14.5) mBeUse = mBe10 + 3.368; // 10Be*(3.368)
        if (Q > -13.5) mBeUse = mBe10;         // 10Be(gs)

        // skip the transition region
        if (mBeUse < 0.0) continue;

        // step 3: build two-body invariant mass
        TVector3 pvAlpha = MeasuredMomentum(ek[0], mHe4, theta[0], phi[0]);

        TVector3 pvBe = MeasuredMomentum(ek[1], mBe10, theta[1], phi[1]);

        TLorentzVector lvAlpha(pvAlpha, ek[0] + mHe4);
        TLorentzVector lvBe   (pvBe,    ek[1] + mBeUse);

        double Minv = (lvAlpha + lvBe).M();
        double Ex   = Minv - mC14;

        hExcal->Fill(Ex, weight_total);
        hExcalQ->Fill(Ex, Q, weight_total);

        selectedWeight += weight_total;
        if (processID == 1) hWeight += weight_total;
        if ((Q < -14.5) == (exBe10 > 1.0)) correctBranchWeight += weight_total;

        // truth for comparison
        hEx->Fill(exC14, weight_total);
        hExQ->Fill(exC14, Q, weight_total);
    }
    if (selectedWeight > 0) {
        std::cout << "H fraction after both gates = " << hWeight/selectedWeight << std::endl;
        std::cout << "Correct branch fraction = " << correctBranchWeight/selectedWeight << std::endl;
    }

    TCanvas *c2 = new TCanvas("c2", "Invariant mass", 800, 800);
    c2->Divide(2, 2);

    c2->cd(1);
    hExcal->Draw("hist");

    c2->cd(2);
    hExcalQ->Draw("colz");
    gPad->SetLogz();

    c2->cd(3);
    hEx->Draw("hist");

    c2->cd(4);
    hExQ->Draw("colz");
    gPad->SetLogz();

    c2->Draw();
}

MassSpectrum_CHn();

H fraction after both gates = 0.999526
Correct branch fraction = 0.998855

比较上、下两行时，注意重建谱和真值谱来自同一批通过 H 门及 Q 分支门的事件，不是与未选的输入全谱比较。峰位接近说明此近似在该样本中有效；剩余偏差和展宽来自测量分辨、分支混叠及未测 γ 反冲的近似处理。

同时检查二维图：经验 Q 门是否切掉某个激发能区间，是否改变分支的相对强度。这里不能仅凭峰位对齐就声称恢复了真实产额。程序额外输出门内 H 比例和分支识别正确比例，用模拟真值评价选择，但真值不用于做选择。

作业：利用 $Q_3$ 条件重建 H 靶三体测量下的不变质量谱

在 4.5 节中，已经利用三体动量守恒方法重建了 H 靶反应的 $Q_3$ 谱。下面进一步以此为基础，完成 H 靶三体全末态测量条件下的不变质量谱重建。

选出模拟 H 靶样本作为已知反应道的练习，使用三个带电粒子构造 Q3。根据 Q3 选择基态与 3.368 MeV 分支，按本节相同近似重建母核激发能；将结果与两体方法及同一选中样本的真值比较。关注分支重叠、峰位和宽度，不要求两种方法的峰宽必然按固定次序排列。