cascade_feeddown/Pythia/Xi_fraction.cc

0001 #include "HepMC3/ReaderAscii.h"
0002 #include "HepMC3/GenEvent.h"
0003 #include "HepMC3/GenParticle.h"
0004 #include "HepMC3/GenVertex.h"
0005 #include "HepMC3/Units.h"
0006
0007 #include "TFile.h"
0008 #include "TTree.h"
0009 #include "TH1D.h"
0010 #include "TH2F.h"
0011
0012 #include <cmath>
0013 #include <iostream>
0014 #include <memory>
0015 #include <string>
0016 #include <vector>
0017
0018 #include "/sphenix/user/hjheng/sPHENIXRepo/analysis/LightFlavorRatios/util/binning.h"
0019
0020 // Hardcoded MVTX global displacement (x,y,z)= (6.2955, -1.06874, -6.66289) mm, and the minimum radius of MVTX from simulation/g4simulation/g4mvtx/PHG4MvtxDefs.h
0021 const double MVTX_global_displacement_x_cm = 6.2955 / 10.0;  // convert mm to cm
0022 const double MVTX_global_displacement_y_cm = -1.06874 / 10.0;
0023 const double MVTX_global_displacement_z_cm = -6.66289 / 10.0;
0024 const double MVTX_rmin_cm = 24.610 / 10.0; // convert mm to cm, innermost layer
0025
0026 // Return true if this particle is the last copy in a same-PID chain
0027 // This helps avoid double counting intermediate copies in HepMC records
0028 bool is_last_copy(const HepMC3::ConstGenParticlePtr &p)
0029 {
0030     if (!p)
0031         return false;
0032
0033     auto end_vtx = p->end_vertex();
0034     if (!end_vtx)
0035         return true;
0036
0037     for (const auto &d : end_vtx->particles_out())
0038     {
0039         if (!d)
0040             continue;
0041         if (d->pid() == p->pid())
0042             return false;
0043     }
0044     return true;
0045 }
0046
0047 // Convert HepMC length-unit value to cm
0048 double length_unit_to_cm(const HepMC3::GenEvent &event)
0049 {
0050     const auto u = event.length_unit();
0051
0052     if (u == HepMC3::Units::CM)
0053         return 1.0;
0054     if (u == HepMC3::Units::MM)
0055         return 0.1;
0056
0057     return 1.0;
0058 }
0059
0060 // decay length in cm, using production and decay vertices of the particle, option for transverse or 3D
0061 double decay_length_cm(const HepMC3::ConstGenParticlePtr &p, const HepMC3::GenEvent &event, bool transverse = true)
0062 {
0063     if (!p)
0064         return -999.0;
0065
0066     auto prod_vtx = p->production_vertex();
0067     auto end_vtx = p->end_vertex();
0068
0069     if (!prod_vtx || !end_vtx)
0070         return -999.0;
0071
0072     const auto &prod = prod_vtx->position();
0073     const auto &dec = end_vtx->position();
0074
0075     const double dx = dec.x() - prod.x();
0076     const double dy = dec.y() - prod.y();
0077     const double dz = dec.z() - prod.z();
0078
0079     if (transverse)
0080     {
0081         const double L = std::sqrt(dx * dx + dy * dy);
0082         return L * length_unit_to_cm(event);
0083     }
0084     else
0085     {
0086         const double L = std::sqrt(dx * dx + dy * dy + dz * dz);
0087         return L * length_unit_to_cm(event);
0088     }
0089 }
0090
0091 // Distance from the origin to the displaced inner MVTX boundary along azimuth phi.
0092 double mvtx_inner_radius_at_phi_cm(const double phi)
0093 {
0094     // convert this phi in HEP definition to [0, 2pi]
0095     double phi_local = phi;
0096     if (phi_local < 0)
0097     {
0098         phi_local += 2 * M_PI;
0099     }
0100
0101     const double cos_phi = std::cos(phi_local);
0102     const double sin_phi = std::sin(phi_local);
0103
0104     // calculate the angle for the displacement vector, make it in the range of [0, 2pi] as well
0105     double disp_phi = std::atan2(MVTX_global_displacement_y_cm, MVTX_global_displacement_x_cm);
0106     if (disp_phi < 0)
0107     {
0108         disp_phi += 2 * M_PI;
0109     }
0110     double d = std::sqrt(MVTX_global_displacement_x_cm * MVTX_global_displacement_x_cm + MVTX_global_displacement_y_cm * MVTX_global_displacement_y_cm);
0111     double rho = d * std::cos(phi_local - disp_phi) + std::sqrt(std::pow(MVTX_rmin_cm, 2) - std::pow(d * std::sin(phi_local - disp_phi), 2));
0112
0113     return rho;
0114 }
0115
0116 // Check whether this Lambda comes directly from a Xi parent at its production vertex
0117 // 3312: from charged Xi (+/-3312)
0118 // 3322: from neutral Xi (+/-3322)
0119 // 0: otherwise
0120 int lambda_parent_xi_type(const HepMC3::ConstGenParticlePtr &lambda)
0121 {
0122     if (!lambda)
0123         return 0;
0124
0125     auto prod_vtx = lambda->production_vertex();
0126     if (!prod_vtx)
0127         return 0;
0128
0129     for (const auto &parent : prod_vtx->particles_in())
0130     {
0131         if (!parent)
0132             continue;
0133
0134         const int apid = std::abs(parent->pid());
0135
0136         if (apid == 3312)
0137             return 3312; // charged Xi
0138         if (apid == 3322)
0139             return 3322; // neutral Xi
0140     }
0141
0142     return 0;
0143 }
0144
0145 int main(int argc, char *argv[])
0146 {
0147     if (argc != 3)
0148     {
0149         std::cerr << "Usage: " << argv[0] << " input.hepmc output.root\n";
0150         return 1;
0151     }
0152
0153     const char *hepmc_file = argv[1];
0154     const char *root_file = argv[2];
0155
0156     HepMC3::ReaderAscii reader(hepmc_file);
0157     HepMC3::GenEvent event;
0158
0159     TFile *fout = new TFile(root_file, "RECREATE");
0160
0161     // Tree for selected Lambdas
0162     TTree *tree = new TTree("lambda_tree", "lambda_tree");
0163
0164     int event_id = -1;
0165     int lambda_pid = 0;
0166     int parent_xi_type = 0; // 0, 3312, 3322
0167     double lambda_pt = -999.;
0168     double lambda_eta = -999.;
0169     double lambda_y = -999.;
0170     double lambda_phi = -999.;
0171     double lambda_mass = -999.;
0172     double lambda_transverse_decay_length_cm = -999.;
0173
0174     tree->Branch("event", &event_id, "event/I");
0175     tree->Branch("lambda_pid", &lambda_pid, "lambda_pid/I");
0176     tree->Branch("parent_xi_type", &parent_xi_type, "parent_xi_type/I");
0177     tree->Branch("lambda_pt", &lambda_pt, "lambda_pt/D");
0178     tree->Branch("lambda_eta", &lambda_eta, "lambda_eta/D");
0179     tree->Branch("lambda_y", &lambda_y, "lambda_y/D");
0180     tree->Branch("lambda_mass", &lambda_mass, "lambda_mass/D");
0181     tree->Branch("lambda_transverse_decay_length_cm", &lambda_transverse_decay_length_cm, "lambda_transverse_decay_length_cm/D");
0182
0183     const auto &lambda_pt_bins = BinInfo::final_pt_bins.bins;
0184     const auto &lambda_rapidity_bins = BinInfo::final_rapidity_bins.bins;
0185     const auto &lambda_eta_bins = BinInfo::final_eta_bins.bins;
0186     const auto &lambda_phi_bins = BinInfo::final_phi_bins.bins;
0187
0188     // Denominator: all selected Lambdas
0189     TH1D *h_lambda_pt_all = new TH1D("h_lambda_pt_all", ";#Lambda^{0}(+c.c) p_{T} [GeV];Counts", lambda_pt_bins.size() - 1, lambda_pt_bins.data());
0190     TH1D *h_lambda_rapidity_all = new TH1D("h_lambda_rapidity_all", ";#Lambda^{0}(+c.c) rapidity;Counts", lambda_rapidity_bins.size() - 1, lambda_rapidity_bins.data());
0191     TH1D *h_lambda_eta_all = new TH1D("h_lambda_eta_all", ";#Lambda^{0}(+c.c) #eta;Counts", lambda_eta_bins.size() - 1, lambda_eta_bins.data());
0192     TH1D *h_lambda_phi_all = new TH1D("h_lambda_phi_all", ";#Lambda^{0}(+c.c) #phi;Counts", lambda_phi_bins.size() - 1, lambda_phi_bins.data());
0193     // Numerators: selected Lambdas from Xi
0194     TH1D *h_lambda_pt_from_xi_all = new TH1D("h_lambda_pt_from_xi_all", ";#Lambda^{0}(+c.c) p_{T} [GeV];Counts", lambda_pt_bins.size() - 1, lambda_pt_bins.data());
0195     TH1D *h_lambda_pt_from_xi_charged = new TH1D("h_lambda_pt_from_xi_charged", ";#Lambda^{0}(+c.c) p_{T} [GeV];Counts", lambda_pt_bins.size() - 1, lambda_pt_bins.data());
0196     TH1D *h_lambda_pt_from_xi_neutral = new TH1D("h_lambda_pt_from_xi_neutral", ";#Lambda^{0}(+c.c) p_{T} [GeV];Counts", lambda_pt_bins.size() - 1, lambda_pt_bins.data());
0197     TH1D *h_lambda_rapidity_from_xi_all = new TH1D("h_lambda_rapidity_from_xi_all", ";#Lambda^{0}(+c.c) rapidity;Counts", lambda_rapidity_bins.size() - 1, lambda_rapidity_bins.data());
0198     TH1D *h_lambda_rapidity_from_xi_charged = new TH1D("h_lambda_rapidity_from_xi_charged", ";#Lambda^{0}(+c.c) rapidity;Counts", lambda_rapidity_bins.size() - 1, lambda_rapidity_bins.data());
0199     TH1D *h_lambda_rapidity_from_xi_neutral = new TH1D("h_lambda_rapidity_from_xi_neutral", ";#Lambda^{0}(+c.c) rapidity;Counts", lambda_rapidity_bins.size() - 1, lambda_rapidity_bins.data());
0200     TH1D *h_lambda_eta_from_xi_all = new TH1D("h_lambda_eta_from_xi_all", ";#Lambda^{0}(+c.c) #eta;Counts", lambda_eta_bins.size() - 1, lambda_eta_bins.data());
0201     TH1D *h_lambda_eta_from_xi_charged = new TH1D("h_lambda_eta_from_xi_charged", ";#Lambda^{0}(+c.c) #eta;Counts", lambda_eta_bins.size() - 1, lambda_eta_bins.data());
0202     TH1D *h_lambda_eta_from_xi_neutral = new TH1D("h_lambda_eta_from_xi_neutral", ";#Lambda^{0}(+c.c) #eta;Counts", lambda_eta_bins.size() - 1, lambda_eta_bins.data());
0203     TH1D *h_lambda_phi_from_xi_all = new TH1D("h_lambda_phi_from_xi_all", ";#Lambda^{0}(+c.c) #phi;Counts", lambda_phi_bins.size() - 1, lambda_phi_bins.data());
0204     TH1D *h_lambda_phi_from_xi_charged = new TH1D("h_lambda_phi_from_xi_charged", ";#Lambda^{0}(+c.c) #phi;Counts", lambda_phi_bins.size() - 1, lambda_phi_bins.data());
0205     TH1D *h_lambda_phi_from_xi_neutral = new TH1D("h_lambda_phi_from_xi_neutral", ";#Lambda^{0}(+c.c) #phi;Counts", lambda_phi_bins.size() - 1, lambda_phi_bins.data());
0206     // Event counter (for check)
0207     TH1D *h_event_counter = new TH1D("h_event_counter", ";dummy;events", 1, 0.5, 1.5);
0208     // Diagnostic histograms
0209     TH1D *h_lambda_transverse_decay_length_cm = new TH1D("h_lambda_transverse_decay_length_cm", ";#Lambda transverse decay length [cm];Counts", 100, 0, 5);
0210     TH2D *h_phi_decaylengthcut = new TH2D("h_phi_decaylengthcut", ";#Lambda #phi [rad];Transverse decay length cut [cm]", 128, -3.2, 3.2, 100, 0, 4);
0211
0212     int iev = 0;
0213
0214     while (!reader.failed())
0215     {
0216         reader.read_event(event);
0217         if (reader.failed())
0218             break;
0219
0220         ++iev;
0221
0222         if (iev % 5000 == 0)
0223         {
0224             std::cout << "Processing event " << iev << std::endl;
0225         }
0226
0227         event_id = iev;
0228         h_event_counter->Fill(1.0);
0229
0230         for (const auto &p : event.particles())
0231         {
0232             if (!p)
0233                 continue;
0234
0235             // Select Lambda0 / anti-Lambda0
0236             if (std::abs(p->pid()) != 3122)
0237                 continue;
0238
0239             // Avoid double counting if same-PID copies exist in a chain
0240             if (!is_last_copy(p))
0241                 continue;
0242
0243             // Need both production and decay vertices for the decay-length cut
0244             if (!p->production_vertex() || !p->end_vertex())
0245                 continue;
0246
0247             lambda_pid = p->pid();
0248             lambda_pt = p->momentum().pt();
0249             lambda_eta = p->momentum().eta();
0250             lambda_y = p->momentum().rap();
0251             lambda_phi = p->momentum().phi();
0252             lambda_mass = p->momentum().m();
0253             lambda_transverse_decay_length_cm = decay_length_cm(p, event, true); // transverse decay length
0254             h_lambda_transverse_decay_length_cm->Fill(lambda_transverse_decay_length_cm);
0255
0256             // Kinematic selection on Lambda
0257             if (std::abs(lambda_eta) >= 1.3)
0258                 continue;
0259
0260             // Topological selection on transverse decay length: require the decay to occur before the displaced inner MVTX boundary at this azimuth.
0261             const double mvtx_transverse_cut_cm = mvtx_inner_radius_at_phi_cm(lambda_phi);
0262             h_phi_decaylengthcut->Fill(lambda_phi, mvtx_transverse_cut_cm);
0263             if (mvtx_transverse_cut_cm < 0 ||
0264                 lambda_transverse_decay_length_cm > mvtx_transverse_cut_cm ||
0265                 lambda_transverse_decay_length_cm < 0)
0266                 continue;
0267
0268             // Denominator
0269             h_lambda_pt_all->Fill(lambda_pt);
0270             h_lambda_rapidity_all->Fill(lambda_y);
0271             h_lambda_eta_all->Fill(lambda_eta);
0272             h_lambda_phi_all->Fill(lambda_phi);
0273
0274             // Parent classification
0275             parent_xi_type = lambda_parent_xi_type(p);
0276
0277             if (parent_xi_type == 3312)
0278             {
0279                 h_lambda_pt_from_xi_all->Fill(lambda_pt);
0280                 h_lambda_pt_from_xi_charged->Fill(lambda_pt);
0281                 h_lambda_rapidity_from_xi_all->Fill(lambda_y);
0282                 h_lambda_rapidity_from_xi_charged->Fill(lambda_y);
0283                 h_lambda_eta_from_xi_all->Fill(lambda_eta);
0284                 h_lambda_eta_from_xi_charged->Fill(lambda_eta);
0285                 h_lambda_phi_from_xi_all->Fill(lambda_phi);
0286                 h_lambda_phi_from_xi_charged->Fill(lambda_phi);
0287             }
0288             else if (parent_xi_type == 3322)
0289             {
0290                 h_lambda_pt_from_xi_all->Fill(lambda_pt);
0291                 h_lambda_pt_from_xi_neutral->Fill(lambda_pt);
0292                 h_lambda_rapidity_from_xi_all->Fill(lambda_y);
0293                 h_lambda_rapidity_from_xi_neutral->Fill(lambda_y);
0294                 h_lambda_eta_from_xi_all->Fill(lambda_eta);
0295                 h_lambda_eta_from_xi_neutral->Fill(lambda_eta);
0296                 h_lambda_phi_from_xi_all->Fill(lambda_phi);
0297                 h_lambda_phi_from_xi_neutral->Fill(lambda_phi);
0298             }
0299
0300             tree->Fill();
0301         }
0302
0303         event.clear();
0304     }
0305
0306     // Feeddown fractions
0307     TH1D *h_feeddown_frac_xi_all = (TH1D *)h_lambda_pt_from_xi_all->Clone("h_feeddown_frac_xi_all");
0308     h_feeddown_frac_xi_all->SetTitle(";#Lambda^{0}(+c.c) p_{T} [GeV];#Lambda from Xi / all #Lambda");
0309
0310     TH1D *h_feeddown_frac_xi_charged = (TH1D *)h_lambda_pt_from_xi_charged->Clone("h_feeddown_frac_xi_charged");
0311     h_feeddown_frac_xi_charged->SetTitle(";#Lambda^{0}(+c.c) p_{T} [GeV];#Lambda from charged Xi / all #Lambda");
0312
0313     TH1D *h_feeddown_frac_xi_neutral = (TH1D *)h_lambda_pt_from_xi_neutral->Clone("h_feeddown_frac_xi_neutral");
0314     h_feeddown_frac_xi_neutral->SetTitle(";#Lambda^{0}(+c.c) p_{T} [GeV];#Lambda from neutral Xi / all #Lambda");
0315
0316     // Binomial errors are appropriate since numerator is a subset of denominator.
0317     h_feeddown_frac_xi_all->Divide(h_lambda_pt_from_xi_all, h_lambda_pt_all, 1.0, 1.0, "B");
0318     h_feeddown_frac_xi_charged->Divide(h_lambda_pt_from_xi_charged, h_lambda_pt_all, 1.0, 1.0, "B");
0319     h_feeddown_frac_xi_neutral->Divide(h_lambda_pt_from_xi_neutral, h_lambda_pt_all, 1.0, 1.0, "B");
0320
0321     fout->Write();
0322     fout->Close();
0323
0324     std::cout << "Processed events = " << iev << "\n";
0325     std::cout << "Saved ROOT file  = " << root_file << "\n";
0326
0327     return 0;
0328 }