// 只要外推到靶上的误差小于 1.5 mm 的事件
TCut good_track_cut = "sigma_tx < 1.5 && sigma_theta < 0.05";
tree->Draw("tx >> h1", good_track_cut);

# 在终端中启动 ROOT 并加载数据文件
root -l f8ppac001.root
root [1] tree->MakeClass("tracking");
root [2] .q

#ifndef tracking_h
#define tracking_h

#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>

class tracking {
public :
   TTree          *fChain;   //!pointer to the analyzed TTree or TChain
   Int_t           fCurrent; //!current Tree number in a TChain

   // ================= 原始数据 Branch 变量 =================
   Float_t         PPACF8[5][5];
   Float_t         F8PPACRawData[5][5];
   Int_t           beamTrig;
   Int_t           must2Trig;
   Float_t         targetX, targetY;

   // ================= 用户自定义物理分析变量 =================
   Double_t        xx[3], xz[3], yy[3], yz[3];    // 参考层 1A, 2A, 3 的位置坐标
   Double_t        xx2b[2], yy2b[2], xz2b, yz2b;  // 待测层 2B: [0]实际测量值, [1]拟合外推值
   Double_t        anode2b;                       // 待测层 2B 阳极信号
   Double_t        dx[3], dy[3];                  // 拟合残差 (Residual)
   
   // 运动学物理量与本征误差
   Double_t        tx, ty;                        // 靶上外推坐标 (Target pos)
   Double_t        theta_x, theta_y;              // 粒子出射角 (单位: rad)
   Double_t        sigma_tx, sigma_ty;            // 靶上坐标的绝对物理误差
   Double_t        sigma_thetax, sigma_thetay;    // 出射角的绝对物理误差
   
   Double_t        c2nx, c2ny;                    // 直线拟合质量评估 chi2/ndf

   // 原始数据 Branch 指针 
   TBranch        *b_PPACF8;          //!
   TBranch        *b_F8PPACRawData;   //!
   TBranch        *b_beamTrig;        //!
   TBranch        *b_must2Trig;       //!
   TBranch        *b_targetX;         //!   
   TBranch        *b_targetY;         //!
   
   // ================= 标准系统函数 =================

   tracking(TTree *tree=0);
   virtual ~tracking();
   virtual Int_t    Cut(Long64_t entry);
   virtual Int_t    GetEntry(Long64_t entry);
   virtual Long64_t LoadTree(Long64_t entry);
   virtual void     Init(TTree *tree);
   virtual void     Loop();
   virtual Bool_t   Notify();
   virtual void     Show(Long64_t entry = -1);

   // ================= 用户自定义成员函数 =================
   virtual void     SetBranch(TTree *tree);       
   virtual void     TrackInit();                  
   virtual void     SetTrace(TH2D *h, Double_t k, Double_t b, Int_t min, Int_t max); 
};

// ... (下方的 tracking_cxx 实现部分保持原样不变) ...
#endif

#define tracking_cxx
#include "tracking.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TF1.h>
#include <TGraphErrors.h>   // 【修改】必须包含带有误差的图类
#include <TFitResult.h>
#include <TMatrixDSym.h>    // 【新增】用于接收协方差矩阵
#include <iostream>
#include <cmath>

using namespace std;

