检查 time-walk 修正

先检查各 cluster 内高能参考及全能区的分布，再检查不同 cluster 的高能时间差。比较相同能量选择下的峰位和宽度；全能区出现尾部时，应回到能量—时间差图查找来源。

%jsroot on

%%cpp -d
#include <iostream>

const int MAXHIT = 1024;
const int NCLUSTER = 12;
const int NSEG = 7;

// high-reference walk:
//   E_i vs t_i-t_j, with E_i > EminTarget and E_j > EminRef
TH2F *h_walk_cluster[NCLUSTER] = {0};

// all-energy walk:
//   E_i vs t_i-t_j, no energy gate
TH2F *h_walk_cluster_allE[NCLUSTER] = {0};

// high-energy time difference:
//   t_i-t_j, with E_i > EminRef and E_j > EminRef
TH1F *h_dt_cluster[NCLUSTER] = {0};
TH1F *h_dt_all = 0;

// all-energy time difference:
//   t_i-t_j, no energy gate
TH1F *h_dt_cluster_allE[NCLUSTER] = {0};
TH1F *h_dt_all_allE = 0;

TF1 *f_dt_cluster[NCLUSTER] = {0};
TF1 *f_dt_cluster_allE[NCLUSTER] = {0};

TF1 *f_dt_all = 0;
TF1 *f_dt_all_allE = 0;

// Make corrected walk plots and time-difference spectra

一次读取数据，同时填充高能参考和全能区的检查谱。二维图使用有方向的 i–j 组合；一维时间差按晶体编号排序，每一对只填一次。

%%cpp -d
void make_corrected_cluster_walk(
    const char *filename = "eurica_time_pair.root",
    double EminTarget = 30.0,
    double EminRef = 600.0)
{
  TH1::AddDirectory(kFALSE);
  for (int c = 0; c < NCLUSTER; c++) {
    if (h_walk_cluster[c]) {
      delete h_walk_cluster[c];
      h_walk_cluster[c] = 0;
    }

    h_walk_cluster[c] =
      new TH2F(Form("h_walk_cluster_%d", c),
               Form("cluster %d: corrected walk, E_{ref}>%.0f; t_i-t_j (ns); E_i",
                    c, EminRef),
               120, -600, 600,
               130, 0, 1300);

    h_walk_cluster[c]->SetDirectory(0);
    if (h_walk_cluster_allE[c]) {
      delete h_walk_cluster_allE[c];
      h_walk_cluster_allE[c] = 0;
    }

    h_walk_cluster_allE[c] =
      new TH2F(Form("h_walk_cluster_allE_%d", c),
               Form("cluster %d: corrected walk, all energy; t_i-t_j (ns); E_i",
                    c),
               120, -600, 600,
               130, 0, 1300);

    h_walk_cluster_allE[c]->SetDirectory(0);
    if (h_dt_cluster[c]) {
      delete h_dt_cluster[c];
      h_dt_cluster[c] = 0;
    }

    h_dt_cluster[c] =
      new TH1F(Form("h_dt_cluster_%d", c),
               Form("cluster %d: high-energy time difference; t_i-t_j (ns); counts",
                    c),
               240, -600, 600);

    h_dt_cluster[c]->SetDirectory(0);
    if (h_dt_cluster_allE[c]) {
      delete h_dt_cluster_allE[c];
      h_dt_cluster_allE[c] = 0;
    }

    h_dt_cluster_allE[c] =
      new TH1F(Form("h_dt_cluster_allE_%d", c),
               Form("cluster %d: all-energy time difference; t_i-t_j (ns); counts",
                    c),
               240, -600, 600);

    h_dt_cluster_allE[c]->SetDirectory(0);
  }
  if (h_dt_all) {
    delete h_dt_all;
    h_dt_all = 0;
  }

  h_dt_all =
    new TH1F("h_dt_all",
             "high-energy time difference: all clusters; t_i-t_j (ns); counts",
             240, -600, 600);

  h_dt_all->SetDirectory(0);
  if (h_dt_all_allE) {
    delete h_dt_all_allE;
    h_dt_all_allE = 0;
  }

  h_dt_all_allE =
    new TH1F("h_dt_all_allE",
             "all-energy time difference: all clusters; t_i-t_j (ns); counts",
             240, -600, 600);

  h_dt_all_allE->SetDirectory(0);

  TFile *fin = new TFile(filename);
  if (!fin || fin->IsZombie()) {
    std::cout << "cannot open " << filename << std::endl;
    return;
  }

  TTree *tree = (TTree*)fin->Get("tree");
  if (!tree) {
    std::cout << "cannot find tree in " << filename << std::endl;
    fin->Close();
    return;
  }

  int ghit;
  int gid[MAXHIT];
  double ge[MAXHIT];
  double gt[MAXHIT];
  if (tree->GetMaximum("ghit") > MAXHIT)
    throw std::runtime_error("Increase MAXHIT before reading branches");
  tree->SetBranchAddress("ghit", &ghit);
  tree->SetBranchAddress("gid", gid);
  tree->SetBranchAddress("ge", ge);
  tree->SetBranchAddress("gt", gt);

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

    tree->GetEntry(ientry);
    for (int a = 0; a < ghit; a++) {

      int cluster_a = gid[a] / 7;
      int seg_a = gid[a] % 7;
      if (cluster_a < 0 || cluster_a >= NCLUSTER) continue;
      if (seg_a < 0 || seg_a >= NSEG) continue;
      for (int b = 0; b < ghit; b++) {
        if (a == b) continue;

        int cluster_b = gid[b] / 7;
        int seg_b = gid[b] % 7;
        if (cluster_b != cluster_a) continue;
        if (seg_b < 0 || seg_b >= NSEG) continue;
        if (seg_b == seg_a) continue;

        double dt = gt[a] - gt[b];

        // Walk plots

        // all-energy walk: no energy gate
        h_walk_cluster_allE[cluster_a]->Fill(dt, ge[a]);

        // high-reference walk:
        // target hit can be low energy; reference hit is high energy
        if (ge[a] >= EminTarget && ge[b] >= EminRef) {
          h_walk_cluster[cluster_a]->Fill(dt, ge[a]);
        }

        // Time-difference spectra
        // Fill each detector pair once.
        // Use gid ordering, not array ordering.
        if (gid[a] >= gid[b]) continue;

        // all-energy time difference
        h_dt_cluster_allE[cluster_a]->Fill(dt);
        h_dt_all_allE->Fill(dt);

        // high-energy time difference
        if (ge[a] >= EminRef && ge[b] >= EminRef) {
          h_dt_cluster[cluster_a]->Fill(dt);
          h_dt_all->Fill(dt);
        }
      }
    }
  }

  fin->Close();
  for (int c = 0; c < NCLUSTER; c++) {
    std::cout << "cluster " << c
              << " high-ref walk entries = "
              << h_walk_cluster[c]->GetEntries()
              << ", all-energy walk entries = "
              << h_walk_cluster_allE[c]->GetEntries()
              << ", high-energy dt entries = "
              << h_dt_cluster[c]->GetEntries()
              << ", all-energy dt entries = "
              << h_dt_cluster_allE[c]->GetEntries()
              << std::endl;
  }

  std::cout << "all high-energy dt entries = "
            << h_dt_all->GetEntries()
            << std::endl;

  std::cout << "all all-energy dt entries = "
            << h_dt_all_allE->GetEntries()
            << std::endl;
}

