irt/scripts/drich-hepmc-loop-writer.C

0001 //

0002 //

0003 //

0004
0005 #define _AEROGEL_
0006
0007 #include "HepMC3/GenEvent.h"
0008 #include "HepMC3/ReaderAscii.h"
0009 #include "HepMC3/WriterAscii.h"
0010 #include "HepMC3/Print.h"
0011
0012 #include <iostream>
0013 #include <random>
0014 #include <cmath>
0015 #include <math.h>
0016 #include <TMath.h>
0017
0018 using namespace HepMC3;
0019
0020 void drich_hepmc_loop_writer(int n_events)
0021 {
0022   auto *DatabasePDG = new TDatabasePDG();
0023
0024   int pdg[] = {211, 321, 2212};
0025   unsigned gdim = sizeof(pdg) / sizeof(pdg[0]);
0026
0027   double eta[] = {1.2, 1.6, 2.0, 2.4, 2.8, 3.0, 3.2, 3.4, 3.5, 3.6};
0028   unsigned edim = sizeof(eta) / sizeof(eta[0]);
0029
0030 #ifdef _AEROGEL_
0031   double mom[] = {0.9, 1.4, 2.9,  4.2,  5.5, 10.0, 15.0, 20.0, 27.0, 40.0};
0032 #else
0033   double mom[] = {4.0, 5.0, 7.0, 14.0, 17.0, 21.0, 27.0, 30.0, 32.0, 35.0, 40.0, 45.0, 50.0, 55.0, 60.0};
0034 #endif
0035   unsigned mdim = sizeof(mom) / sizeof(mom[0]);
0036
0037   char fname[128-1];
0038 #ifdef _AEROGEL_
0039   const char *format = "HEPMC.A/drich-data.%04d..%04.2f-%04.2f..%04.1f-%04.1f..hepmc";
0040 #else
0041   const char *format = "HEPMC.G/drich-data.%04d..%04.2f-%04.2f..%04.1f-%04.1f..hepmc";
0042 #endif
0043
0044   unsigned int seed = 0x12345678;
0045   std::cout << "init seed for random generator is " << seed << std::endl;
0046
0047   TRandom *rdmn_gen = new TRandom(seed);
0048
0049   for(unsigned ig=0; ig<gdim; ig++) {
0050     auto *particle = DatabasePDG->GetParticle(pdg[ig]);
0051
0052     for(unsigned ie=0; ie<edim-1; ie++)
0053       for(unsigned im=0; im<mdim-1; im++) {
0054 #ifdef _AEROGEL_
0055     if (pdg[ig] == 321  && mom[im] <  2.9) continue;
0056     if (pdg[ig] == 2212 && mom[im] <  5.5) continue;
0057 #else
0058     if (pdg[ig] == 321  && mom[im] < 14.0) continue;
0059     if (pdg[ig] == 2212 && mom[im] < 27.0) continue;
0060 #endif
0061
0062     snprintf(fname, 128-1, format, pdg[ig], eta[ie], eta[ie+1], mom[im], mom[im+1]);
0063
0064     WriterAscii hepmc_output(fname);
0065     int events_parsed = 0;
0066     GenEvent evt(Units::GEV, Units::MM);
0067
0068     for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
0069       GenParticlePtr p1 =
0070         std::make_shared<GenParticle>(FourVector(0.0, 0.0,  12.0,  12.0),     11, 4);
0071       GenParticlePtr p2 =
0072         std::make_shared<GenParticle>(FourVector(0.0, 0.0, 100.0, 100.004), 2212, 4);
0073
0074       GenVertexPtr v1 = std::make_shared<GenVertex>();
0075       v1->add_particle_in(p1);
0076       v1->add_particle_in(p2);
0077
0078       // type 1 is final state;

0079       Double_t qeta  = rdmn_gen->Uniform(eta[ie], eta[ie+1]);
0080       Double_t th    = 2*std::atan(exp(-qeta));
0081       Double_t p     = rdmn_gen->Uniform(mom[im], mom[im+1]);
0082       Double_t phi   = rdmn_gen->Uniform(0.0, 2*M_PI);
0083
0084       Double_t px    = p * std::cos(phi) * std::sin(th);
0085       Double_t py    = p * std::sin(phi) * std::sin(th);
0086       Double_t pz    = p * std::cos(th);
0087
0088       GenParticlePtr pq = std::make_shared<GenParticle>(FourVector(
0089                                        px, py, pz,
0090                                        sqrt(p*p + pow(particle->Mass(), 2))),
0091                                 pdg[ig], 1);
0092       v1->add_particle_out(pq);
0093
0094       evt.add_vertex(v1);
0095
0096       if (events_parsed == 0) {
0097         std::cout << "First event: " << std::endl;
0098         Print::listing(evt);
0099       }
0100
0101       hepmc_output.write_event(evt);
0102       if (events_parsed % 10000 == 0) {
0103         std::cout << "Event: " << events_parsed << std::endl;
0104       }
0105       evt.clear();
0106     } //for ev

0107
0108     hepmc_output.close();
0109     std::cout << "Events parsed and written: " << events_parsed << std::endl;
0110       } //for ie..im

0111   } //for ig

0112
0113   exit(0);
0114 }