HepMC3 event record library
ReaderLHEF.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // This file is part of HepMC
4 // Copyright (C) 2014-2021 The HepMC collaboration (see AUTHORS for details)
5 //
6 /**
7  * @file ReaderLHEF.cc
8  * @brief Implementation of \b class ReaderLHEF
9  *
10  */
11 #include "HepMC3/ReaderLHEF.h"
12 namespace HepMC3
13 {
14 ReaderLHEF::ReaderLHEF(const std::string& filename)
15 {
16  m_reader = std::make_shared<LHEF::Reader>(filename);
17  init();
18 }
19 ReaderLHEF::ReaderLHEF(std::istream & stream)
20 {
21  m_reader = std::make_shared<LHEF::Reader>(stream);
22  init();
23 }
24 
25 ReaderLHEF::ReaderLHEF(std::shared_ptr<std::istream> s_stream)
26  : m_shared_stream(s_stream)
27 {
28  m_reader = std::make_shared<LHEF::Reader>(*(m_shared_stream.get()));
29  init();
30 }
31 
32 bool ReaderLHEF::skip(const int n)
33 {
34  GenEvent evt;
35  for (int nn = n; nn > 0; --nn)
36  {
37  if (!read_event(evt)) return false;
38  evt.clear();
39  }
40  return !failed();
41 }
42 
43 
45 {
46  m_neve = 0;
47  m_failed = false;
48  // Create a HEPRUP attribute and initialize it from the reader.
49  m_hepr = std::make_shared<HEPRUPAttribute>();
50  m_hepr->heprup = m_reader->heprup;
51  // There may be some XML tags in the LHE file which are
52  // non-standard, but we can save them as well.
53  m_hepr->tags = LHEF::XMLTag::findXMLTags(m_reader->headerBlock + m_reader->initComments);
54  // This code is ugly and should be replaced.
55  size_t nweights = 0;
56  for (auto t1: m_hepr->tags) {
57  if (t1->name != "header") continue;
58  for (auto t2: t1->tags) {
59  if (t2->name != "initrwgt") continue;
60  for (auto t3: t2->tags) {
61  if (t3->name != "weightgroup") continue;
62  for (auto t4: t3->tags) if (t4->name == "weight") nweights++;
63  break;
64  }
65  break;
66  }
67  break;
68  }
69  //
70  // Now we want to create a GenRunInfo object for the HepMC file, and
71  // we add the LHEF attribute to that.
72  set_run_info(std::make_shared<GenRunInfo>());
73  run_info()->add_attribute("HEPRUP", m_hepr);
74 
75  // This is just a test to make sure we can add other attributes as
76  // well.
77  run_info()->add_attribute("NPRUP",
78  std::make_shared<FloatAttribute>(m_hepr->heprup.NPRUP));
79 
80  // We want to be able to convey the different event weights to
81  // HepMC. In particular we need to add the names of the weights to
82  // the GenRunInfo object.
83 
84  std::vector<std::string> weightnames;
85  for ( int i = 0, N = m_hepr->heprup.weightinfo.size(); i < N; ++i ) weightnames.push_back(m_hepr->heprup.weightNameHepMC(i));
86  if (nweights == 0) nweights=1;
87  for ( size_t i = weightnames.size(); i < nweights; ++i ) weightnames.push_back(std::to_string(i));
88  run_info()->set_weight_names(weightnames);
89 
90  // We also want to convey the information about which generators was
91  // used.
92  for ( int i = 0, N = m_hepr->heprup.generators.size(); i < N; ++i )
93  {
95  tool.name = m_hepr->heprup.generators[i].name;
96  tool.version = m_hepr->heprup.generators[i].version;
97  tool.description = m_hepr->heprup.generators[i].contents;
98  run_info()->tools().push_back(tool);
99  }
100 }
101 /// @brief Destructor
103 
105 {
106  if (m_storage.size() > 0)
107  {
108  ev = m_storage.front();
109  m_storage.pop_front();
110  return m_failed;
111  }
112  m_failed = !(m_reader->readEvent());
113  if (m_failed) return m_failed;
114  // To each GenEvent we want to add an attribute corresponding to
115  // the HEPEUP. Also here there may be additional non-standard
116  // information outside the LHEF <event> tags, which we may want to
117  // add.
118  std::shared_ptr<HEPEUPAttribute> hepe = std::make_shared<HEPEUPAttribute>();
119  if ( m_reader->outsideBlock.length() )
120  hepe->tags = LHEF::XMLTag::findXMLTags(m_reader->outsideBlock);
121 
122  hepe->hepeup = m_reader->hepeup;
123  std::vector<LHEF::HEPEUP*> input;
124  if (m_reader->hepeup.subevents.size() > 0) input.insert(input.end(), hepe->hepeup.subevents.begin(), hepe->hepeup.subevents.end());
125  else { input.push_back(&m_reader->hepeup);}
126  int first_group_event = m_neve;
127  m_neve++;
128  for (auto ahepeup: input)
129  {
130  GenEvent evt;
131  evt.set_event_number(first_group_event);
132  evt.add_attribute("AlphaQCD", std::make_shared<DoubleAttribute>(ahepeup->AQCDUP));
133  evt.add_attribute("AlphaEM", std::make_shared<DoubleAttribute>(ahepeup->AQEDUP));
134  evt.add_attribute("NUP", std::make_shared<IntAttribute>(ahepeup->NUP));
135  evt.add_attribute("IDPRUP", std::make_shared<LongAttribute>(ahepeup->IDPRUP));
136  // Now add the Particles from the LHE event to HepMC
137  std::vector<GenParticlePtr> particles;
138  std::map< std::pair<int, int>, GenVertexPtr> vertices;
139  for ( int i = 0; i < ahepeup->NUP; ++i )
140  {
141  FourVector mom((ahepeup->PUP)[i][0], (ahepeup->PUP)[i][1], (ahepeup->PUP)[i][2], (ahepeup->PUP)[i][3]);
142  particles.push_back(std::make_shared<GenParticle>(mom, ahepeup->IDUP[i], ahepeup->ISTUP[i]));
143  if ( i < 2 ) continue;
144  std::pair<int, int> vertex_index(ahepeup->MOTHUP[i].first, ahepeup->MOTHUP[i].second);
145  if (vertices.find(vertex_index) == vertices.end()) vertices[vertex_index] = std::make_shared<GenVertex>();
146  vertices[vertex_index]->add_particle_out(particles.back());
147  }
148  for ( auto v: vertices )
149  {
150  std::pair<int, int> vertex_index = v.first;
151  GenVertexPtr vertex = v.second;
152  for (int i = vertex_index.first-1; i < vertex_index.second; ++i)
153  if ( i >= 0 && i < (int)particles.size())
154  vertex->add_particle_in(particles[i]);
155  }
156  std::pair<int, int> vertex_index(0, 0);
157  if (vertices.find(vertex_index) == vertices.end()) vertices[vertex_index] = std::make_shared<GenVertex>();
158  for (size_t i = 0; i < particles.size(); ++i)
159  if (!particles[i]->end_vertex() && !particles[i]->production_vertex())
160  {
161  if ( i < 2 ) vertices[vertex_index]->add_particle_in(particles[i]);
162  else vertices[vertex_index]->add_particle_out(particles[i]);
163  }
164  for ( auto v: vertices ) evt.add_vertex(v.second);
165  if (particles.size() > 1)
166  {
167  particles[0]->set_status(4);
168  particles[1]->set_status(4);
169  evt.set_beam_particles(particles[0], particles[1]);
170  }
171  // And we also want to add the weights.
172 
173 
174  std::vector<double> wts;
175  for ( int i = 0, N = ahepeup->weights.size(); i < N; ++i )
176  {
177  wts.push_back(ahepeup->weights[i].first);
178  }
179  evt.weights() = wts;
180  m_storage.push_back(evt);
181  }
182  ev = m_storage.front();
183  m_storage.pop_front();
184  return m_failed;
185 }
186 /// @brief Return status of the stream
187 bool ReaderLHEF::failed() { return m_failed;}
188 
189 /// @brief Close file stream
191 } // namespace HepMC3
void init()
Init helper.
Definition: ReaderLHEF.cc:44
std::string version
The version of the tool.
Definition: GenRunInfo.h:44
void close() override
Close.
Definition: ReaderLHEF.cc:190
static std::vector< XMLTag * > findXMLTags(std::string str, std::string *leftover=0)
Definition: LHEF.h:198
void add_vertex(GenVertexPtr v)
Add vertex.
Definition: GenEvent.cc:96
bool read_event(GenEvent &ev) override
Reading event.
Definition: ReaderLHEF.cc:104
std::shared_ptr< LHEF::Reader > m_reader
The actual reader.
Definition: ReaderLHEF.h:56
int m_neve
Event counter.
Definition: ReaderLHEF.h:58
ReaderLHEF(std::istream &)
The ctor to read from stream.
Definition: ReaderLHEF.cc:19
~ReaderLHEF()
Destructor.
Definition: ReaderLHEF.cc:102
std::string description
Other information about how the tool was used in the run.
Definition: GenRunInfo.h:48
std::deque< GenEvent > m_storage
storage used for subevents.
Definition: ReaderLHEF.h:60
bool skip(const int) override
skip events
Definition: ReaderLHEF.cc:32
bool failed() override
State.
Definition: ReaderLHEF.cc:187
void add_attribute(const std::string &name, const std::shared_ptr< Attribute > &att, const int &id=0)
Definition: GenEvent.cc:806
std::shared_ptr< std::istream > m_shared_stream
Holds temporary stream.
Definition: ReaderLHEF.h:55
Stores event-related information.
Definition: GenEvent.h:41
Generic 4-vector.
Definition: FourVector.h:36
void set_beam_particles(GenParticlePtr p1, GenParticlePtr p2)
Set incoming beam particles.
Definition: GenEvent.cc:760
std::shared_ptr< GenRunInfo > run_info() const
Get the global GenRunInfo object.
Definition: Reader.h:44
Interrnal struct for keeping track of tools.
Definition: GenRunInfo.h:38
void set_run_info(std::shared_ptr< GenRunInfo > run)
Set the global GenRunInfo object.
Definition: Reader.h:64
const std::vector< double > & weights() const
Get event weight values as a vector.
Definition: GenEvent.h:98
std::string name
The name of the tool.
Definition: GenRunInfo.h:41
bool m_failed
State of reader.
Definition: ReaderLHEF.h:59
void clear()
Remove contents of this event.
Definition: GenEvent.cc:599
Definition of class ReaderLHEF.
void set_event_number(const int &num)
Set event number.
Definition: GenEvent.h:150
std::shared_ptr< HEPRUPAttribute > m_hepr
Holder of attributes.
Definition: ReaderLHEF.h:57