In [1]:
%%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
// ------------------------------------------------------------
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];
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;
}
// ------------------------------------------------------------
// Draw corrected walk plots: high-reference version
// ------------------------------------------------------------
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();
}
// ------------------------------------------------------------
// Draw corrected walk plots: all-energy version
// ------------------------------------------------------------
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();
}
// ------------------------------------------------------------
// Draw high-energy time-difference spectra by cluster
// ------------------------------------------------------------
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();
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();
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();
}
// ------------------------------------------------------------
// Draw all-energy time-difference spectra by cluster
// ------------------------------------------------------------
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();
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();
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();
}
// ------------------------------------------------------------
// Draw cumulative high-energy time-difference spectrum
// ------------------------------------------------------------
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();
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();
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;
}
// ------------------------------------------------------------
// Draw cumulative all-energy time-difference spectrum
// ------------------------------------------------------------
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();
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();
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¶
In [3]:
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¶
In [5]:
draw_corrected_dt_by_cluster();
draw_corrected_dt_all();
high-energy cluster 0: mean = 4.08921 ns, sigma = 23.3062 ns high-energy cluster 1: mean = 5.2245 ns, sigma = 25.145 ns high-energy cluster 2: mean = 3.75998 ns, sigma = 25.726 ns high-energy cluster 3: mean = 0.638946 ns, sigma = 23.0765 ns high-energy cluster 4: mean = 0.728055 ns, sigma = 24.1417 ns high-energy cluster 5: mean = 2.3464 ns, sigma = 23.7511 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.2949 ns, sigma = 23.3773 ns high-energy cluster 10: mean = -3.91272 ns, sigma = 25.485 ns high-energy cluster 11: mean = 1.36527 ns, sigma = 25.4948 ns all clusters high-energy: mean = 1.20075 ns, sigma = 26.9197 ns
Corrected walk plots ： all-energy¶
In [7]:
draw_corrected_cluster_walk_allE();
Time-difference spectra ： all-energy¶
In [9]:
draw_corrected_dt_by_cluster_allE();
draw_corrected_dt_all_allE();
all-energy cluster 0: mean = 1.69707 ns, sigma = 42.9055 ns all-energy cluster 1: mean = 5.74527 ns, sigma = 43.0695 ns all-energy cluster 2: mean = 5.91469 ns, sigma = 45.4182 ns all-energy cluster 3: mean = -1.22497 ns, sigma = 44.9568 ns all-energy cluster 4: mean = 0.792093 ns, sigma = 43.5485 ns all-energy cluster 5: mean = 2.01992 ns, sigma = 45.3595 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.20442 ns, sigma = 44.6229 ns all-energy cluster 10: mean = -1.00156 ns, sigma = 44.2354 ns all-energy cluster 11: mean = 1.29496 ns, sigma = 48.3454 ns all clusters all-energy: mean = 1.02185 ns, sigma = 50.3818 ns
In [ ]: