CaloDebugTools/UtilityMacros/MakeHepmcFileWithParticleGun.cxx

0001 //////////////////////////////////////////////////////////////
0002 // 'MakeHepmcFileWithParticleGun.cxx'
0003 // Dhevan Gangadharan, Nico Schmidt
0004 //
0005 // Derived from emcal_barrel_particles_gen.cxx. Generates
0006 // HepMC files for events of 1 or more particles.
0007 //////////////////////////////////////////////////////////////
0008 #include <cmath>
0009 #include <string>
0010 #include <iostream>
0011 #include <math.h>
0012 #include <random>
0013 #include <cctype>
0014 #include <fmt/core.h>
0015
0016 #include "TFile.h"
0017 #include "TMath.h"
0018 #include "TRandom3.h"
0019 #include "TF1.h"
0020 #include "TF2.h"
0021 #include "TH2D.h"
0022 #include "TLorentzVector.h"
0023
0024 #include "HepMC3/GenEvent.h"
0025 #include "HepMC3/Print.h"
0026 #include "HepMC3/ReaderAscii.h"
0027 #include "HepMC3/WriterAscii.h"
0028
0029 using namespace std;
0030 using namespace HepMC3;
0031
0032 std::tuple <int, double> extract_particle_parameters(std::string particle_name);
0033
0034 // converter center: -55609.5, full thickness = 1 mm
0035 // end of spectrometer magnet: -56390
0036 double Vz = 0; // Primary vertex location in mm
0037
0038 double Z = 1; // ion beam particle charge
0039 double electronPz = -18;
0040 double hadronPz = 275;
0041
0042 double prefactor = 2.3179; // 4 alpha r_e^2 (mb)
0043 double protonMass = 0.938272;
0044 double electronMass = 0.51099895e-3;
0045 double photonMass = 0;
0046 double muonMass = 0.1056583745;
0047 double pionZeroMass = 0.1349768;
0048 double pionMass = 0.13957039;
0049
0050 // output name [Derek, 02.22.2024]
0051 const string sOutBase = "forBHCalMultiPartSim.e10h11pipXpim.d8m3y2024";
0052
0053 void MakeHepmcFileWithParticleGun(
0054   int n_events = 1e4,
0055   double E_start = 10.0,
0056   double E_end = 10.0,
0057   double eta_min = -1.1,
0058   double eta_max = 1.1,
0059   string out_fname= sOutBase + ".hepmc"
0060 ) {
0061
0062   // make output root file name [Derek, 02.22.2024]
0063   const string sOutRootName = sOutBase + ".hepmc.root";
0064
0065   TFile *fout = new TFile(sOutRootName.data(), "RECREATE");
0066
0067   cout << "Generating " << n_events << " events" << endl;
0068   cout << "particle energy range: " << E_start << " - " << E_end << " GeV" << endl;
0069   cout << "particle eta range: " << eta_min << " - " << eta_max << endl;
0070
0071
0072   WriterAscii hepmc_output(out_fname);
0073
0074   int events_parsed = 0;
0075
0076   GenEvent evt(Units::GEV, Units::MM);
0077
0078   // Random number generator
0079   TRandom* r1 = new TRandom3();
0080
0081   // Create events
0082   for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
0083     // FourVector( px, py, pz, e, pdgid, status )
0084
0085     // Create a vertex for the beam particles
0086     GenVertexPtr v1 = std::make_shared<GenVertex>();
0087
0088     // Generate beam particles
0089     // status 1 (final state particle),
0090     // status 2 (decayed physical particle)
0091     // status 3 (Documentation line, can be used to indicate in/out particles of hard process)
0092     // status 4 (beam particle)
0093     GenParticlePtr p1 = std::make_shared<GenParticle>(
0094         FourVector(0.0, 0.0, electronPz, sqrt( pow(electronPz, 2) + pow(electronMass, 2) ) ), 11, 4);
0095     GenParticlePtr p2 = std::make_shared<GenParticle>(
0096         FourVector(0.0, 0.0, hadronPz, sqrt( pow(hadronPz, 2) + pow(protonMass, 2) ) ), 2212, 4);
0097     GenParticlePtr p2_2 = std::make_shared<GenParticle>(
0098         FourVector(0.0, 0.0, hadronPz, sqrt( pow(hadronPz, 2) + pow(protonMass, 2) ) ), 2212, 1);
0099
0100     v1->add_particle_in( p1 );
0101     v1->add_particle_in( p2 );
0102     v1->add_particle_out( p2_2 ); // make the output particle a beam proton (not to be in analysis)
0103
0104     /////////////////////////////////////////////////////////////////
0105     // Now add particles of interest at separate vertices
0106     // vertices must have at least 1 outgoing and 1 incoming particle
0107
0108     double E = r1->Uniform(E_start, E_end);
0109     double p = sqrt( pow(E,2) - pow(pionMass,2) );
0110     double eta1 = r1->Uniform(eta_min, eta_max);
0111     double eta2 = r1->Uniform(eta_min, eta_max);
0112     double theta1 = 2*atan(exp(-eta1));
0113     double theta2 = 2*atan(exp(-eta2));
0114     // double theta1 = TMath::Pi();
0115     // double theta2 = TMath::Pi();
0116     double phi1 = r1->Uniform(0, 2*TMath::Pi());
0117     double phi2 = r1->Uniform(0, 2*TMath::Pi());
0118     auto [id, mass] = extract_particle_parameters( "piplus" );
0119
0120     // pion 1
0121     GenVertexPtr v2 = std::make_shared<GenVertex>();
0122     v2->set_position( FourVector(0.0, 0.0, 0.0, 0) ); // in mm
0123     // v2->set_position( FourVector(0.0, 120.0, -64000, 0) ); // in mm
0124     GenParticlePtr p3_in = std::make_shared<GenParticle>(
0125         FourVector(p*sin(theta1)*cos(phi1), p*sin(theta1)*sin(phi1), p*cos(theta1), E), id, 3);
0126     GenParticlePtr p3_out = std::make_shared<GenParticle>(
0127         FourVector(p*sin(theta1)*cos(phi1), p*sin(theta1)*sin(phi1), p*cos(theta1), E), id, 1);
0128
0129     v2->add_particle_in( p3_in );
0130     v1->add_particle_out( p3_out );
0131
0132
0133     E = r1->Uniform(E_start, E_end);
0134     p = sqrt( pow(E,2) - pow(pionMass,2) );
0135
0136     auto [id2, mass2] = extract_particle_parameters( "piminus" );
0137     // pion 2
0138     GenVertexPtr v3 = std::make_shared<GenVertex>();
0139     v3->set_position( FourVector(0.0, 0.0, 0.0, 0) ); // in mm
0140     GenParticlePtr p4_in = std::make_shared<GenParticle>(
0141         FourVector(p*sin(theta2)*cos(phi2), p*sin(theta2)*sin(phi2), p*cos(theta2), E), id2, 3);
0142     GenParticlePtr p4_out = std::make_shared<GenParticle>(
0143         FourVector(p*sin(theta2)*cos(phi2), p*sin(theta2)*sin(phi2), p*cos(theta2), E), id2, 1);
0144
0145     v3->add_particle_in( p4_in );
0146     v1->add_particle_out( p4_out );
0147
0148     evt.add_vertex( v1 );
0149     //evt.add_vertex( v2 );
0150     //evt.add_vertex( v3 );
0151
0152
0153     // Shift the whole event to specific point
0154     evt.shift_position_to( FourVector(0,0, Vz, 0) );
0155
0156     if (events_parsed == 0) {
0157       std::cout << "First event: " << std::endl;
0158       Print::listing(evt);
0159     }
0160
0161     hepmc_output.write_event(evt);
0162
0163     if (events_parsed % 10000 == 0) {
0164       std::cout << "Event: " << events_parsed << std::endl;
0165     }
0166
0167     evt.clear();
0168   }
0169   hepmc_output.close();
0170   std::cout << "Events parsed and written: " << events_parsed << std::endl;
0171
0172   fout->Close();
0173
0174 }
0175
0176 //----------------------------------------------------------------------------
0177 // Returns particle pdgID and mass in [GeV]
0178 std::tuple <int, double> extract_particle_parameters(std::string particle_name) {
0179   if (particle_name == "electron") return std::make_tuple(11,    electronMass);
0180   if (particle_name == "photon")   return std::make_tuple(22,    photonMass);
0181   if (particle_name == "positron") return std::make_tuple(-11,   electronMass);
0182   if (particle_name == "proton")   return std::make_tuple(2212,  protonMass);
0183   if (particle_name == "muon")     return std::make_tuple(13,    muonMass);
0184   if (particle_name == "pi0")      return std::make_tuple(111,   pionZeroMass);
0185   if (particle_name == "piplus")   return std::make_tuple(211,   pionMass);
0186   if (particle_name == "piminus")  return std::make_tuple(-211,  pionMass);
0187
0188   std::cout << "wrong particle name" << std::endl;
0189   abort();
0190 }
0191