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