%%cpp -d
void draw_corrected_cluster_walk()
{
  TCanvas *c_old = (TCanvas*)gROOT->FindObject("c_corr_cluster_walk");
  if (c_old) delete c_old;

  TCanvas *c = new TCanvas("c_corr_cluster_walk",
                           "corrected walk by cluster: high-reference",
                           900, 600);

  c->Divide(4, 3);
  gStyle->SetOptStat(0);
  for (int ic = 0; ic < NCLUSTER; ic++) {
    c->cd(ic + 1);
    gPad->SetLogz(0);
    if (h_walk_cluster[ic]) {
      h_walk_cluster[ic]->Draw("colz");
    }
  }

  c->Draw();
}

%%cpp -d
void draw_corrected_cluster_walk_allE()
{
  TCanvas *c_old = (TCanvas*)gROOT->FindObject("c_corr_cluster_walk_allE");
  if (c_old) delete c_old;

  TCanvas *c = new TCanvas("c_corr_cluster_walk_allE",
                           "corrected walk by cluster: all energy",
                           900, 600);

  c->Divide(4, 3);
  gStyle->SetOptStat(0);
  for (int ic = 0; ic < NCLUSTER; ic++) {
    c->cd(ic + 1);
    gPad->SetLogz(0);
    if (h_walk_cluster_allE[ic]) {
      h_walk_cluster_allE[ic]->Draw("colz");
    }
  }

  c->Draw();
}

逐 cluster 拟合高能时间差的中心峰；初值由最高 bin 给出。拟合窗口只描述 prompt 核心。