void tracking::SetBranch(TTree *tree)
{
    tree->Branch("xx", xx, "xx[3]/D");
    tree->Branch("xz", xz, "xz[3]/D");
    tree->Branch("yy", yy, "yy[3]/D");
    tree->Branch("yz", yz, "yz[3]/D");
    tree->Branch("dx", dx, "dx[3]/D");
    tree->Branch("dy", dy, "dy[3]/D");
    tree->Branch("xx2b", xx2b, "xx2b[2]/D");
    tree->Branch("yy2b", yy2b, "yy2b[2]/D");
    tree->Branch("anode2b", &anode2b, "anode2b/D");  
  
    // 【新增】注册所有运动学中心值及其物理误差
    tree->Branch("tx", &tx, "tx/D");
    tree->Branch("ty", &ty, "ty/D");
    tree->Branch("theta_x", &theta_x, "theta_x/D");
    tree->Branch("theta_y", &theta_y, "theta_y/D");
    tree->Branch("sigma_tx", &sigma_tx, "sigma_tx/D");
    tree->Branch("sigma_ty", &sigma_ty, "sigma_ty/D");
    tree->Branch("sigma_thetax", &sigma_thetax, "sigma_thetax/D");
    tree->Branch("sigma_thetay", &sigma_thetay, "sigma_thetay/D");
  
    tree->Branch("c2nx", &c2nx, "c2nx/D");
    tree->Branch("c2ny", &c2ny, "c2ny/D");
    tree->Branch("beamTrig", &beamTrig, "beamTrig/I");
    tree->Branch("must2Trig", &must2Trig, "must2Trig/I");
    tree->Branch("targetX", &targetX, "targetX");
    tree->Branch("targetY", &targetY, "targetY");  
}

void tracking::TrackInit()
{
    // 初始化所有计算变量为无效值，防止上一个事件的数据污染
    tx = -999; ty = -999;
    theta_x = -999; theta_y = -999;
    sigma_tx = -1; sigma_ty = -1;             // 误差初始化为负数代表无效
    sigma_thetax = -1; sigma_thetay = -1;
  
    xx[0] = PPACF8[0][0];  yy[0] = PPACF8[0][1];  xz[0] = PPACF8[0][2];  yz[0] = PPACF8[0][3];
    xx[1] = PPACF8[2][0];  yy[1] = PPACF8[2][1];  xz[1] = PPACF8[2][2];  yz[1] = PPACF8[2][3];
    xx[2] = PPACF8[4][0];  yy[2] = PPACF8[4][1];  xz[2] = PPACF8[4][2];  yz[2] = PPACF8[4][3];

    xx2b[0] = PPACF8[3][0]; yy2b[0] = PPACF8[3][1];
    xz2b    = PPACF8[3][2]; yz2b    = PPACF8[3][3];
    anode2b = PPACF8[3][4];
  
    xx2b[1] = -1000; yy2b[1] = -1000;
}

void tracking::SetTrace(TH2D *h, Double_t k, Double_t b, Int_t min, Int_t max){
    if(h == 0 || min >= max) return;
    for(int i = min; i < max; i++){
        h->Fill(i, (Int_t)(i * k + b));
    }
}

