zdc/scripts/gen_zdc_particles.cxx

0001 //////////////////////////////////////////////////////////////
0002 // ZDC detector
0003 // Single Neutron dataset
0004 // J.KIM 05/07/2021
0005 //////////////////////////////////////////////////////////////
0006 #include "HepMC3/GenEvent.h"
0007 #include "HepMC3/Print.h"
0008 #include "HepMC3/ReaderAscii.h"
0009 #include "HepMC3/WriterAscii.h"
0010
0011 #include <cmath>
0012 #include <iostream>
0013 #include <math.h>
0014 #include <random>
0015
0016 #include "TMath.h"
0017 #include "TRandom.h"
0018
0019 #include <Math/Vector3D.h>
0020 #include <Math/RotationY.h>
0021
0022 using namespace HepMC3;
0023
0024 void gen_zdc_particles(int n_events = 1e6, const std::string& particle = "neutron", double p_start = 125.0, double p_end = 145.0, const std::string& out_fname = "./data/zdc_neutrons.hepmc") {
0025   WriterAscii hepmc_output(out_fname);
0026   int events_parsed = 0;
0027   GenEvent evt(Units::GEV, Units::MM);
0028
0029   // Random number generator
0030   TRandom* r1 = new TRandom();
0031
0032   // Set crossing angle [rad]
0033   double crossing_angle = -0.025;
0034
0035   // Constraining the solid angle, but larger than that subtended by the
0036   // detector
0037   // https://indico.bnl.gov/event/7449/contributions/35966/attachments/27177/41430/EIC-DWG-Calo-03192020.pdf
0038   // See a figure on slide 26
0039   double cos_theta_min = std::cos(0.0);
0040   double cos_theta_max = std::cos(0.005);
0041
0042   for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
0043     // FourVector(px,py,pz,e,pdgid,status)
0044     // type 4 is beam
0045     // pdgid 11 - electron
0046     // pdgid 111 - pi0
0047     // pdgid 2212 - proton
0048     // pdgid 2112 - neutron
0049     GenParticlePtr p1 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
0050     GenParticlePtr p2 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
0051
0052     // Define momentum
0053     Double_t p_mag    = r1->Uniform(p_start, p_end);
0054     Double_t phi      = r1->Uniform(0.0, 2.0 * M_PI);
0055     Double_t costheta = r1->Uniform(cos_theta_min, cos_theta_max);
0056     Double_t theta    = std::acos(costheta);
0057     Double_t p0_x     = p_mag * std::cos(phi) * std::sin(theta);
0058     Double_t p0_y     = p_mag * std::sin(phi) * std::sin(theta);
0059     Double_t p0_z     = p_mag * std::cos(theta);
0060
0061     // Rotate the vector in Y by crossing angle when particles are being generated
0062     ROOT::Math::XYZVector p0{p0_x, p0_y, p0_z};
0063     ROOT::Math::RotationY r(crossing_angle);
0064     auto p = r * p0;
0065     auto px = p.X();
0066     auto py = p.Y();
0067     auto pz = p.Z();
0068
0069     // type 1 is final state
0070     // pdgid 11 - electron 0.510 MeV/c^2
0071     // pdgid 22 - photon massless
0072     // pdgid 2112 - neutron 939.565 MeV/c^2
0073     GenParticlePtr p3;
0074     if (particle == "neutron") {
0075       p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, std::hypot(p_mag, 0.939565)), 2112, 1);
0076     } else if (particle == "photon") {
0077       p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, p_mag), 22, 1);
0078     } else {
0079       std::cout << "Particle type " << particle << " not recognized!" << std::endl;
0080       exit(-1);
0081     }
0082
0083     GenVertexPtr v1 = std::make_shared<GenVertex>();
0084     v1->add_particle_in(p1);
0085     v1->add_particle_in(p2);
0086
0087     v1->add_particle_out(p3);
0088     evt.add_vertex(v1);
0089
0090     if (events_parsed == 0) {
0091       std::cout << "First event: " << std::endl;
0092       Print::listing(evt);
0093     }
0094
0095     hepmc_output.write_event(evt);
0096     if (events_parsed % 10000 == 0) {
0097       std::cout << "Event: " << events_parsed << std::endl;
0098     }
0099     evt.clear();
0100   }
0101   hepmc_output.close();
0102   std::cout << "Events parsed and written: " << events_parsed << std::endl;
0103 }