InttSeedTrackPerformance-Jingyu/macro/Pt_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 Pt_reso_vspt()
0012 {
0013     // 输入 ROOT file 名
0014     // TFile *infile = TFile::Open("/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/InttSeedTrackPerformance/output/Ptfunc_Performance_positron.root");
0015     // TFile *infile = TFile::Open("/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/InttSeedTrackPerformance/output/Ptfunc_Performance_positron_wTa.root");
0016     // electron
0017     // TFile *infile = TFile::Open("/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/InttSeedTrackPerformance/output/Ptfunc_Performance.root");
0018     // TFile *infile = TFile::Open("/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/InttSeedTrackPerformance/output/Ptfunc_Performance_electron_wTa.root");
0019
0020     // ML electron
0021     // TFile *infile = TFile::Open("/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/ML4Reco/version2/outputFile/pt_relative_error_INTT_CaloI.root");
0022     // TFile *infile = TFile::Open("/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/ML4Reco/version4/outputFile/pt_relative_error_INTT_CaloI.root");
0023     // TFile *infile = TFile::Open("/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/ML4Reco/combine2/outputFile/pt_relative_error_combined.root");
0024     TFile *infile = TFile::Open("/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/ML4Reco/combine2_gate/outputFile/pt_relative_error_combined.root");
0025
0026     // 获取 TH2D
0027     // TH2D *h2 = (TH2D *)infile->Get("h2_ptreso");
0028     // TH2D *h2 = (TH2D *)infile->Get("h2_pt_vs_relerr_all");
0029     // TH2D *h2 = (TH2D *)infile->Get("h2_pt_vs_relerr_combined");
0030     TH2D *h2 = (TH2D *)infile->Get("h2_pt_vs_relerr_0.0_10.0");
0031     if (!h2)
0032     {
0033         std::cerr << "Error: Cannot find h2_ptreso!" << std::endl;
0034         infile->Close();
0035         return;
0036     }
0037
0038     // 输出 ROOT file
0039     // TFile *outfile = new TFile("../output/projected_p_ptreso_TA.root", "RECREATE");
0040     // TFile *outfile = new TFile("../output/projected_e_ptreso.root", "RECREATE");
0041
0042     // TFile *outfile = new TFile("../output/ML/projected_e_ptreso_ML_E_woDropout.root", "RECREATE");
0043     // TFile *outfile = new TFile("../output/ML/projected_e_ptreso_ML_Phi.root", "RECREATE");
0044     TFile *outfile = new TFile("../output/ML/projected_e_ptreso_ML_Combined.root", "RECREATE");
0045
0046     // 准备循环
0047     int nxbins = h2->GetNbinsX();
0048     int group_size = 1;
0049     int n_groups = 100 / group_size;
0050     if (nxbins % group_size != 0) n_groups += 1;  // 处理非整除情况
0051
0052     // 准备 TGraphErrors 存 mean 和 sigma
0053     std::vector<double> vec_pt_center;
0054     std::vector<double> vec_mean, vec_mean_err;
0055     std::vector<double> vec_sigma, vec_sigma_err;
0056
0057     for (int i = 5; i < n_groups; ++i)
0058     {
0059         int first_bin = i * group_size + 1; // bin 从 1 开始
0060         int last_bin = std::min((i + 1) * group_size, nxbins);
0061
0062         // 投影
0063         TString proj_name = Form("proj_y_binGroup%d", i);
0064         TString proj_title = Form("Phi_truth - Phi_reco %d to %d GeV", i, i+1);
0065         TH1D *proj = h2->ProjectionY(proj_name, first_bin, last_bin);
0066         proj->SetTitle(proj_title);
0067
0068         // 寻峰
0069         int peak_bin = proj->GetMaximumBin();
0070         double peak_center = proj->GetBinCenter(peak_bin);
0071
0072         // 找出 content 最高的 3 个 bin，并计算它们的 bin center 平均值
0073         int nBins = proj->GetNbinsX();
0074         std::vector<std::pair<double,int>> bins;
0075         bins.reserve(nBins);
0076
0077         for (int i = 1; i <= nBins; ++i) {
0078             bins.emplace_back(proj->GetBinContent(i), i);
0079         }
0080
0081         // 按 content 从大到小排序
0082         std::sort(bins.begin(), bins.end(),
0083                   [](auto &a, auto &b){ return a.first > b.first; });
0084
0085         // 取前 3 个峰，累加它们的 bin center
0086         int nPeaks = 3;
0087         double sumCenters = 0;
0088         for (int i = 0; i < nPeaks && i < (int)bins.size(); ++i) {
0089             int bin = bins[i].second;
0090             sumCenters += proj->GetBinCenter(bin);
0091         }
0092
0093         // 计算平均 peak center
0094         double peak_center_avg = sumCenters / nPeaks;
0095
0096         // 拟合
0097         // TF1 *fit_func = new TF1(Form("gaus_fit_%d", i), "gaus", peak_center - 0.10, peak_center + 0.10);
0098         TF1 *fit_func = new TF1(Form("gaus_fit_%d", i), "gaus", peak_center_avg - 0.10, peak_center_avg + 0.10);
0099         proj->Fit(fit_func, "RQ");
0100
0101         // 取拟合参数
0102         double fit_mean = fit_func->GetParameter(1);
0103         double fit_mean_err = fit_func->GetParError(1);
0104         double fit_sigma = fit_func->GetParameter(2);
0105         double fit_sigma_err = fit_func->GetParError(2);
0106
0107         // 获取 pt bin 中心
0108         double pt_low = h2->GetXaxis()->GetBinLowEdge(first_bin);
0109         double pt_high = h2->GetXaxis()->GetBinUpEdge(last_bin);
0110         double pt_center = 0.5 * (pt_low + pt_high);
0111
0112         // 存入 vector
0113         vec_pt_center.push_back(pt_center);
0114         vec_mean.push_back(fit_mean);
0115         vec_mean_err.push_back(fit_mean_err);
0116         vec_sigma.push_back(fit_sigma);
0117         vec_sigma_err.push_back(fit_sigma_err);
0118
0119         // 保存 proj 和拟合
0120         proj->Write();
0121         // fit_func->Write();
0122     }
0123
0124     // 生成 TGraphErrors
0125     int npoints = vec_pt_center.size();
0126     TGraphErrors *gr_mean = new TGraphErrors(npoints);
0127     TGraphErrors *gr_sigma = new TGraphErrors(npoints);
0128
0129     for (int i = 0; i < npoints; ++i)
0130     {
0131         gr_mean->SetPoint(i, vec_pt_center[i], vec_mean[i]);
0132         gr_mean->SetPointError(i, 0.0, vec_mean_err[i]);
0133
0134         gr_sigma->SetPoint(i, vec_pt_center[i], vec_sigma[i]);
0135         gr_sigma->SetPointError(i, 0.0, vec_sigma_err[i]);
0136     }
0137
0138     gr_mean->SetTitle("Peak value vs Pt;Pt (GeV);Peak Mean");
0139     gr_mean->SetMarkerStyle(20);
0140     gr_mean->SetMarkerColor(kBlue);
0141
0142     gr_sigma->SetMinimum(0.0);
0143     gr_sigma->SetMaximum(0.2);
0144     gr_sigma->SetTitle("Resolution (Sigma) vs Pt;Pt (GeV);Sigma");
0145     gr_sigma->SetMarkerStyle(21);
0146     gr_sigma->SetMarkerColor(kRed);
0147
0148     // 保存 TGraphErrors
0149     gr_mean->Write("gr_mean_vs_pt");
0150     gr_sigma->Write("gr_sigma_vs_pt");
0151
0152     // 可以画出来
0153     TCanvas *c1 = new TCanvas("c1", "Peak Mean vs Pt", 800, 600);
0154     gr_mean->Draw("AP");
0155     c1->SaveAs("../output/mean_vs_pt.pdf");
0156
0157     TCanvas *c2 = new TCanvas("c2", "Sigma vs Pt", 800, 600);
0158     gr_sigma->Draw("AP");
0159     c2->SaveAs("../output/sigma_vs_pt.pdf");
0160
0161     // 关闭文件
0162     outfile->Close();
0163     infile->Close();
0164
0165     std::cout << "Done. Projected TH1D saved to root file" << std::endl;
0166 }