void tracking::Loop()
{
    if (fChain == 0) return;

    TFile *opf = new TFile("tracking.root", "RECREATE");
    TTree *tree = new TTree("tree", "PPAC Tracking Data");
    SetBranch(tree);

    TH2D *htf8xz = new TH2D("htf8xz", "X-Z Plane Trace; Z (mm); X (mm)", 2200, -2000, 200, 300, -150, 150);
    TH2D *htf8yz = new TH2D("htf8yz", "Y-Z Plane Trace; Z (mm); Y (mm)", 2200, -2000, 200, 300, -150, 150);

    // 【核心修改】使用 TGraphErrors 替代 TGraph，赋予探测器真实的本征分辨率
    TGraphErrors *grx = new TGraphErrors(3);
    TGraphErrors *gry = new TGraphErrors(3);
    TF1 *fx = new TF1("fx", "pol1", -2000, 0); 
    TF1 *fy = new TF1("fy", "pol1", -2000, 0);

    // 假设：所有PPAC每层的本征位置分辨率为 1.0 mm
    const double det_resolution = 1.0; 
    const double z_target = 0.0; // 物理靶所在Z坐标位置

    Long64_t nentries = fChain->GetEntriesFast();
    Long64_t nbytes = 0, nb = 0;
    
    for (Long64_t jentry = 0; jentry < nentries; jentry++) {
        Long64_t ientry = LoadTree(jentry);
        if (ientry < 0) break;
        nb = fChain->GetEntry(jentry);   nbytes += nb;
       
        TrackInit();  
       
        bool b1a = abs(xx[0]) < 150 && abs(yy[0]) < 150; 
        bool b2a = abs(xx[1]) < 150 && abs(yy[1]) < 150;
        bool b3  = abs(xx[2]) < 100 && abs(yy[2]) < 100;
        if(!b1a || !b2a || !b3) continue; 
      
        // ================= X-Z 平面径迹拟合与误差计算 =================
        for(int i=0; i<3; i++) {
            grx->SetPoint(i, xz[i], xx[i]); 
            grx->SetPointError(i, 0.0, det_resolution); // 关键：输入Z和X的误差
        }
        
        // 【核心修改】去除 "W" 选项。S=保存结果(以获取矩阵), Q=静默模式
        TFitResultPtr rx = grx->Fit(fx, "SQ"); 
        
        if (rx->IsValid()) {
            double p0_x = fx->GetParameter(0);
            double p1_x = fx->GetParameter(1);
            
            // 提取中心值
            xx2b[1] = fx->Eval(xz2b);   
            tx      = p0_x + p1_x * z_target;
            theta_x = atan(p1_x); // 物理出射角 (rad)
            
            // 提取协方差矩阵并计算严谨物理误差
            TMatrixDSym cov_x = rx->GetCovarianceMatrix();
            double var_p0 = cov_x(0, 0);
            double var_p1 = cov_x(1, 1);
            double cov_p0_p1 = cov_x(0, 1);
            
            // 外推位置误差传递公式
            double err2_tx = var_p0 + (z_target * z_target * var_p1) + (2.0 * z_target * cov_p0_p1);
            sigma_tx = sqrt(err2_tx);
            
            // 角度非线性误差传递公式
            sigma_thetax = sqrt(var_p1) / (1.0 + p1_x * p1_x);
            
            c2nx = rx->Chi2() / rx->Ndf();
            SetTrace(htf8xz, p1_x, p0_x, -1800, 0);
            for(int i=0; i<3; i++) dx[i] = xx[i] - fx->Eval(xz[i]); 
        }

        // ================= Y-Z 平面径迹拟合与误差计算 =================
        for(int i=0; i<3; i++) {
            gry->SetPoint(i, yz[i], yy[i]);
            gry->SetPointError(i, 0.0, det_resolution); 
        }
        
        TFitResultPtr ry = gry->Fit(fy, "SQ");
        
        if (ry->IsValid()) {
            double p0_y = fy->GetParameter(0);
            double p1_y = fy->GetParameter(1);
            
            yy2b[1] = fy->Eval(yz2b);
            ty      = p0_y + p1_y * z_target;
            theta_y = atan(p1_y);
            
            TMatrixDSym cov_y = ry->GetCovarianceMatrix();
            double var_p0 = cov_y(0, 0);
            double var_p1 = cov_y(1, 1);
            double cov_p0_p1 = cov_y(0, 1);
            
            double err2_ty = var_p0 + (z_target * z_target * var_p1) + (2.0 * z_target * cov_p0_p1);
            sigma_ty = sqrt(err2_ty);
            
            sigma_thetay = sqrt(var_p1) / (1.0 + p1_y * p1_y);
            
            c2ny = ry->Chi2() / ry->Ndf(); 
            SetTrace(htf8yz, p1_y, p0_y, -1800, 0);
            for(int i=0; i<3; i++) dy[i] = yy[i] - fy->Eval(yz[i]);
        }
     
        // 将本事件结果写入 Tree (包括新算出的误差)
        tree->Fill();
        
        if(jentry % 10000 == 0) cout << "Processing Event: " << jentry << " / " << nentries << endl;
    }
    
    // 释放内存并保存结果
    delete grx; delete gry;
    delete fx;  delete fy;

    htf8xz->Write();
    htf8yz->Write();
    tree->Write();
    opf->Close();
    
    cout << "Tracking finished. Output saved to tracking.root" << endl;
}

