File indexing completed on 2024-09-27 07:02:37
0001
0002
0003
0004
0005
0006
0007
0008 #include "HepMC3/GenEvent.h"
0009 #include "HepMC3/Print.h"
0010 #include "HepMC3/ReaderAscii.h"
0011 #include "HepMC3/WriterAscii.h"
0012
0013 #include "TMath.h"
0014 #include "TRandom.h"
0015
0016 #include <cmath>
0017 #include <iostream>
0018 #include <fstream>
0019 #include <math.h>
0020 #include <random>
0021 #include <cstdlib>
0022 #include <fmt/core.h>
0023
0024 #include "hcal_barrel_common_functions.h"
0025
0026 using namespace HepMC3;
0027
0028 void hcal_barrel_particles_gen(int n_events = 1e6, double e_start = 0.0, double e_end = 20.0, std::string particle_name = "muon") {
0029 std::string out_fname;
0030 auto env_fname = getenv("JUGGLER_GEN_FILE");
0031 if (env_fname != nullptr)
0032 out_fname = env_fname;
0033 else
0034 out_fname = fmt::format("./data/hcal_barrel_{}.hepmc", particle_name);
0035 WriterAscii hepmc_output(out_fname);
0036 int events_parsed = 0;
0037 GenEvent evt(Units::GEV, Units::MM);
0038
0039
0040 TRandom* r1 = new TRandom();
0041
0042
0043
0044
0045
0046 double cos_theta_min = std::cos(M_PI * (45.0 / 180.0));
0047 double cos_theta_max = std::cos(M_PI * (135.0 / 180.0));
0048
0049 for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
0050
0051
0052
0053
0054 GenParticlePtr p1 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
0055 GenParticlePtr p2 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
0056
0057
0058 Double_t p = r1->Uniform(e_start, e_end);
0059 Double_t phi = r1->Uniform(0.0, 2.0 * M_PI);
0060 Double_t costheta = r1->Uniform(cos_theta_min, cos_theta_max);
0061 Double_t theta = std::acos(costheta);
0062 Double_t px = p * std::cos(phi) * std::sin(theta);
0063 Double_t py = p * std::sin(phi) * std::sin(theta);
0064 Double_t pz = p * std::cos(theta);
0065
0066
0067 auto [id, mass] = extract_particle_parameters(particle_name);
0068 GenParticlePtr p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (mass * mass))), id, 1);
0069
0070 GenVertexPtr v1 = std::make_shared<GenVertex>();
0071 v1->add_particle_in(p1);
0072 v1->add_particle_in(p2);
0073 v1->add_particle_out(p3);
0074 evt.add_vertex(v1);
0075
0076 if (events_parsed == 0) {
0077 std::cout << "First event: " << std::endl;
0078 Print::listing(evt);
0079 }
0080
0081 hepmc_output.write_event(evt);
0082 if (events_parsed % 10000 == 0) {
0083 std::cout << "Event: " << events_parsed << std::endl;
0084 }
0085 evt.clear();
0086 }
0087 hepmc_output.close();
0088 std::cout << "Events parsed and written: " << events_parsed << std::endl;
0089 }
0090