gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
getHepMCEvent.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
36 
37 #include "gambit/cmake/cmake_variables.hpp"
38 
39 #ifndef EXCLUDE_HEPMC
40 
43 #include "HepMC3/ReaderAsciiHepMC2.h"
45 #include "HepMC3/GenEvent.h"
46 #include "HepMC3/GenParticle.h"
47 #include "HepMC3/ReaderAscii.h"
48 
49 #define DEBUG_PREFIX "DEBUG: OMP thread " << omp_get_thread_num() << ": "
50 //#define COLLIDERBIT_DEBUG
51 
52 namespace Gambit
53 {
54 
55  namespace ColliderBit
56  {
57 
59  void getHepMCEvent(HEPUtils::Event& result)
60  {
61  using namespace Pipes::getHepMCEvent;
62 
63  result.clear();
64 
65  // Get yaml options and initialise the HepMC reader
66  const static double antiktR = runOptions->getValueOrDef<double>(0.4, "antiktR");
67  const static double jet_pt_min = runOptions->getValueOrDef<double>(10.0, "jet_pt_min");
68  const static str HepMC_filename = runOptions->getValueOrDef<str>("", "hepmc_filename");
69  static int HepMC_file_version = -1;
70 
71  static bool first = true;
72  if (first)
73  {
74  if (not Utils::file_exists(HepMC_filename)) throw std::runtime_error("HepMC event file "+HepMC_filename+" not found. Quitting...");
75 
76  // Figure out if the file is HepMC2 or HepMC3
77  std::ifstream infile(HepMC_filename);
78  if (infile.good())
79  {
80  std::string line;
81  while(std::getline(infile, line))
82  {
83  // Skip blank lines
84  if(line == "") continue;
85 
86  // We look for "HepMC::Version 2" or "HepMC::Version 3",
87  // so we only need the first 16 characters of the line
88  std::string short_line = line.substr(0,16);
89 
90  if (short_line == "HepMC::Version 2")
91  {
92  HepMC_file_version = 2;
93  break;
94  }
95  else if (short_line == "HepMC::Version 3")
96  {
97  // Check the text format
98  std::getline(infile, line);
99  std::string text_format = line.substr(0,14);
100  if (text_format == "HepMC::Asciiv3")
101  {
102  HepMC_file_version = 3;
103  break;
104  }else if (text_format == "HepMC::IO_GenE"){
105  HepMC_file_version = 2;
106  break;
107  } else
108  {
109  throw std::runtime_error("Could not determine HepMC version from the string '"+text_format+"' extracted from the line '"+line+"'. Quitting...");
110  }
111  }
112  else
113  {
114  throw std::runtime_error("Could not determine HepMC version from the string '"+short_line+"' extracted from the line '"+line+"'. Quitting...");
115  }
116  }
117  }
118  first = false;
119  }
120 
121  if(HepMC_file_version != 2 and HepMC_file_version != 3)
122  {
123  throw std::runtime_error("Failed to determine HepMC version for input file "+HepMC_filename+". Quitting...");
124  }
125 
126  static HepMC3::Reader *HepMCio;
127 
128  // Initialize the reader on the first iteration
129  if (*Loop::iteration == BASE_INIT)
130  {
131  if (HepMC_file_version == 2)
132  HepMCio = new HepMC3::ReaderAsciiHepMC2(HepMC_filename);
133  else
134  HepMCio = new HepMC3::ReaderAscii(HepMC_filename);
135  }
136 
137  // Delete the reader in the last iteration
138  if (*Loop::iteration == BASE_FINALIZE)
139  delete HepMCio;
140 
141  // Don't do anything else during special iterations
142  if (*Loop::iteration < 0) return;
143 
144  #ifdef COLLIDERBIT_DEBUG
145  cout << DEBUG_PREFIX << "Event number: " << *Loop::iteration << endl;
146  #endif
147 
148  // Attempt to read the next HepMC event as a HEPUtils event. If there are no more events, wrap up the loop and skip the rest of this iteration.
149  HepMC3::GenEvent ge(HepMC3::Units::GEV, HepMC3::Units::MM);
150  bool event_retrieved = true;
151  #pragma omp critical (reading_HepMCEvent)
152  {
153  event_retrieved = HepMCio->read_event(ge);
154 
155  // FIXME This is a temp solution to ensure that the event reading
156  // stops when there are no more events in the HepMC file.
157  // Remove this once bugfix is implemented in HepMC.
158  if ((ge.particles().size() == 0) && (ge.vertices().size() == 0)) event_retrieved = false;
159  }
160  if (not event_retrieved)
161  {
162  // Tell the MCLoopInfo instance that we have reached the end of the file
163  Dep::RunMC->report_end_of_event_file();
164  Loop::halt();
165  }
166 
167  //Set the weight
168  result.set_weight(ge.weight());
169 
170  //Translate to HEPUtils event by calling the unified HEPMC/Pythia event converter:
171  Gambit::ColliderBit::convertParticleEvent(ge.particles(), result, antiktR, jet_pt_min);
172  }
173  }
174 
175 }
176 
177 #endif
Declarations common to all ColliderBit event loop functions.
General small utility functions.
EXPORT_SYMBOLS bool file_exists(const std::string &filename)
Check if a file exists.
#define DEBUG_PREFIX
std::string infile("FlavBit/data/example.slha")
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
void convertParticleEvent(const EventT &pevt, HEPUtils::Event &result, double antiktR, double jet_pt_min)
Convert a hadron-level EventT into an unsmeared HEPUtils::Event.
TODO: see if we can use this one:
Definition: Analysis.hpp:33
Helper functions for converting between Pythia8 events and other event types.
void getHepMCEvent(HEPUtils::Event &result)
A nested function that reads in HepMC event files and converts them to HEPUtils::Event format...