%%cpp -d
void draw_corrected_dt_by_cluster()
{
  TCanvas *c_old = (TCanvas*)gROOT->FindObject("c_corr_dt_cluster");
  if (c_old) delete c_old;

  TCanvas *c = new TCanvas("c_corr_dt_cluster",
                           "high-energy time difference by cluster",
                           900, 600);

  c->Divide(4, 3);
  gStyle->SetOptStat(0);
  for (int ic = 0; ic < NCLUSTER; ic++) {

    c->cd(ic + 1);
    gPad->SetLogy(0);
    if (!h_dt_cluster[ic]) continue;
    if (f_dt_cluster[ic]) {
      delete f_dt_cluster[ic];
      f_dt_cluster[ic] = 0;
    }
    if (h_dt_cluster[ic]->GetEntries() < 20) {
      h_dt_cluster[ic]->Draw("hist");
      continue;
    }

    int maxbin = h_dt_cluster[ic]->GetMaximumBin();
    double x0 = h_dt_cluster[ic]->GetBinCenter(maxbin);

    f_dt_cluster[ic] =
      new TF1(Form("f_dt_cluster_%d", ic),
              "gaus",
              x0 - 50.0,
              x0 + 50.0);

    f_dt_cluster[ic]->SetParameter(0, h_dt_cluster[ic]->GetMaximum());
    f_dt_cluster[ic]->SetParameter(1, x0);
    f_dt_cluster[ic]->SetParameter(2, 60.0);

    h_dt_cluster[ic]->Fit(f_dt_cluster[ic], "RQ0");

    h_dt_cluster[ic]->Draw("hist");
    f_dt_cluster[ic]->SetLineColor(kRed);
    f_dt_cluster[ic]->Draw("same");

    std::cout << "high-energy cluster " << ic
              << ": mean = " << f_dt_cluster[ic]->GetParameter(1)
              << " ns, sigma = " << f_dt_cluster[ic]->GetParameter(2)
              << " ns"
              << std::endl;
  }

  c->Draw();
}

用相同流程检查全能区时间差，与高能结果比较低能 time walk 的影响。

%%cpp -d
void draw_corrected_dt_by_cluster_allE()
{
  TCanvas *c_old = (TCanvas*)gROOT->FindObject("c_corr_dt_cluster_allE");
  if (c_old) delete c_old;

  TCanvas *c = new TCanvas("c_corr_dt_cluster_allE",
                           "all-energy time difference by cluster",
                           900, 600);

  c->Divide(4, 3);
  gStyle->SetOptStat(0);
  for (int ic = 0; ic < NCLUSTER; ic++) {

    c->cd(ic + 1);
    gPad->SetLogy(0);
    if (!h_dt_cluster_allE[ic]) continue;
    if (f_dt_cluster_allE[ic]) {
      delete f_dt_cluster_allE[ic];
      f_dt_cluster_allE[ic] = 0;
    }
    if (h_dt_cluster_allE[ic]->GetEntries() < 20) {
      h_dt_cluster_allE[ic]->Draw("hist");
      continue;
    }

    int maxbin = h_dt_cluster_allE[ic]->GetMaximumBin();
    double x0 = h_dt_cluster_allE[ic]->GetBinCenter(maxbin);

    f_dt_cluster_allE[ic] =
      new TF1(Form("f_dt_cluster_allE_%d", ic),
              "gaus",
              x0 - 100.0,
              x0 + 100.0);

    f_dt_cluster_allE[ic]->SetParameter(0, h_dt_cluster_allE[ic]->GetMaximum());
    f_dt_cluster_allE[ic]->SetParameter(1, x0);
    f_dt_cluster_allE[ic]->SetParameter(2, 60.0);

    h_dt_cluster_allE[ic]->Fit(f_dt_cluster_allE[ic], "RQ0");

    h_dt_cluster_allE[ic]->Draw("hist");
    f_dt_cluster_allE[ic]->SetLineColor(kRed);
    f_dt_cluster_allE[ic]->Draw("same");

    std::cout << "all-energy cluster " << ic
              << ": mean = " << f_dt_cluster_allE[ic]->GetParameter(1)
              << " ns, sigma = " << f_dt_cluster_allE[ic]->GetParameter(2)
              << " ns"
              << std::endl;
  }

  c->Draw();
}

将各 cluster 的高能时间差合并，显示拟合参数。

