benchmarks/zdc_neutron/gen_forward_neutrons.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 <cmath>
0015 #include <math.h>
0016 #include <TMath.h>
0017
0018 using namespace HepMC3;
0019
0020 //Generate single neutrons in the forward detector region.
0021 void gen_forward_neutrons(int n_events = 10000, UInt_t seed = 0, const char* out_fname = "fwd_neutrons.hepmc")
0022 {
0023
0024   double theta_min = 0.0; //in mRad
0025   double theta_max = 6.0; //in mRad
0026   double cost_min = std::cos(theta_max/1000.) ; //Minimum value of cos(theta)
0027   double cost_max = std::cos(theta_min/1000.) ; //Maximum value of cos(theta)
0028
0029   double p_min = 100.; //in GeV/c
0030   double p_max = 100.; //in GeV/c
0031
0032   WriterAscii hepmc_output(out_fname);
0033   int events_parsed = 0;
0034   GenEvent evt(Units::GEV, Units::MM);
0035
0036   // Random number generator
0037   TRandom3 *r1 = new TRandom3(seed); //Default = 0, which uses clock to set seed
0038   cout<<"Random number seed is "<<r1->GetSeed()<<"!"<<endl;
0039
0040   // Getting generated particle information
0041   TDatabasePDG *pdg = new TDatabasePDG();
0042   TParticlePDG *particle = pdg->GetParticle("neutron");
0043   const double mass = particle->Mass();
0044   const int pdgID = particle->PdgCode();
0045
0046   for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
0047
0048     //Set the event number
0049     evt.set_event_number(events_parsed);
0050
0051     // FourVector(px,py,pz,e,pdgid,status)
0052     // type 4 is beam
0053     // pdgid 11 - electron
0054     // pdgid 2212 - proton
0055     GenParticlePtr p1 =
0056         std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
0057     GenParticlePtr p2 = std::make_shared<GenParticle>(
0058         FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
0059
0060     // Define momentum with respect to EIC proton beam direction
0061     Double_t p     = r1->Uniform(p_min,p_max); //Uniform momentum generation
0062     Double_t phi   = r1->Uniform(0.0, 2.0 * M_PI); //Uniform phi generation
0063
0064     Double_t cost = r1->Uniform(cost_min,cost_max); //Uniform cos(theta) generation
0065     Double_t th    = std::acos(cost);
0066
0067     Double_t px    = p * std::cos(phi) * std::sin(th);
0068     Double_t py    = p * std::sin(phi) * std::sin(th);
0069     Double_t pz    = p * std::cos(th);
0070     Double_t E     = sqrt(p*p + mass*mass);
0071
0072     //Rotate to lab coordinate system
0073     TVector3 pvec(px,py,pz);
0074     double cross_angle = -25./1000.; //in Rad
0075     TVector3 pbeam_dir(sin(cross_angle),0,cos(cross_angle)); //proton beam direction
0076     pvec.RotateY(-pbeam_dir.Theta()); // Theta is returned positive, beam in negative X
0077
0078     // type 1 is final state
0079     GenParticlePtr p3 = std::make_shared<GenParticle>(
0080         FourVector(pvec.X(), pvec.Y(), pvec.Z(), E), pdgID, 1 );
0081
0082     GenVertexPtr v1 = std::make_shared<GenVertex>();
0083     v1->add_particle_in(p1);
0084     v1->add_particle_in(p2);
0085
0086     v1->add_particle_out(p3);
0087     evt.add_vertex(v1);
0088
0089     if (events_parsed == 0) {
0090       std::cout << "First event: " << std::endl;
0091       Print::listing(evt);
0092     }
0093
0094     hepmc_output.write_event(evt);
0095     if (events_parsed % 10000 == 0) {
0096       std::cout << "Event: " << events_parsed << std::endl;
0097     }
0098     evt.clear();
0099   }
0100   hepmc_output.close();
0101   std::cout << "Events parsed and written: " << events_parsed << std::endl;
0102 }