NewCode2024/KshortMomentum/valid.C

0001 int Nevt = 100;
0002
0003 void valid()
0004 {
0005     // Read in two ntuple files into Rdataframe
0006     // ROOT::EnableImplicitMT();
0007     // ROOT::RDataFrame df_ori("EventTree", "/sphenix/tg/tg01/hf/hjheng/ppg02/dst/Sim_HIJING_ana457_20250110/ntuple_00000.root");
0008     ROOT::RDataFrame df_ori("EventTree", "/sphenix/tg/tg01/hf/hjheng/ppg02/dst/Sim_HIJING_MDC2_ana467_20250225/ntuple_00000.root");
0009     ROOT::RDataFrame df_add40perc("EventTree", "/sphenix/tg/tg01/hf/hjheng/ppg02/dst/Sim_HIJING_addStrange40perc_ana457_20250113/ntuple_00000.root");
0010     // ROOT::RDataFrame df_add40perc("EventTree", "/sphenix/tg/tg01/hf/hjheng/ppg02/dst/Sim_HIJING_strangeness_MDC2_ana467_20250226/ntuple_00000.root");
0011     ROOT::RDataFrame df_add100perc("EventTree", "/sphenix/tg/tg01/hf/hjheng/ppg02/dst/Sim_HIJING_addStrange100perc_ana457_20250114/ntuple_00000.root");
0012
0013     std::vector<int> list_strangePID = {310, 3122, -3122}; // list of PID for strange particles to check
0014     // From the HepMCFSPrtl_PID column, count the number of strange particles (PID is in list_strangePID) and create a new column for an integer, Nstrange_HepMC, for the count
0015     // And from the PrimaryG4P_PID column, count the number of strange particles and add a new column, Nstrange_PHG4, for the count
0016     auto df_ori_Nstrange = df_ori
0017                                .Define("Nstrange_HepMC",
0018                                        [list_strangePID](const std::vector<int> &HepMCFSPrtl_PID)
0019                                        {
0020                                            // Count occurrences of PIDs in list_strangePID within HepMCFSPrtl_PID
0021                                            return static_cast<int>(std::count_if(HepMCFSPrtl_PID.begin(), HepMCFSPrtl_PID.end(), [&list_strangePID](int pid) { return std::find(list_strangePID.begin(), list_strangePID.end(), pid) != list_strangePID.end(); }));
0022                                        },
0023                                        {"HepMCFSPrtl_PID"})
0024                                .Define("Nstrange_PHG4",
0025                                        [list_strangePID](const std::vector<int> &PrimaryG4P_PID)
0026                                        {
0027                                            // Count occurrences of PIDs in list_strangePID within PrimaryG4P_PID
0028                                            return static_cast<int>(std::count_if(PrimaryG4P_PID.begin(), PrimaryG4P_PID.end(), [&list_strangePID](int pid) { return std::find(list_strangePID.begin(), list_strangePID.end(), pid) != list_strangePID.end(); }));
0029                                        },
0030                                        {"PrimaryG4P_PID"});
0031
0032     auto df_add40perc_Nstrange = df_add40perc
0033                                      .Define("Nstrange_HepMC",
0034                                              [list_strangePID](const std::vector<int> &HepMCFSPrtl_PID)
0035                                              {
0036                                                  // Count occurrences of PIDs in list_strangePID within HepMCFSPrtl_PID
0037                                                  return static_cast<int>(std::count_if(HepMCFSPrtl_PID.begin(), HepMCFSPrtl_PID.end(), [&list_strangePID](int pid) { return std::find(list_strangePID.begin(), list_strangePID.end(), pid) != list_strangePID.end(); }));
0038                                              },
0039                                              {"HepMCFSPrtl_PID"})
0040                                      .Define("Nstrange_PHG4",
0041                                              [list_strangePID](const std::vector<int> &PrimaryG4P_PID)
0042                                              {
0043                                                  // Count occurrences of PIDs in list_strangePID within PrimaryG4P_PID
0044                                                  return static_cast<int>(std::count_if(PrimaryG4P_PID.begin(), PrimaryG4P_PID.end(), [&list_strangePID](int pid) { return std::find(list_strangePID.begin(), list_strangePID.end(), pid) != list_strangePID.end(); }));
0045                                              },
0046                                              {"PrimaryG4P_PID"});
0047
0048     auto df_add100perc_Nstrange = df_add100perc
0049                                       .Define("Nstrange_HepMC",
0050                                               [list_strangePID](const std::vector<int> &HepMCFSPrtl_PID)
0051                                               {
0052                                                   // Count occurrences of PIDs in list_strangePID within HepMCFSPrtl_PID
0053                                                   return static_cast<int>(std::count_if(HepMCFSPrtl_PID.begin(), HepMCFSPrtl_PID.end(), [&list_strangePID](int pid) { return std::find(list_strangePID.begin(), list_strangePID.end(), pid) != list_strangePID.end(); }));
0054                                               },
0055                                               {"HepMCFSPrtl_PID"})
0056                                       .Define("Nstrange_PHG4",
0057                                               [list_strangePID](const std::vector<int> &PrimaryG4P_PID)
0058                                               {
0059                                                   // Count occurrences of PIDs in list_strangePID within PrimaryG4P_PID
0060                                                   return static_cast<int>(std::count_if(PrimaryG4P_PID.begin(), PrimaryG4P_PID.end(), [&list_strangePID](int pid) { return std::find(list_strangePID.begin(), list_strangePID.end(), pid) != list_strangePID.end(); }));
0061                                               },
0062                                               {"PrimaryG4P_PID"});
0063
0064     // Histograms of the PrimaryG4P_Pt and PrimaryG4P_Eta of the 40% and 100% additional strangeness samples
0065     auto h_ori_PHG4_Pt = df_ori.Histo1D({"h_ori_PHG4_Pt", ";PHG4Particle p_{T} [GeV]; Normalized entries", 100, 0, 10}, "PrimaryG4P_Pt");
0066     auto h_ori_PHG4_Eta = df_ori.Histo1D({"h_ori_PHG4_Eta", ";PHG4Particle #eta; Normalized entries", 40, -1, 1}, "PrimaryG4P_Eta");
0067     auto h_add40perc_PHG4_Pt = df_add40perc.Histo1D({"h_add40perc_PHG4_Pt", ";PHG4Particle p_{T} [GeV]; Normalized entries", 100, 0, 10}, "PrimaryG4P_Pt");
0068     auto h_add40perc_PHG4_Eta = df_add40perc.Histo1D({"h_add40perc_PHG4_Eta", ";PHG4Particle #eta; Normalized entries", 40, -1, 1}, "PrimaryG4P_Eta");
0069     auto h_add100perc_PHG4_Pt = df_add100perc.Histo1D({"h_add100perc_PHG4_Pt", ";PHG4Particle p_{T} [GeV]; Normalized entries", 100, 0, 10}, "PrimaryG4P_Pt");
0070     auto h_add100perc_PHG4_Eta = df_add100perc.Histo1D({"h_add100perc_PHG4_Eta", ";PHG4Particle #eta; Normalized entries", 40, -1, 1}, "PrimaryG4P_Eta");
0071
0072     // Make histogram/graphs for the number of strange particles in HepMC and PHG4 as a function of event number
0073     // First, get the Nstrange_HepMC/Nstrange_PHG4 columns as vectors
0074     auto Nstrange_HepMC_ori = df_ori_Nstrange.Take<int>("Nstrange_HepMC");
0075     auto Nstrange_PHG4_ori = df_ori_Nstrange.Take<int>("Nstrange_PHG4");
0076     auto Nstrange_HepMC_add40perc = df_add40perc_Nstrange.Take<int>("Nstrange_HepMC");
0077     auto Nstrange_PHG4_add40perc = df_add40perc_Nstrange.Take<int>("Nstrange_PHG4");
0078     auto Nstrange_HepMC_add100perc = df_add100perc_Nstrange.Take<int>("Nstrange_HepMC");
0079     auto Nstrange_PHG4_add100perc = df_add100perc_Nstrange.Take<int>("Nstrange_PHG4");
0080
0081     TH1F *h_Nstrange_HepMC_ori = new TH1F("h_Nstrange_HepMC_ori", ";Event index;N_{strange (K_{s}^{0},#Lambda)}", Nevt+1, -0.5, Nevt+0.5);
0082     TH1F *h_Nstrange_PHG4_ori = new TH1F("h_Nstrange_PHG4_ori", ";Event index;N_{strange (K_{s}^{0},#Lambda)}", Nevt+1, -0.5, Nevt+0.5);
0083     TH1F *h_Nstrange_HepMC_add40perc = new TH1F("h_Nstrange_HepMC_add40perc", ";Event index;N_{strange (K_{s}^{0},#Lambda)}", Nevt+1, -0.5, Nevt+0.5);
0084     TH1F *h_Nstrange_PHG4_add40perc = new TH1F("h_Nstrange_PHG4_add40perc", ";Event index;N_{strange (K_{s}^{0},#Lambda)}", Nevt+1, -0.5, Nevt+0.5);
0085     TH1F *h_Nstrange_HepMC_add100perc = new TH1F("h_Nstrange_HepMC_add100perc", ";Event index;N_{strange (K_{s}^{0},#Lambda)}", Nevt+1, -0.5, Nevt+0.5);
0086     TH1F *h_Nstrange_PHG4_add100perc = new TH1F("h_Nstrange_PHG4_add100perc", ";Event index;N_{strange (K_{s}^{0},#Lambda)}", Nevt+1, -0.5, Nevt+0.5);
0087
0088     for (int i = 0; i < Nevt; i++)
0089     {
0090         cout << "Event " << i << ": Nstrange HepMC (ori,40%,100%): " << (*Nstrange_HepMC_ori)[i] << ", " << (*Nstrange_HepMC_add40perc)[i] << ", " << (*Nstrange_HepMC_add100perc)[i] << "; Nstrange PHG4 (ori,40%,100%): " << (*Nstrange_PHG4_ori)[i] << ", " << (*Nstrange_PHG4_add40perc)[i] << ", " << (*Nstrange_PHG4_add100perc)[i] << endl;
0091         h_Nstrange_HepMC_ori->SetBinContent(i + 1, (*Nstrange_HepMC_ori)[i]);
0092         h_Nstrange_PHG4_ori->SetBinContent(i + 1, (*Nstrange_PHG4_ori)[i]);
0093         h_Nstrange_HepMC_add40perc->SetBinContent(i + 1, (*Nstrange_HepMC_add40perc)[i]);
0094         h_Nstrange_PHG4_add40perc->SetBinContent(i + 1, (*Nstrange_PHG4_add40perc)[i]);
0095         h_Nstrange_HepMC_add100perc->SetBinContent(i + 1, (*Nstrange_HepMC_add100perc)[i]);
0096         h_Nstrange_PHG4_add100perc->SetBinContent(i + 1, (*Nstrange_PHG4_add100perc)[i]);
0097     }
0098
0099     TH1F *h_NstrangeRatio_PHG4_add40perc = (TH1F*) h_Nstrange_PHG4_add40perc->Clone("h_NstrangeRatio_PHG4_add40perc");
0100     h_NstrangeRatio_PHG4_add40perc->Divide(h_Nstrange_HepMC_ori);
0101     TH1F *h_NstrangeRatio_PHG4_add100perc = (TH1F*) h_Nstrange_PHG4_add100perc->Clone("h_NstrangeRatio_PHG4_add100perc");
0102     h_NstrangeRatio_PHG4_add100perc->Divide(h_Nstrange_HepMC_ori);
0103
0104     // get the maximum value among all histograms
0105     double maxval = std::max({h_Nstrange_HepMC_ori->GetMaximum(), h_Nstrange_PHG4_ori->GetMaximum(), h_Nstrange_HepMC_add40perc->GetMaximum(), h_Nstrange_PHG4_add40perc->GetMaximum(), h_Nstrange_HepMC_add100perc->GetMaximum(), h_Nstrange_PHG4_add100perc->GetMaximum()});
0106
0107     TCanvas *c1 = new TCanvas("c1", "c1", 800, 700);
0108     gPad->SetTopMargin(0.08);
0109     c1->cd();
0110     h_Nstrange_HepMC_ori->SetLineColor(kBlack);
0111     h_Nstrange_PHG4_ori->SetLineColor(4);
0112     h_Nstrange_HepMC_add40perc->SetLineColor(7);
0113     h_Nstrange_PHG4_add40perc->SetLineColor(3);
0114     h_Nstrange_HepMC_add100perc->SetLineColor(5);
0115     h_Nstrange_PHG4_add100perc->SetLineColor(2);
0116     h_Nstrange_HepMC_ori->GetYaxis()->SetRangeUser(0, maxval * 1.8);
0117     h_Nstrange_HepMC_ori->GetYaxis()->SetTitleOffset(1.6);
0118     h_Nstrange_HepMC_ori->Draw("L");
0119     h_Nstrange_PHG4_ori->Draw("L same");
0120     h_Nstrange_HepMC_add40perc->Draw("L same");
0121     h_Nstrange_PHG4_add40perc->Draw("L same");
0122     h_Nstrange_HepMC_add100perc->Draw("L same");
0123     h_Nstrange_PHG4_add100perc->Draw("L same");
0124     TLegend *leg = new TLegend(0.2, 0.67, 0.4, 0.88);
0125     leg->AddEntry(h_Nstrange_HepMC_ori, "HepMC Particle - original", "l");
0126     leg->AddEntry(h_Nstrange_PHG4_ori, "PHG4Particle - original", "l");
0127     leg->AddEntry(h_Nstrange_HepMC_add40perc, "HepMC Particle - add 40% strange particles", "l");
0128     leg->AddEntry(h_Nstrange_PHG4_add40perc, "PHG4Particle - add 40% strange particles", "l");
0129     leg->AddEntry(h_Nstrange_HepMC_add100perc, "HepMC Particle - add 100% strange particles", "l");
0130     leg->AddEntry(h_Nstrange_PHG4_add100perc, "PHG4Particle - add 100% strange particles", "l");
0131     leg->SetTextSize(0.035);
0132     leg->SetBorderSize(0);
0133     leg->SetFillStyle(0);
0134     leg->Draw();
0135     c1->SaveAs("validation_NstrangeEvt.png");
0136     c1->SaveAs("validation_NstrangeEvt.pdf");
0137
0138     h_ori_PHG4_Pt->Scale(1.0 / h_ori_PHG4_Pt->Integral(-1, -1));
0139     h_add40perc_PHG4_Pt->Scale(1.0 / h_add40perc_PHG4_Pt->Integral(-1, -1));
0140     h_add100perc_PHG4_Pt->Scale(1.0 / h_add100perc_PHG4_Pt->Integral(-1, -1));
0141     float maxval_PHG4_Pt = std::max({h_ori_PHG4_Pt->GetMaximum(), h_add40perc_PHG4_Pt->GetMaximum(), h_add100perc_PHG4_Pt->GetMaximum()});
0142     float minval_PHG4_Pt = std::min({h_ori_PHG4_Pt->GetMinimum(0), h_add40perc_PHG4_Pt->GetMinimum(0), h_add100perc_PHG4_Pt->GetMinimum(0)});
0143     c1->Clear();
0144     c1->cd();
0145     c1->SetLogy(1);
0146     h_ori_PHG4_Pt->SetLineColor(kBlack);
0147     h_add40perc_PHG4_Pt->SetLineColor(4);
0148     h_add100perc_PHG4_Pt->SetLineColor(2);
0149     h_ori_PHG4_Pt->GetYaxis()->SetRangeUser(minval_PHG4_Pt*0.7, maxval_PHG4_Pt*1.5);
0150     h_ori_PHG4_Pt->Draw("hist");
0151     h_add40perc_PHG4_Pt->Draw("hist same");
0152     h_add100perc_PHG4_Pt->Draw("hist same");
0153     TLegend *leg1_pt = new TLegend(0.45, 0.73, 0.8, 0.88);
0154     leg1_pt->AddEntry(h_ori_PHG4_Pt.GetPtr(), "Original", "l");
0155     leg1_pt->AddEntry(h_add40perc_PHG4_Pt.GetPtr(), "Add 40% strange particles", "l");
0156     leg1_pt->AddEntry(h_add100perc_PHG4_Pt.GetPtr(), "Add 100% strange particles", "l");
0157     leg1_pt->SetTextSize(0.04);
0158     leg1_pt->SetBorderSize(0);
0159     leg1_pt->SetFillStyle(0);
0160     leg1_pt->Draw();
0161     c1->SaveAs("validation_PHG4_Pt.png");
0162     c1->SaveAs("validation_PHG4_Pt.pdf");
0163
0164     h_ori_PHG4_Eta->Scale(1.0 / h_ori_PHG4_Eta->Integral());
0165     h_add40perc_PHG4_Eta->Scale(1.0 / h_add40perc_PHG4_Eta->Integral());
0166     h_add100perc_PHG4_Eta->Scale(1.0 / h_add100perc_PHG4_Eta->Integral());
0167     float maxval_PHG4_Eta = std::max({h_ori_PHG4_Eta->GetMaximum(), h_add40perc_PHG4_Eta->GetMaximum(), h_add100perc_PHG4_Eta->GetMaximum()});
0168     c1->Clear();
0169     c1->cd();
0170     c1->SetLogy(0);
0171     h_ori_PHG4_Eta->SetLineColor(kBlack);
0172     h_add40perc_PHG4_Eta->SetLineColor(4);
0173     h_add100perc_PHG4_Eta->SetLineColor(2);
0174     h_ori_PHG4_Eta->GetYaxis()->SetRangeUser(0, maxval_PHG4_Eta*1.5);
0175     h_ori_PHG4_Eta->GetYaxis()->SetTitleOffset(1.6);
0176     h_ori_PHG4_Eta->Draw("hist");
0177     h_add40perc_PHG4_Eta->Draw("ep same");
0178     h_add100perc_PHG4_Eta->Draw("hist same");
0179     TLegend *leg1_eta = new TLegend(0.45, 0.73, 0.8, 0.88);
0180     leg1_eta->AddEntry(h_ori_PHG4_Eta.GetPtr(), "Original", "l");
0181     leg1_eta->AddEntry(h_add40perc_PHG4_Eta.GetPtr(), "Add 40% strange particles", "l");
0182     leg1_eta->AddEntry(h_add100perc_PHG4_Eta.GetPtr(), "Add 100% strange particles", "l");
0183     leg1_eta->SetTextSize(0.04);
0184     leg1_eta->SetBorderSize(0);
0185     leg1_eta->SetFillStyle(0);
0186     leg1_eta->Draw();
0187     c1->SaveAs("validation_PHG4_Eta.png");
0188     c1->SaveAs("validation_PHG4_Eta.pdf");
0189
0190     c1->Clear();
0191     c1->cd();
0192     h_NstrangeRatio_PHG4_add40perc->SetLineColor(3);
0193     h_NstrangeRatio_PHG4_add100perc->SetLineColor(2);
0194     h_NstrangeRatio_PHG4_add40perc->GetYaxis()->SetRangeUser(0, 5);
0195     h_NstrangeRatio_PHG4_add40perc->Draw("L");
0196     h_NstrangeRatio_PHG4_add100perc->Draw("L same");
0197     TLegend *leg2 = new TLegend(0.2, 0.7, 0.4, 0.9);
0198     leg2->AddEntry(h_NstrangeRatio_PHG4_add40perc, "PHG4Particle - add 40% strange particles", "l");
0199     leg2->AddEntry(h_NstrangeRatio_PHG4_add100perc, "PHG4Particle - add 100% strange particles", "l");
0200     leg2->SetTextSize(0.035);
0201     leg2->SetBorderSize(0);
0202     leg2->SetFillStyle(0);
0203     leg2->Draw();
0204     c1->SaveAs("validation_NstrangeRatioEvt.png");
0205     c1->SaveAs("validation_NstrangeRatioEvt.pdf");
0206 }