barrel_ecal/scripts/emcal_barrel_pi0.cxx

0001 //////////////////////////////////////////////////////////////
0002 // EMCAL Barrel detector
0003 // Single Pi0 dataset
0004 // M. Scott 05/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 using namespace HepMC3;
0020
0021 void emcal_barrel_pi0(int n_events = 1e6, double e_start = 0.0, double e_end = 30.0, const char* out_fname = "./data/emcal_barrel_pi0.hepmc") {
0022   WriterAscii hepmc_output(out_fname);
0023   int events_parsed = 0;
0024   GenEvent evt(Units::GEV, Units::MM);
0025
0026   // Random number generator
0027   TRandom* r1 = new TRandom();
0028
0029   // Constraining the solid angle, but larger than that subtended by the
0030   // detector
0031   // https://indico.bnl.gov/event/7449/contributions/35966/attachments/27177/41430/EIC-DWG-Calo-03192020.pdf
0032   // See a figure on slide 26
0033   double cos_theta_min = std::cos(M_PI * (45.0 / 180.0));
0034   double cos_theta_max = std::cos(M_PI * (135.0 / 180.0));
0035
0036   for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
0037     // FourVector(px,py,pz,e,pdgid,status)
0038     // type 4 is beam
0039     // pdgid 11 - electron
0040     // pdgid 111 - pi0
0041     // pdgid 2212 - proton
0042     GenParticlePtr p1 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
0043     GenParticlePtr p2 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
0044
0045     // Define momentum
0046     Double_t p        = r1->Uniform(e_start, e_end);
0047     Double_t phi      = r1->Uniform(0.0, 2.0 * M_PI);
0048     Double_t costheta = r1->Uniform(cos_theta_min, cos_theta_max);
0049     Double_t theta    = std::acos(costheta);
0050     Double_t px       = p * std::cos(phi) * std::sin(theta);
0051     Double_t py       = p * std::sin(phi) * std::sin(theta);
0052     Double_t pz       = p * std::cos(theta);
0053
0054     // Generates random vectors, uniformly distributed over the surface of a
0055     // sphere of given radius, in this case momentum.
0056     // r1->Sphere(px, py, pz, p);
0057
0058     // type 1 is final state
0059     // pdgid 211 - pion+ 139.570 MeV/c^2
0060     // pdgid 111 - pion0 134.977 MeV/c^2
0061     //GenParticlePtr p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.139570 * 0.139570))), 211, 1);
0062     GenParticlePtr p4 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.134977 * 0.134977))), 111, 1);
0063
0064     GenVertexPtr v1 = std::make_shared<GenVertex>();
0065     v1->add_particle_in(p1);
0066     v1->add_particle_in(p2);
0067
0068     //v1->add_particle_out(p3);
0069     v1->add_particle_out(p4);
0070     evt.add_vertex(v1);
0071
0072     if (events_parsed == 0) {
0073       std::cout << "First event: " << std::endl;
0074       Print::listing(evt);
0075     }
0076
0077     hepmc_output.write_event(evt);
0078     if (events_parsed % 10000 == 0) {
0079       std::cout << "Event: " << events_parsed << std::endl;
0080     }
0081     evt.clear();
0082   }
0083   hepmc_output.close();
0084   std::cout << "Events parsed and written: " << events_parsed << std::endl;
0085 }