zdc_pi0/analysis/gen_pi0_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 pi0 mesons and decay them to 2 photons
0035 void gen_pi0_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname = "pi0_decay.hepmc",
0036               double p_min = 60., // 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 = 4; // 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 pi0_info = GetParticleInfo(pdg, "pi0");
0057   double pi0_mass = std::get<0>(pi0_info);
0058   int pi0_pdgID = std::get<1>(pi0_info);
0059   double pi0_lifetime = std::get<2>(pi0_info);
0060
0061   auto photon_info = GetParticleInfo(pdg, "gamma");
0062   double photon_mass = std::get<0>(photon_info);
0063   int photon_pdgID = std::get<1>(photon_info);
0064
0065   for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
0066
0067     //Set the event number
0068     evt.set_event_number(events_parsed);
0069
0070     // FourVector(px,py,pz,e,pdgid,status)
0071     // type 4 is beam
0072     // pdgid 11 - electron
0073     // pdgid 2212 - proton
0074     GenParticlePtr p1 =
0075         std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
0076     GenParticlePtr p2 = std::make_shared<GenParticle>(
0077         FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
0078
0079     // Define momentum with respect to EIC proton beam direction
0080     Double_t pi0_p     = r1->Uniform(p_min, p_max);
0081     Double_t pi0_phi   = r1->Uniform(0.0, 2.0 * M_PI);
0082     Double_t pi0_th    = r1->Uniform(theta_min/1000., theta_max/1000.); // Divide by 1000 for radians
0083     Double_t pi0_px    = pi0_p * TMath::Cos(pi0_phi) * TMath::Sin(pi0_th);
0084     Double_t pi0_py    = pi0_p * TMath::Sin(pi0_phi) * TMath::Sin(pi0_th);
0085     Double_t pi0_pz    = pi0_p * TMath::Cos(pi0_th);
0086     Double_t pi0_E     = TMath::Sqrt(pi0_p*pi0_p + pi0_mass*pi0_mass);
0087
0088     // Rotate to lab coordinate system
0089     TVector3 pi0_pvec(pi0_px, pi0_py, pi0_pz);
0090     double cross_angle = -25./1000.; // in Rad
0091     TVector3 pbeam_dir(TMath::Sin(cross_angle), 0, TMath::Cos(cross_angle)); //proton beam direction
0092     pi0_pvec.RotateY(cross_angle); // Theta is returned positive, beam in negative X
0093
0094     // type 2 is state that will decay
0095     GenParticlePtr p_pi0 = std::make_shared<GenParticle>(
0096         FourVector(pi0_pvec.X(), pi0_pvec.Y(), pi0_pvec.Z(), pi0_E), pi0_pdgID, 2 );
0097
0098     // Generating pi0 particle, will be generated at origin
0099     // Must have input electron + proton for vertex
0100     GenVertexPtr pi0_initial_vertex = std::make_shared<GenVertex>();
0101     pi0_initial_vertex->add_particle_in(p1);
0102     pi0_initial_vertex->add_particle_in(p2);
0103     pi0_initial_vertex->add_particle_out(p_pi0);
0104     evt.add_vertex(pi0_initial_vertex);
0105
0106     // Generate gamma1 + gamma1 in pi0 rest frame
0107     TLorentzVector gamma1_rest, gamma2_rest;
0108
0109     // Generating uniformly along a sphere
0110     double cost_gamma1_rest = r1->Uniform(-1,1);
0111     double th_gamma1_rest = TMath::ACos(cost_gamma1_rest);
0112     double sint_gamma1_rest = TMath::Sin(th_gamma1_rest);
0113
0114     double phi_gamma1_rest = r1->Uniform(-1.*TMath::Pi(),1.*TMath::Pi());
0115     double cosp_gamma1_rest = TMath::Cos(phi_gamma1_rest);
0116     double sinp_gamma1_rest = TMath::Sin(phi_gamma1_rest);
0117
0118     // Calculate energy of each particle in the pi0 rest frame
0119     // This is half the pi0 mass
0120     double E_gamma1_rest = pi0_mass/2;
0121     double E_gamma2_rest = pi0_mass/2;
0122
0123     // Both particles will have the same momentum, and are massless
0124     double momentum_rest = pi0_mass/2;
0125
0126     gamma1_rest.SetE(E_gamma1_rest);
0127     gamma1_rest.SetPx( momentum_rest * sint_gamma1_rest * cosp_gamma1_rest );
0128     gamma1_rest.SetPy( momentum_rest * sint_gamma1_rest * sinp_gamma1_rest );
0129     gamma1_rest.SetPz( momentum_rest * cost_gamma1_rest );
0130
0131     gamma2_rest.SetE(E_gamma2_rest);
0132     gamma2_rest.SetPx( -gamma1_rest.Px() );
0133     gamma2_rest.SetPy( -gamma1_rest.Py() );
0134     gamma2_rest.SetPz( -gamma1_rest.Pz() );
0135
0136     // Boost both gammas to lab frame
0137     TLorentzVector pi0_lab(pi0_pvec.X(), pi0_pvec.Y(), pi0_pvec.Z(), pi0_E);
0138     TVector3 pi0_boost = pi0_lab.BoostVector();
0139     TLorentzVector gamma1_lab, gamma2_lab;
0140     gamma1_lab = gamma1_rest;
0141     gamma1_lab.Boost(pi0_boost);
0142     gamma2_lab = gamma2_rest;
0143     gamma2_lab.Boost(pi0_boost);
0144
0145     // Calculating position for pi0 decay
0146     TVector3 pi0_unit = pi0_lab.Vect().Unit();
0147     double pi0_decay_length = GetDecayLength(r1, pi0_lifetime, pi0_mass, pi0_lab.P());
0148     TVector3 pi0_decay_position = pi0_unit * pi0_decay_length;
0149     double pi0_decay_time = pi0_decay_length / pi0_lab.Beta() ; // Decay time in lab frame in length units (mm)
0150
0151     // Generating vertex for pi0 decay
0152     GenParticlePtr p_gamma1 = std::make_shared<GenParticle>(
0153       FourVector(gamma1_lab.Px(), gamma1_lab.Py(), gamma1_lab.Pz(), gamma1_lab.E()), photon_pdgID, 1 );
0154
0155     GenParticlePtr p_gamma2 = std::make_shared<GenParticle>(
0156       FourVector(gamma2_lab.Px(), gamma2_lab.Py(), gamma2_lab.Pz(), gamma2_lab.E()), photon_pdgID, 1 );
0157
0158     GenVertexPtr v_pi0_decay = std::make_shared<GenVertex>(FourVector(pi0_decay_position.X(), pi0_decay_position.Y(), pi0_decay_position.Z(), pi0_decay_time));
0159     v_pi0_decay->add_particle_in(p_pi0);
0160     v_pi0_decay->add_particle_out(p_gamma1);
0161     v_pi0_decay->add_particle_out(p_gamma2);
0162
0163     evt.add_vertex(v_pi0_decay);
0164
0165     if (events_parsed == 0) {
0166       std::cout << "First event: " << std::endl;
0167       Print::listing(evt);
0168     }
0169     double zdc_z=35800;
0170     TVector3 extrap_gamma1=pi0_decay_position+gamma1_lab.Vect()*((zdc_z-pbeam_dir.Dot(pi0_decay_position))/(pbeam_dir.Dot(gamma1_lab.Vect())));
0171     TVector3 extrap_gamma2=pi0_decay_position+gamma2_lab.Vect()*((zdc_z-pbeam_dir.Dot(pi0_decay_position))/(pbeam_dir.Dot(gamma2_lab.Vect())));
0172     if (extrap_gamma1.Angle(pbeam_dir)<0.004 && extrap_gamma2.Angle(pbeam_dir)<0.004 && pi0_decay_position.Dot(pbeam_dir)<zdc_z)
0173       hepmc_output.write_event(evt);
0174     if (events_parsed % 1000 == 0) {
0175       std::cout << "Event: " << events_parsed << std::endl;
0176     }
0177     evt.clear();
0178   }
0179   hepmc_output.close();
0180
0181   std::cout << "Events parsed and written: " << events_parsed << std::endl;
0182 }