//%jsroot on
TFile *ipf = new TFile("tracking.root");
TTree *tree = (TTree*) ipf->Get("tree");
TCanvas *c1 = new TCanvas("c1","c1");

TH2D *hxz = (TH2D*) ipf->Get("htf8xz");
hxz->Draw("colz");
c1->Draw();

TH2D *hyz = (TH2D*) ipf->Get("htf8yz");
hyz->Draw("colz");
c1->Draw();

tree->Draw("ty:tx>>htx(120,-60,60,120,-60,60)","must2Trig","colz");
c1->Draw();

tree->Draw("tx:ty>>htx(120,-60,60,120,-60,60)","beamTrig","colz");
c1->Draw();

TF1 *total2 = new TF1("total2", "gaus(0)+gaus(3)");

TF1 *g1 = new TF1("g1", "gaus");
TF1 *g2 = new TF1("g2", "gaus");
TF1 *total = new TF1("total", "gaus(0) + gaus(3)");
TH1F *hdx;
Double_t sigma;

tree->Draw("dx[0]>>hdx(200,-5,5)");
hdx = (TH1F*)gROOT->FindObject("hdx");
hdx->Fit("g1","","",-0.5,0.5);
sigma = g1->GetParameter(2);
total->SetParameter(2,sigma);
total->SetParameter(5,7*sigma);//初始化，估计半宽为2*sigma；
hdx->Fit("total");
gPad->SetLogy(0);
c1->Draw();//residual

****************************************
Minimizer is Minuit2 / Migrad
Chi2                      =      1890.97
NDf                       =           17
Edm                       =  5.74935e-07
NCalls                    =           70
Constant                  =      16429.2   +/-   53.2958     
Mean                      =   -0.0644376   +/-   0.000586192 
Sigma                     =     0.225004   +/-   0.000603074  	 (limited)
****************************************
Minimizer is Minuit2 / Migrad
Chi2                      =      4054.38
NDf                       =          194
Edm                       =  8.03563e-06
NCalls                    =          327
p0                        =      15539.3   +/-   54.0504     
p1                        =   -0.0690443   +/-   0.000592113 
p2                        =     0.216952   +/-   0.000650672 
p3                        =      852.983   +/-   8.11293     
p4                        =   -0.0868271   +/-   0.00572421  
p5                        =      1.36888   +/-   0.00718881

tree->Draw("c2ny:dy[0]>>hh(40,-10,10,200,0,1000)","","colz");
c1->SetLogy(0);
c1->Draw();//从chi2/ndf图上可看出，部分事件的径迹拟合误差很大，这部分要在后续数据处理中去掉。

tree->Draw("c2ny>>hh(200,0,1000)","","");
gPad->SetLogy();
c1->Draw();//从chi2/ndf图上可看出，部分事件的径迹拟合误差很大，这部分要在后续数据处理中去掉。

gPad->SetLogy(0);
tree->Draw("tx:ty>>(120,-60,60)","c2nx<10 && c2ny<10 && beamTrig ","colz");
c1->Draw();//

tree->Draw("tx:ty>>(120,-60,60,120,-60,60)","(c2nx>20 || c2ny>20) && beamTrig ","colz");
c1->Draw();//

// ================= 入射粒子 =================
// 1. 预期穿过：拟合径迹外推至 2B 平面，且落在 200x200 mm 的灵敏面积内
TCut c2btrack = "abs(xx2b[1])<100 && abs(yy2b[1])<100"; 

// 2. 2B 探测器的阳极 (Anode) 确实给出了有效的时间信号
TCut c2ba     = "anode2b>-1000"; 

// ================= 探测数目 =================
// 测得有效位置：实际测量值落在正常范围内 (排除 -1000 等无效值)
TCut c2bx     = "abs(xx2b[0])<100";  // X 有信号
TCut c2by     = "abs(yy2b[0])<100";  // Y 信号

