InttSeedTrackPerformance-Jingyu/macro/Phi_reso_vspt.C

0001 #include <iostream>
0002 #include <string>
0003 #include <vector>
0004 #include <TChain.h>
0005 #include <TTree.h>
0006 #include <TH1F.h>
0007 #include <TFile.h>
0008 #include <TMath.h>
0009 #include <TProfile.h>
0010
0011 void Phi_reso_vspt()
0012 {
0013     // 输入 ROOT file 名
0014     // TFile *infile = TFile::Open("/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/ParticleGen/output/calo_position_electron_dphi_ptbin.root");
0015     // TFile *infile = TFile::Open("/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/ParticleGen/output/calo_positron_dphi_ptbin_wC.root");
0016     TFile *infile = TFile::Open("/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/ParticleGen/output/calo_electron_dphi_ptbin_wC.root");
0017
0018     // 获取 TH2D
0019     TH2D *h2 = (TH2D *)infile->Get("h2_ptreso");
0020     // TH2D *h2 = (TH2D *)infile->Get("h2_dphi_distri_test");
0021     if (!h2)
0022     {
0023         std::cerr << "Error: Cannot find h2_dphi_distri_test!" << std::endl;
0024         infile->Close();
0025         return;
0026     }
0027
0028     // 输出 ROOT file
0029     // TFile *outfile = new TFile("../output/projected_P.root", "RECREATE");
0030     TFile *outfile = new TFile("../output/projected_E.root", "RECREATE");
0031
0032     // 准备循环
0033     int nxbins = h2->GetNbinsX();
0034     int group_size = 5;
0035     int n_groups = 50 / group_size;
0036     if (nxbins % group_size != 0) n_groups += 1;  // 处理非整除情况
0037
0038     // 准备 TGraphErrors 存 mean 和 sigma
0039     std::vector<double> vec_pt_center;
0040     std::vector<double> vec_mean, vec_mean_err;
0041     std::vector<double> vec_sigma, vec_sigma_err;
0042
0043     for (int i = 0; i < n_groups; ++i)
0044     {
0045         int first_bin = i * group_size + 1; // bin 从 1 开始
0046         int last_bin = std::min((i + 1) * group_size, nxbins);
0047
0048         // 投影
0049         TString proj_name = Form("proj_y_binGroup%d", i);
0050         TString proj_title = Form("Phi_truth - Phi_reco %d to %d GeV", i, i+1);
0051         TH1D *proj = h2->ProjectionY(proj_name, first_bin, last_bin);
0052         proj->SetTitle(proj_title);
0053
0054         // 寻峰
0055         int peak_bin = proj->GetMaximumBin();
0056         double peak_center = proj->GetBinCenter(peak_bin);
0057
0058         // 拟合
0059         TF1 *fit_func = new TF1(Form("gaus_fit_%d", i), "gaus", peak_center - 0.005, peak_center + 0.005);
0060         proj->Fit(fit_func, "RQ");
0061
0062         // 取拟合参数
0063         double fit_mean = fit_func->GetParameter(1);
0064         double fit_mean_err = fit_func->GetParError(1);
0065         double fit_sigma = fit_func->GetParameter(2);
0066         double fit_sigma_err = fit_func->GetParError(2);
0067
0068         // 获取 pt bin 中心
0069         double pt_low = h2->GetXaxis()->GetBinLowEdge(first_bin);
0070         double pt_high = h2->GetXaxis()->GetBinUpEdge(last_bin);
0071         double pt_center = 0.5 * (pt_low + pt_high);
0072
0073         // 存入 vector
0074         vec_pt_center.push_back(pt_center);
0075         vec_mean.push_back(fit_mean);
0076         vec_mean_err.push_back(fit_mean_err);
0077         vec_sigma.push_back(fit_sigma);
0078         vec_sigma_err.push_back(fit_sigma_err);
0079
0080         // 保存 proj 和拟合
0081         proj->Write();
0082         // fit_func->Write();
0083     }
0084
0085     // 生成 TGraphErrors
0086     int npoints = vec_pt_center.size();
0087     TGraphErrors *gr_mean = new TGraphErrors(npoints);
0088     TGraphErrors *gr_sigma = new TGraphErrors(npoints);
0089
0090     for (int i = 0; i < npoints; ++i)
0091     {
0092         gr_mean->SetPoint(i, vec_pt_center[i], vec_mean[i]);
0093         gr_mean->SetPointError(i, 0.0, vec_mean_err[i]);
0094
0095         gr_sigma->SetPoint(i, vec_pt_center[i], vec_sigma[i]);
0096         gr_sigma->SetPointError(i, 0.0, vec_sigma_err[i]);
0097     }
0098
0099     gr_mean->SetTitle("Peak value vs Pt;Pt (GeV);Peak Mean");
0100     gr_mean->SetMarkerStyle(20);
0101     gr_mean->SetMarkerColor(kBlue);
0102
0103     gr_sigma->SetTitle("Resolution (Sigma) vs Pt;Pt (GeV);Sigma");
0104     gr_sigma->SetMarkerStyle(21);
0105     gr_sigma->SetMarkerColor(kRed);
0106
0107     // 保存 TGraphErrors
0108     gr_mean->Write("gr_mean_vs_pt");
0109     gr_sigma->Write("gr_sigma_vs_pt");
0110
0111     // 可以画出来
0112     TCanvas *c1 = new TCanvas("c1", "Peak Mean vs Pt", 800, 600);
0113     gr_mean->Draw("AP");
0114     c1->SaveAs("../output/mean_vs_pt.pdf");
0115
0116     TCanvas *c2 = new TCanvas("c2", "Sigma vs Pt", 800, 600);
0117     gr_sigma->Draw("AP");
0118     c2->SaveAs("../output/sigma_vs_pt.pdf");
0119
0120     // 关闭文件
0121     outfile->Close();
0122     infile->Close();
0123
0124     std::cout << "Done. Projected TH1D saved to root file" << std::endl;
0125 }