zdc_lambda/analysis/gen_lambda_decay.cxx

0001 #include "HepMC3/GenEvent.h"
0002 #include "HepMC3/ReaderAscii.h"
0003 #include "HepMC3/WriterAscii.h"
0004 #include "HepMC3/Print.h"
0005
0006 #include "TRandom3.h"
0007 #include "TVector3.h"
0008
0009 #include <TDatabasePDG.h>
0010 #include <TParticlePDG.h>
0011
0012 #include <iostream>
0013 #include <random>
0014 #include <TMath.h>
0015
0016 using namespace HepMC3;
0017
0018 std::tuple<double, int, double> GetParticleInfo(TDatabasePDG* pdg, TString particle_name)
0019 {
0020   TParticlePDG *particle = pdg->GetParticle(particle_name);
0021   const double mass = particle->Mass();
0022   const int pdgID = particle->PdgCode();
0023   const double lifetime = particle->Lifetime();
0024   return std::make_tuple(mass, pdgID, lifetime);
0025 }
0026 // Calculates the decay length of a particle. Samples from an exponential decay.
0027 double GetDecayLength(TRandom3* r1, double lifetime, double mass, double momentum_magnitude)
0028 {
0029   double c_speed = TMath::C() * 1000.; // speed of light im mm/sec
0030   double average_decay_length = (momentum_magnitude/mass) * lifetime * c_speed;
0031   return r1->Exp(average_decay_length);
0032 }
0033
0034 // Generate single lambda mesons and decay them to a neutron + 2 photons
0035 void gen_lambda_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname = "lambda_decay.hepmc",
0036               double p_min = 100., // in GeV/c
0037               double p_max = 275.) // in GeV/c
0038 {
0039
0040   const double theta_min = 0.0; // in mRad
0041   const double theta_max = 3.0; // in mRad
0042   //const double p_min = 100.; // in GeV/c
0043   //const double p_max = 275.; // in GeV/c
0044
0045   WriterAscii hepmc_output(out_fname);
0046   int events_parsed = 0;
0047   GenEvent evt(Units::GEV, Units::MM);
0048
0049   // Random number generator
0050   TRandom3 *r1 = new TRandom3(seed); //Default = 0, which uses clock to set seed
0051   cout<<"Random number seed is "<<r1->GetSeed()<<"!"<<endl;
0052
0053   // Getting generated particle information
0054   TDatabasePDG *pdg = new TDatabasePDG();
0055
0056   auto lambda_info = GetParticleInfo(pdg, "Lambda0");
0057   double lambda_mass = std::get<0>(lambda_info);
0058   int lambda_pdgID = std::get<1>(lambda_info);
0059   double lambda_lifetime = std::get<2>(lambda_info);
0060
0061   auto neutron_info = GetParticleInfo(pdg, "neutron");
0062   double neutron_mass = std::get<0>(neutron_info);
0063   int neutron_pdgID = std::get<1>(neutron_info);
0064
0065   auto pi0_info = GetParticleInfo(pdg, "pi0");
0066   double pi0_mass = std::get<0>(pi0_info);
0067   int pi0_pdgID = std::get<1>(pi0_info);
0068   double pi0_lifetime = std::get<2>(pi0_info);
0069
0070   auto photon_info = GetParticleInfo(pdg, "gamma");
0071   double photon_mass = std::get<0>(photon_info);
0072   int photon_pdgID = std::get<1>(photon_info);
0073
0074   for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
0075
0076     //Set the event number
0077     evt.set_event_number(events_parsed);
0078
0079     // FourVector(px,py,pz,e,pdgid,status)
0080     // type 4 is beam
0081     // pdgid 11 - electron
0082     // pdgid 2212 - proton
0083     GenParticlePtr p1 =
0084         std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
0085     GenParticlePtr p2 = std::make_shared<GenParticle>(
0086         FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
0087
0088     // Define momentum with respect to EIC proton beam direction
0089     Double_t lambda_p     = r1->Uniform(p_min, p_max);
0090     Double_t lambda_phi   = r1->Uniform(0.0, 2.0 * M_PI);
0091     Double_t lambda_th    = r1->Uniform(theta_min/1000., theta_max/1000.); // Divide by 1000 for radians
0092     Double_t lambda_px    = lambda_p * TMath::Cos(lambda_phi) * TMath::Sin(lambda_th);
0093     Double_t lambda_py    = lambda_p * TMath::Sin(lambda_phi) * TMath::Sin(lambda_th);
0094     Double_t lambda_pz    = lambda_p * TMath::Cos(lambda_th);
0095     Double_t lambda_E     = TMath::Sqrt(lambda_p*lambda_p + lambda_mass*lambda_mass);
0096
0097     // Rotate to lab coordinate system
0098     TVector3 lambda_pvec(lambda_px, lambda_py, lambda_pz);
0099     double cross_angle = -25./1000.; // in Rad
0100     TVector3 pbeam_dir(TMath::Sin(cross_angle), 0, TMath::Cos(cross_angle)); //proton beam direction
0101     lambda_pvec.RotateY(cross_angle); // Theta is returned positive, beam in negative X
0102
0103     // type 2 is state that will decay
0104     GenParticlePtr p_lambda = std::make_shared<GenParticle>(
0105         FourVector(lambda_pvec.X(), lambda_pvec.Y(), lambda_pvec.Z(), lambda_E), lambda_pdgID, 2 );
0106
0107     // Generating lambda particle, will be generated at origin
0108     // Must have input electron + proton for vertex
0109     GenVertexPtr lambda_initial_vertex = std::make_shared<GenVertex>();
0110     lambda_initial_vertex->add_particle_in(p1);
0111     lambda_initial_vertex->add_particle_in(p2);
0112     lambda_initial_vertex->add_particle_out(p_lambda);
0113     evt.add_vertex(lambda_initial_vertex);
0114
0115     // Generate neutron + pi0 in lambda rest frame
0116     TLorentzVector neutron_rest, pi0_rest;
0117
0118     // Generating uniformly along a sphere
0119     double cost_neutron_rest = r1->Uniform(-1,1);
0120     double th_neutron_rest = TMath::ACos(cost_neutron_rest);
0121     double sint_neutron_rest = TMath::Sin(th_neutron_rest);
0122
0123     double phi_neutron_rest = r1->Uniform(-1.*TMath::Pi(),1.*TMath::Pi());
0124     double cosp_neutron_rest = TMath::Cos(phi_neutron_rest);
0125     double sinp_neutron_rest = TMath::Sin(phi_neutron_rest);
0126
0127     // Calculate energy of each particle in the lambda rest frame
0128     // See problem 3.19 in Introduction to Elementary Particles, 2nd edition by D. Griffiths
0129     double E_neutron_rest = (-TMath::Power(pi0_mass, 2.) + TMath::Power(lambda_mass, 2.) + TMath::Power(neutron_mass, 2.) ) / (2. * lambda_mass) ;
0130     double E_pi0_rest = (-TMath::Power(neutron_mass, 2.) + TMath::Power(lambda_mass, 2.) + TMath::Power(pi0_mass, 2.) ) / (2. * lambda_mass) ;
0131
0132     // Both particles will have the same momentum, so just use neutron variables
0133     double momentum_rest = TMath::Sqrt( E_neutron_rest*E_neutron_rest - neutron_mass*neutron_mass );
0134
0135     neutron_rest.SetE(E_neutron_rest);
0136     neutron_rest.SetPx( momentum_rest * sint_neutron_rest * cosp_neutron_rest );
0137     neutron_rest.SetPy( momentum_rest * sint_neutron_rest * sinp_neutron_rest );
0138     neutron_rest.SetPz( momentum_rest * cost_neutron_rest );
0139
0140     pi0_rest.SetE(E_pi0_rest);
0141     pi0_rest.SetPx( -neutron_rest.Px() );
0142     pi0_rest.SetPy( -neutron_rest.Py() );
0143     pi0_rest.SetPz( -neutron_rest.Pz() );
0144
0145     // Boost neutron & pion to lab frame
0146     TLorentzVector lambda_lab(lambda_pvec.X(), lambda_pvec.Y(), lambda_pvec.Z(), lambda_E);
0147     TVector3 lambda_boost = lambda_lab.BoostVector();
0148     TLorentzVector neutron_lab, pi0_lab;
0149     neutron_lab = neutron_rest;
0150     neutron_lab.Boost(lambda_boost);
0151     pi0_lab = pi0_rest;
0152     pi0_lab.Boost(lambda_boost);
0153
0154     // Calculating position for lambda decay
0155     TVector3 lambda_unit = lambda_lab.Vect().Unit();
0156     double lambda_decay_length = GetDecayLength(r1, lambda_lifetime, lambda_mass, lambda_lab.P());
0157     TVector3 lambda_decay_position = lambda_unit * lambda_decay_length;
0158     double lambda_decay_time = lambda_decay_length / lambda_lab.Beta() ; // Decay time in lab frame in length units (mm)
0159
0160     // Generating vertex for lambda decay
0161     GenParticlePtr p_neutron = std::make_shared<GenParticle>(
0162       FourVector(neutron_lab.Px(), neutron_lab.Py(), neutron_lab.Pz(), neutron_lab.E()), neutron_pdgID, 1 );
0163
0164     GenParticlePtr p_pi0 = std::make_shared<GenParticle>(
0165       FourVector(pi0_lab.Px(), pi0_lab.Py(), pi0_lab.Pz(), pi0_lab.E()), pi0_pdgID, 2 );
0166
0167     GenVertexPtr v_lambda_decay = std::make_shared<GenVertex>(FourVector(lambda_decay_position.X(), lambda_decay_position.Y(), lambda_decay_position.Z(), lambda_decay_time));
0168     v_lambda_decay->add_particle_in(p_lambda);
0169     v_lambda_decay->add_particle_out(p_neutron);
0170     v_lambda_decay->add_particle_out(p_pi0);
0171
0172     evt.add_vertex(v_lambda_decay);
0173
0174     // Generate two photons from pi0 decay
0175     TLorentzVector gamma1_rest, gamma2_rest;
0176
0177     // Generating uniformly along a sphere
0178     double cost_gamma1_rest = r1->Uniform(-1,1);
0179     double th_gamma1_rest = TMath::ACos(cost_gamma1_rest);
0180     double sint_gamma1_rest = TMath::Sin(th_gamma1_rest);
0181
0182     double phi_gamma1_rest = r1->Uniform(-1.*TMath::Pi(),1.*TMath::Pi());
0183     double cosp_gamma1_rest = TMath::Cos(phi_gamma1_rest);
0184     double sinp_gamma1_rest = TMath::Sin(phi_gamma1_rest);
0185
0186     // Photons are massless so they each get equal energies
0187     gamma1_rest.SetE(pi0_mass/2.);
0188     gamma1_rest.SetPx( (pi0_mass/2.)*sint_gamma1_rest*cosp_gamma1_rest );
0189     gamma1_rest.SetPy( (pi0_mass/2.)*sint_gamma1_rest*sinp_gamma1_rest );
0190     gamma1_rest.SetPz( (pi0_mass/2.)*cost_gamma1_rest );
0191
0192     gamma2_rest.SetE(pi0_mass/2.);
0193     gamma2_rest.SetPx( -gamma1_rest.Px() );
0194     gamma2_rest.SetPy( -gamma1_rest.Py() );
0195     gamma2_rest.SetPz( -gamma1_rest.Pz() );
0196
0197     // Boost neutron & pion to lab frame
0198     TVector3 pi0_boost = pi0_lab.BoostVector();
0199     TLorentzVector gamma1_lab, gamma2_lab;
0200     gamma1_lab = gamma1_rest;
0201     gamma1_lab.Boost(pi0_boost);
0202     gamma2_lab = gamma2_rest;
0203     gamma2_lab.Boost(pi0_boost);
0204
0205     GenParticlePtr p_gamma1 = std::make_shared<GenParticle>(
0206       FourVector(gamma1_lab.Px(), gamma1_lab.Py(), gamma1_lab.Pz(), gamma1_lab.E()), photon_pdgID, 1 );
0207
0208     GenParticlePtr p_gamma2 = std::make_shared<GenParticle>(
0209       FourVector(gamma2_lab.Px(), gamma2_lab.Py(), gamma2_lab.Pz(), gamma2_lab.E()), photon_pdgID, 1 );
0210
0211     // Generate pi0 at same position as the lambda. Approximating pi0 decay as instantaneous
0212     GenVertexPtr v_pi0_decay = std::make_shared<GenVertex>(FourVector(lambda_decay_position.X(), lambda_decay_position.Y(), lambda_decay_position.Z(), lambda_decay_time));
0213     v_pi0_decay->add_particle_in(p_pi0);
0214     v_pi0_decay->add_particle_out(p_gamma1);
0215     v_pi0_decay->add_particle_out(p_gamma2);
0216
0217     //std::cout<<  lambda_pvec.Angle(pbeam_dir) << " " << neutron_lab.Angle(pbeam_dir) << " " << gamma1_lab.Angle(pbeam_dir) << " " << gamma2_lab.Angle(pbeam_dir) << std::endl;
0218
0219     evt.add_vertex(v_pi0_decay);
0220
0221     if (events_parsed == 0) {
0222       std::cout << "First event: " << std::endl;
0223       Print::listing(evt);
0224     }
0225     double zdc_z=35800;
0226     TVector3 extrap_neutron=lambda_decay_position+neutron_lab.Vect()*((zdc_z-pbeam_dir.Dot(lambda_decay_position))/(pbeam_dir.Dot(neutron_lab.Vect())));
0227     TVector3 extrap_gamma1=lambda_decay_position+gamma1_lab.Vect()*((zdc_z-pbeam_dir.Dot(lambda_decay_position))/(pbeam_dir.Dot(gamma1_lab.Vect())));
0228     TVector3 extrap_gamma2=lambda_decay_position+gamma2_lab.Vect()*((zdc_z-pbeam_dir.Dot(lambda_decay_position))/(pbeam_dir.Dot(gamma2_lab.Vect())));
0229     if (extrap_neutron.Angle(pbeam_dir)<0.004 && extrap_gamma1.Angle(pbeam_dir)<0.004 && extrap_gamma2.Angle(pbeam_dir)<0.004 && lambda_decay_position.Dot(pbeam_dir)<zdc_z)
0230       hepmc_output.write_event(evt);
0231     if (events_parsed % 1000 == 0) {
0232       std::cout << "Event: " << events_parsed << std::endl;
0233     }
0234     evt.clear();
0235   }
0236   hepmc_output.close();
0237
0238   std::cout << "Events parsed and written: " << events_parsed << std::endl;
0239 }