InttSeedTrackPerformance-Jingyu/macro/comparePNcharge.C

0001 void comparePNcharge()
0002 {
0003     // 输入文件路径
0004     TString file1P_name = "/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/ParticleGen/output/calo_position_positron.root";
0005     TString file2N_name = "/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/ParticleGen/output/calo_position_electron.root";
0006
0007     // 打开输入文件
0008     TFile *file1_Pcharge = TFile::Open(file1P_name, "READ");
0009     TFile *file2_Ncharge = TFile::Open(file2N_name, "READ");
0010
0011     if (!file1_Pcharge || file1_Pcharge->IsZombie())
0012     {
0013       std::cerr << "无法打开 " << file1P_name << std::endl;
0014       return;
0015     }
0016     if (!file2_Ncharge || file2_Ncharge->IsZombie())
0017     {
0018       std::cerr << "无法打开 " << file2N_name << std::endl;
0019       return;
0020     }
0021
0022     // 从两个文件中读取 TH1D
0023     TH1F *h1_DR_pr1_clustinnr_P = nullptr;
0024     TH1F *h1_DR_pr1_clustinnr_N = nullptr;
0025     file1_Pcharge->GetObject("h1_DR_pr1_clustinnr",h1_DR_pr1_clustinnr_P);
0026     file2_Ncharge->GetObject("h1_DR_pr1_clustinnr",h1_DR_pr1_clustinnr_N);
0027
0028     // // 从两个文件中读取 TH1D
0029     // TH1* h_temp = nullptr;
0030     // file1_Pcharge->GetObject("h1_dphi_distri_pr1_clustinnr", h_temp);
0031
0032     // 获取两个 TProfile
0033     TProfile *TP_DR_pr1_clustinnr_pt_P = nullptr;
0034     TProfile *TP_DR_pr1_clustinnr_pt_N = nullptr;
0035     file1_Pcharge->GetObject("TP_DR_pr1_clustinnr_pt", TP_DR_pr1_clustinnr_pt_P);
0036     file2_Ncharge->GetObject("TP_DR_pr1_clustinnr_pt", TP_DR_pr1_clustinnr_pt_N);
0037
0038     // 准备新 profile，结构与 prof1 相同
0039     auto TP_DR_pr1_clustinnr_pt_PN  = new TProfile("TP_DR_pr1_clustinnr_pt_PN", "TP_DR_pr1_clustinnr_pt_PN;charge(-e)*pt(GeV);R(cm)", 22, -11, 11, 0, 20);
0040
0041     // 合并 P charge
0042     for (int i = 1; i <= TP_DR_pr1_clustinnr_pt_P->GetNbinsX(); ++i)
0043     {
0044         double x = TP_DR_pr1_clustinnr_pt_P->GetBinCenter(i);
0045         double y = TP_DR_pr1_clustinnr_pt_P->GetBinContent(i);
0046         double ey = TP_DR_pr1_clustinnr_pt_P->GetBinError(i);
0047         double n  = TP_DR_pr1_clustinnr_pt_P->GetBinEntries(i);
0048
0049         if (n <= 0) continue;
0050
0051         // 反推 RMS²，模拟填充
0052         double rms2 = ey * ey * n;
0053         double rms = std::sqrt(rms2 / n); // RMS
0054
0055         // 模拟填入 n 次，带有 RMS 的偏移
0056         for (int j = 0; j < (int)n; ++j)
0057         {
0058             double y_j = gRandom->Gaus(y, rms);  // 模拟原始测量值
0059             TP_DR_pr1_clustinnr_pt_PN->Fill(x, y);
0060         }
0061     }
0062
0063     TProfile *flipped = new TProfile("flipped_profile", "Flipped profile", 11, 0, 11, 0, 20);
0064     // 合并 N charge（反转X轴）
0065     for (int i = 1; i <= TP_DR_pr1_clustinnr_pt_N->GetNbinsX(); ++i)
0066     {
0067         // 翻转 X
0068         int flipped_bin = TP_DR_pr1_clustinnr_pt_N->GetNbinsX() - i + 1;
0069
0070         double x = TP_DR_pr1_clustinnr_pt_N->GetBinCenter(i);
0071         double flipped_x = -x;
0072         double y = TP_DR_pr1_clustinnr_pt_N->GetBinContent(i);
0073         double ey = TP_DR_pr1_clustinnr_pt_N->GetBinError(i);
0074         double n  = TP_DR_pr1_clustinnr_pt_N->GetBinEntries(i);
0075
0076         if (n <= 0) continue;
0077
0078         double rms2 = ey * ey * n;
0079         // double rms = std::sqrt(rms2 / n);
0080         double rms = ey * ey * n;
0081
0082
0083         for (int j = 0; j < (int)n; ++j)
0084         {
0085             double y_j = gRandom->Gaus(y, rms);  // 模拟原始测量值
0086             TP_DR_pr1_clustinnr_pt_PN->Fill(flipped_x, y);
0087         }
0088     }
0089     TP_DR_pr1_clustinnr_pt_PN->SetMarkerStyle(29);
0090
0091
0092     // 创建画布并画图 ---------------------------------------------------------------------------------------------------
0093     TCanvas *c1 = new TCanvas("h1_DR_pr1_clustinnr_PN", "P vs N", 800, 600);
0094     c1->cd();
0095
0096     // 设置样式
0097     h1_DR_pr1_clustinnr_P->SetLineColor(kRed);
0098     h1_DR_pr1_clustinnr_P->SetLineWidth(2);
0099     h1_DR_pr1_clustinnr_P->SetTitle("Comparison of e+ and e-");
0100     h1_DR_pr1_clustinnr_P->Draw("");
0101
0102     h1_DR_pr1_clustinnr_N->SetLineColor(kBlue);
0103     h1_DR_pr1_clustinnr_N->SetLineWidth(2);
0104     h1_DR_pr1_clustinnr_N->Draw("SAME");
0105
0106     // 添加图例
0107     TLegend *leg1 = new TLegend(0.65, 0.75, 0.88, 0.88);
0108     leg1->AddEntry(h1_DR_pr1_clustinnr_P, "Positron", "l");
0109     leg1->AddEntry(h1_DR_pr1_clustinnr_N, "Electron", "l");
0110     leg1->Draw();
0111
0112     // 保存图像
0113     c1->SaveAs("comparePNcharge_R.pdf");
0114
0115     TCanvas *c2 = new TCanvas("h1_dphi_distri_pr1_clustinnr_PN", "P vs N", 800, 600);
0116     c2->cd();
0117     // // 设置样式
0118     // h1_dphi_distri_pr1_clustinnr_P->SetLineColor(kRed);
0119     // h1_dphi_distri_pr1_clustinnr_P->SetLineWidth(2);
0120     // h1_dphi_distri_pr1_clustinnr_P->SetTitle("Comparison of e+ and e-");
0121     // h1_dphi_distri_pr1_clustinnr_P->Draw("");
0122
0123     // h1_dphi_distri_pr1_clustinnr_N->SetLineColor(kBlue);
0124     // h1_dphi_distri_pr1_clustinnr_N->SetLineWidth(2);
0125     // h1_dphi_distri_pr1_clustinnr_N->Draw("SAME");
0126
0127     // // 添加图例
0128     // TLegend *leg2 = new TLegend(0.65, 0.75, 0.88, 0.88);
0129     // leg2->AddEntry(h1_dphi_distri_pr1_clustinnr_P, "Positron", "l");
0130     // leg2->AddEntry(h1_dphi_distri_pr1_clustinnr_N, "Electron", "l");
0131     // leg2->Draw();
0132
0133     // 保存图像
0134     c2->SaveAs("comparePNcharge_dphi.pdf");
0135
0136     // 保存到新的文件
0137     TFile *output = new TFile("PNchargeComp.root", "RECREATE");
0138     c1->Write();
0139     c2->Write();
0140     TP_DR_pr1_clustinnr_pt_PN->Write();
0141     output->Close();
0142
0143     std::cout << "完成！新直方图已保存到 output.root 中。" << std::endl;
0144
0145     // 清理
0146     file1_Pcharge->Close();
0147     file2_Ncharge->Close();
0148   }
0149