1.2 读取ROOT的tree数据，进行逐事件分析¶

  1. 学习进行逐事件分析的方法
  2. 掌握中子探测器位置刻度、飞行时间刻度的方法
  • 从实验测量值 $t_u, t_d, q_u, q_d$ 计算得到物理量
  1. 学习数据分析的思路

tree.root 内 TTree格式¶

Image

代码1 : 逐事件读取文件内TTree变量信息。¶

  • 将下列代码保存到 readTree.cc
  • ROOT环境下运行: .x readTree.cc
//readTree.C

TH1D *hTOF = new TH1D("hTOF", "Time of flight", 1000, 0, 100);
void readTree()
{
// 1.打开文件，得到TTree指针
  TFile *ipf = new TFile("tree.root");//打开ROOT文件
  if (ipf->IsZombie()) {
   cout << "Error opening file!" << endl;
   exit(-1);
  }
  ipf->cd();
  TTree *tree=(TTree*)ipf->Get("tree");//得到名字为“tree”的TTree指针

//2. 声明tree的Branch变量

  Double_t x;
  Double_t e;
  int pid;
  Double_t tof, ctof;
  Double_t tu, td;
  Double_t qu, qd;

//3. 将变量指向对应Branch的地址
  tree->SetBranchAddress("ctof", &ctof);//将ROOT文件内tree内名为"ctof"的branch的数据的指针指向ctof的变量。
  tree->SetBranchAddress("tof", &tof);  
  tree->SetBranchAddress("pid", &pid);
  tree->SetBranchAddress("tu", &tu);   
  tree->SetBranchAddress("td", &td);
  tree->SetBranchAddress("qu", &qu);   
  tree->SetBranchAddress("qd", &qd);

//4. 逐事件读取tree的branch数据
//以下三行为固定格式
  Long64_t nentries = tree->GetEntries();  //得到tree的事件总数
  for(Long64_t jentry = 0; jentry < nentries; jentry++) {//对事件进行遍历
    tree->GetEntry(jentry);//将第jentry个事件数据填入步骤3.中指向的变量。

    cout <<"jentry: " << jentry << " tof: " << tof <<" pid: " << pid << endl;

  }   
   ipf->Close();  

}

代码2 : 逐事件读取已有文件的TTree变量，计算生成新变量，并将其存入新的root文件。¶

//以下代码中标 ... 的部分需要用户填写.

//重复readTree.C 中的1.-3.部分的代码
...

//4. 与新变量相关变量声明

  //将数据写入新ROOT文件
  //calibration parameters
  Double_t kt = 4.0;
  Double_t bt = 3.1;
  ...
  //new tree parameters
  Double_t tx;
  ...

//5. 新ROOT文件声明，内部TTree数据结构声明

  TFile *opf = new TFile("tree2.root","recreate");  //输出文件名
  opt->cd(); 
  TTree *opt = new TTree("tree","tree");            //输出文件内TTree
  opt->SetBrach("tx", &tx, "tx/D");
  ...

//5. 逐事件读取原文件tree的数据
//以下三行为固定格式
  Long64_t nentries = tree->GetEntries();  
  for(Long64_t jentry = 0; jentry < nentries; jentry++) {
    tree->GetEntry(jentry);

   // calculate new parameters
   tx = kx * (tu - td) + kb;
   ...

   opt->Fill();

   if(jentry % 100000 == 0) 
        cout << "process " << jentry << " of " << nentries << endl;
  }

  ipf->Close();

//6. 写入新ROOT文件    
  opt->Write();
  opf->Close();
}