%%cpp -d
void draw_corrected_dt_all()
{
  TCanvas *c_old = (TCanvas*)gROOT->FindObject("c_corr_dt_all");
  if (c_old) delete c_old;

  TCanvas *c = new TCanvas("c_corr_dt_all",
                           "cumulative high-energy time difference",
                           800, 600);

  c->SetLogy();
  if (!h_dt_all) {
    std::cout << "h_dt_all does not exist. Run make_corrected_cluster_walk() first."
              << std::endl;
    return;
  }
  if (f_dt_all) {
    delete f_dt_all;
    f_dt_all = 0;
  }
  if (h_dt_all->GetEntries() < 20) {
    h_dt_all->Draw("hist");
    c->Draw();
    return;
  }

  int maxbin = h_dt_all->GetMaximumBin();
  double x0 = h_dt_all->GetBinCenter(maxbin);

  f_dt_all =
    new TF1("f_dt_all",
            "gaus",
            x0 - 100.0,
            x0 + 100.0);

  f_dt_all->SetParameter(0, h_dt_all->GetMaximum());
  f_dt_all->SetParameter(1, x0);
  f_dt_all->SetParameter(2, 60.0);

  h_dt_all->Fit(f_dt_all, "RQ0");

  h_dt_all->Draw("hist");
  f_dt_all->SetLineColor(kRed);
  f_dt_all->Draw("same");

  c->Draw();

  std::cout << "all clusters high-energy:"
            << " mean = " << f_dt_all->GetParameter(1)
            << " ns, sigma = " << f_dt_all->GetParameter(2)
            << " ns"
            << std::endl;
}

将全能区时间差合并，与高能参考的结果比较。

%%cpp -d
void draw_corrected_dt_all_allE()
{
  TCanvas *c_old = (TCanvas*)gROOT->FindObject("c_corr_dt_all_allE");
  if (c_old) delete c_old;

  TCanvas *c = new TCanvas("c_corr_dt_all_allE",
                           "cumulative all-energy time difference",
                           800, 600);

  c->SetLogy();
  if (!h_dt_all_allE) {
    std::cout << "h_dt_all_allE does not exist. Run make_corrected_cluster_walk() first."
              << std::endl;
    return;
  }
  if (f_dt_all_allE) {
    delete f_dt_all_allE;
    f_dt_all_allE = 0;
  }
  if (h_dt_all_allE->GetEntries() < 20) {
    h_dt_all_allE->Draw("hist");
    c->Draw();
    return;
  }

  int maxbin = h_dt_all_allE->GetMaximumBin();
  double x0 = h_dt_all_allE->GetBinCenter(maxbin);

  f_dt_all_allE =
    new TF1("f_dt_all_allE",
            "gaus",
            x0 - 200.0,
            x0 + 200.0);

  f_dt_all_allE->SetParameter(0, h_dt_all_allE->GetMaximum());
  f_dt_all_allE->SetParameter(1, x0);
  f_dt_all_allE->SetParameter(2, 60.0);

  h_dt_all_allE->Fit(f_dt_all_allE, "RQ0");

  h_dt_all_allE->Draw("hist");
  f_dt_all_allE->SetLineColor(kRed);

  f_dt_all_allE->Draw("same");

  c->Draw();

  std::cout << "all clusters all-energy:"
            << " mean = " << f_dt_all_allE->GetParameter(1)
            << " ns, sigma = " << f_dt_all_allE->GetParameter(2)
            << " ns"
            << std::endl;
}

Corrected walk plots ： high-energy reference

make_corrected_cluster_walk("eurica_time_pair.root");
draw_corrected_cluster_walk();

cluster 0 high-ref walk entries = 39878, all-energy walk entries = 257318, high-energy dt entries = 1175, all-energy dt entries = 128659
cluster 1 high-ref walk entries = 16848, all-energy walk entries = 119404, high-energy dt entries = 600, all-energy dt entries = 59702
cluster 2 high-ref walk entries = 38320, all-energy walk entries = 305710, high-energy dt entries = 1203, all-energy dt entries = 152855
cluster 3 high-ref walk entries = 33064, all-energy walk entries = 228046, high-energy dt entries = 1005, all-energy dt entries = 114023
cluster 4 high-ref walk entries = 39754, all-energy walk entries = 283122, high-energy dt entries = 1266, all-energy dt entries = 141561
cluster 5 high-ref walk entries = 44449, all-energy walk entries = 328816, high-energy dt entries = 1418, all-energy dt entries = 164408
cluster 6 high-ref walk entries = 66762, all-energy walk entries = 447384, high-energy dt entries = 1915, all-energy dt entries = 223692
cluster 7 high-ref walk entries = 0, all-energy walk entries = 0, high-energy dt entries = 0, all-energy dt entries = 0
cluster 8 high-ref walk entries = 13768, all-energy walk entries = 94514, high-energy dt entries = 393, all-energy dt entries = 47257
cluster 9 high-ref walk entries = 37841, all-energy walk entries = 272640, high-energy dt entries = 1418, all-energy dt entries = 136320
cluster 10 high-ref walk entries = 38395, all-energy walk entries = 257078, high-energy dt entries = 1281, all-energy dt entries = 128539
cluster 11 high-ref walk entries = 50814, all-energy walk entries = 356328, high-energy dt entries = 1558, all-energy dt entries = 178164
all high-energy dt entries = 13232
all all-energy dt entries = 1.47518e+06

