HepMC3 event record library
LHEF_example_cat.cc
1 /**
2  * @example LHEF_example_cat.cc
3  * @brief Basic example of use of LHEF for reading and writing LHE files
4  */
6 #include "HepMC3/ReaderAscii.h"
7 #include "HepMC3/WriterAscii.h"
8 #include "HepMC3/GenEvent.h"
9 #include "HepMC3/GenParticle.h"
10 #include "HepMC3/GenVertex.h"
12 #include <iomanip>
13 
14 using namespace HepMC3;
15 
16 int main(int /*argc*/, char ** /*argv*/) {
17 
18  // Create Reader to read the example LHE file.
19  LHEF::Reader reader("LHEF_example.lhe");
20 
21  // Create a HEPRUP attribute and initialize it from the reader.
22  std::shared_ptr<HEPRUPAttribute> hepr = std::make_shared<HEPRUPAttribute>();
23  hepr->heprup = reader.heprup;
24 
25  // There may be some XML tags in the LHE file which are
26  // non-standard, but we can save them as well.
27  hepr->tags = LHEF::XMLTag::findXMLTags(reader.headerBlock + reader.initComments);
28 
29  // Nowwe want to create a GenRunInfo object for the HepMC file, and
30  // we add the LHEF attribute to that.
31  std::shared_ptr<GenRunInfo> runinfo = std::make_shared<GenRunInfo>();
32  runinfo->add_attribute("HEPRUP", hepr);
33 
34  // This is just a test to make sure we can add other attributes as
35  // well.
36  runinfo->add_attribute("NPRUP",
37  std::make_shared<FloatAttribute>(hepr->heprup.NPRUP));
38 
39  // We want to be able to convey the different event weights to
40  // HepMC. In particular we need to add the names of the weights to
41  // the GenRunInfo object.
42  std::vector<std::string> weightnames;
43  weightnames.push_back("0"); // The first weight is always the
44  // default weight with name "0".
45  for ( int i = 0, N = hepr->heprup.weightinfo.size(); i < N; ++i )
46  weightnames.push_back(hepr->heprup.weightNameHepMC(i));
47  runinfo->set_weight_names(weightnames);
48 
49  // We also want to convey the information about which generators was
50  // used to HepMC.
51  for ( int i = 0, N = hepr->heprup.generators.size(); i < N; ++i ) {
53  tool.name = hepr->heprup.generators[i].name;
54  tool.version = hepr->heprup.generators[i].version;
55  tool.description = hepr->heprup.generators[i].contents;
56  runinfo->tools().push_back(tool);
57  }
58 
59  // Now we want to start reading events from the LHE file and
60  // translate them to HepMC.
61  WriterAscii output("LHEF_example.hepmc3", runinfo);
62  int neve = 0;
63  while ( reader.readEvent() ) {
64  ++neve;
65 
66  // To each GenEvent we want to add an attribute corresponding to
67  // the HEPEUP. Also here there may be additional non-standard
68  // information outside the LHEF <event> tags, which we may want to
69  // add.
70  std::shared_ptr<HEPEUPAttribute> hepe = std::make_shared<HEPEUPAttribute>();
71  if ( reader.outsideBlock.length() )
72  hepe->tags = LHEF:: XMLTag::findXMLTags(reader.outsideBlock);
73  hepe->hepeup = reader.hepeup;
74  GenEvent ev(runinfo, Units::GEV, Units::MM);
75  ev.set_event_number(neve);
76 
77  // This is just a text to check that we can add additional
78  // attributes to each event.
79  ev.add_attribute("HEPEUP", hepe);
80  ev.add_attribute("AlphaQCD",
81  std:: make_shared<DoubleAttribute>(hepe->hepeup.AQCDUP));
82  ev.add_attribute("AlphaEM",
83  std::make_shared<DoubleAttribute>(hepe->hepeup.AQEDUP));
84  ev.add_attribute("NUP",
85  std::make_shared<IntAttribute>(hepe->hepeup.NUP));
86  ev.add_attribute("IDPRUP",
87  std::make_shared<LongAttribute>(hepe->hepeup.IDPRUP));
88 
89  // Now add the Particles from the LHE event to HepMC
90  std::vector<GenParticlePtr> particles;
91  std::map< std::pair<int,int>, GenVertexPtr> vertices;
92  for ( int i = 0; i < hepe->hepeup.NUP; ++i )
93  {
94  particles.push_back(std::make_shared<GenParticle>(hepe->momentum(i),hepe->hepeup.IDUP[i],hepe->hepeup.ISTUP[i]));
95  if (i<2) continue;
96  std::pair<int,int> vertex_index(hepe->hepeup.MOTHUP[i].first,hepe->hepeup.MOTHUP[i].second);
97  if (vertices.find(vertex_index)==vertices.end())vertices[vertex_index]=std::make_shared<GenVertex>();
98  vertices[vertex_index]->add_particle_out(particles.back());
99  }
100  for ( auto v: vertices )
101  {
102  std::pair<int,int> vertex_index=v.first;
103  GenVertexPtr vertex=v.second;
104  for (int i=vertex_index.first-1; i<vertex_index.second; i++) if (i>=0&&i<(int)particles.size()) vertex->add_particle_in(particles[i]);
105  }
106  for ( auto v: vertices ) ev.add_vertex(v.second);
107 
108  // And we also want to add the weights.
109  std::vector<double> wts;
110  for ( int i = 0, N = hepe->hepeup.weights.size(); i < N; ++i )
111  wts.push_back(hepe->hepeup.weights[i].first);
112  ev.weights() = wts;
113 
114  // Let's see if we can associate p1 and p2.
115  ev.add_attribute("OtherIncoming",
116  std::make_shared<AssociatedParticle>(particles[1]), particles[0]->id());
117 
118 
119  // And then we are done and can write out the GenEvent.
120  output.write_event(ev);
121 
122  }
123 
124  output.close();
125 
126  // Now we wnat to make sure we can read in the HepMC file and
127  // recreate the same info. To check this we try to reconstruct the
128  // LHC file we read in above.
129  ReaderAscii input("LHEF_example.hepmc3");
130  LHEF::Writer writer("LHEF_example_out.lhe");
131  hepr = std::shared_ptr<HEPRUPAttribute>();
132 
133  // The loop over all events in the HepMC file.
134  while ( true ) {
135 
136  // Read in the next event.
137  GenEvent ev(Units::GEV, Units::MM);
138  if ( !input.read_event(ev) || ev.event_number() == 0 ) break;
139 
140  // Check that the first incoming particle still refers to the second.
141  std::shared_ptr<AssociatedParticle> assoc =
142  ev.attribute<AssociatedParticle>("OtherIncoming", 1);
143  if ( !assoc || !assoc->associated() ||
144  assoc->associated() != ev.particles()[1] ) return 3;
145 
146  // Make sure the weight names are the same.
147  if ( input.run_info()->weight_names() != weightnames ) return 2;
148 
149  // For the first event we also go in and reconstruct the HEPRUP
150  // information, and write it out to the new LHE file.
151  if ( !hepr ) {
152  hepr = ev.attribute<HEPRUPAttribute>("HEPRUP");
153 
154  // Here we also keep track of the additional non-standard info
155  // we found in the original LHE file.
156  for ( int i = 0, N = hepr->tags.size(); i < N; ++i )
157  if ( hepr->tags[i]->name != "init" )
158  hepr->tags[i]->print(writer.headerBlock());
159 
160  // This is just a test that we can access other attributes
161  // included in the GenRunInfo.
162  hepr->heprup.NPRUP =
163  int(input.run_info()->
164  attribute<FloatAttribute>("NPRUP")->value());
165 
166  // Then we write out the HEPRUP object.
167  writer.heprup = hepr->heprup;
168  if ( writer.heprup.eventfiles.size() >= 2 ) {
169  writer.heprup.eventfiles[0].filename = "LHEF_example_1_out.plhe";
170  writer.heprup.eventfiles[1].filename = "LHEF_example_2_out.plhe";
171  }
172  writer.init();
173 
174  }
175 
176  // Now we can access the HEPEUP attribute of the current event.
177  std::shared_ptr<HEPEUPAttribute> hepe =
178  ev.attribute<HEPEUPAttribute>("HEPEUP");
179 
180  // Again, there may be addisional non-standard information we want
181  // to keep.
182  for ( int i = 0, N = hepe->tags.size(); i < N; ++i )
183  if ( hepe->tags[i]->name != "event" &&
184  hepe->tags[i]->name != "eventgroup" )
185  hepe->tags[i]->print(writer.eventComments());
186 
187  // This is just a test that we can access other attributes
188  // included in the GenRunInfo.
189  hepe->hepeup.AQCDUP =
190  ev.attribute<DoubleAttribute>("AlphaQCD")->value();
191  hepe->hepeup.AQEDUP =
192  ev.attribute<DoubleAttribute>("AlphaEM")->value();
193  hepe->hepeup.NUP =
194  ev.attribute<IntAttribute>("NUP")->value();
195  hepe->hepeup.IDPRUP =
196  ev.attribute<LongAttribute>("IDPRUP")->value();
197 
198  // And then we can write out the HEPEUP object.
199  writer.hepeup = hepe->hepeup;
200  writer.hepeup.heprup = &writer.heprup;
201  writer.writeEvent();
202 
203  }
204 
205 }
Definition of class AssociatedParticle,.
Definition of class GenEvent.
Definition of class GenParticle.
Definition of class GenVertex.
Definition of class HEPRUPAttribute and class HEPEUAttribute.
Definition of class ReaderAscii.
Definition of class WriterAscii.
Attribute class allowing eg. a GenParticle to refer to another GenParticle.
Attribute that holds a real number as a double.
Definition: Attribute.h:241
Stores event-related information.
Definition: GenEvent.h:41
Class for storing data for LHEF run information.
std::vector< LHEF::XMLTag * > tags
The parsed XML-tags.
Class for storing data for LHEF run information.
std::vector< LHEF::XMLTag * > tags
The parsed XML-tags.
Attribute that holds an Integer implemented as an int.
Definition: Attribute.h:157
Attribute that holds an Integer implemented as an int.
Definition: Attribute.h:198
GenEvent I/O parsing for structured text files.
Definition: ReaderAscii.h:29
GenEvent I/O serialization for structured text files.
Definition: WriterAscii.h:25
HepMC3 main namespace.
int main(int argc, char **argv)
Interrnal struct for keeping track of tools.
Definition: GenRunInfo.h:38
std::string description
Other information about how the tool was used in the run.
Definition: GenRunInfo.h:48
std::string version
The version of the tool.
Definition: GenRunInfo.h:44
std::string name
The name of the tool.
Definition: GenRunInfo.h:41
static std::vector< XMLTag * > findXMLTags(std::string str, std::string *leftover=0)
Definition: LHEF.h:198