LightFlavorRatios/lfratio_pythia/lfratio_pythia.cc

0001 #include "HepMC3/GenEvent.h"
0002 #include "HepMC3/GenParticle.h"
0003 #include "HepMC3/GenVertex.h"
0004 #include "HepMC3/ReaderAscii.h"
0005
0006 #include "TFile.h"
0007 #include "TH1.h"
0008 #include "TH1D.h"
0009 #include "TH2F.h"
0010 #include "TTree.h"
0011
0012 #include <cmath>
0013 #include <cstdlib>
0014 #include <iostream>
0015 #include <string>
0016 #include <vector>
0017
0018 #include "/sphenix/user/hjheng/sPHENIXRepo/analysis/LightFlavorRatios/util/binning.h"
0019
0020 // Return true if this particle is the last copy in a same-PID chain
0021 // This avoids double counting intermediate copies in HepMC records
0022 bool is_last_copy(const HepMC3::ConstGenParticlePtr &p)
0023 {
0024     if (!p)
0025         return false;
0026
0027     auto end_vtx = p->end_vertex();
0028     if (!end_vtx)
0029         return true;
0030
0031     for (const auto &d : end_vtx->particles_out())
0032     {
0033         if (!d)
0034             continue;
0035
0036         if (d->pid() == p->pid())
0037             return false;
0038     }
0039
0040     return true;
0041 }
0042
0043 bool has_charm_or_bottom_digit(const int apid)
0044 {
0045     int x = apid / 10;
0046
0047     for (int i = 0; i < 3; ++i)
0048     {
0049         const int digit = x % 10;
0050
0051         if (digit == 4 || digit == 5)
0052             return true;
0053
0054         x /= 10;
0055     }
0056
0057     return false;
0058 }
0059
0060 // trace back through same-PID, single-step copies to the earliest copy
0061 HepMC3::ConstGenParticlePtr trace_to_first_copy(HepMC3::ConstGenParticlePtr p)
0062 {
0063     while (p)
0064     {
0065         auto prod = p->production_vertex();
0066         if (!prod)
0067             break;
0068
0069         HepMC3::ConstGenParticlePtr same_pid_parent = nullptr;
0070         for (const auto &in : prod->particles_in())
0071         {
0072             if (in && in->pid() == p->pid())
0073             {
0074                 same_pid_parent = in;
0075                 break;
0076             }
0077         }
0078
0079         if (!same_pid_parent)
0080             break;
0081
0082         p = same_pid_parent;
0083     }
0084
0085     return p;
0086 }
0087
0088 // Prompt definition used here:
0089 //
0090 // Keep K0S / Lambda / anti-Lambda if it is not produced from
0091 // weak strange-hadron feed-down or charm/bottom feed-down
0092 //
0093 // For K0S, do NOT reject K0/K0bar parents, since the HepMC record may
0094 // represent neutral-kaon mixing as K0/K0bar -> K0S
0095 bool is_from_feeddown_parent(const HepMC3::ConstGenParticlePtr &p)
0096 {
0097     if (!p)
0098         return false;
0099
0100     auto first = trace_to_first_copy(p);
0101     if (!first)
0102         return false;
0103
0104     auto prod_vtx = first->production_vertex();
0105     if (!prod_vtx)
0106         return false;
0107
0108     for (const auto &parent : prod_vtx->particles_in())
0109     {
0110         if (!parent)
0111             continue;
0112
0113         const int apid = std::abs(parent->pid());
0114
0115         // Xi-, Xi0, Omega-
0116         if (apid == 3312 || apid == 3322 || apid == 3334)
0117             return true;
0118
0119         // charm/bottom hadron feed-down
0120         if (has_charm_or_bottom_digit(apid))
0121             return true;
0122     }
0123
0124     return false;
0125 }
0126
0127 int rapidity_region(const double y)
0128 {
0129     const double absy = std::abs(y);
0130
0131     if (absy <= 0.5)
0132         return 0;
0133
0134     if (absy <= 1.0)
0135         return 1;
0136
0137     return 2;
0138 }
0139
0140 int main(int argc, char *argv[])
0141 {
0142     if (argc != 3)
0143     {
0144         std::cerr << "Usage: " << argv[0] << " input.hepmc output.root\n";
0145         return 1;
0146     }
0147
0148     const char *hepmc_file = argv[1];
0149     const char *root_file = argv[2];
0150
0151     HepMC3::ReaderAscii reader(hepmc_file);
0152     HepMC3::GenEvent event;
0153
0154     TFile *fout = new TFile(root_file, "RECREATE");
0155
0156     if (!fout || fout->IsZombie())
0157     {
0158         std::cerr << "Error: cannot create output ROOT file: " << root_file << std::endl;
0159         return 1;
0160     }
0161
0162     const auto &pt_bins = BinInfo::final_pt_bins.bins;
0163     const auto &eta_bins = BinInfo::final_eta_bins.bins;
0164     const auto &rapidity_bins = BinInfo::final_rapidity_bins.bins;
0165     const auto &phi_bins = BinInfo::final_phi_bins.bins;
0166
0167     TH1D *h_event_counter = new TH1D("h_event_counter", ";dummy;Events", 1, 0.5, 1.5);
0168     TH1D *h_prompt_kshort_count = new TH1D("h_prompt_kshort_count", ";dummy;Prompt K^{0}_{S}", 1, 0.5, 1.5);
0169     TH1D *h_prompt_lambda_count = new TH1D("h_prompt_lambda_count", ";dummy;Prompt #Lambda^{0}(+c.c.)", 1, 0.5, 1.5);
0170
0171     TH1D *h_prompt_kshort_pt = new TH1D("h_prompt_kshort_pt", ";K^{0}_{S} p_{T} [GeV];Counts", pt_bins.size() - 1, pt_bins.data());
0172     TH1D *h_prompt_lambda_pt = new TH1D("h_prompt_lambda_pt", ";#Lambda^{0}(+c.c.) p_{T} [GeV];Counts", pt_bins.size() - 1, pt_bins.data());
0173
0174     TH1D *h_prompt_kshort_eta = new TH1D("h_prompt_kshort_eta", ";K^{0}_{S} #eta;Counts", eta_bins.size() - 1, eta_bins.data());
0175     TH1D *h_prompt_lambda_eta = new TH1D("h_prompt_lambda_eta", ";#Lambda^{0}(+c.c.) #eta;Counts", eta_bins.size() - 1, eta_bins.data());
0176
0177     TH1D *h_prompt_kshort_y = new TH1D("h_prompt_kshort_y", ";K^{0}_{S} rapidity;Counts", rapidity_bins.size() - 1, rapidity_bins.data());
0178     TH1D *h_prompt_lambda_y = new TH1D("h_prompt_lambda_y", ";#Lambda^{0}(+c.c.) rapidity;Counts", rapidity_bins.size() - 1, rapidity_bins.data());
0179
0180     TH1D *h_prompt_kshort_phi = new TH1D("h_prompt_kshort_phi", ";K^{0}_{S} #phi;Counts", phi_bins.size() - 1, phi_bins.data());
0181     TH1D *h_prompt_lambda_phi = new TH1D("h_prompt_lambda_phi", ";#Lambda^{0}(+c.c.) #phi;Counts", phi_bins.size() - 1, phi_bins.data());
0182
0183     TH1D *h_prompt_kshort_pt_y_le_0p5 = new TH1D("h_prompt_kshort_pt_y_le_0p5", ";K^{0}_{S} p_{T} [GeV];Counts, |y| #leq 0.5", pt_bins.size() - 1, pt_bins.data());
0184
0185     TH1D *h_prompt_lambda_pt_y_le_0p5 = new TH1D("h_prompt_lambda_pt_y_le_0p5", ";#Lambda^{0}(+c.c.) p_{T} [GeV];Counts, |y| #leq 0.5", pt_bins.size() - 1, pt_bins.data());
0186
0187     TH1D *h_prompt_kshort_pt_y_0p5_1p0 = new TH1D("h_prompt_kshort_pt_y_0p5_1p0", ";K^{0}_{S} p_{T} [GeV];Counts, 0.5 < |y| #leq 1.0", pt_bins.size() - 1, pt_bins.data());
0188
0189     TH1D *h_prompt_lambda_pt_y_0p5_1p0 = new TH1D("h_prompt_lambda_pt_y_0p5_1p0", ";#Lambda^{0}(+c.c.) p_{T} [GeV];Counts, 0.5 < |y| #leq 1.0", pt_bins.size() - 1, pt_bins.data());
0190
0191     TH1D *h_prompt_kshort_pt_y_gt_1p0 = new TH1D("h_prompt_kshort_pt_y_gt_1p0", ";K^{0}_{S} p_{T} [GeV];Counts, |y| > 1.0", pt_bins.size() - 1, pt_bins.data());
0192
0193     TH1D *h_prompt_lambda_pt_y_gt_1p0 = new TH1D("h_prompt_lambda_pt_y_gt_1p0", ";#Lambda^{0}(+c.c.) p_{T} [GeV];Counts, |y| > 1.0", pt_bins.size() - 1, pt_bins.data());
0194
0195     TTree *tree = new TTree("prompt_strange_tree", "prompt_strange_tree");
0196
0197     int event_id = -1;
0198     int particle_pid = 0;
0199     int particle_y_region = -1;
0200     double particle_pt = -999.;
0201     double particle_eta = -999.;
0202     double particle_y = -999.;
0203     double particle_phi = -999.;
0204     double particle_mass = -999.;
0205
0206     tree->Branch("event", &event_id, "event/I");
0207     tree->Branch("pid", &particle_pid, "pid/I");
0208     tree->Branch("y_region", &particle_y_region, "y_region/I");
0209     tree->Branch("pt", &particle_pt, "pt/D");
0210     tree->Branch("eta", &particle_eta, "eta/D");
0211     tree->Branch("y", &particle_y, "y/D");
0212     tree->Branch("phi", &particle_phi, "phi/D");
0213     tree->Branch("mass", &particle_mass, "mass/D");
0214
0215     int iev = 0;
0216
0217     while (!reader.failed())
0218     {
0219         reader.read_event(event);
0220
0221         if (reader.failed())
0222             break;
0223
0224         ++iev;
0225
0226         if (iev % 5000 == 0)
0227         {
0228             std::cout << "Processing event " << iev << std::endl;
0229         }
0230
0231         event_id = iev;
0232         h_event_counter->Fill(1.0);
0233
0234         for (const auto &p : event.particles())
0235         {
0236             if (!p)
0237                 continue;
0238
0239             const int pid = p->pid();
0240             const int apid = std::abs(pid);
0241
0242             // select K0S, Lambda, and anti-Lambda
0243             if (pid != 310 && apid != 3122)
0244                 continue;
0245
0246             // avoid double counting same-PID copies
0247             if (!is_last_copy(p))
0248                 continue;
0249
0250             if (apid == 3122 && is_from_feeddown_parent(p))
0251                 continue;
0252
0253             particle_pid = pid;
0254             particle_pt = p->momentum().pt();
0255             particle_eta = p->momentum().eta();
0256             particle_y = p->momentum().rap();
0257             particle_phi = p->momentum().phi();
0258             particle_mass = p->momentum().m();
0259             particle_y_region = rapidity_region(particle_y);
0260
0261             if (pid == 310)
0262             {
0263                 h_prompt_kshort_count->Fill(1.0);
0264                 h_prompt_kshort_pt->Fill(particle_pt);
0265                 h_prompt_kshort_eta->Fill(particle_eta);
0266                 h_prompt_kshort_y->Fill(particle_y);
0267                 h_prompt_kshort_phi->Fill(particle_phi);
0268
0269                 if (particle_y_region == 0)
0270                     h_prompt_kshort_pt_y_le_0p5->Fill(particle_pt);
0271                 else if (particle_y_region == 1)
0272                     h_prompt_kshort_pt_y_0p5_1p0->Fill(particle_pt);
0273                 else
0274                     h_prompt_kshort_pt_y_gt_1p0->Fill(particle_pt);
0275             }
0276             else if (apid == 3122)
0277             {
0278                 h_prompt_lambda_count->Fill(1.0);
0279                 h_prompt_lambda_pt->Fill(particle_pt);
0280                 h_prompt_lambda_eta->Fill(particle_eta);
0281                 h_prompt_lambda_y->Fill(particle_y);
0282                 h_prompt_lambda_phi->Fill(particle_phi);
0283
0284                 if (particle_y_region == 0)
0285                     h_prompt_lambda_pt_y_le_0p5->Fill(particle_pt);
0286                 else if (particle_y_region == 1)
0287                     h_prompt_lambda_pt_y_0p5_1p0->Fill(particle_pt);
0288                 else
0289                     h_prompt_lambda_pt_y_gt_1p0->Fill(particle_pt);
0290             }
0291
0292             tree->Fill();
0293         }
0294
0295         event.clear();
0296     }
0297
0298     TH1D *h_ratio_lambda_over_kshort_pt = static_cast<TH1D *>(h_prompt_lambda_pt->Clone("h_ratio_lambda_over_kshort_pt"));
0299     h_ratio_lambda_over_kshort_pt->SetTitle(";p_{T} [GeV];#Lambda^{0}(+c.c.) / K^{0}_{S}");
0300     h_ratio_lambda_over_kshort_pt->Divide(h_prompt_lambda_pt, h_prompt_kshort_pt, 1.0, 1.0);
0301
0302     TH1D *h_ratio_lambda_over_2kshort_pt = static_cast<TH1D *>(h_prompt_lambda_pt->Clone("h_ratio_lambda_over_2kshort_pt"));
0303     h_ratio_lambda_over_2kshort_pt->SetTitle(";p_{T} [GeV];#Lambda^{0}(+c.c.) / 2K^{0}_{S}");
0304     h_ratio_lambda_over_2kshort_pt->Divide(h_prompt_lambda_pt, h_prompt_kshort_pt, 1.0, 2.0);
0305
0306     TH1D *h_ratio_lambda_over_kshort_eta = static_cast<TH1D *>(h_prompt_lambda_eta->Clone("h_ratio_lambda_over_kshort_eta"));
0307     h_ratio_lambda_over_kshort_eta->SetTitle(";#eta;#Lambda^{0}(+c.c.) / K^{0}_{S}");
0308     h_ratio_lambda_over_kshort_eta->Divide(h_prompt_lambda_eta, h_prompt_kshort_eta, 1.0, 1.0);
0309
0310     TH1D *h_ratio_lambda_over_2kshort_eta = static_cast<TH1D *>(h_prompt_lambda_eta->Clone("h_ratio_lambda_over_2kshort_eta"));
0311     h_ratio_lambda_over_2kshort_eta->SetTitle(";#eta;#Lambda^{0}(+c.c.) / 2K^{0}_{S}");
0312     h_ratio_lambda_over_2kshort_eta->Divide(h_prompt_lambda_eta, h_prompt_kshort_eta, 1.0, 2.0);
0313
0314     TH1D *h_ratio_lambda_over_kshort_y = static_cast<TH1D *>(h_prompt_lambda_y->Clone("h_ratio_lambda_over_kshort_y"));
0315     h_ratio_lambda_over_kshort_y->SetTitle(";rapidity;#Lambda^{0}(+c.c.) / K^{0}_{S}");
0316     h_ratio_lambda_over_kshort_y->Divide(h_prompt_lambda_y, h_prompt_kshort_y, 1.0, 1.0);
0317
0318     TH1D *h_ratio_lambda_over_2kshort_y = static_cast<TH1D *>(h_prompt_lambda_y->Clone("h_ratio_lambda_over_2kshort_y"));
0319     h_ratio_lambda_over_2kshort_y->SetTitle(";rapidity;#Lambda^{0}(+c.c.) / 2K^{0}_{S}");
0320     h_ratio_lambda_over_2kshort_y->Divide(h_prompt_lambda_y, h_prompt_kshort_y, 1.0, 2.0);
0321
0322     TH1D *h_ratio_lambda_over_kshort_pt_y_le_0p5 = static_cast<TH1D *>(h_prompt_lambda_pt_y_le_0p5->Clone("h_ratio_lambda_over_kshort_pt_y_le_0p5"));
0323     h_ratio_lambda_over_kshort_pt_y_le_0p5->SetTitle(";p_{T} [GeV];#Lambda^{0}(+c.c.) / K^{0}_{S}, |y| #leq 0.5");
0324     h_ratio_lambda_over_kshort_pt_y_le_0p5->Divide(h_prompt_lambda_pt_y_le_0p5, h_prompt_kshort_pt_y_le_0p5, 1.0, 1.0);
0325
0326     TH1D *h_ratio_lambda_over_2kshort_pt_y_le_0p5 = static_cast<TH1D *>(h_prompt_lambda_pt_y_le_0p5->Clone("h_ratio_lambda_over_2kshort_pt_y_le_0p5"));
0327     h_ratio_lambda_over_2kshort_pt_y_le_0p5->SetTitle(";p_{T} [GeV];#Lambda^{0}(+c.c.) / 2K^{0}_{S}, |y| #leq 0.5");
0328     h_ratio_lambda_over_2kshort_pt_y_le_0p5->Divide(h_prompt_lambda_pt_y_le_0p5, h_prompt_kshort_pt_y_le_0p5, 1.0, 2.0);
0329
0330     TH1D *h_ratio_lambda_over_kshort_pt_y_0p5_1p0 = static_cast<TH1D *>(h_prompt_lambda_pt_y_0p5_1p0->Clone("h_ratio_lambda_over_kshort_pt_y_0p5_1p0"));
0331     h_ratio_lambda_over_kshort_pt_y_0p5_1p0->SetTitle(";p_{T} [GeV];#Lambda^{0}(+c.c.) / K^{0}_{S}, 0.5 < |y| #leq 1.0");
0332     h_ratio_lambda_over_kshort_pt_y_0p5_1p0->Divide(h_prompt_lambda_pt_y_0p5_1p0, h_prompt_kshort_pt_y_0p5_1p0, 1.0, 1.0);
0333
0334     TH1D *h_ratio_lambda_over_2kshort_pt_y_0p5_1p0 = static_cast<TH1D *>(h_prompt_lambda_pt_y_0p5_1p0->Clone("h_ratio_lambda_over_2kshort_pt_y_0p5_1p0"));
0335     h_ratio_lambda_over_2kshort_pt_y_0p5_1p0->SetTitle(";p_{T} [GeV];#Lambda^{0}(+c.c.) / 2K^{0}_{S}, 0.5 < |y| #leq 1.0");
0336     h_ratio_lambda_over_2kshort_pt_y_0p5_1p0->Divide(h_prompt_lambda_pt_y_0p5_1p0, h_prompt_kshort_pt_y_0p5_1p0, 1.0, 2.0);
0337
0338     TH1D *h_ratio_lambda_over_kshort_pt_y_gt_1p0 = static_cast<TH1D *>(h_prompt_lambda_pt_y_gt_1p0->Clone("h_ratio_lambda_over_kshort_pt_y_gt_1p0"));
0339     h_ratio_lambda_over_kshort_pt_y_gt_1p0->SetTitle(";p_{T} [GeV];#Lambda^{0}(+c.c.) / K^{0}_{S}, |y| > 1.0");
0340     h_ratio_lambda_over_kshort_pt_y_gt_1p0->Divide(h_prompt_lambda_pt_y_gt_1p0, h_prompt_kshort_pt_y_gt_1p0, 1.0, 1.0);
0341
0342     TH1D *h_ratio_lambda_over_2kshort_pt_y_gt_1p0 = static_cast<TH1D *>(h_prompt_lambda_pt_y_gt_1p0->Clone("h_ratio_lambda_over_2kshort_pt_y_gt_1p0"));
0343     h_ratio_lambda_over_2kshort_pt_y_gt_1p0->SetTitle(";p_{T} [GeV];#Lambda^{0}(+c.c.) / 2K^{0}_{S}, |y| > 1.0");
0344     h_ratio_lambda_over_2kshort_pt_y_gt_1p0->Divide(h_prompt_lambda_pt_y_gt_1p0, h_prompt_kshort_pt_y_gt_1p0, 1.0, 2.0);
0345
0346     const double n_prompt_kshort = h_prompt_kshort_count->GetBinContent(1);
0347     const double n_prompt_lambda = h_prompt_lambda_count->GetBinContent(1);
0348
0349     const double ratio_lambda_over_kshort = (n_prompt_kshort > 0.0) ? n_prompt_lambda / n_prompt_kshort : 0.0;
0350     const double ratio_lambda_over_2kshort = (n_prompt_kshort > 0.0) ? n_prompt_lambda / (2.0 * n_prompt_kshort) : 0.0;
0351
0352     std::cout << "Processed events              = " << iev << "\n";
0353     std::cout << "Prompt K0S                    = " << n_prompt_kshort << "\n";
0354     std::cout << "Prompt Lambda(+c.c.)          = " << n_prompt_lambda << "\n";
0355     std::cout << "Lambda(+c.c.) / K0S           = " << ratio_lambda_over_kshort << "\n";
0356     std::cout << "Lambda(+c.c.) / 2K0S          = " << ratio_lambda_over_2kshort << "\n";
0357     std::cout << "Saved ROOT file               = " << root_file << "\n";
0358
0359     fout->Write();
0360     fout->Close();
0361
0362     delete fout;
0363
0364     return 0;
0365 }