练习：¶

  • 利用代码1生成ROOT文件，逐事件读出 tu + td, ln(qu * qd）数据。

分析按上述步骤生成的ROOT文件¶

In [1]:
//%jsroot on
TH1D *tdiff = new TH1D("tdiff", "td-tu", 140, -20, 50);  
TCanvas *c1 = new TCanvas("c1", "c1");

Double_t tu, td;
Double_t ctof,tof,x;
int pid;
TFile *ipf = new TFile("tree.root");//打开ROOT文件
TTree *tree = (TTree*)ipf->Get("tree");//得到tree的指针

tree->SetBranchAddress("tu", &tu);   
tree->SetBranchAddress("td", &td);
tree->SetBranchAddress("ctof", &ctof);
tree->SetBranchAddress("tof", &tof);
tree->SetBranchAddress("x", &x);
tree->SetBranchAddress("pid", &pid);
In [2]:
Long64_t nentries = tree->GetEntries();//得到事件总数
for(Long64_t jentry = 0; jentry < nentries; jentry++) {//对每个事件进行遍历
    tree->GetEntry(jentry);
    tdiff->Fill(td - tu);  // if(ng == 1) tx->Fill(tu-td), 只写入满足给定条件的事件      
  }
tdiff->Draw();
c1->Draw();

确定均匀分布边界的方法-微分法¶

  • 原理见参考文献2.
In [3]:
TH1D *dtd = new TH1D("dtd", "dN/dx", 141, -20.25, 50.25);

for(int i = 1;i <= tdiff->GetNbinsX();i++) {
    Double_t df = tdiff->GetBinContent(i+1) - tdiff->GetBinContent(i);
    dtd->Fill(tdiff->GetBinCenter(i+1), df);
}
dtd->Sumw2(0);//不显示传递误差, 等价于 dtd->Draw("hist");
dtd->Fit("gaus", "", "", -14, -8);//txl
c1->Draw();

 FCN=95.674 FROM MIGRAD    STATUS=CONVERGED      81 CALLS          82 TOTAL
                     EDM=9.09603e-07    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  Constant     2.92771e+02   1.19978e+01   4.62970e-02   1.11585e-04
   2  Mean        -1.12749e+01   2.10422e-02   1.00580e-04   3.96569e-04
   3  Sigma        6.30115e-01   1.56904e-02   2.60032e-05   2.38449e-01
In [4]:
dtd->Fit("gaus", "", "", 39.5, 43);//为什么拟合有问题？
c1->Draw();

 FCN=1055 FROM HESSE     STATUS=FAILED         11 CALLS         131 TOTAL
                     EDM=0    STRATEGY= 1  ERROR MATRIX UNCERTAINTY 100.0 per cent
  EXT PARAMETER                APPROXIMATE        STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  Constant    -9.18837e+05   1.41421e+00   0.00000e+00   0.00000e+00
   2  Mean        -1.12749e+01   1.41421e+00   0.00000e+00   0.00000e+00
   3  Sigma        6.30115e-01   5.11833e+00   0.00000e+00   0.00000e+00
In [5]:
TF1 *f1 = new TF1("f1","[0]*TMath::Exp(-0.5*((x-[1])/[2])^2)",39.5,45);//定义函数的方法 TF1
//gaus:f(x) = p0*exp(-0.5*((x-p1)/p2)^2)

设定拟合函数的初值¶

In [6]:
//进行参数拟合时，设置合理的初始参数至关重要！
c1->SetLogy(0);
f1->SetParameter(0,-350);//大致估计范围
f1->SetParameter(1,41.5);
f1->SetParameter(2,0.5);
dtd->Fit("f1","R");
dtd->Draw();
c1->Draw();

 FCN=200.49 FROM MIGRAD    STATUS=CONVERGED     128 CALLS         129 TOTAL
                     EDM=9.39984e-07    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  p0          -3.37168e+02   1.47578e+01   7.88452e-02  -2.79834e-05
   2  p1           4.21311e+01   1.98002e-02   1.21753e-04   7.59427e-02
   3  p2          -5.19832e-01   1.45390e-02   7.21177e-05   2.88199e-02

探测器时间差-位置的刻度方法¶

  • $t_x=t_d-t_u$
  • $t_{xl}=-11.76,$ $\space t_{xr}=41.57$ - 左右边界
  • $t_{xoff}=(t_{xl}+t_{xr})/2=14.90$. - offset
$$ x=\frac{2L}{t_{xr}-t_{xl}}*(t_x-t_{xoff})=3.750*(t_x-14.90) $$
In [7]:
TH1D *htx = new TH1D("htx", "htx", 500, -120, 120);
TH2D *hdx = new TH2D("hdx", "htx-hx:hx", 100, -20, 20, 500, -120, 120);

for(Long64_t jentry = 0; jentry < nentries; jentry++) {//对事件进行遍历
    tree->GetEntry(jentry);
    Double_t tx = 3.750 * (td - tu - 14.90);
    htx->Fill(tx);
    hdx->Fill(tx - x,x);//difference
  }
htx->Draw();
c1->Draw();

评估计算结果¶

  • 实际参数 x,计算参数tx
  • 观察 tx-x 是否与x有关联性：二维图
  • 观察 tx-x的中心值是否为零：一维图
In [8]:
hdx->Draw("colz");//为一条竖线，没有关联
c1->Draw();
In [9]:
TH1D *hdx1 = hdx->ProjectionX("projx");
hdx1->Draw();
hdx1->Fit("gaus");//如果中心值 mean x偏离0, 说明确定td-tu的边界时可能有些问题。[3],[6]
c1->Draw();

 FCN=35.5366 FROM MIGRAD    STATUS=CONVERGED      57 CALLS          58 TOTAL
                     EDM=9.30043e-07    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  Constant     7.05945e+03   2.73424e+01   6.58995e-02   3.18615e-06
   2  Mean        -1.05085e-02   7.14929e-03   2.11001e-05  -1.88113e-01
   3  Sigma        2.25974e+00   5.05437e-03   1.79614e-06  -2.07212e-01

对TOF进行飞行距离修正¶

In [10]:
tree->Draw("ctof>>(1000,20,120)");
c1->Draw();
In [11]:
TH2D *hgtofx=new TH2D("hgtofx", "hgtofx", 100, -120, 120, 100, 39, 48);
TH1D *hgctof=new TH1D("hgctof", "hgctof", 100, 39, 48);
In [12]:
for(Long64_t jentry = 0; jentry < nentries; jentry++) {//对事件进行遍历
    tree->GetEntry(jentry);
    Double_t tx = 3.745 * (td - tu - 14.75);
    if(ctof > 42 && ctof < 44.5) { //选择gamma
        Double_t d = TMath::Sqrt(502.5 * 502.5 + tx * tx);
        Double_t ctofa = (ctof)/d * 500.;//normalized to 500cm
        hgtofx->Fill(tx,ctofa);
        if(abs(tx) < 5) hgctof->Fill(ctof);//gamma hits the center of the det.
    }
  }
hgtofx->Draw("colz");
c1->Draw();

TOF已进行飞行距离修正，修正后的gamma的TOF与应与x无关联性。上图说明当前做法存在问题。¶

TOF绝对刻度的方法¶

实际测得的时间信号还有一个额外附加的常数项，$t_{Off}$：来源包括PMT度越时间，和电缆，电子学的传输时间。

  • 这个常数项一般是无法直接测定的，需要根据实验数据计算。

方法¶

  • $t_u=TOF+(L-x)/v_{sc}+t_{uoff}$
  • $t_d=TOF+(L+x)/v_{sc}+t_{doff}$
  • $TOF=(t_u+t_d)/2 -L/v_{sc}-(t_{uoff}+t_{doff})/2$
  • 将上式常数项设为C，$TOF=(t_u+t_d)/2 +C$

对于$\gamma$ 有,¶

  • $ TOF=(t_u+t_d)/2 +C=c*d $

  • 可得， $ C=c*d-(t_u+t_d)/2 $

In [13]:
TH1D *hC = new TH1D("hC", "", 300, -40, -10);
TH2D *hCx = new TH2D("hCx", "", 100, -105, 105, 300, -30, -23);
In [14]:
for(Long64_t jentry = 0; jentry < nentries; jentry++) {
    tree->GetEntry(jentry);
    Double_t tx = 3.745 * (td - tu - 14.75);
    
    if(tx > -100 && tx < 100 && ctof > 42 && ctof < 44.5) {//gamma
        Double_t d = TMath::Sqrt(502.5 * 502.5 + tx * tx);
        hC->Fill(3.333 * d * 0.01 - ctof);
        hCx->Fill(tx, 3.333 * d * 0.01 - ctof);
    }
  }
hC->Draw();
c1->Draw();

从上图可得 C=-26.2 ns¶

检验C与x是否有关联¶

In [15]:
hCx->Draw("colz");//检验C与x是否有关联
c1->Draw();
In [16]:
TH1D *htof = new TH1D("tof", "", 500, 0, 100);//real TOF, normalized to 500cm
TH1D *htofc = new TH1D("tofc", "", 500, 0, 100);//calculated TOF, normalized to 500cm
TH2D *htofcx = new TH2D("htofcx", "htofcx", 100, -120, 120, 100, 12, 22);
In [17]:
for(Long64_t jentry = 0; jentry < nentries; jentry++) {
    tree->GetEntry(jentry);
    
    Double_t tx = 3.745 * (td - tu - 14.75);
    Double_t d = TMath::Sqrt(502.5 * 502.5 + tx * tx);    
    htof->Fill(tof/d*500);
    htofc->Fill((ctof - 26.2)/d * 500);
    if(ctof > 42 && ctof < 44.5)
        htofcx->Fill(tx,(ctof - 26.2)/d * 500);
  }
In [18]:
c1->SetLogy();
htofc->Draw();//修正后的飞行时间谱
htof->SetLineColor(kRed);//kGreen,kBlack
htof->Draw("same");
c1->Draw();

检验tofc与x是否有关联¶

In [19]:
c1->SetLogy(0);
htofcx->Draw("colz");
c1->Draw();