Time-difference spectra ： high-energy reference

以下谱只统计同一 cluster 内的 pair；“all”表示把这些谱合并，不是跨 cluster 的时间差。Gaussian 拟合描述 prompt 峰核心，不能代表全部尾部的效率。

draw_corrected_dt_by_cluster();
draw_corrected_dt_all();

high-energy cluster 0: mean = 4.00458 ns, sigma = 22.4145 ns
high-energy cluster 1: mean = -2.98153 ns, sigma = 26.1133 ns
high-energy cluster 2: mean = 1.86119 ns, sigma = 21.7397 ns
high-energy cluster 3: mean = -2.30874 ns, sigma = 24.7694 ns
high-energy cluster 4: mean = 3.52715 ns, sigma = 24.2324 ns
high-energy cluster 5: mean = 1.45389 ns, sigma = 24.155 ns
high-energy cluster 6: mean = -1.15531 ns, sigma = 23.2892 ns
high-energy cluster 8: mean = -2.49556 ns, sigma = 22.8304 ns
high-energy cluster 9: mean = 3.77452 ns, sigma = 23.9822 ns
high-energy cluster 10: mean = -3.91272 ns, sigma = 25.485 ns
high-energy cluster 11: mean = -3.9237 ns, sigma = 25.3389 ns
all clusters high-energy: mean = 0.0868766 ns, sigma = 26.5682 ns

Corrected walk plots ： all-energy

draw_corrected_cluster_walk_allE();

Time-difference spectra ： all-energy

draw_corrected_dt_by_cluster_allE();
draw_corrected_dt_all_allE();

all-energy cluster 0: mean = 2.04605 ns, sigma = 42.8448 ns
all-energy cluster 1: mean = -2.16476 ns, sigma = 43.941 ns
all-energy cluster 2: mean = 4.45549 ns, sigma = 43.5069 ns
all-energy cluster 3: mean = -2.79493 ns, sigma = 45.3416 ns
all-energy cluster 4: mean = 3.68107 ns, sigma = 43.1379 ns
all-energy cluster 5: mean = 1.05046 ns, sigma = 45.4297 ns
all-energy cluster 6: mean = -0.384232 ns, sigma = 44.037 ns
all-energy cluster 8: mean = 1.51517 ns, sigma = 43.4519 ns
all-energy cluster 9: mean = -2.03786 ns, sigma = 44.4928 ns
all-energy cluster 10: mean = -1.00156 ns, sigma = 44.2354 ns
all-energy cluster 11: mean = -4.23082 ns, sigma = 48.3964 ns
all clusters all-energy: mean = 0.0134643 ns, sigma = 50.1657 ns

跨 cluster 的高能时间差

固定 cluster 0 为参考，画 t(cluster)−t(cluster 0)，两条能量均取 600–3000 keV。每个无序 pair 只填一次，不补相反符号。若某个 cluster 的 prompt ridge 系统性偏离零，需要先校正该常数 offset，再应用共同的 prompt gate。

TFile *fcross=TFile::Open("eurica_time_pair.root");
TTree *tcross=(TTree*)fcross->Get("tree");
int nh, id[MAXHIT]; double energy[MAXHIT], time[MAXHIT];
tcross->SetBranchAddress("ghit",&nh);
tcross->SetBranchAddress("gid",id);
tcross->SetBranchAddress("ge",energy);
tcross->SetBranchAddress("gt",time);
TH2F *hcross=new TH2F("hcross","High-energy cross-cluster timing;Cluster;#Deltat (ns)",11,0.5,11.5,160,-800,800);
for (Long64_t n=0;n<tcross->GetEntries();++n) {
    tcross->GetEntry(n);
    for (int i=0;i<nh;++i) {
        if (id[i]/7!=0 || energy[i]<600 || energy[i]>3000) continue;
        for (int j=0;j<nh;++j) {
            if (id[j]/7==0 || energy[j]<600 || energy[j]>3000) continue;
            hcross->Fill(id[j]/7,time[j]-time[i]);
        }
    }
}
TCanvas *ccross=new TCanvas("ccross","Cross-cluster timing",800,450);
hcross->Draw("colz");
ccross->Draw();