// 绘制外推位置的 2D 分布 (预期击中位置)
tree->Draw("yy2b[1]:xx2b[1]>>h_expected(200,-100,100,200,-100,100)", c2btrack, "colz");
c1->Draw();
// 统计预期穿过灵敏区域的总数
Long64_t N_track = tree->GetEntries(c2btrack);
cout << "预期穿过 PPAC2B 灵敏区的粒子总数 N_track = " << N_track << endl;

预期穿过 PPAC2B 灵敏区的粒子总数 N_track = 232296

Long64_t Na;
Long64_t Nx, Ny, Nxy;
Double_t effa;
Double_t effx1, effy1, effxy1;
Double_t effx2, effy2, effxy2;

Long64_t Na, Nx, Ny, Nxy;

// 统计各项实际响应的事件数
Na  = tree->GetEntries(c2ba && c2btrack);       // 阳极触发数
Nx  = tree->GetEntries(c2bx && c2btrack);       // 仅要求 X 
Ny  = tree->GetEntries(c2by && c2btrack);       // 仅要求 Y 
Nxy = tree->GetEntries(c2bx && c2by && c2btrack); // 要求 X 和 Y 

// ==========================================================
// 方法一：径迹外推的“绝对探测效率” 
// ==========================================================
Double_t eff_a  = Double_t(Na)  / N_track; // 阳极本征效率
Double_t eff_x1 = Double_t(Nx)  / N_track; // X 绝对效率
Double_t eff_y1 = Double_t(Ny)  / N_track; // Y 绝对效率
Double_t eff_xy1 = Double_t(Nxy) / N_track; // (X,Y) 二维绝对效率


// ==========================================================
// 方法二：以探测器阳极信号为基准
// ==========================================================
Double_t eff_x2  = Double_t(Nx)  / Na; 
Double_t eff_y2  = Double_t(Ny)  / Na; 
Double_t eff_xy2 = Double_t(Nxy) / Na; 

// 格式化输出结果
printf("\n--- PPAC2B 绝对探测效率 (基于 Tracking) ---\n");
printf("阳极触发效率   eff_a  = %5.2f%%\n", eff_a * 100);
printf(" X 绝对效率   eff_x1 = %5.2f%%\n", eff_x1 * 100);
printf(" Y 绝对效率   eff_y1 = %5.2f%%\n", eff_y1 * 100);
printf(" X-Y 联合效率 eff_xy1 = %5.2f%%\n", eff_xy1 * 100);

printf("\n--- PPAC2B 位置读出效率 (基于 Anode) ---\n");
printf(" X 读出效率   eff_x2 = %5.2f%%\n", eff_x2 * 100);
printf(" Y 读出效率   eff_y2 = %5.2f%%\n", eff_y2 * 100);
printf(" X-Y 联合效率 eff_xy2 = %5.2f%%\n", eff_xy2 * 100);
c1->Draw();

--- PPAC2B 绝对探测效率 (基于 Tracking) ---
阳极触发效率   eff_a  = 99.23%
 X 绝对效率   eff_x1 = 93.11%
 Y 绝对效率   eff_y1 = 92.98%
 X-Y 联合效率 eff_xy1 = 86.87%

--- PPAC2B 位置读出效率 (基于 Anode) ---
 X 读出效率   eff_x2 = 93.84%
 Y 读出效率   eff_y2 = 93.70%
 X-Y 联合效率 eff_xy2 = 87.54%

探测器层级	原始 Branch	物理意义 (已刻度)	讲义代码对应变量
PPAC 1 Layer A	`PPACF8[0][0]` `PPACF8[0][1]` `PPACF8[0][2]` `PPACF8[0][3]`	$X$ 坐标 (mm) $Y$ 坐标 (mm) $X$ 平面对应的 $Z$ 坐标 $Y$ 平面对应的 $Z$ 坐标	`xx[0]` `yy[0]` `xz[0]` `yz[0]`
PPAC 1 Layer B	`PPACF8[1][0~4]`	(同上结构，含 Anode 时间)	(本示例暂未使用)
PPAC 2 Layer A	`PPACF8[2][0~4]`	(同上结构)	`xx[1], yy[1]` `xz[1], yz[1]`
PPAC 2 Layer B (作为待测层 DUT)	`PPACF8[3][0~4]`	(同上结构)	`xx2b[0], yy2b[0]` `xz2b, yz2b` `anode2b`
PPAC 3	`PPACF8[4][0~4]`	(同上结构)	`xx[2], yy[2]` `xz[2], yz[2]`

