obsolete/ak/example_hepmc_reader.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 example_hepmc_reader.C+("out.hepmc", 100)

0008 //

0009
0010 #include "HepMC3/GenEvent.h"
0011 #include "HepMC3/Data/GenEventData.h"
0012 #include "HepMC3/ReaderAscii.h"
0013 //#include "HepMC3/WriterAscii.h"

0014 //#include "HepMC3/Print.h"

0015
0016 //#include <iostream>

0017 //#include <random>

0018 //#include <cmath>

0019 //#include <math.h>

0020
0021 //#include <TMath.h>

0022 //#include <TRandom.h>

0023 //#include <TDatabasePDG.h>

0024
0025 using namespace HepMC3;
0026
0027 void example_hepmc_reader(const char* in_fname, int n_events)
0028 {
0029   ReaderAscii hepmc_input(in_fname);
0030
0031   int events_parsed = 0;
0032
0033   while(!hepmc_input.failed()) {
0034     GenEvent evt(Units::GEV,Units::MM);
0035
0036     // Read event from input file

0037
0038     hepmc_input.read_event(evt);
0039
0040     // If reading failed - exit loop

0041
0042     if( hepmc_input.failed() ) break;
0043
0044     {
0045       printf("%d\n", evt.particles_size());
0046
0047       GenEventData data;
0048       evt.write_data(data);
0049       for(auto const &part: data.particles) {
0050     printf("   %8d %8d -> %f %f %f\n", part.pid, part.status, part.momentum.x(), part.momentum.y(), part.momentum.z());
0051       } //for part

0052       //printf("%d %d %d\n", data.particles.size(), data.links1.size(), data.links2.size());

0053     }
0054
0055     ++events_parsed;
0056     if (events_parsed == n_events) break;
0057   } //while

0058
0059   hepmc_input.close();
0060
0061 #if _TODAY_
0062   auto *DatabasePDG = new TDatabasePDG();
0063   auto *pion = DatabasePDG->GetParticle(211), *kaon = DatabasePDG->GetParticle(321);
0064
0065   WriterAscii hepmc_output(out_fname);
0066   int events_parsed = 0;
0067   GenEvent evt(Units::GEV, Units::MM);
0068
0069   // Random number generator

0070   TRandom *rdmn_gen = new TRandom(0x12345678);
0071
0072   for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
0073     // type 4 is beam

0074     GenParticlePtr p1 =
0075         std::make_shared<GenParticle>(FourVector(0.0, 0.0, 12.0, 12.0), 11, 4);
0076     GenParticlePtr p2 = std::make_shared<GenParticle>(
0077         FourVector(0.0, 0.0, 100.0, 100.004), 2212, 4);
0078
0079     GenVertexPtr v1 = std::make_shared<GenVertex>();
0080     v1->add_particle_in(p1);
0081     v1->add_particle_in(p2);
0082
0083     {
0084       Double_t eta   = rdmn_gen->Uniform(2.50, 2.51);
0085       Double_t th    = 2*std::atan(exp(-eta));
0086       Double_t p     = rdmn_gen->Uniform(7.0, 7.1);
0087       Double_t phi   = rdmn_gen->Uniform(0.0, 2*M_PI);
0088
0089       Double_t px    = p * std::cos(phi) * std::sin(th);
0090       Double_t py    = p * std::sin(phi) * std::sin(th);
0091       Double_t pz    = p * std::cos(th);
0092
0093       auto particle = events_parsed%2 ? pion : kaon;
0094       GenParticlePtr pq = std::make_shared<GenParticle>(FourVector(
0095                                    px, py, pz,
0096                                    sqrt(p*p + pow(particle->Mass(), 2))),
0097                             particle->PdgCode(), 1);
0098       v1->add_particle_out(pq);
0099     }
0100
0101     evt.add_vertex(v1);
0102
0103     if (events_parsed == 0) {
0104       std::cout << "First event: " << std::endl;
0105       Print::listing(evt);
0106     }
0107
0108     hepmc_output.write_event(evt);
0109     if (events_parsed % 10000 == 0) {
0110       std::cout << "Event: " << events_parsed << std::endl;
0111     }
0112     evt.clear();
0113   }
0114   hepmc_output.close();
0115   std::cout << "Events parsed and written: " << events_parsed << std::endl;
0116 #endif
0117   exit(0);
0118 }