source/pfRICH/hepmc_writer.C

0001 //

0002 // export LD_LIBRARY_PATH=/home/eic/lib64:${LD_LIBRARY_PATH}

0003 //

0004 // root -l

0005 //

0006 // root [0] gSystem->AddIncludePath("-I/home/eic/include");

0007 // root [1] .x hepmc_writer_two_particles.C+("out.hepmc", 1000)

0008 //

0009
0010 #include "HepMC3/GenEvent.h"
0011 #include "HepMC3/ReaderAscii.h"
0012 #include "HepMC3/WriterAscii.h"
0013 #include "HepMC3/Print.h"
0014
0015 #include <iostream>
0016 #include <random>
0017 #include <cmath>
0018 #include <math.h>
0019
0020 #include <TMath.h>
0021 #include <TRandom3.h>
0022 #include <TDatabasePDG.h>
0023
0024 using namespace HepMC3;
0025
0026 void hepmc_writer(const char *out_fname, int part, Double_t thMin, Double_t thMax, Double_t pMin, Double_t pMax, int n_events)
0027 {
0028   auto *DatabasePDG = new TDatabasePDG();
0029   auto *electron = DatabasePDG->GetParticle(11);
0030   auto *pion = DatabasePDG->GetParticle(211); //

0031   auto *kaon = DatabasePDG->GetParticle(321);
0032   auto *proton = DatabasePDG->GetParticle(2212);
0033   //const char * out_fname= "out.hepmc";

0034   WriterAscii hepmc_output(out_fname);
0035   int events_parsed = 0;
0036   GenEvent evt(Units::GEV, Units::MM);
0037
0038   // Random number generator

0039   //TRandom *rdmn_gen = new TRandom(0x12345678);

0040   TRandom3 *rdmn_gen = new TRandom3(0);
0041
0042   for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
0043     // type 4 is beam

0044     GenParticlePtr p1 =
0045         std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
0046     GenParticlePtr p2 = std::make_shared<GenParticle>(
0047         FourVector(0.0, 0.0, 100.0, 100.004), 2212, 4);
0048
0049     GenVertexPtr v1 = std::make_shared<GenVertex>();
0050     v1->add_particle_in(p1);
0051     v1->add_particle_in(p2);
0052
0053     // Set Kinematics

0054     //Double_t eta   = rdmn_gen->Uniform(etaMin, etaMax);

0055     //Double_t th    = 2*std::atan(exp(-eta));

0056     Double_t th    = rdmn_gen->Uniform(thMin, thMax);
0057     Double_t p     = rdmn_gen->Uniform(pMin, pMax);
0058     Double_t phi   = rdmn_gen->Uniform(0.0*M_PI/180.0, 360.0*M_PI/180.);
0059     //Double_t phi   = rdmn_gen->Uniform(85.0*M_PI/180.0, 95.0*M_PI/180.);

0060
0061     Double_t px    = p * std::cos(phi) * std::sin(th);
0062     Double_t py    = p * std::sin(phi) * std::sin(th);
0063     Double_t pz    = p * std::cos(th);
0064
0065     /*

0066     if(events_parsed < 50)

0067       {

0068     cout << events_parsed << " " << p << " " << eta << " " << th << " " << phi << " " << px << " " << py << " " << pz << endl;

0069       }

0070     */
0071
0072     TParticlePDG *particle;
0073     switch(part)
0074       {
0075       case 0:
0076     particle = electron;
0077     break;
0078       case 1:
0079     particle = pion;
0080     break;
0081       case 2:
0082     particle = kaon;
0083     break;
0084       case 3:
0085     particle = proton;
0086     break;
0087       default:
0088     cout << "Invalid Particle Selection - Default to Pion" << endl;
0089     particle = pion;
0090     break;
0091       }
0092
0093     //auto particle = pion;

0094     //if(part == 0) particle = electron;

0095     //if(part == 1) particle = pion;

0096     //if(part == 2) particle = kaon;

0097     //if(part == 3) particle = proton;

0098
0099     GenParticlePtr pq = std::make_shared<GenParticle>(FourVector(
0100                                  px, py, pz,
0101                                  sqrt(p*p + pow(particle->Mass(), 2))),
0102                               particle->PdgCode(), 1);
0103     v1->add_particle_out(pq);
0104
0105     evt.add_vertex(v1);
0106
0107     if (events_parsed == 0) {
0108       std::cout << "First event: " << std::endl;
0109       Print::listing(evt);
0110     }
0111
0112     hepmc_output.write_event(evt);
0113     if (events_parsed % 10000 == 0) {
0114       std::cout << "Event: " << events_parsed << std::endl;
0115     }
0116     evt.clear();
0117   }
0118   hepmc_output.close();
0119   std::cout << "Events parsed and written: " << events_parsed << std::endl;
0120   exit(0);
0121 }