2.2 PPAC 信号处理-II¶

1. Tracking¶

本实验中采用的重建策略¶

束流径迹重建结果（所有PPAC的信息都有x/y位置信息）¶

2. ROOT 文件中 TTree 的 Branch 说明¶

3. 径迹重建 (Tracking) 事例代码¶

靶点坐标 $(t_x, t_y)$ 与出射角 $\theta$ 的误差评估¶

1. 物理动机：为什么要逐事件评估追踪误差？¶

用法一：Quality Cut¶

用法二：Inverse-Variance Weighting¶

2. 协方差矩阵 (Covariance Matrix)¶

3. 原理推导：外推靶点 $(t_x)$ 与出射角 $(\theta)$ 的误差传递¶

A. 靶上位置 $t_x$ 的误差推导（线性传递）¶

B. 物理出射角 $\theta$ 的误差推导（非线性传递）¶

4. 位置分辨率¶

2.1 准直束刻度方法（离线测试 / 放射源）¶

2.2 多探测器径迹残差法¶

5. PPAC 探测效率¶

6. 分析代码框架¶

`tracking.h`¶

`tracking.C`¶

束流径迹¶

束流在靶上投影¶

残差，chi2/NDF¶

PPAC2B x,y,x-y的探测效率¶

0. 定义物理筛选条件 (TCut)¶

1. 计算穿过探测器灵敏面积的预期粒子数¶

2. 计算并对比两类探测效率¶

讨论：¶

作业¶

Task A: 物理靶点外推与重建信息记录¶

1. 径迹拟合与状态记录¶

2. 靶位坐标¶

3. 数据分析¶

Task B: PPAC 探测效率的多重验证¶

2.2 PPAC 信号处理-II¶

1. Tracking¶

本实验中采用的重建策略¶

束流径迹重建结果（所有PPAC的信息都有x/y位置信息）¶

2. ROOT 文件中 TTree 的 Branch 说明¶

3. 径迹重建 (Tracking) 事例代码¶

靶点坐标 $(t_x, t_y)$ 与出射角 $\theta$ 的误差评估¶

1. 物理动机：为什么要逐事件评估追踪误差？¶

用法一：Quality Cut¶

用法二：Inverse-Variance Weighting¶

2. 协方差矩阵 (Covariance Matrix)¶

3. 原理推导：外推靶点 $(t_x)$ 与出射角 $(\theta)$ 的误差传递¶

A. 靶上位置 $t_x$ 的误差推导（线性传递）¶

B. 物理出射角 $\theta$ 的误差推导（非线性传递）¶

4. 位置分辨率¶

2.1 准直束刻度方法（离线测试 / 放射源）¶

2.2 多探测器径迹残差法¶

5. PPAC 探测效率¶

6. 分析代码框架¶

tracking.h¶

tracking.C¶

束流径迹¶

束流在靶上投影¶

残差，chi2/NDF¶

PPAC2B x,y,x-y的探测效率¶

0. 定义物理筛选条件 (TCut)¶

1. 计算穿过探测器灵敏面积的预期粒子数¶

2. 计算并对比两类探测效率¶

讨论：¶

作业¶

Task A: 物理靶点外推与重建信息记录¶

1. 径迹拟合与状态记录¶

2. 靶位坐标¶

3. 数据分析¶

Task B: PPAC 探测效率的多重验证¶

`tracking.h`¶

`tracking.C`¶