tracking_detectors/scripts/gen_track_hits.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 "TRandom.h"
0013
0014 using namespace HepMC3;
0015
0016 /** Generate multiple electrons/positron tracks in the central region.
0017  *  This is for testing detectors in the "barrel" region.
0018  */
0019 void gen_track_hits(int n_events = 100,
0020                          const char* out_fname = "track_hits.hepmc",
0021                          int n_parts = 2)
0022 {
0023   double cos_theta_min = std::cos( 1.0*(M_PI/180.0));
0024   double cos_theta_max = std::cos(189.0*(M_PI/180.0));
0025
0026   WriterAscii hepmc_output(out_fname);
0027   int events_parsed = 0;
0028   GenEvent evt(Units::GEV, Units::MM);
0029
0030   // Random number generator
0031   TRandom *r1 = new TRandom();
0032
0033   for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
0034     // FourVector(px,py,pz,e,pdgid,status)
0035     // type 4 is beam
0036     // pdgid 11 - electron
0037     // pdgid 111 - pi0
0038     // pdgid 2212 - proton
0039     for (int ip = 0; ip < n_parts; ip++) {
0040       GenParticlePtr p1 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
0041       GenParticlePtr p2 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
0042
0043       // Define momentum
0044       Double_t p     = r1->Uniform(1.0, 10.0);
0045       Double_t phi   = r1->Uniform(0.0, 2.0 * M_PI);
0046       Double_t costh = r1->Uniform(cos_theta_min, cos_theta_max);
0047       Double_t th    = std::acos(costh);
0048       Double_t px    = p * std::cos(phi) * std::sin(th);
0049       Double_t py    = p * std::sin(phi) * std::sin(th);
0050       Double_t pz    = p * std::cos(th);
0051       // Generates random vectors, uniformly distributed over the surface of a
0052       // sphere of given radius, in this case momentum.
0053       // r1->Sphere(px, py, pz, p);
0054
0055       // std::cout << std::sqrt(px*px + py*py + pz*pz) - p << " is zero? \n";
0056
0057       // type 1 is final state
0058       // pdgid 11 - electron 0.510 MeV/c^2
0059       GenParticlePtr p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.000511 * 0.000511))),
0060                                                         ((ip % 2 == 0) ? 11 : -11), 1);
0061
0062       GenVertexPtr v1 = std::make_shared<GenVertex>();
0063       v1->add_particle_in(p1);
0064       v1->add_particle_in(p2);
0065
0066       v1->add_particle_out(p3);
0067       evt.add_vertex(v1);
0068     }
0069
0070     if (events_parsed == 0) {
0071       std::cout << "First event: " << std::endl;
0072       Print::listing(evt);
0073     }
0074
0075     hepmc_output.write_event(evt);
0076     if (events_parsed % 10000 == 0) {
0077       std::cout << "Event: " << events_parsed << std::endl;
0078     }
0079     evt.clear();
0080   }
0081   hepmc_output.close();
0082   std::cout << "Events parsed and written: " << events_parsed << std::endl;
0083 }