zdc_sigma/analysis/gen_sigma_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 sigma baryons and decay them to a neutron + 2 photons
0035 void gen_sigma_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname = "sigma_decay.hepmc",
0036               double p_min = 100., // in GeV/c
0037               double p_max = 275.) // in GeV/c
0038 {
0039   const double theta_min = 0.0; // in mRad
0040   const double theta_max = 3.0; // in mRad
0041   //const double p_min = 100.; // in GeV/c
0042   //const double p_max = 275.; // in GeV/c
0043
0044   WriterAscii hepmc_output(out_fname);
0045   GenEvent evt(Units::GEV, Units::MM);
0046
0047   // Random number generator
0048   TRandom3 *r1 = new TRandom3(seed); //Default = 0, which uses clock to set seed
0049   cout<<"Random number seed is "<<r1->GetSeed()<<"!"<<endl;
0050
0051   // Getting generated particle information
0052   TDatabasePDG *pdg = new TDatabasePDG();
0053
0054   auto sigma_info = GetParticleInfo(pdg, "Sigma0");
0055   double sigma_mass = std::get<0>(sigma_info);
0056   int sigma_pdgID = std::get<1>(sigma_info);
0057   double sigma_lifetime = std::get<2>(sigma_info);
0058
0059   auto lambda_info = GetParticleInfo(pdg, "Lambda0");
0060   double lambda_mass = std::get<0>(lambda_info);
0061   int lambda_pdgID = std::get<1>(lambda_info);
0062   double lambda_lifetime = std::get<2>(lambda_info);
0063
0064   auto neutron_info = GetParticleInfo(pdg, "neutron");
0065   double neutron_mass = std::get<0>(neutron_info);
0066   int neutron_pdgID = std::get<1>(neutron_info);
0067
0068   auto pi0_info = GetParticleInfo(pdg, "pi0");
0069   double pi0_mass = std::get<0>(pi0_info);
0070   int pi0_pdgID = std::get<1>(pi0_info);
0071   double pi0_lifetime = std::get<2>(pi0_info);
0072
0073   auto photon_info = GetParticleInfo(pdg, "gamma");
0074   double photon_mass = std::get<0>(photon_info);
0075   int photon_pdgID = std::get<1>(photon_info);
0076
0077   int accepted_events = 0;
0078   while (accepted_events < n_events) {
0079
0080     //Set the event number
0081     evt.set_event_number(accepted_events);
0082
0083     // FourVector(px,py,pz,e,pdgid,status)
0084     // type 4 is beam
0085     // pdgid 11 - electron
0086     // pdgid 2212 - proton
0087     GenParticlePtr p1 =
0088         std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
0089     GenParticlePtr p2 = std::make_shared<GenParticle>(
0090         FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
0091
0092     // Define momentum with respect to EIC proton beam direction
0093     Double_t sigma_p     = r1->Uniform(p_min, p_max);
0094     Double_t sigma_phi   = r1->Uniform(0.0, 2.0 * M_PI);
0095     Double_t sigma_th    = r1->Uniform(theta_min/1000., theta_max/1000.); // Divide by 1000 for radians
0096     Double_t sigma_px    = sigma_p * TMath::Cos(sigma_phi) * TMath::Sin(sigma_th);
0097     Double_t sigma_py    = sigma_p * TMath::Sin(sigma_phi) * TMath::Sin(sigma_th);
0098     Double_t sigma_pz    = sigma_p * TMath::Cos(sigma_th);
0099     Double_t sigma_E     = TMath::Sqrt(sigma_p*sigma_p + sigma_mass*sigma_mass);
0100
0101
0102     TVector3 sigma_pvec(sigma_px, sigma_py, sigma_pz);
0103
0104     double cross_angle = -25./1000.; // in Rad
0105     TVector3 pbeam_dir(TMath::Sin(cross_angle), 0, TMath::Cos(cross_angle)); //proton beam direction
0106     sigma_pvec.RotateY(cross_angle); // Theta is returned positive, beam in negative X
0107
0108     // type 2 is state that will decay
0109     GenParticlePtr p_sigma = std::make_shared<GenParticle>(
0110         FourVector(sigma_pvec.X(), sigma_pvec.Y(), sigma_pvec.Z(), sigma_E), sigma_pdgID, 2 );
0111     // Generating sigma particle, will be generated at origin
0112         // Must have input electron + proton for vertex
0113         GenVertexPtr sigma_initial_vertex = std::make_shared<GenVertex>();
0114         sigma_initial_vertex->add_particle_in(p1);
0115         sigma_initial_vertex->add_particle_in(p2);
0116         sigma_initial_vertex->add_particle_out(p_sigma);
0117         evt.add_vertex(sigma_initial_vertex);
0118
0119         // Generate lambda + gamma in sigma rest frame
0120         TLorentzVector lambda_rest, gamma_rest;
0121
0122         // Generating uniformly along a sphere
0123         double cost_lambda_rest = r1->Uniform(-1,1);
0124         double th_lambda_rest = TMath::ACos(cost_lambda_rest);
0125         double sint_lambda_rest = TMath::Sin(th_lambda_rest);
0126
0127         double phi_lambda_rest = r1->Uniform(-1.*TMath::Pi(),1.*TMath::Pi());
0128         double cosp_lambda_rest = TMath::Cos(phi_lambda_rest);
0129         double sinp_lambda_rest = TMath::Sin(phi_lambda_rest);
0130
0131         // Calculate energy of each particle in the sigma rest frame
0132         // See problem 3.19 in Introduction to Elementary Particles, 2nd edition by D. Griffiths
0133         double E_lambda_rest = (-TMath::Power(photon_mass, 2.) + TMath::Power(sigma_mass, 2.) + TMath::Power(lambda_mass, 2.) ) / (2. * sigma_mass) ;
0134         double E_gamma_rest = (-TMath::Power(lambda_mass, 2.) + TMath::Power(sigma_mass, 2.) + TMath::Power(photon_mass, 2.) ) / (2. * sigma_mass) ;
0135
0136         // Both particles will have the same momentum, so just use lambda variables
0137         double momentum_rest = TMath::Sqrt( E_lambda_rest*E_lambda_rest - lambda_mass*lambda_mass );
0138
0139         lambda_rest.SetE(E_lambda_rest);
0140         lambda_rest.SetPx( momentum_rest * sint_lambda_rest * cosp_lambda_rest );
0141         lambda_rest.SetPy( momentum_rest * sint_lambda_rest * sinp_lambda_rest );
0142         lambda_rest.SetPz( momentum_rest * cost_lambda_rest );
0143
0144         gamma_rest.SetE(E_gamma_rest);
0145         gamma_rest.SetPx( -lambda_rest.Px() );
0146         gamma_rest.SetPy( -lambda_rest.Py() );
0147         gamma_rest.SetPz( -lambda_rest.Pz() );
0148
0149         // Boost lambda & pion to lab frame
0150         TLorentzVector sigma_lab(sigma_pvec.X(), sigma_pvec.Y(), sigma_pvec.Z(), sigma_E);
0151         TVector3 sigma_boost = sigma_lab.BoostVector();
0152         TLorentzVector lambda_lab, gamma_lab;
0153         lambda_lab = lambda_rest;
0154         lambda_lab.Boost(sigma_boost);
0155         gamma_lab = gamma_rest;
0156         gamma_lab.Boost(sigma_boost);
0157
0158     // Calculating position for sigma decay
0159     TVector3 sigma_unit = sigma_lab.Vect().Unit();
0160     double sigma_decay_length = GetDecayLength(r1, sigma_lifetime, sigma_mass, sigma_lab.P());
0161     TVector3 sigma_decay_position = sigma_unit * sigma_decay_length;
0162     double sigma_decay_time = sigma_decay_length / sigma_lab.Beta() ; // Decay time in lab frame in length units (mm)
0163
0164     // Generating vertex for sigma decay
0165     GenParticlePtr p_lambda = std::make_shared<GenParticle>(
0166       FourVector(lambda_lab.Px(), lambda_lab.Py(), lambda_lab.Pz(), lambda_lab.E()), lambda_pdgID, 2 );
0167
0168     GenParticlePtr p_gamma = std::make_shared<GenParticle>(
0169       FourVector(gamma_lab.Px(), gamma_lab.Py(), gamma_lab.Pz(), gamma_lab.E()), photon_pdgID, 1 );
0170
0171     GenVertexPtr v_sigma_decay = std::make_shared<GenVertex>(FourVector(sigma_decay_position.X(), sigma_decay_position.Y(), sigma_decay_position.Z(), sigma_decay_time));
0172     v_sigma_decay->add_particle_in(p_sigma);
0173     v_sigma_decay->add_particle_out(p_lambda);
0174     v_sigma_decay->add_particle_out(p_gamma);
0175
0176     evt.add_vertex(v_sigma_decay);
0177
0178     // Generate neutron + pi0 in lambda rest frame
0179     TLorentzVector neutron_rest, pi0_rest;
0180
0181     // Generating uniformly along a sphere
0182     double cost_neutron_rest = r1->Uniform(-1,1);
0183     double th_neutron_rest = TMath::ACos(cost_neutron_rest);
0184     double sint_neutron_rest = TMath::Sin(th_neutron_rest);
0185
0186     double phi_neutron_rest = r1->Uniform(-1.*TMath::Pi(),1.*TMath::Pi());
0187     double cosp_neutron_rest = TMath::Cos(phi_neutron_rest);
0188     double sinp_neutron_rest = TMath::Sin(phi_neutron_rest);
0189
0190     // Calculate energy of each particle in the lambda rest frame
0191     // See problem 3.19 in Introduction to Elementary Particles, 2nd edition by D. Griffiths
0192     double E_neutron_rest = (-TMath::Power(pi0_mass, 2.) + TMath::Power(lambda_mass, 2.) + TMath::Power(neutron_mass, 2.) ) / (2. * lambda_mass) ;
0193     double E_pi0_rest = (-TMath::Power(neutron_mass, 2.) + TMath::Power(lambda_mass, 2.) + TMath::Power(pi0_mass, 2.) ) / (2. * lambda_mass) ;
0194
0195     // Both particles will have the same momentum, so just use neutron variables
0196     momentum_rest = TMath::Sqrt( E_neutron_rest*E_neutron_rest - neutron_mass*neutron_mass );
0197
0198     neutron_rest.SetE(E_neutron_rest);
0199     neutron_rest.SetPx( momentum_rest * sint_neutron_rest * cosp_neutron_rest );
0200     neutron_rest.SetPy( momentum_rest * sint_neutron_rest * sinp_neutron_rest );
0201     neutron_rest.SetPz( momentum_rest * cost_neutron_rest );
0202
0203     pi0_rest.SetE(E_pi0_rest);
0204     pi0_rest.SetPx( -neutron_rest.Px() );
0205     pi0_rest.SetPy( -neutron_rest.Py() );
0206     pi0_rest.SetPz( -neutron_rest.Pz() );
0207
0208     // Boost neutron & pion to lab frame
0209     TVector3 lambda_boost = lambda_lab.BoostVector();
0210     TLorentzVector neutron_lab, pi0_lab;
0211     neutron_lab = neutron_rest;
0212     neutron_lab.Boost(lambda_boost);
0213     pi0_lab = pi0_rest;
0214     pi0_lab.Boost(lambda_boost);
0215
0216     // Calculating position for lambda decay
0217     TVector3 lambda_unit = lambda_lab.Vect().Unit();
0218     double lambda_decay_length = GetDecayLength(r1, lambda_lifetime, lambda_mass, lambda_lab.P());
0219     TVector3 lambda_decay_position = lambda_unit * lambda_decay_length;
0220     double lambda_decay_time = lambda_decay_length / lambda_lab.Beta() ; // Decay time in lab frame in length units (mm)
0221
0222     // Generating vertex for lambda decay
0223     GenParticlePtr p_neutron = std::make_shared<GenParticle>(
0224       FourVector(neutron_lab.Px(), neutron_lab.Py(), neutron_lab.Pz(), neutron_lab.E()), neutron_pdgID, 1 );
0225
0226     GenParticlePtr p_pi0 = std::make_shared<GenParticle>(
0227       FourVector(pi0_lab.Px(), pi0_lab.Py(), pi0_lab.Pz(), pi0_lab.E()), pi0_pdgID, 2 );
0228
0229     GenVertexPtr v_lambda_decay = std::make_shared<GenVertex>(FourVector(lambda_decay_position.X(), lambda_decay_position.Y(), lambda_decay_position.Z(), lambda_decay_time));
0230     v_lambda_decay->add_particle_in(p_lambda);
0231     v_lambda_decay->add_particle_out(p_neutron);
0232     v_lambda_decay->add_particle_out(p_pi0);
0233
0234     evt.add_vertex(v_lambda_decay);
0235
0236     // Generate two photons from pi0 decay
0237     TLorentzVector gamma1_rest, gamma2_rest;
0238
0239     // Generating uniformly along a sphere
0240     double cost_gamma1_rest = r1->Uniform(-1,1);
0241     double th_gamma1_rest = TMath::ACos(cost_gamma1_rest);
0242     double sint_gamma1_rest = TMath::Sin(th_gamma1_rest);
0243
0244     double phi_gamma1_rest = r1->Uniform(-1.*TMath::Pi(),1.*TMath::Pi());
0245     double cosp_gamma1_rest = TMath::Cos(phi_gamma1_rest);
0246     double sinp_gamma1_rest = TMath::Sin(phi_gamma1_rest);
0247
0248     // Photons are massless so they each get equal energies
0249     gamma1_rest.SetE(pi0_mass/2.);
0250     gamma1_rest.SetPx( (pi0_mass/2.)*sint_gamma1_rest*cosp_gamma1_rest );
0251     gamma1_rest.SetPy( (pi0_mass/2.)*sint_gamma1_rest*sinp_gamma1_rest );
0252     gamma1_rest.SetPz( (pi0_mass/2.)*cost_gamma1_rest );
0253
0254     gamma2_rest.SetE(pi0_mass/2.);
0255     gamma2_rest.SetPx( -gamma1_rest.Px() );
0256     gamma2_rest.SetPy( -gamma1_rest.Py() );
0257     gamma2_rest.SetPz( -gamma1_rest.Pz() );
0258
0259     // Boost neutron & pion to lab frame
0260     TVector3 pi0_boost = pi0_lab.BoostVector();
0261     TLorentzVector gamma1_lab, gamma2_lab;
0262     gamma1_lab = gamma1_rest;
0263     gamma1_lab.Boost(pi0_boost);
0264     gamma2_lab = gamma2_rest;
0265     gamma2_lab.Boost(pi0_boost);
0266
0267     GenParticlePtr p_gamma1 = std::make_shared<GenParticle>(
0268       FourVector(gamma1_lab.Px(), gamma1_lab.Py(), gamma1_lab.Pz(), gamma1_lab.E()), photon_pdgID, 1 );
0269
0270     GenParticlePtr p_gamma2 = std::make_shared<GenParticle>(
0271       FourVector(gamma2_lab.Px(), gamma2_lab.Py(), gamma2_lab.Pz(), gamma2_lab.E()), photon_pdgID, 1 );
0272
0273     // Generate pi0 at same position as the lambda. Approximating pi0 decay as instantaneous
0274     GenVertexPtr v_pi0_decay = std::make_shared<GenVertex>(FourVector(lambda_decay_position.X(), lambda_decay_position.Y(), lambda_decay_position.Z(), lambda_decay_time));
0275     v_pi0_decay->add_particle_in(p_pi0);
0276     v_pi0_decay->add_particle_out(p_gamma1);
0277     v_pi0_decay->add_particle_out(p_gamma2);
0278
0279     //std::cout<<  lambda_pvec.Angle(pbeam_dir) << " " << neutron_lab.Angle(pbeam_dir) << " " << gamma1_lab.Angle(pbeam_dir) << " " << gamma2_lab.Angle(pbeam_dir) << std::endl;
0280
0281     evt.add_vertex(v_pi0_decay);
0282
0283     if (accepted_events == 0) {
0284       std::cout << "First event: " << std::endl;
0285       Print::listing(evt);
0286     }
0287     double zdc_z=35800;
0288     TVector3 extrap_gamma=sigma_decay_position+gamma_lab.Vect()*((zdc_z-pbeam_dir.Dot(sigma_decay_position))/(pbeam_dir.Dot(gamma_lab.Vect())));
0289     TVector3 extrap_neutron=lambda_decay_position+neutron_lab.Vect()*((zdc_z-pbeam_dir.Dot(lambda_decay_position))/(pbeam_dir.Dot(neutron_lab.Vect())));
0290     TVector3 extrap_gamma1=lambda_decay_position+gamma1_lab.Vect()*((zdc_z-pbeam_dir.Dot(lambda_decay_position))/(pbeam_dir.Dot(gamma1_lab.Vect())));
0291     TVector3 extrap_gamma2=lambda_decay_position+gamma2_lab.Vect()*((zdc_z-pbeam_dir.Dot(lambda_decay_position))/(pbeam_dir.Dot(gamma2_lab.Vect())));
0292     if (extrap_neutron.Angle(pbeam_dir)<0.004 && extrap_gamma1.Angle(pbeam_dir)<0.004 && extrap_gamma2.Angle(pbeam_dir)<0.004 && extrap_gamma.Angle(pbeam_dir)<0.004 && lambda_decay_position.Dot(pbeam_dir)<zdc_z){
0293       hepmc_output.write_event(evt);
0294       if (accepted_events % 1000 == 0) {
0295         std::cout << "Event: " << accepted_events << std::endl;
0296       }
0297       accepted_events ++;
0298     }
0299     evt.clear();
0300   }
0301   hepmc_output.close();
0302
0303   std::cout << "Events parsed and written: " << accepted_events << std::endl;
0304 }