irt/scripts/pfrich-hepmc-writer.C

0001 //

0002 // root -l 'rich-hepmc-writer.cxx("out.hepmc", 100)'

0003 //

0004
0005 #include "HepMC3/GenEvent.h"
0006 #include "HepMC3/ReaderAscii.h"
0007 #include "HepMC3/WriterAscii.h"
0008 #include "HepMC3/Print.h"
0009
0010 #include <iostream>
0011 #include <random>
0012 #include <cmath>
0013 #include <math.h>
0014 #include <TMath.h>
0015
0016 using namespace HepMC3;
0017
0018 /** Generate single muon event with fixed three momentum **/
0019 void pfrich_hepmc_writer(const char* out_fname, int n_events)
0020 {
0021   auto *DatabasePDG = new TDatabasePDG();
0022   int pdg = 211;
0023   auto *particle = DatabasePDG->GetParticle(pdg);
0024
0025   WriterAscii hepmc_output(out_fname);
0026   int events_parsed = 0;
0027   GenEvent evt(Units::GEV, Units::MM);
0028
0029   //std::random_device rd;

0030   unsigned int seed = 0x12345678;//(unsigned int)abs(rd());

0031   std::cout << "init seed for random generator is " << seed << std::endl;
0032   // Random number generator

0033   TRandom *rdmn_gen = new TRandom(seed);
0034
0035   for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
0036     //FourVector(px,py,pz,e,pdgid,status)

0037     // type 4 is beam

0038     // pdgid 2212 - proton

0039     GenParticlePtr p1 =
0040         std::make_shared<GenParticle>(FourVector(0.0, 0.0, 12.0, 12.0), 11, 4);
0041     GenParticlePtr p2 = std::make_shared<GenParticle>(
0042         FourVector(0.0, 0.0, 100.0, 100.004), 2212, 4);
0043
0044     GenVertexPtr v1 = std::make_shared<GenVertex>();//FourVector(0,0,30,0));

0045     v1->add_particle_in(p1);
0046     v1->add_particle_in(p2);
0047
0048     // type 1 is final state; 211: pion; FIXME: give a proper mass;

0049     for(int iq=0; iq</*2*/1; iq++){
0050       Double_t eta   = rdmn_gen->Uniform(-2.0, -1.9);
0051       //Double_t eta   = rdmn_gen->Uniform(1.9, 2.0);//-2.0, -1.9);

0052       //Double_t eta   = rdmn_gen->Uniform(1.5, 1.6);//2.9, 3.0);//2.0, 2.1);

0053       Double_t th    = 2*std::atan(exp(-eta));
0054       //Double_t th    = rdmn_gen->Uniform(180.0-16.0, 180.0-14.0)*M_PI/180;

0055       Double_t p     = rdmn_gen->Uniform(10.0, 10.1);
0056       Double_t phi   = rdmn_gen->Uniform(0.0, 2*M_PI);
0057       //Double_t phi   = rdmn_gen->Uniform(-20.0, 20.0)*M_PI/180;

0058
0059       Double_t px    = p * std::cos(phi) * std::sin(th);
0060       Double_t py    = p * std::sin(phi) * std::sin(th);
0061       Double_t pz    = p * std::cos(th);
0062
0063       //cout<<"px,py,pz: "<<px<<" "<<py<<" "<<pz<<endl;

0064       GenParticlePtr pq = std::make_shared<GenParticle>(FourVector(
0065                                    px, py, pz,
0066                                    sqrt(p*p + pow(particle->Mass(), 2))),
0067                             pdg, 1);
0068       v1->add_particle_out(pq);
0069     }//iq

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   exit(0);
0086 }