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