b0_tracker/scripts/gen_forward_protons.cxx

0001 #include "HepMC3/GenEvent.h"
0002 #include "HepMC3/ReaderAscii.h"
0003 #include "HepMC3/WriterAscii.h"
0004 #include "HepMC3/Print.h"
0005
0006 #include <iostream>
0007 #include <random>
0008 #include <cmath>
0009 #include <math.h>
0010
0011 #include "TMath.h"
0012 #include "Math/Vector3D.h"
0013 #include "Math/Rotation3D.h"
0014 #include "Math/RotationY.h"
0015 #include "TRandom.h"
0016
0017 #include "common_bench/particles.h"
0018
0019 using namespace HepMC3;
0020
0021 /** Generate electrons in the central region.
0022  *  This is for testing detectors in the "barrel" region.
0023  */
0024 void gen_forward_protons(int n_events = 100,
0025                      const char* out_fname = "forward_protons.hepmc")
0026 {
0027   double cos_theta_min = std::cos(0.5*(M_PI/180.0));
0028   double cos_theta_max = std::cos(0.0*(M_PI/180.0));
0029
0030   const double M_p = common_bench::particleMap.at(2212).mass;
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   TRandom *r1 = new TRandom();
0038
0039   for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
0040     // FourVector(px,py,pz,e,pdgid,status)
0041     // type 4 is beam
0042     // pdgid 11 - electron
0043     // pdgid 111 - pi0
0044     // pdgid 2212 - proton
0045     GenParticlePtr p1 =
0046         std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 2212, 4);
0047     GenParticlePtr p2 = std::make_shared<GenParticle>(
0048         FourVector(0.0, 0.0, 0.0, M_p), 2212, 4);
0049
0050     // Define momentum
0051     Double_t p     = r1->Uniform(200.0, 275.0);
0052     Double_t phi   = r1->Uniform(0.0, 2.0 * M_PI);
0053     Double_t costh = r1->Uniform(cos_theta_min, cos_theta_max);
0054     Double_t th    = std::acos(costh);
0055     Double_t px    = p * std::cos(phi) * std::sin(th);
0056     Double_t py    = p * std::sin(phi) * std::sin(th);
0057     Double_t pz    = p * std::cos(th);
0058
0059
0060     ROOT::Math::XYZVector p0 = {px,py,pz};
0061
0062     //ROOT::Math::Rotation3D r = (-0.025);
0063     ROOT::Math::RotationY r(-0.025);
0064     auto p_rot = r*p0;
0065
0066
0067     // Generates random vectors, uniformly distributed over the surface of a
0068     // sphere of given radius, in this case momentum.
0069     // r1->Sphere(px, py, pz, p);
0070
0071     //std::cout << std::sqrt(px*px + py*py + pz*pz) - p << " is zero? \n";
0072
0073     // type 1 is final state
0074     // pdgid 11 - electron 0.510 MeV/c^2
0075     GenParticlePtr p3 = std::make_shared<GenParticle>(
0076         FourVector(p_rot.x(), p_rot.y(), p_rot.z(), sqrt(p * p + (M_p * M_p))), 2212, 1);
0077
0078     GenVertexPtr v1 = std::make_shared<GenVertex>();
0079     v1->add_particle_in(p1);
0080     v1->add_particle_in(p2);
0081
0082     v1->add_particle_out(p3);
0083     evt.add_vertex(v1);
0084
0085     if (events_parsed == 0) {
0086       std::cout << "First event: " << std::endl;
0087       Print::listing(evt);
0088     }
0089
0090     hepmc_output.write_event(evt);
0091     if (events_parsed % 10000 == 0) {
0092       std::cout << "Event: " << events_parsed << std::endl;
0093     }
0094     evt.clear();
0095   }
0096   hepmc_output.close();
0097   std::cout << "Events parsed and written: " << events_parsed << std::endl;
0098 }