HepMC3 event record library
LHEF.h
1 // -*- C++ -*-
2 #ifndef HEPMC3_LHEF_H
3 #define HEPMC3_LHEF_H
4 //
5 // This is the declaration of the Les Houches Event File classes,
6 // implementing a simple C++ parser/writer for Les Houches Event files.
7 // Copyright (C) 2009-2013 Leif Lonnblad
8 //
9 // The code is licenced under version 2 of the GPL, see COPYING for details.
10 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
11 //
12 
13 #include <iostream>
14 #include <iomanip>
15 #include <sstream>
16 #include <fstream>
17 #include <string>
18 #include <vector>
19 #include <map>
20 #include <set>
21 #include <utility>
22 #include <stdexcept>
23 #include <cstdlib>
24 #include <cmath>
25 #include <limits>
26 #ifndef M_PI
27 #define M_PI 3.14159265358979323846264338327950288
28 #endif
29 
30 
31 /**
32  * @brief Les Houches event file classes.
33  *
34  * The namespace containing helper classes and Reader and Writer
35  * classes for handling Les Houches event files.
36  *
37  * @ingroup LHEF
38  */
39 namespace LHEF {
40 
41 /**
42  * Helper struct for output of attributes.
43  */
44 template <typename T>
45 struct OAttr {
46 
47  /**
48  * Constructor
49  */
50  OAttr(std::string n, const T & v): name(n), val(v) {}
51 
52  /**
53  * The name of the attribute being printed.
54  */
55  std::string name;
56 
57  /**
58  * The value of the attribute being printed.
59  */
60  T val;
61 
62 };
63 
64 /**
65  * Output manipulator for writing attributes.
66  */
67 template <typename T>
68 OAttr<T> oattr(std::string name, const T & value) {
69  return OAttr<T>(name, value);
70 }
71 
72 /**
73  * Output operator for attributes.
74  */
75 template <typename T>
76 std::ostream & operator<<(std::ostream & os, const OAttr<T> & oa) {
77  os << " " << oa.name << "=\"" << oa.val << "\"";
78  return os;
79 }
80 
81 /**
82  * The XMLTag struct is used to represent all information within an
83  * XML tag. It contains the attributes as a map, any sub-tags as a
84  * vector of pointers to other XMLTag objects, and any other
85  * information as a single string.
86  */
87 struct XMLTag {
88 
89  /**
90  * Convenient typdef.
91  */
92  typedef std::string::size_type pos_t;
93 
94  /**
95  * Convenient typdef.
96  */
97  typedef std::map<std::string,std::string> AttributeMap;
98 
99  /**
100  * Convenient alias for npos.
101  */
102  static const pos_t end = std::string::npos;
103 
104  /**
105  * Default constructor.
106  */
107  XMLTag() {}
108 
109  /**
110  * The destructor also destroys any sub-tags.
111  */
113  for ( int i = 0, N = tags.size(); i < N; ++i ) delete tags[i];
114  }
115 
116  /**
117  * The name of this tag.
118  */
119  std::string name;
120 
121  /**
122  * The attributes of this tag.
123  */
125 
126  /**
127  * A vector of sub-tags.
128  */
129  std::vector<XMLTag*> tags;
130 
131  /**
132  * The contents of this tag.
133  */
134  std::string contents;
135 
136  /**
137  * Find an attribute named \a n and set the double variable \a v to
138  * the corresponding value. @return false if no attribute was found.
139  */
140  bool getattr(std::string n, double & v) const {
141  AttributeMap::const_iterator it = attr.find(n);
142  if ( it == attr.end() ) return false;
143  v = std::atof(it->second.c_str());
144  return true;
145  }
146 
147  /**
148  * Find an attribute named \a n and set the bool variable \a v to
149  * true if the corresponding value is "yes". @return false if no
150  * attribute was found.
151  */
152  bool getattr(std::string n, bool & v) const {
153  AttributeMap::const_iterator it = attr.find(n);
154  if ( it == attr.end() ) return false;
155  if ( it->second == "yes" ) v = true;
156  return true;
157  }
158 
159  /**
160  * Find an attribute named \a n and set the long variable \a v to
161  * the corresponding value. @return false if no attribute was found.
162  */
163  bool getattr(std::string n, long & v) const {
164  AttributeMap::const_iterator it = attr.find(n);
165  if ( it == attr.end() ) return false;
166  v = std::atoi(it->second.c_str());
167  return true;
168  }
169 
170  /**
171  * Find an attribute named \a n and set the long variable \a v to
172  * the corresponding value. @return false if no attribute was found.
173  */
174  bool getattr(std::string n, int & v) const {
175  AttributeMap::const_iterator it = attr.find(n);
176  if ( it == attr.end() ) return false;
177  v = int(std::atoi(it->second.c_str()));
178  return true;
179  }
180 
181  /**
182  * Find an attribute named \a n and set the string variable \a v to
183  * the corresponding value. @return false if no attribute was found.
184  */
185  bool getattr(std::string n, std::string & v) const {
186  AttributeMap::const_iterator it = attr.find(n);
187  if ( it == attr.end() ) return false;
188  v = it->second;
189  return true;
190  }
191 
192  /**
193  * Scan the given string and return all XML tags found as a vector
194  * of pointers to XMLTag objects. Text which does not belong to any
195  * tag is stored in tags without name and in the string pointed to
196  * by leftover (if not null).
197  */
198  static std::vector<XMLTag*> findXMLTags(std::string str,
199  std::string * leftover = 0) {
200  std::vector<XMLTag*> tags;
201  pos_t curr = 0;
202 
203  while ( curr != end ) {
204 
205  // Find the first tag
206  pos_t begin = str.find("<", curr);
207 
208  // Check for comments
209  if ( begin != end && str.find("<!--", curr) == begin ) {
210  pos_t endcom = str.find("-->", begin);
211  tags.push_back(new XMLTag());
212  if ( endcom == end ) {
213  tags.back()->contents = str.substr(curr);
214  if ( leftover ) *leftover += str.substr(curr);
215  return tags;
216  }
217  tags.back()->contents = str.substr(curr, endcom - curr);
218  if ( leftover ) *leftover += str.substr(curr, endcom - curr);
219  curr = endcom;
220  continue;
221  }
222 
223  if ( begin != curr ) {
224  tags.push_back(new XMLTag());
225  tags.back()->contents = str.substr(curr, begin - curr);
226  if ( leftover ) *leftover += str.substr(curr, begin - curr);
227  }
228  if ( begin == end || begin > str.length() - 3 || str[begin + 1] == '/' )
229  return tags;
230 
231  pos_t close = str.find(">", curr);
232  if ( close == end ) return tags;
233 
234  // find the tag name.
235  curr = str.find_first_of(" \t\n/>", begin);
236  tags.push_back(new XMLTag());
237  tags.back()->name = str.substr(begin + 1, curr - begin - 1);
238 
239  while ( true ) {
240 
241  // Now skip some white space to see if we can find an attribute.
242  curr = str.find_first_not_of(" \t\n", curr);
243  if ( curr == end || curr >= close ) break;
244 
245  pos_t tend = str.find_first_of("= \t\n", curr);
246  if ( tend == end || tend >= close ) break;
247 
248  std::string namex = str.substr(curr, tend - curr);
249  curr = str.find("=", curr) + 1;
250 
251  // OK now find the beginning and end of the atribute.
252  curr = str.find_first_of("\"'", curr);
253  if ( curr == end || curr >= close ) break;
254  char quote = str[curr];
255  pos_t bega = ++curr;
256  curr = str.find(quote, curr);
257  while ( curr != end && str[curr - 1] == '\\' )
258  curr = str.find(quote, curr + 1);
259 
260  std::string value = str.substr(bega, curr == end? end: curr - bega);
261 
262  tags.back()->attr[namex] = value;
263 
264  ++curr;
265 
266  }
267 
268  curr = close + 1;
269  if ( str[close - 1] == '/' ) continue;
270 
271  pos_t endtag = str.find("</" + tags.back()->name + ">", curr);
272  if ( endtag == end ) {
273  tags.back()->contents = str.substr(curr);
274  curr = endtag;
275  } else {
276  tags.back()->contents = str.substr(curr, endtag - curr);
277  curr = endtag + tags.back()->name.length() + 3;
278  }
279 
280  std::string leftovers;
281  tags.back()->tags = findXMLTags(tags.back()->contents, &leftovers);
282  if ( leftovers.find_first_not_of(" \t\n") == end ) leftovers="";
283  tags.back()->contents = leftovers;
284 
285  }
286  return tags;
287 
288  }
289 
290  /**
291  * Delete all tags in a vector.
292  */
293  static void deleteAll(std::vector<XMLTag*> & tags) {
294  while ( tags.size() && tags.back() ) {
295  delete tags.back();
296  tags.pop_back();
297  }
298  }
299  /**
300  * Print out this tag to a stream.
301  */
302  void print(std::ostream & os) const {
303  if ( name.empty() ) {
304  os << contents;
305  return;
306  }
307  os << "<" << name;
308  for ( AttributeMap::const_iterator it = attr.begin();
309  it != attr.end(); ++it )
310  os << oattr(it->first, it->second);
311  if ( contents.empty() && tags.empty() ) {
312  os << "/>" << std::endl;
313  return;
314  }
315  os << ">";
316  for ( int i = 0, N = tags.size(); i < N; ++i )
317  tags[i]->print(os);
318 
319  os << contents << "</" << name << ">" << std::endl;
320  }
321 
322 };
323 
324 /**
325  * Helper function to make sure that each line in the string \a s starts with a
326  * #-character and that the string ends with a new-line.
327  */
328 inline std::string hashline(std::string s) {
329  std::string ret;
330  std::istringstream is(s);
331  std::string ss;
332  while ( getline(is, ss) ) {
333  if ( ss.empty() ) continue;
334  if ( ss.find_first_not_of(" \t") == std::string::npos ) continue;
335  if ( ss.find('#') == std::string::npos ||
336  ss.find('#') != ss.find_first_not_of(" \t") ) ss = "# " + ss;
337  ret += ss + '\n';
338  }
339  return ret;
340 }
341 
342 /**
343  * This is the base class of all classes representing xml tags.
344  */
345 struct TagBase {
346 
347  /**
348  * Convenient typedef.
349  */
351 
352  /**
353  * Default constructor does nothing.
354  */
355  TagBase() {}
356 
357  /**
358  * Main constructor stores the attributes and contents of a tag.
359  */
360  TagBase(const AttributeMap & attr, std::string conts = std::string()): attributes(attr), contents(conts) {}
361 
362  /**
363  * Find an attribute named \a n and set the double variable \a v to
364  * the corresponding value. Remove the correspondig attribute from
365  * the list if found and \a erase is true. @return false if no
366  * attribute was found.
367  */
368  bool getattr(std::string n, double & v, bool erase = true) {
369  AttributeMap::iterator it = attributes.find(n);
370  if ( it == attributes.end() ) return false;
371  v = std::atof(it->second.c_str());
372  if ( erase) attributes.erase(it);
373  return true;
374  }
375 
376  /**
377  * Find an attribute named \a n and set the bool variable \a v to
378  * true if the corresponding value is "yes". Remove the correspondig
379  * attribute from the list if found and \a erase is true. @return
380  * false if no attribute was found.
381  */
382  bool getattr(std::string n, bool & v, bool erase = true) {
383  AttributeMap::iterator it = attributes.find(n);
384  if ( it == attributes.end() ) return false;
385  if ( it->second == "yes" ) v = true;
386  if ( erase) attributes.erase(it);
387  return true;
388  }
389 
390  /**
391  * Find an attribute named \a n and set the long variable \a v to
392  * the corresponding value. Remove the correspondig attribute from
393  * the list if found and \a erase is true. @return false if no
394  * attribute was found.
395  */
396  bool getattr(std::string n, long & v, bool erase = true) {
397  AttributeMap::iterator it = attributes.find(n);
398  if ( it == attributes.end() ) return false;
399  v = std::atoi(it->second.c_str());
400  if ( erase) attributes.erase(it);
401  return true;
402  }
403 
404  /**
405  * Find an attribute named \a n and set the long variable \a v to
406  * the corresponding value. Remove the correspondig attribute from
407  * the list if found and \a erase is true. @return false if no
408  * attribute was found.
409  */
410  bool getattr(std::string n, int & v, bool erase = true) {
411  AttributeMap::iterator it = attributes.find(n);
412  if ( it == attributes.end() ) return false;
413  v = int(std::atoi(it->second.c_str()));
414  if ( erase) attributes.erase(it);
415  return true;
416  }
417 
418  /**
419  * Find an attribute named \a n and set the string variable \a v to
420  * the corresponding value. Remove the correspondig attribute from
421  * the list if found and \a erase is true. @return false if no
422  * attribute was found.
423  */
424  bool getattr(std::string n, std::string & v, bool erase = true) {
425  AttributeMap::iterator it = attributes.find(n);
426  if ( it == attributes.end() ) return false;
427  v = it->second;
428  if ( erase) attributes.erase(it);
429  return true;
430  }
431 
432  /**
433  * print out ' name="value"' for all unparsed attributes.
434  */
435  void printattrs(std::ostream & file) const {
436  for ( AttributeMap::const_iterator it = attributes.begin();
437  it != attributes.end(); ++ it )
438  file << oattr(it->first, it->second);
439  }
440 
441  /**
442  * Print out end of tag marker. Print contents if not empty else
443  * print simple close tag.
444  */
445  void closetag(std::ostream & file, std::string tag) const {
446  if ( contents.empty() )
447  file << "/>\n";
448  else if ( contents.find('\n') != std::string::npos )
449  file << ">\n" << contents << "\n</" << tag << ">\n";
450  else
451  file << ">" << contents << "</" << tag << ">\n";
452  }
453 
454  /**
455  * The attributes of this tag;
456  */
458 
459  /**
460  * The contents of this tag.
461  */
462  mutable std::string contents;
463 
464  /**
465  * Static string token for truth values.
466  */
467  static std::string yes() { return "yes"; }
468 
469 };
470 
471 /**
472  * The Generator class contains information about a generator used in a run.
473  */
474 struct Generator : public TagBase {
475 
476  /**
477  * Construct from XML tag.
478  */
479  Generator(const XMLTag & tag)
480  : TagBase(tag.attr, tag.contents) {
481  getattr("name", name);
482  getattr("version", version);
483  }
484 
485  /**
486  * Print the object as a generator tag.
487  */
488  void print(std::ostream & file) const {
489  file << "<generator";
490  if ( !name.empty() ) file << oattr("name", name);
491  if ( !version.empty() ) file << oattr("version", version);
492  printattrs(file);
493  closetag(file, "generator");
494  }
495 
496  /**
497  * The name of the generator.
498  */
499  std::string name;
500 
501  /**
502  * The version of the generator.
503  */
504  std::string version;
505 
506 };
507 
508 /**
509  * The XSecInfo class contains information given in the xsecinfo tag.
510  */
511 struct XSecInfo : public TagBase {
512 
513  /**
514  * Intitialize default values.
515  */
516  XSecInfo(): neve(-1), ntries(-1), totxsec(0.0), xsecerr(0.0), maxweight(1.0),
517  meanweight(1.0), negweights(false), varweights(false) {}
518 
519  /**
520  * Create from XML tag
521  */
522  XSecInfo(const XMLTag & tag)
523  : TagBase(tag.attr, tag.contents), neve(-1), ntries(-1),
524  totxsec(0.0), xsecerr(0.0), maxweight(1.0), meanweight(1.0),
525  negweights(false), varweights(false) {
526  if ( !getattr("neve", neve) )
527  throw std::runtime_error("Found xsecinfo tag without neve attribute "
528  "in Les Houches Event File.");
529  ntries = neve;
530  getattr("ntries", ntries);
531  if ( !getattr("totxsec", totxsec) )
532  throw std::runtime_error("Found xsecinfo tag without totxsec "
533  "attribute in Les Houches Event File.");
534  getattr("xsecerr", xsecerr);
535  getattr("weightname", weightname);
536  getattr("maxweight", maxweight);
537  getattr("meanweight", meanweight);
538  getattr("negweights", negweights);
539  getattr("varweights", varweights);
540 
541  }
542 
543  /**
544  * Print out an XML tag.
545  */
546  void print(std::ostream & file) const {
547  file << "<xsecinfo" << oattr("neve", neve)
548  << oattr("totxsec", totxsec);
549  if ( maxweight != 1.0 )
550  file << oattr("maxweight", maxweight)
551  << oattr("meanweight", meanweight);
552  if ( ntries > neve ) file << oattr("ntries", ntries);
553  if ( xsecerr > 0.0 ) file << oattr("xsecerr", xsecerr);
554  if ( !weightname.empty() ) file << oattr("weightname", weightname);
555  if ( negweights ) file << oattr("negweights", yes());
556  if ( varweights ) file << oattr("varweights", yes());
557  printattrs(file);
558  closetag(file, "xsecinfo");
559  }
560 
561  /**
562  * The number of events.
563  */
564  long neve;
565 
566  /**
567  * The number of attempte that was needed to produce the neve events.
568  */
569  long ntries;
570 
571  /**
572  * The total cross section in pb.
573  */
574  double totxsec;
575 
576  /**
577  * The estimated statistical error on totxsec.
578  */
579  double xsecerr;
580 
581  /**
582  * The maximum weight.
583  */
584  double maxweight;
585 
586  /**
587  * The average weight.
588  */
589  double meanweight;
590 
591  /**
592  * Does the file contain negative weights?
593  */
595 
596  /**
597  * Does the file contain varying weights?
598  */
600 
601  /**
602  * The named weight to which this object belongs.
603  */
604  std::string weightname;
605 
606 };
607 
608 /**
609  * Convenient Alias.
610  */
611 typedef std::map<std::string,XSecInfo> XSecInfos;
612 
613 /**
614  * Simple struct to store information about separate eventfiles to be
615  * loaded.
616  */
617 struct EventFile : public TagBase {
618 
619  /**
620  * Intitialize default values.
621  */
622  EventFile(): filename(""), neve(-1), ntries(-1) {}
623 
624  /**
625  * Create from XML tag
626  */
627  EventFile(const XMLTag & tag)
628  : TagBase(tag.attr, tag.contents), filename(""), neve(-1), ntries(-1) {
629  if ( !getattr("name", filename) )
630  throw std::runtime_error("Found eventfile tag without name attribute "
631  "in Les Houches Event File.");
632  getattr("neve", neve);
633  ntries = neve;
634  getattr("ntries", ntries);
635  }
636 
637  /**
638  * Print out an XML tag.
639  */
640  void print(std::ostream & file) const {
641  if ( filename.empty() ) return;
642  file << " <eventfile" << oattr("name", filename);
643  if ( neve > 0 ) file << oattr("neve", neve);
644  if ( ntries > neve ) file << oattr("ntries", ntries);
645  printattrs(file);
646  closetag(file, "eventfile");
647  }
648 
649  /**
650  * The name of the file.
651  */
652  std::string filename;
653 
654  /**
655  * The number of events.
656  */
657  long neve;
658 
659  /**
660  * The number of attempte that was needed to produce the neve events.
661  */
662  long ntries;
663 
664 };
665 
666 /**
667  * The Cut class represents a cut used by the Matrix Element generator.
668  */
669 struct Cut : public TagBase {
670 
671  /**
672  * Intitialize default values.
673  */
674  Cut(): min(-0.99*std::numeric_limits<double>::max()),
675  max(0.99*std::numeric_limits<double>::max()) {}
676 
677  /**
678  * Create from XML tag.
679  */
680  Cut(const XMLTag & tag,
681  const std::map<std::string,std::set<long> >& ptypes)
682  : TagBase(tag.attr),
683  min(-0.99*std::numeric_limits<double>::max()),
684  max(0.99*std::numeric_limits<double>::max()) {
685  if ( !getattr("type", type) )
686  throw std::runtime_error("Found cut tag without type attribute "
687  "in Les Houches file");
688  long tmp;
689  if ( tag.getattr("p1", np1) ) {
690  if ( ptypes.find(np1) != ptypes.end() ) {
691  p1 = ptypes.find(np1)->second;
692  attributes.erase("p1");
693  } else {
694  getattr("p1", tmp);
695  p1.insert(tmp);
696  np1 = "";
697  }
698  }
699  if ( tag.getattr("p2", np2) ) {
700  if ( ptypes.find(np2) != ptypes.end() ) {
701  p2 = ptypes.find(np2)->second;
702  attributes.erase("p2");
703  } else {
704  getattr("p2", tmp);
705  p2.insert(tmp);
706  np2 = "";
707  }
708  }
709 
710  std::istringstream iss(tag.contents);
711  iss >> min;
712  if ( iss >> max ) {
713  if ( min >= max )
714  min = -0.99*std::numeric_limits<double>::max();
715  } else
716  max = 0.99*std::numeric_limits<double>::max();
717  }
718 
719  /**
720  * Print out an XML tag.
721  */
722  void print(std::ostream & file) const {
723  file << "<cut" << oattr("type", type);
724  if ( !np1.empty() )
725  file << oattr("p1", np1);
726  else
727  if ( p1.size() == 1 ) file << oattr("p1", *p1.begin());
728  if ( !np2.empty() )
729  file << oattr("p2", np2);
730  else
731  if ( p2.size() == 1 ) file << oattr("p2", *p2.begin());
732  printattrs(file);
733 
734  file << ">";
735  if ( min > -0.9*std::numeric_limits<double>::max() )
736  file << min;
737  else
738  file << max;
739  if ( max < 0.9*std::numeric_limits<double>::max() )
740  file << " " << max;
741  if ( !contents.empty() ) file << std::endl << contents << std::endl;
742  file << "</cut>" << std::endl;
743  }
744 
745  /**
746  * Check if a \a id1 matches p1 and \a id2 matches p2. Only non-zero
747  * values are considered.
748  */
749  bool match(long id1, long id2 = 0) const {
750  std::pair<bool,bool> ret(false, false);
751  if ( !id2 ) ret.second = true;
752  if ( !id1 ) ret.first = true;
753  if ( p1.find(0) != p1.end() ) ret.first = true;
754  if ( p1.find(id1) != p1.end() ) ret.first = true;
755  if ( p2.find(0) != p2.end() ) ret.second = true;
756  if ( p2.find(id2) != p2.end() ) ret.second = true;
757  return ret.first && ret.second;
758  }
759 
760  /**
761  * Check if the particles given as a vector of PDG \a id numbers,
762  * and a vector of vectors of momentum components, \a p, will pass
763  * the cut defined in this event.
764  */
765  bool passCuts(const std::vector<long> & id,
766  const std::vector< std::vector<double> >& p ) const {
767  if ( ( type == "m" && !p2.size() ) || type == "kt" || type == "eta" ||
768  type == "y" || type == "E" ) {
769  for ( int i = 0, N = id.size(); i < N; ++i )
770  if ( match(id[i]) ) {
771  if ( type == "m" ) {
772  double v = p[i][4]*p[i][4] - p[i][3]*p[i][3] - p[i][2]*p[i][2]
773  - p[i][1]*p[i][1];
774  v = v >= 0.0? std::sqrt(v): -std::sqrt(-v);
775  if ( outside(v) ) return false;
776  }
777  else if ( type == "kt" ) {
778  if ( outside(std::sqrt(p[i][2]*p[i][2] + p[i][1]*p[i][1])) )
779  return false;
780  }
781  else if ( type == "E" ) {
782  if ( outside(p[i][4]) ) return false;
783  }
784  else if ( type == "eta" ) {
785  if ( outside(eta(p[i])) ) return false;
786  }
787  else if ( type == "y" ) {
788  if ( outside(rap(p[i])) ) return false;
789  }
790  }
791  }
792  else if ( type == "m" || type == "deltaR" ) {
793  for ( int i = 1, N = id.size(); i < N; ++i )
794  for ( int j = 0; j < i; ++j )
795  if ( match(id[i], id[j]) || match(id[j], id[i]) ) {
796  if ( type == "m" ) {
797  double v = (p[i][4] + p[j][4])*(p[i][4] + p[j][4])
798  - (p[i][3] + p[j][3])*(p[i][3] + p[j][3])
799  - (p[i][2] + p[j][2])*(p[i][2] + p[j][2])
800  - (p[i][1] + p[j][1])*(p[i][1] + p[j][1]);
801  v = v >= 0.0? std::sqrt(v): -std::sqrt(-v);
802  if ( outside(v) ) return false;
803  }
804  else if ( type == "deltaR" ) {
805  if ( outside(deltaR(p[i], p[j])) ) return false;
806  }
807  }
808  }
809  else if ( type == "ETmiss" ) {
810  double x = 0.0;
811  double y = 0.0;
812  for ( int i = 0, N = id.size(); i < N; ++i )
813  if ( match(id[i]) && !match(0, id[i]) ) {
814  x += p[i][1];
815  y += p[i][2];
816  }
817  if ( outside(std::sqrt(x*x + y*y)) ) return false;
818  }
819  else if ( type == "HT" ) {
820  double pt = 0.0;
821  for ( int i = 0, N = id.size(); i < N; ++i )
822  if ( match(id[i]) && !match(0, id[i]) )
823  pt += std::sqrt(p[i][1]*p[i][1] + p[i][2]*p[i][2]);
824  if ( outside(pt) ) return false;
825  }
826  return true;
827  }
828 
829  /**
830  * Return the pseudorapidity of a particle with momentum \a p.
831  */
832  static double eta(const std::vector<double> & p) {
833  double pt2 = p[2]*p[2] + p[1]*p[1];
834  if ( pt2 != 0.0 ) {
835  double dum = std::sqrt(pt2 + p[3]*p[3]) + p[3];
836  if ( dum != 0.0 )
837  return std::log(dum/std::sqrt(pt2));
838  }
839  return p[3] < 0.0? -std::numeric_limits<double>::max():
840  std::numeric_limits<double>::max();
841  }
842 
843  /**
844  * Return the true rapidity of a particle with momentum \a p.
845  */
846  static double rap(const std::vector<double> & p) {
847  double pt2 = p[5]*p[5] + p[2]*p[2] + p[1]*p[1];
848  if ( pt2 != 0.0 ) {
849  double dum = std::sqrt(pt2 + p[3]*p[3]) + p[3];
850  if ( dum != 0.0 )
851  return std::log(dum/std::sqrt(pt2));
852  }
853  return p[3] < 0.0? -std::numeric_limits<double>::max():
854  std::numeric_limits<double>::max();
855  }
856 
857  /**
858  * Return the delta-R of a particle pair with momenta \a p1 and \a p2.
859  */
860  static double deltaR(const std::vector<double> & p1,
861  const std::vector<double> & p2) {
862  double deta = eta(p1) - eta(p2);
863  double dphi = std::atan2(p1[1], p1[2]) - std::atan2(p2[1], p2[2]);
864  if ( dphi > M_PI ) dphi -= 2.0*M_PI;
865  if ( dphi < -M_PI ) dphi += 2.0*M_PI;
866  return std::sqrt(dphi*dphi + deta*deta);
867  }
868 
869  /**
870  * Return true if the given \a value is outside limits.
871  */
872  bool outside(double value) const {
873  return value < min || value >= max;
874  }
875 
876  /**
877  * The variable in which to cut.
878  */
879  std::string type;
880 
881  /**
882  * The first types particle types for which this cut applies.
883  */
884  std::set<long> p1;
885 
886  /**
887  * Symbolic name for p1.
888  */
889  std::string np1;
890 
891  /**
892  * The second type of particles for which this cut applies.
893  */
894  std::set<long> p2;
895 
896  /**
897  * Symbolic name for p1.
898  */
899  std::string np2;
900 
901  /**
902  * The minimum value of the variable
903  */
904  double min;
905  /**
906  * The maximum value of the variable
907  */
908  double max;
909 
910 };
911 
912 /**
913  * The ProcInfo class represents the information in a procinfo tag.
914  */
915 struct ProcInfo : public TagBase {
916 
917  /**
918  * Intitialize default values.
919  */
920  ProcInfo(): iproc(0), loops(0), qcdorder(-1), eworder(-1) {}
921 
922  /**
923  * Create from XML tag.
924  */
925  ProcInfo(const XMLTag & tag)
926  : TagBase(tag.attr, tag.contents),
927  iproc(0), loops(0), qcdorder(-1), eworder(-1) {
928  getattr("iproc", iproc);
929  getattr("loops", loops);
930  getattr("qcdorder", qcdorder);
931  getattr("eworder", eworder);
932  getattr("rscheme", rscheme);
933  getattr("fscheme", fscheme);
934  getattr("scheme", scheme);
935  }
936 
937  /**
938  * Print out an XML tag.
939  */
940  void print(std::ostream & file) const {
941  file << "<procinfo" << oattr("iproc", iproc);
942  if ( loops >= 0 ) file << oattr("loops", loops);
943  if ( qcdorder >= 0 ) file << oattr("qcdorder", qcdorder);
944  if ( eworder >= 0 ) file<< oattr("eworder", eworder);
945  if ( !rscheme.empty() ) file << oattr("rscheme", rscheme);
946  if ( !fscheme.empty() ) file << oattr("fscheme", fscheme);
947  if ( !scheme.empty() ) file << oattr("scheme", scheme);
948  printattrs(file);
949  closetag(file, "procinfo");
950  }
951 
952  /**
953  * The id number for the process.
954  */
955  int iproc;
956 
957  /**
958  * The number of loops
959  */
960  int loops;
961 
962  /**
963  * The number of QCD vertices.
964  */
965  int qcdorder;
966 
967  /**
968  * The number of electro-weak vertices.
969  */
970  int eworder;
971 
972  /**
973  * The factorization scheme used.
974  */
975  std::string fscheme;
976 
977  /**
978  * The renormalization scheme used.
979  */
980  std::string rscheme;
981 
982  /**
983  * The NLO scheme used.
984  */
985  std::string scheme;
986 
987 };
988 
989 /**
990  * The MergeInfo class represents the information in a mergeinfo tag.
991  */
992 struct MergeInfo : public TagBase {
993 
994  /**
995  * Intitialize default values.
996  */
997  MergeInfo(): iproc(0), mergingscale(0.0), maxmult(false) {}
998 
999  /**
1000  * Create from XML tag.
1001  */
1002  MergeInfo(const XMLTag & tag)
1003  : TagBase(tag.attr, tag.contents),
1004  iproc(0), mergingscale(0.0), maxmult(false) {
1005  getattr("iproc", iproc);
1006  getattr("mergingscale", mergingscale);
1007  getattr("maxmult", maxmult);
1008  }
1009 
1010  /**
1011  * Print out an XML tag.
1012  */
1013  void print(std::ostream & file) const {
1014  file << "<mergeinfo" << oattr("iproc", iproc);
1015  if ( mergingscale > 0.0 ) file << oattr("mergingscale", mergingscale);
1016  if ( maxmult ) file << oattr("maxmult", yes());
1017  printattrs(file);
1018  closetag(file, "mergeinfo");
1019  }
1020 
1021  /**
1022  * The id number for the process.
1023  */
1024  int iproc;
1025 
1026  /**
1027  * The merging scale used if different from the cut definitions.
1028  */
1030 
1031  /**
1032  * Is this event reweighted as if it was the maximum multiplicity.
1033  */
1034  bool maxmult;
1035 
1036 };
1037 
1038 /**
1039  * The WeightInfo class encodes the description of a given weight
1040  * present for all events.
1041  */
1042 struct WeightInfo : public TagBase {
1043 
1044  /**
1045  * Constructors
1046  */
1047  WeightInfo(): inGroup(-1), isrwgt(false),
1048  muf(1.0), mur(1.0), pdf(0), pdf2(0) {}
1049 
1050  /**
1051  * Construct from the XML tag
1052  */
1053  WeightInfo(const XMLTag & tag)
1054  : TagBase(tag.attr, tag.contents),
1055  inGroup(-1), isrwgt(tag.name == "weight"),
1056  muf(1.0), mur(1.0), pdf(0), pdf2(0) {
1057  getattr("mur", mur);
1058  getattr("muf", muf);
1059  getattr("pdf", pdf);
1060  getattr("pdf2", pdf2);
1061  if ( isrwgt )
1062  getattr("id", name);
1063  else
1064  getattr("name", name);
1065  }
1066 
1067  /**
1068  * Print out an XML tag.
1069  */
1070  void print(std::ostream & file) const {
1071 
1072  if ( isrwgt )
1073  file << "<weight" << oattr("id", name);
1074  else
1075  file << "<weightinfo" << oattr("name", name);
1076  if ( mur != 1.0 ) file << oattr("mur", mur);
1077  if ( muf != 1.0 ) file << oattr("muf", muf);
1078  if ( pdf != 0 ) file << oattr("pdf", pdf);
1079  if ( pdf2 != 0 ) file << oattr("pdf2", pdf2);
1080  printattrs(file);
1081  if ( isrwgt )
1082  closetag(file, "weight");
1083  else
1084  closetag(file, "weightinfo");
1085  }
1086 
1087  /**
1088  * If inside a group, this is the index of that group.
1089  */
1090  int inGroup;
1091 
1092  /**
1093  * Is this a weightinfo or an rwgt tag?
1094  */
1095  bool isrwgt;
1096 
1097  /**
1098  * The name.
1099  */
1100  std::string name;
1101 
1102  /**
1103  * Factor multiplying the nominal factorization scale for this weight.
1104  */
1105  double muf;
1106 
1107  /**
1108  * Factor multiplying the nominal renormalization scale for this weight.
1109  */
1110  double mur;
1111 
1112  /**
1113  * The LHAPDF set relevant for this weight
1114  */
1115  long pdf;
1116 
1117  /**
1118  * The LHAPDF set for the second beam relevant for this weight if
1119  * different from pdf.
1120  */
1121  long pdf2;
1122 
1123 };
1124 
1125 /**
1126  * The WeightGroup assigns a group-name to a set of WeightInfo objects.
1127  */
1128 struct WeightGroup : public TagBase {
1129 
1130  /**
1131  * Default constructor;
1132  */
1134 
1135  /**
1136  * Construct a group of WeightInfo objects from an XML tag and
1137  * insert them in the given vector.
1138  */
1139  WeightGroup(const XMLTag & tag, int groupIndex, std::vector<WeightInfo> & wiv)
1140  : TagBase(tag.attr) {
1141  getattr("type", type);
1142  getattr("combine", combine);
1143  for ( int i = 0, N = tag.tags.size(); i < N; ++i ) {
1144  if ( tag.tags[i]->name == "weight" ||
1145  tag.tags[i]->name == "weightinfo" ) {
1146  WeightInfo wi(*tag.tags[i]);
1147  wi.inGroup = groupIndex;
1148  wiv.push_back(wi);
1149  }
1150  }
1151  }
1152 
1153  /**
1154  * The type.
1155  */
1156  std::string type;
1157 
1158  /**
1159  * The way in which these weights should be combined.
1160  */
1161  std::string combine;
1162 
1163 };
1164 
1165 
1166 /**
1167  * The Weight class represents the information in a weight tag.
1168  */
1169 struct Weight : public TagBase {
1170 
1171  /**
1172  * Initialize default values.
1173  */
1174  Weight(): iswgt(false), born(0.0), sudakov(0.0) {}
1175 
1176  /**
1177  * Create from XML tag
1178  */
1179  Weight(const XMLTag & tag)
1180  : TagBase(tag.attr, tag.contents),
1181  iswgt(tag.name == "wgt"), born(0.0), sudakov(0.0) {
1182  if ( iswgt )
1183  getattr("id", name);
1184  else
1185  getattr("name", name);
1186  getattr("born", born);
1187  getattr("sudakov", sudakov);
1188  std::istringstream iss(tag.contents);
1189  double w;
1190  while ( iss >> w ) weights.push_back(w);
1191  indices.resize(weights.size(), 0.0);
1192  }
1193 
1194  /**
1195  * Print out an XML tag.
1196  */
1197  void print(std::ostream & file) const {
1198  if ( iswgt )
1199  file << "<wgt" << oattr("id", name);
1200  else {
1201  file << "<weight";
1202  if ( !name.empty() ) file << oattr("name", name);
1203  }
1204  if ( born != 0.0 ) file << oattr("born", born);
1205  if ( sudakov != 0.0 ) file << oattr("sudakov", sudakov);
1206  file << ">";
1207  for ( int j = 0, M = weights.size(); j < M; ++j ) file << " " << weights[j];
1208  if ( iswgt )
1209  file << "</wgt>" << std::endl;
1210  else
1211  file << "</weight>" << std::endl;
1212  }
1213 
1214  /**
1215  * The identifyer for this set of weights.
1216  */
1217  std::string name;
1218 
1219  /**
1220  * Is this a wgt or a weight tag
1221  */
1222  bool iswgt;
1223 
1224  /**
1225  * The relative size of the born cross section of this event.
1226  */
1227  double born;
1228 
1229  /**
1230  * The relative size of the sudakov applied to this event.
1231  */
1232  double sudakov;
1233 
1234  /**
1235  * The weights of this event.
1236  */
1237  mutable std::vector<double> weights;
1238 
1239  /**
1240  * The indices where the weights are stored.
1241  */
1242  std::vector<int> indices;
1243 
1244 };
1245 
1246 /**
1247  * The Clus class represents a clustering of two particle entries into
1248  * one as defined in a clustering tag.
1249  */
1250 struct Clus : public TagBase {
1251 
1252  /**
1253  * Initialize default values.
1254  */
1255  Clus(): p1(0), p2(0), p0(0), scale(-1.0), alphas(-1.0) {}
1256 
1257  /**
1258  * Initialize default values.
1259  */
1260  Clus(const XMLTag & tag)
1261  : TagBase(tag.attr, tag.contents), scale(-1.0), alphas(-1.0) {
1262  getattr("scale", scale);
1263  getattr("alphas", alphas);
1264  std::istringstream iss(tag.contents);
1265  iss >> p1 >> p2;
1266  if ( !( iss >> p0 ) ) p0 = p1;
1267  }
1268 
1269  /**
1270  * Print out an XML tag.
1271  */
1272  void print(std::ostream & file) const {
1273  file << "<clus";
1274  if ( scale > 0.0 ) file << oattr("scale", scale);
1275  if ( alphas > 0.0 ) file << oattr("alphas", alphas);
1276  file << ">" << p1 << " " << p2;
1277  if ( p1 != p0 ) file << " " << p0;
1278  file << "</clus>" << std::endl;
1279  }
1280 
1281  /**
1282  * The first particle entry that has been clustered.
1283  */
1284  int p1;
1285 
1286  /**
1287  * The second particle entry that has been clustered.
1288  */
1289  int p2;
1290 
1291  /**
1292  * The particle entry corresponding to the clustered particles.
1293  */
1294  int p0;
1295 
1296  /**
1297  * The scale in GeV associated with the clustering.
1298  */
1299  double scale;
1300 
1301  /**
1302  * The alpha_s used in the corresponding vertex, if this was used in
1303  * the cross section.
1304  */
1305  double alphas;
1306 
1307 };
1308 
1309 /**
1310  * Store special scales from within a scales tag.
1311  */
1312 
1313 struct Scale : public TagBase {
1314 
1315  /**
1316  * Empty constructor
1317  */
1318  Scale(std::string st = "veto", int emtr = 0, double sc = 0.0)
1319  : stype(st), emitter(emtr), scale(sc) {}
1320 
1321  /**
1322  * Construct from an XML-tag.
1323  */
1324  Scale(const XMLTag & tag)
1325  : TagBase(tag.attr, tag.contents),stype("veto"), emitter(0) {
1326  if ( !getattr("stype", stype) )
1327  throw std::runtime_error("Found scale tag without stype attribute "
1328  "in Les Houches Event File.");
1329  std::string pattr;
1330  if ( getattr("pos", pattr) ) {
1331  std::istringstream pis(pattr);
1332  if ( !(pis >> emitter) ) emitter = 0;
1333  else {
1334  int rec = 0;
1335  while ( pis >> rec ) recoilers.insert(rec);
1336  }
1337  }
1338 
1339  std::string eattr;
1340  if ( getattr("etype", eattr) ) {
1341  if ( eattr == "QCD" ) eattr = "-5 -4 -3 -2 -1 1 2 3 4 5 21";
1342  if ( eattr == "EW" ) eattr = "-13 -12 -11 11 12 13 22 23 24";
1343  std::istringstream eis(eattr);
1344  int pdg = 0;
1345  while ( eis >> pdg ) emitted.insert(pdg);
1346  }
1347  std::istringstream cis(tag.contents);
1348  cis >> scale;
1349 
1350  }
1351 
1352  /**
1353  * Print out an XML tag.
1354  */
1355  void print(std::ostream & file) const {
1356  file << "<scale" << oattr("stype", stype);
1357  if ( emitter > 0 ) {
1358  std::ostringstream pos;
1359  pos << emitter;
1360  for ( std::set<int>::iterator it = recoilers.begin();
1361  it != recoilers.end(); ++it )
1362  pos << " " << *it;
1363  file << oattr("pos", pos.str());
1364  }
1365  if ( emitted.size() > 0 ) {
1366  std::set<int>::iterator it = emitted.begin();
1367  std::ostringstream eos;
1368  eos << *it;
1369  while ( ++it != emitted.end() ) eos << " " << *it;
1370  if ( eos.str() == "-5 -4 -3 -2 -1 1 2 3 4 5 21" )
1371  file << oattr("etype", std::string("QCD"));
1372  else if ( eos.str() == "-13 -12 -11 11 12 13 22 23 24" )
1373  file << oattr("etype", std::string("EW"));
1374  else
1375  file << oattr("etype", eos.str());
1376  }
1377  std::ostringstream os;
1378  os << scale;
1379  contents = os.str();
1380  closetag(file, "scale");
1381  }
1382 
1383  /**
1384  * The type of scale this represents. Predefined values are "veto"
1385  * and "start".
1386  */
1387  std::string stype;
1388 
1389  /**
1390  * The emitter this scale applies to. This is the index of a
1391  * particle in HEPEUP (starting at 1). Zero corresponds to any
1392  * particle in HEPEUP.
1393  */
1394  int emitter;
1395 
1396  /**
1397  * The set of recoilers for which this scale applies.
1398  */
1399  std::set<int> recoilers;
1400 
1401  /**
1402  * The set of emitted particles (PDG id) this applies to.
1403  */
1404  std::set<int> emitted;
1405 
1406  /**
1407  * The actual scale given.
1408  */
1409  double scale;
1410 
1411 };
1412 
1413 /**
1414  * Collect different scales relevant for an event.
1415  */
1416 struct Scales : public TagBase {
1417 
1418  /**
1419  * Empty constructor.
1420  */
1421  Scales(double defscale = -1.0, int npart = 0)
1422  : muf(defscale), mur(defscale), mups(defscale), SCALUP(defscale) {
1423  (void) npart; // avoid "unused variable" compiler warning
1424  }
1425 
1426  /**
1427  * Construct from an XML-tag
1428  */
1429  Scales(const XMLTag & tag, double defscale = -1.0, int npart = 0)
1430  : TagBase(tag.attr, tag.contents),
1431  muf(defscale), mur(defscale), mups(defscale),
1432  SCALUP(defscale) {
1433  getattr("muf", muf);
1434  getattr("mur", mur);
1435  getattr("mups", mups);
1436  for ( int i = 0, N = tag.tags.size(); i < N; ++i )
1437  if ( tag.tags[i]->name == "scale" )
1438  scales.push_back(Scale(*tag.tags[i]));
1439  for ( int i = 0; i < npart; ++i ) {
1440  std::ostringstream pttag;
1441  pttag << "pt_start_" << i + 1;
1442  double sc = 0.0;
1443  if ( getattr(pttag.str(), sc) )
1444  scales.push_back(Scale("start", i + 1, sc));
1445  }
1446 
1447  }
1448 
1449  /**
1450  * Check if this object contains useful information besides SCALUP.
1451  */
1452  bool hasInfo() const {
1453  return muf != SCALUP || mur != SCALUP || mups != SCALUP ||
1454  !scales.empty();
1455  }
1456 
1457  /**
1458  * Print out the corresponding XML-tag.
1459  */
1460  void print(std::ostream & file) const {
1461  if ( !hasInfo() ) return;
1462  file << "<scales";
1463  if ( muf != SCALUP ) file << oattr("muf", muf);
1464  if ( mur != SCALUP ) file << oattr("mur", mur);
1465  if ( mups != SCALUP ) file << oattr("mups", mups);
1466  printattrs(file);
1467 
1468  if ( !scales.empty() ) {
1469  std::ostringstream os;
1470  for ( int i = 0, N = scales.size(); i < N; ++i )
1471  scales[i].print(os);
1472  contents = os.str();
1473  }
1474  closetag(file, "scales");
1475  }
1476 
1477  /**
1478  * Return the scale of type st for a given emission of particle type
1479  * pdgem from the emitter with number emr and a recoiler rec. (Note
1480  * that the indices for emr and rec starts at 1 and 0 is interpreted
1481  * as any particle.) First it will check for Scale object with an
1482  * exact match. If not found, it will search for an exact match for
1483  * the emitter and recoiler with an undefined emitted particle. If
1484  * not found, it will look for a match for only emitter and emitted,
1485  * of if not found, a match for only the emitter. Finally a general
1486  * Scale object will be used, or if nothing matches, the mups will
1487  * be returned.
1488  */
1489  double getScale(std::string st, int pdgem, int emr, int rec) const {
1490  for ( int i = 0, N = scales.size(); i < N; ++i ) {
1491  if ( scales[i].emitter == emr && st == scales[i].stype &&
1492  ( emr == rec ||
1493  scales[i].recoilers.find(rec) != scales[i].recoilers.end() ) &&
1494  scales[i].emitted.find(pdgem) != scales[i].emitted.end() )
1495  return scales[i].scale;
1496  }
1497  for ( int i = 0, N = scales.size(); i < N; ++i ) {
1498  if ( scales[i].emitter == emr && st == scales[i].stype &&
1499  ( emr == rec ||
1500  scales[i].recoilers.find(rec) != scales[i].recoilers.end() ) &&
1501  scales[i].emitted.empty() )
1502  return scales[i].scale;
1503  }
1504  if ( emr != rec ) return getScale(st, pdgem, emr, emr);
1505  if ( emr == rec ) return getScale(st, pdgem, 0, 0);
1506  return mups;
1507  }
1508 
1509  /**
1510  * The factorization scale used for this event.
1511  */
1512  double muf;
1513 
1514  /**
1515  * The renormalization scale used for this event.
1516  */
1517  double mur;
1518 
1519  /**
1520  * The starting scale for the parton shower as suggested by the
1521  * matrix element generator.
1522  */
1523  double mups;
1524 
1525  /**
1526  * The default scale in this event.
1527  */
1528  double SCALUP;
1529 
1530  /**
1531  * The list of special scales.
1532  */
1533  std::vector<Scale> scales;
1534 
1535 };
1536 
1537 /**
1538  * The PDFInfo class represents the information in a pdfinto tag.
1539  */
1540 struct PDFInfo : public TagBase {
1541 
1542  /**
1543  * Initialize default values.
1544  */
1545  PDFInfo(double defscale = -1.0): p1(0), p2(0), x1(-1.0), x2(-1.0),
1546  xf1(-1.0), xf2(-1.0), scale(defscale), SCALUP(defscale) {}
1547 
1548  /**
1549  * Create from XML tag.
1550  */
1551  PDFInfo(const XMLTag & tag, double defscale = -1.0)
1552  : TagBase(tag.attr, tag.contents),
1553  p1(0), p2(0), x1(-1.0), x2(-1.0), xf1(-1.0), xf2(-1.0),
1554  scale(defscale), SCALUP(defscale) {
1555  getattr("scale", scale);
1556  getattr("p1", p1);
1557  getattr("p2", p2);
1558  getattr("x1", x1);
1559  getattr("x2", x2);
1560  }
1561 
1562  /**
1563  * Print out an XML tag.
1564  */
1565  void print(std::ostream & file) const {
1566  if ( xf1 <= 0 ) return;
1567  file << "<pdfinfo";
1568  if ( p1 != 0 ) file << oattr("p1", p1);
1569  if ( p2 != 0 ) file << oattr("p2", p2);
1570  if ( x1 > 0 ) file << oattr("x1", x1);
1571  if ( x2 > 0 ) file << oattr("x2", x2);
1572  if ( scale != SCALUP ) file << oattr("scale", scale);
1573  printattrs(file);
1574  file << ">" << xf1 << " " << xf2 << "</pdfinfo>" << std::endl;
1575  }
1576 
1577  /**
1578  * The type of the incoming particle 1.
1579  */
1580  long p1;
1581 
1582  /**
1583  * The type of the incoming particle 2.
1584  */
1585  long p2;
1586 
1587  /**
1588  * The x-value used for the incoming particle 1.
1589  */
1590  double x1;
1591 
1592  /**
1593  * The x-value used for the incoming particle 2.
1594  */
1595  double x2;
1596 
1597  /**
1598  * The value of the pdf for the incoming particle 1.
1599  */
1600  double xf1;
1601 
1602  /**
1603  * The value of the pdf for the incoming particle 2.
1604  */
1605  double xf2;
1606 
1607  /**
1608  * The scale used in the PDF:s
1609  */
1610  double scale;
1611 
1612  /**
1613  * THe default scale in the event.
1614  */
1615  double SCALUP;
1616 
1617 };
1618 
1619 /**
1620  * The HEPRUP class is a simple container corresponding to the Les Houches
1621  * accord (<A HREF="http://arxiv.org/abs/hep-ph/0109068">hep-ph/0109068</A>)
1622  * common block with the same name. The members are named in the same
1623  * way as in the common block. However, fortran arrays are represented
1624  * by vectors, except for the arrays of length two which are
1625  * represented by pair objects.
1626  */
1627 class HEPRUP : public TagBase {
1628 
1629 public:
1630 
1631  /** @name Standard constructors and destructors. */
1632  //@{
1633  /**
1634  * Default constructor.
1635  */
1637  : IDWTUP(0), NPRUP(0), version(3),
1638  dprec(std::numeric_limits<double>::digits10) {}
1639 
1640  /**
1641  * Copy constructor
1642  */
1643  HEPRUP(const HEPRUP &) = default;
1644 
1645  /**
1646  * Assignment operator.
1647  */
1648  HEPRUP & operator=(const HEPRUP & x) = default;
1649 
1650  /**
1651  * Construct from a given init tag.
1652  */
1653  HEPRUP(const XMLTag & tagin, int versin)
1654  : TagBase(tagin.attr, tagin.contents), version(versin),
1655  dprec(std::numeric_limits<double>::digits10) {
1656 
1657  std::vector<XMLTag*> tags = tagin.tags;
1658 
1659  // The first (anonymous) tag should just be the init block.
1660  std::istringstream iss(tags[0]->contents);
1661  if ( !( iss >> IDBMUP.first >> IDBMUP.second >> EBMUP.first >> EBMUP.second
1662  >> PDFGUP.first >> PDFGUP.second >> PDFSUP.first >> PDFSUP.second
1663  >> IDWTUP >> NPRUP ) ) {
1664  throw std::runtime_error("Could not parse init block "
1665  "in Les Houches Event File.");
1666  }
1667  resize();
1668 
1669  for ( int i = 0; i < NPRUP; ++i ) {
1670  if ( !( iss >> XSECUP[i] >> XERRUP[i] >> XMAXUP[i] >> LPRUP[i] ) ) {
1671  throw std::runtime_error("Could not parse processes in init block "
1672  "in Les Houches Event File.");
1673  }
1674  }
1675 
1676  for ( int i = 1, N = tags.size(); i < N; ++i ) {
1677  const XMLTag & tag = *tags[i];
1678 
1679  if ( tag.name.empty() ) junk += tag.contents;
1680 
1681  if ( tag.name == "initrwgt" ) {
1682  for ( int j = 0, M = tag.tags.size(); j < M; ++j ) {
1683  if ( tag.tags[j]->name == "weightgroup" )
1684  weightgroup.push_back(WeightGroup(*tag.tags[j], weightgroup.size(),
1685  weightinfo));
1686  if ( tag.tags[j]->name == "weight" )
1687  weightinfo.push_back(WeightInfo(*tag.tags[j]));
1688 
1689  }
1690  }
1691  if ( tag.name == "weightinfo" ) {
1692  weightinfo.push_back(WeightInfo(tag));
1693  }
1694  if ( tag.name == "weightgroup" ) {
1695  weightgroup.push_back(WeightGroup(tag, weightgroup.size(),
1696  weightinfo));
1697  }
1698  if ( tag.name == "eventfiles" ) {
1699  for ( int j = 0, M = tag.tags.size(); j < M; ++j ) {
1700  XMLTag & eftag = *tag.tags[j];
1701  if ( eftag.name == "eventfile" )
1702  eventfiles.push_back(EventFile(eftag));
1703  }
1704  }
1705  if ( tag.name == "xsecinfo" ) {
1706  XSecInfo xsecinfo = XSecInfo(tag);
1707  xsecinfos[xsecinfo.weightname] = xsecinfo;
1708  }
1709  if ( tag.name == "generator" ) {
1710  generators.push_back(Generator(tag));
1711  }
1712  else if ( tag.name == "cutsinfo" ) {
1713  for ( int j = 0, M = tag.tags.size(); j < M; ++j ) {
1714  XMLTag & ctag = *tag.tags[j];
1715 
1716  if ( ctag.name == "ptype" ) {
1717  std::string tname = ctag.attr["name"];
1718  long id;
1719  std::istringstream isss(ctag.contents);
1720  while ( isss >> id ) ptypes[tname].insert(id);
1721  }
1722  else if ( ctag.name == "cut" )
1723  cuts.push_back(Cut(ctag, ptypes));
1724  }
1725  }
1726  else if ( tag.name == "procinfo" ) {
1727  ProcInfo proc(tag);
1728  procinfo[proc.iproc] = proc;
1729  }
1730  else if ( tag.name == "mergeinfo" ) {
1731  MergeInfo merge(tag);
1732  mergeinfo[merge.iproc] = merge;
1733  }
1734 
1735  }
1736 
1737  weightmap.clear();
1738  for ( int i = 0, N = weightinfo.size(); i < N; ++i )
1739  weightmap[weightinfo[i].name] = i + 1;
1740 
1741  }
1742 
1743  //@}
1744 
1745 public:
1746 
1747  /**
1748  * Return the name of the weight with given index suitable to ne
1749  * used for HepMC3 output.
1750  */
1751  std::string weightNameHepMC(int i) const {
1752  std::string name;
1753  if ( i < 0 || i >= (int)weightinfo.size() ) return name;
1754  if ( weightinfo[i].inGroup >= 0 )
1755  name = weightgroup[weightinfo[i].inGroup].type + "/"
1756  + weightgroup[weightinfo[i].inGroup].combine + "/";
1757  name += weightinfo[i].name;
1758  return name;
1759  }
1760 
1761 
1762  /**
1763  * Print out the corresponding XML tag to a stream.
1764  */
1765  void print(std::ostream & file) const {
1766 
1767  file << std::setprecision(dprec);
1768 
1769  file << "<init>\n"
1770  << " " << std::setw(8) << IDBMUP.first
1771  << " " << std::setw(8) << IDBMUP.second
1772  << " " << std::setw(14) << EBMUP.first
1773  << " " << std::setw(14) << EBMUP.second
1774  << " " << std::setw(4) << PDFGUP.first
1775  << " " << std::setw(4) << PDFGUP.second
1776  << " " << std::setw(4) << PDFSUP.first
1777  << " " << std::setw(4) << PDFSUP.second
1778  << " " << std::setw(4) << IDWTUP
1779  << " " << std::setw(4) << NPRUP << std::endl;
1780 
1781  for ( int i = 0; i < NPRUP; ++i )
1782  file << " " << std::setw(14) << XSECUP[i]
1783  << " " << std::setw(14) << XERRUP[i]
1784  << " " << std::setw(14) << XMAXUP[i]
1785  << " " << std::setw(6) << LPRUP[i] << std::endl;
1786 
1787  for ( int i = 0, N = generators.size(); i < N; ++i )
1788  generators[i].print(file);
1789 
1790  if ( !eventfiles.empty() ) {
1791  file << "<eventfiles>\n";
1792  for ( int i = 0, N = eventfiles.size(); i < N; ++i )
1793  eventfiles[i].print(file);
1794  file << "</eventfiles>\n";
1795  }
1796  //AV if ( !xsecinfos.empty() > 0 )
1797  if ( !xsecinfos.empty())
1798  for ( XSecInfos::const_iterator it = xsecinfos.begin();
1799  it != xsecinfos.end(); ++it )
1800  if ( it->second.neve > 0 ) it->second.print(file);
1801 
1802  if ( cuts.size() > 0 ) {
1803  file << "<cutsinfo>" << std::endl;
1804 
1805  for ( std::map<std::string, std::set<long> >::const_iterator ptit =
1806  ptypes.begin(); ptit != ptypes.end(); ++ptit ) {
1807  file << "<ptype" << oattr("name", ptit->first) << ">";
1808  for ( std::set<long>::const_iterator it = ptit->second.begin();
1809  it != ptit->second.end(); ++it )
1810  file << " " << *it;
1811  file << "</ptype>" << std::endl;
1812  }
1813 
1814  for ( int i = 0, N = cuts.size(); i < N; ++i )
1815  cuts[i].print(file);
1816  file << "</cutsinfo>" << std::endl;
1817  }
1818 
1819  for ( std::map<long,ProcInfo>::const_iterator it = procinfo.begin();
1820  it != procinfo.end(); ++it )
1821  it->second.print(file);
1822 
1823  for ( std::map<long,MergeInfo>::const_iterator it = mergeinfo.begin();
1824  it != mergeinfo.end(); ++it )
1825  it->second.print(file);
1826 
1827  bool isrwgt = false;
1828  int ingroup = -1;
1829  for ( int i = 0, N = weightinfo.size(); i < N; ++i ) {
1830  if ( weightinfo[i].isrwgt ) {
1831  if ( !isrwgt ) file << "<initrwgt>\n";
1832  isrwgt = true;
1833  } else {
1834  if ( isrwgt ) file << "</initrwgt>\n";
1835  isrwgt = false;
1836  }
1837  int group = weightinfo[i].inGroup;
1838  if ( group != ingroup ) {
1839  if ( ingroup != -1 ) file << "</weightgroup>\n";
1840  if ( group != -1 ) {
1841  file << "<weightgroup"
1842  << oattr("type", weightgroup[group].type);
1843  if ( !weightgroup[group].combine.empty() )
1844  file << oattr("combine", weightgroup[group].combine);
1845  file << ">\n";
1846  }
1847  ingroup = group;
1848  }
1849  weightinfo[i].print(file);
1850  }
1851  if ( ingroup != -1 ) file << "</weightgroup>\n";
1852  if ( isrwgt ) file << "</initrwgt>\n";
1853 
1854 
1855  file << hashline(junk) << "</init>" << std::endl;
1856 
1857  }
1858 
1859  /**
1860  * Clear all information.
1861  */
1862  void clear() {
1863  procinfo.clear();
1864  mergeinfo.clear();
1865  weightinfo.clear();
1866  weightgroup.clear();
1867  cuts.clear();
1868  ptypes.clear();
1869  junk.clear();
1870  }
1871 
1872  /**
1873  * Set the NPRUP variable, corresponding to the number of
1874  * sub-processes, to \a nrup, and resize all relevant vectors
1875  * accordingly.
1876  */
1877  void resize(int nrup) {
1878  NPRUP = nrup;
1879  resize();
1880  }
1881 
1882  /**
1883  * Assuming the NPRUP variable, corresponding to the number of
1884  * sub-processes, is correctly set, resize the relevant vectors
1885  * accordingly.
1886  */
1887  void resize() {
1888  XSECUP.resize(NPRUP);
1889  XERRUP.resize(NPRUP);
1890  XMAXUP.resize(NPRUP);
1891  LPRUP.resize(NPRUP);
1892  }
1893 
1894  /**
1895  * @return the index of the weight with the given \a name
1896  */
1897  int weightIndex(std::string name) const {
1898  std::map<std::string, int>::const_iterator it = weightmap.find(name);
1899  if ( it != weightmap.end() ) return it->second;
1900  return 0;
1901  }
1902 
1903  /**
1904  * @return the number of weights (including the nominial one).
1905  */
1906  int nWeights() const {
1907  return weightmap.size() + 1;
1908  }
1909 
1910  /**
1911  * @return the XSecInfo object corresponding to the named weight \a
1912  * weithname. If no such object exists, it will be created.
1913  */
1914  XSecInfo & getXSecInfo(std::string weightname = "") {
1915  XSecInfo & xi = xsecinfos[weightname];
1916  xi.weightname = weightname;
1917  return xi;
1918  }
1919 
1920  /**
1921  * @return the XSecInfo object corresponding to the named weight \a
1922  * weithname. If no such object exists, an empty XSecInfo will be
1923  * returned..
1924  */
1925  const XSecInfo & getXSecInfo(std::string weightname = "") const {
1926  static XSecInfo noinfo;
1927  XSecInfos::const_iterator it = xsecinfos.find(weightname);
1928  if ( it != xsecinfos.end() ) return it->second;
1929  else return noinfo;
1930  }
1931 
1932 
1933 public:
1934 
1935  /**
1936  * PDG id's of beam particles. (first/second is in +/-z direction).
1937  */
1938  std::pair<long,long> IDBMUP;
1939 
1940  /**
1941  * Energy of beam particles given in GeV.
1942  */
1943  std::pair<double,double> EBMUP;
1944 
1945  /**
1946  * The author group for the PDF used for the beams according to the
1947  * PDFLib specification.
1948  */
1949  std::pair<int,int> PDFGUP;
1950 
1951  /**
1952  * The id number the PDF used for the beams according to the
1953  * PDFLib specification.
1954  */
1955  std::pair<int,int> PDFSUP;
1956 
1957  /**
1958  * Master switch indicating how the ME generator envisages the
1959  * events weights should be interpreted according to the Les Houches
1960  * accord.
1961  */
1962  int IDWTUP;
1963 
1964  /**
1965  * The number of different subprocesses in this file.
1966  */
1967  int NPRUP;
1968 
1969  /**
1970  * The cross sections for the different subprocesses in pb.
1971  */
1972  std::vector<double> XSECUP;
1973 
1974  /**
1975  * The statistical error in the cross sections for the different
1976  * subprocesses in pb.
1977  */
1978  std::vector<double> XERRUP;
1979 
1980  /**
1981  * The maximum event weights (in HEPEUP::XWGTUP) for different
1982  * subprocesses.
1983  */
1984  std::vector<double> XMAXUP;
1985 
1986  /**
1987  * The subprocess code for the different subprocesses.
1988  */
1989  std::vector<int> LPRUP;
1990 
1991  /**
1992  * Contents of the xsecinfo tags.
1993  */
1995 
1996  /**
1997  * A vector of EventFiles where the events are stored separate fron
1998  * the init block.
1999  */
2000  std::vector<EventFile> eventfiles;
2001 
2002  /**
2003  * Contents of the cuts tag.
2004  */
2005  std::vector<Cut> cuts;
2006 
2007  /**
2008  * A map of codes for different particle types.
2009  */
2010  std::map<std::string, std::set<long> > ptypes;
2011 
2012  /**
2013  * Contents of the procinfo tags
2014  */
2015  std::map<long,ProcInfo> procinfo;
2016 
2017  /**
2018  * Contents of the mergeinfo tags
2019  */
2020  std::map<long,MergeInfo> mergeinfo;
2021 
2022  /**
2023  * The names of the programs and their version information used to
2024  * create this file.
2025  */
2026  std::vector<Generator> generators;
2027 
2028  /**
2029  * The vector of WeightInfo objects for this file.
2030  */
2031  std::vector<WeightInfo> weightinfo;
2032 
2033  /**
2034  * A map relating names of weights to indices of the weightinfo vector.
2035  */
2036  std::map<std::string,int> weightmap;
2037 
2038  /**
2039  * The vector of WeightGroup objects in this file.
2040  */
2041  std::vector<WeightGroup> weightgroup;
2042 
2043  /**
2044  * Just to be on the safe side we save any junk inside the init-tag.
2045  */
2046  std::string junk;
2047 
2048  /**
2049  * The main version of the information stored.
2050  */
2051  int version;
2052 
2053  /**
2054  * The precision used for outputing real numbers.
2055  */
2056  int dprec;
2057 
2058 };
2059 
2060 /**
2061  * Forward declaration of the HEPEUP class.
2062  */
2063 class HEPEUP;
2064 
2065 /**
2066  * The EventGroup represents a set of events which are to be
2067  * considered together.
2068  */
2069 struct EventGroup: public std::vector<HEPEUP*> {
2070 
2071  /**
2072  * Initialize default values.
2073  */
2074  inline EventGroup(): nreal(-1), ncounter(-1) {}
2075 
2076  /**
2077  * The copy constructor also copies the included HEPEUP object.
2078  */
2079  inline EventGroup(const EventGroup &);
2080 
2081  /**
2082  * The assignment also copies the included HEPEUP object.
2083  */
2084  inline EventGroup & operator=(const EventGroup &);
2085 
2086  /**
2087  * Remove all subevents.
2088  */
2089  inline void clear();
2090 
2091  /**
2092  * The destructor deletes the included HEPEUP objects.
2093  */
2094  inline ~EventGroup();
2095 
2096  /**
2097  * The number of real events in this event group.
2098  */
2099  int nreal;
2100 
2101  /**
2102  * The number of counter events in this event group.
2103  */
2105 
2106 };
2107 
2108 
2109 /**
2110  * The HEPEUP class is a simple container corresponding to the Les Houches accord
2111  * (<A HREF="http://arxiv.org/abs/hep-ph/0109068">hep-ph/0109068</A>)
2112  * common block with the same name. The members are named in the same
2113  * way as in the common block. However, fortran arrays are represented
2114  * by vectors, except for the arrays of length two which are
2115  * represented by pair objects.
2116  */
2117 class HEPEUP : public TagBase {
2118 
2119 public:
2120 
2121  /** @name Standard constructors and destructors. */
2122  //@{
2123  /**
2124  * Default constructor.
2125  */
2127  : NUP(0), IDPRUP(0), XWGTUP(0.0), XPDWUP(0.0, 0.0),
2128  SCALUP(0.0), AQEDUP(0.0), AQCDUP(0.0), heprup(0), currentWeight(0),
2129  ntries(1), isGroup(false) {}
2130 
2131  /**
2132  * Copy constructor
2133  */
2134  HEPEUP(const HEPEUP & x)
2135  : TagBase(x), isGroup(false) {
2136  operator=(x);
2137  }
2138 
2139  /**
2140  * Copy information from the given HEPEUP. Sub event information is
2141  * left untouched.
2142  */
2143  HEPEUP & setEvent(const HEPEUP & x) {
2144  NUP = x.NUP;
2145  IDPRUP = x.IDPRUP;
2146  XWGTUP = x.XWGTUP;
2147  XPDWUP = x.XPDWUP;
2148  SCALUP = x.SCALUP;
2149  AQEDUP = x.AQEDUP;
2150  AQCDUP = x.AQCDUP;
2151  IDUP = x.IDUP;
2152  ISTUP = x.ISTUP;
2153  MOTHUP = x.MOTHUP;
2154  ICOLUP = x.ICOLUP;
2155  PUP = x.PUP;
2156  VTIMUP = x.VTIMUP;
2157  SPINUP = x.SPINUP;
2158  heprup = x.heprup;
2160  weights = x.weights;
2161  pdfinfo = x.pdfinfo;
2162  PDFGUPsave = x.PDFGUPsave;
2163  PDFSUPsave = x.PDFSUPsave;
2164  clustering = x.clustering;
2165  scales = x.scales;
2166  junk = x.junk;
2168  ntries = x.ntries;
2169  return *this;
2170  }
2171 
2172  /**
2173  * Assignment operator.
2174  */
2175  HEPEUP & operator=(const HEPEUP & x) {
2176  if ( &x == this ) return *this;
2177  TagBase::operator=(x);
2178  clear();
2179  setEvent(x);
2180  subevents = x.subevents;
2181  isGroup = x.isGroup;
2182  return *this;
2183  }
2184 
2185  /**
2186  * Destructor.
2187  */
2189  clear();
2190  };
2191  //@}
2192 
2193 public:
2194 
2195 
2196  /**
2197  * Constructor from an event or eventgroup tag.
2198  */
2199  HEPEUP(const XMLTag & tagin, HEPRUP & heprupin)
2200  : TagBase(tagin.attr), NUP(0), IDPRUP(0), XWGTUP(0.0), XPDWUP(0.0, 0.0),
2201  SCALUP(0.0), AQEDUP(0.0), AQCDUP(0.0), heprup(&heprupin),
2202  currentWeight(0), ntries(1), isGroup(tagin.name == "eventgroup") {
2203 
2204  if ( heprup->NPRUP < 0 )
2205  throw std::runtime_error("Tried to read events but no processes defined "
2206  "in init block of Les Houches file.");
2207 
2208  std::vector<XMLTag*> tags = tagin.tags;
2209 
2210  if ( isGroup ) {
2211  getattr("nreal", subevents.nreal);
2212  getattr("ncounter", subevents.ncounter);
2213  for ( int i = 0, N = tags.size(); i < N; ++i )
2214  if ( tags[i]->name == "event" )
2215  subevents.push_back(new HEPEUP(*tags[i], heprupin));
2216  return;
2217  } else
2218  getattr("ntries", ntries);
2219 
2220 
2221 
2222  // The event information should be in the first (anonymous) tag
2223  std::istringstream iss(tags[0]->contents);
2224  if ( !( iss >> NUP >> IDPRUP >> XWGTUP >> SCALUP >> AQEDUP >> AQCDUP ) )
2225  throw std::runtime_error("Failed to parse event in Les Houches file.");
2226 
2227  resize();
2228 
2229  // Read all particle lines.
2230  for ( int i = 0; i < NUP; ++i ) {
2231  if ( !( iss >> IDUP[i] >> ISTUP[i] >> MOTHUP[i].first >> MOTHUP[i].second
2232  >> ICOLUP[i].first >> ICOLUP[i].second
2233  >> PUP[i][0] >> PUP[i][1] >> PUP[i][2]
2234  >> PUP[i][3] >> PUP[i][4]
2235  >> VTIMUP[i] >> SPINUP[i] ) )
2236  throw std::runtime_error("Failed to parse event in Les Houches file.");
2237  }
2238 
2239  junk.clear();
2240  std::string ss;
2241  while ( getline(iss, ss) ) junk += ss + '\n';
2242 
2243  scales = Scales(SCALUP, NUP);
2244  pdfinfo = PDFInfo(SCALUP);
2245  namedweights.clear();
2246  weights.clear();
2247  weights.resize(heprup->nWeights(),
2248  std::make_pair(XWGTUP, (WeightInfo*)(0)));
2249  weights.front().first = XWGTUP;
2250  for ( int i = 1, N = weights.size(); i < N; ++i )
2251  weights[i].second = &heprup->weightinfo[i - 1];
2252 
2253  for ( int i = 1, N = tags.size(); i < N; ++i ) {
2254  XMLTag & tag = *tags[i];
2255 
2256  if ( tag.name.empty() ) junk += tag.contents;
2257 
2258  if ( tag.name == "weights" ) {
2259  weights.resize(heprup->nWeights(),
2260  std::make_pair(XWGTUP, (WeightInfo*)(0)));
2261  weights.front().first = XWGTUP;
2262  for ( int ii = 1, NN = weights.size(); ii < NN; ++ii )
2263  weights[ii].second = &heprup->weightinfo[ii - 1];
2264  double w = 0.0;
2265  int iii = 0;
2266  std::istringstream isss(tag.contents);
2267  while ( isss >> w )
2268  if ( ++iii < int(weights.size()) )
2269  weights[iii].first = w;
2270  else
2271  weights.push_back(std::make_pair(w, (WeightInfo*)(0)));
2272  }
2273  if ( tag.name == "weight" ) {
2274  namedweights.push_back(Weight(tag));
2275  }
2276  if ( tag.name == "rwgt" ) {
2277  for ( int j = 0, M = tag.tags.size(); j < M; ++j ) {
2278  if ( tag.tags[j]->name == "wgt" ) {
2279  namedweights.push_back(Weight(*tag.tags[j]));
2280  }
2281  }
2282  }
2283  else if ( tag.name == "clustering" ) {
2284  for ( int j = 0, M= tag.tags.size(); j < M; ++j ) {
2285  if ( tag.tags[j]->name == "clus" )
2286  clustering.push_back(Clus(*tag.tags[j]));
2287  }
2288  }
2289  else if ( tag.name == "pdfinfo" ) {
2290  pdfinfo = PDFInfo(tag, SCALUP);
2291  }
2292  else if ( tag.name == "scales" ) {
2293  scales = Scales(tag, SCALUP, NUP);
2294  }
2295 
2296  }
2297 
2298  for ( int i = 0, N = namedweights.size(); i < N; ++i ) {
2299  int indx = heprup->weightIndex(namedweights[i].name);
2300  if ( indx > 0 ) {
2301  weights[indx].first = namedweights[i].weights[0];
2302  namedweights[i].indices[0] = indx;
2303  } else {
2304  weights.push_back(std::make_pair(namedweights[i].weights[0],
2305  (WeightInfo*)(0)));
2306  namedweights[i].indices[0] = weights.size() - 1;
2307  }
2308  for ( int j = 1, M = namedweights[i].weights.size(); j < M; ++j ) {
2309  weights.push_back(std::make_pair(namedweights[i].weights[j],
2310  (WeightInfo*)(0)));
2311  namedweights[i].indices[j] = weights.size() - 1;
2312  }
2313  }
2314 
2315  }
2316 
2317  /**
2318  * Print out the event (group) as an XML tag.
2319  */
2320  void print(std::ostream & file) const {
2321 
2322  file << std::setprecision(heprup->dprec);
2323 
2324 
2325  if ( isGroup ) {
2326  file << "<eventgroup";
2327  if ( subevents.nreal > 0 )
2328  file << oattr("nreal", subevents.nreal);
2329  if ( subevents.ncounter > 0 )
2330  file << oattr("ncounter", subevents.ncounter);
2331  printattrs(file);
2332  file << ">\n";
2333  for ( int i = 0, N = subevents.size(); i < N; ++i )
2334  subevents[i]->print(file);
2335  file << "</eventgroup>\n";
2336  return;
2337  }
2338 
2339  file << "<event";
2340  if ( ntries > 1 ) file << oattr("ntries", ntries);
2341  printattrs(file);
2342  file << ">\n";
2343  file << " " << std::setw(4) << NUP
2344  << " " << std::setw(6) << IDPRUP
2345  << " " << std::setw(14) << XWGTUP
2346  << " " << std::setw(14) << SCALUP
2347  << " " << std::setw(14) << AQEDUP
2348  << " " << std::setw(14) << AQCDUP << "\n";
2349 
2350  for ( int i = 0; i < NUP; ++i )
2351  file << " " << std::setw(8) << IDUP[i]
2352  << " " << std::setw(2) << ISTUP[i]
2353  << " " << std::setw(4) << MOTHUP[i].first
2354  << " " << std::setw(4) << MOTHUP[i].second
2355  << " " << std::setw(4) << ICOLUP[i].first
2356  << " " << std::setw(4) << ICOLUP[i].second
2357  << " " << std::setw(14) << PUP[i][0]
2358  << " " << std::setw(14) << PUP[i][1]
2359  << " " << std::setw(14) << PUP[i][2]
2360  << " " << std::setw(14) << PUP[i][3]
2361  << " " << std::setw(14) << PUP[i][4]
2362  << " " << std::setw(1) << VTIMUP[i]
2363  << " " << std::setw(1) << SPINUP[i] << std::endl;
2364 
2365  if ( weights.size() > 0 ) {
2366  file << "<weights>";
2367  for ( int i = 1, N = weights.size(); i < N; ++i )
2368  file << " " << weights[i].first;
2369  file << "</weights>\n";
2370  }
2371 
2372  bool iswgt = false;
2373  for ( int i = 0, N = namedweights.size(); i < N; ++i ) {
2374  if ( namedweights[i].iswgt ) {
2375  if ( !iswgt ) file << "<rwgt>\n";
2376  iswgt = true;
2377  } else {
2378  if ( iswgt ) file << "</rwgt>\n";
2379  iswgt = false;
2380  }
2381  for ( int j = 0, M = namedweights[i].indices.size(); j < M; ++j )
2382  namedweights[i].weights[j] = weight(namedweights[i].indices[j]);
2383  namedweights[i].print(file);
2384  }
2385  if ( iswgt ) file << "</rwgt>\n";
2386 
2387  if ( !clustering.empty() ) {
2388  file << "<clustering>" << std::endl;
2389  for ( int i = 0, N = clustering.size(); i < N; ++i )
2390  clustering[i].print(file);
2391  file << "</clustering>" << std::endl;
2392  }
2393 
2394  pdfinfo.print(file);
2395  scales.print(file);
2396 
2397  file << hashline(junk) << "</event>\n";
2398 
2399  }
2400 
2401  /**
2402  * Reset the HEPEUP object (does not touch the sub events).
2403  */
2404  void reset() {
2405  setWeightInfo(0);
2406  NUP = 0;
2407  clustering.clear();
2408  weights.clear();
2409  }
2410 
2411  /**
2412  * Clear the HEPEUP object.
2413  */
2414  void clear() {
2415  reset();
2416  subevents.clear();
2417  }
2418 
2419  /**
2420  * Set the NUP variable, corresponding to the number of particles in
2421  * the current event, to \a nup, and resize all relevant vectors
2422  * accordingly.
2423  */
2424  void resize(int nup) {
2425  NUP = nup;
2426  resize();
2427  }
2428 
2429  /**
2430  * Return the total weight for this event (including all sub
2431  * evenets) for the given index.
2432  */
2433  double totalWeight(int i = 0) const {
2434  if ( subevents.empty() ) return weight(i);
2435  double w = 0.0;
2436  for ( int ii = 0, N = subevents.size(); ii < N; ++ii )
2437  w += subevents[ii]->weight(i);
2438  return w;
2439  }
2440 
2441  /**
2442  * Return the total weight for this event (including all sub
2443  * evenets) for the given weight name.
2444  */
2445  double totalWeight(std::string name) const {
2446  return totalWeight(heprup->weightIndex(name));
2447  }
2448 
2449  /**
2450  * Return the weight for the given index.
2451  */
2452  double weight(int i = 0) const {
2453  return weights[i].first;
2454  }
2455 
2456  /**
2457  * Return the weight for the given weight name.
2458  */
2459  double weight(std::string name) const {
2460  return weight(heprup->weightIndex(name));
2461  }
2462 
2463  /**
2464  * Set the weight with the given index.
2465  */
2466  void setWeight(int i, double w) {
2467  weights[i].first = w;
2468  }
2469  /**
2470  * Set the weight with the given name.
2471  */
2472  bool setWeight(std::string name, double w) {
2473  int i = heprup->weightIndex(name);
2474  if ( i >= int(weights.size()) ) return false;
2475  setWeight(i, w);
2476  return true;
2477  }
2478 
2479 
2480  /**
2481  * Assuming the NUP variable, corresponding to the number of
2482  * particles in the current event, is correctly set, resize the
2483  * relevant vectors accordingly.
2484  */
2485  void resize() {
2486  IDUP.resize(NUP);
2487  ISTUP.resize(NUP);
2488  MOTHUP.resize(NUP);
2489  ICOLUP.resize(NUP);
2490  PUP.resize(NUP, std::vector<double>(5));
2491  VTIMUP.resize(NUP);
2492  SPINUP.resize(NUP);
2493  }
2494 
2495  /**
2496  * Setup the current event to use weight i. If zero, the default
2497  * weight will be used.
2498  */
2499  bool setWeightInfo(unsigned int i) {
2500  if ( i >= weights.size() ) return false;
2501  if ( currentWeight ) {
2506  }
2507  XWGTUP = weights[i].first;
2508  currentWeight = weights[i].second;
2509  if ( currentWeight ) {
2514  if ( currentWeight->pdf ) {
2515  heprup->PDFGUP.first = heprup->PDFGUP.second = 0;
2516  heprup->PDFSUP.first = heprup->PDFSUP.second = currentWeight->pdf;
2517  }
2518  if ( currentWeight->pdf2 ) {
2519  heprup->PDFSUP.second = currentWeight->pdf2;
2520  }
2521 
2522  }
2523  return true;
2524  }
2525 
2526  /**
2527  * Setup the current event to use sub event i. If zero, no sub event
2528  * will be chsen.
2529  */
2530  bool setSubEvent(unsigned int i) {
2531  if ( i > subevents.size() || subevents.empty() ) return false;
2532  if ( i == 0 ) {
2533  reset();
2534  weights = subevents[0]->weights;
2535  for ( int ii = 1, N = subevents.size(); ii < N; ++ii )
2536  for ( int j = 0, M = weights.size(); j < M; ++j )
2537  weights[j].first += subevents[ii]->weights[j].first;
2538  currentWeight = 0;
2539  } else {
2540  setEvent(*subevents[i - 1]);
2541  }
2542  return true;
2543  }
2544 
2545 public:
2546 
2547  /**
2548  * The number of particle entries in the current event.
2549  */
2550  int NUP;
2551 
2552  /**
2553  * The subprocess code for this event (as given in LPRUP).
2554  */
2555  int IDPRUP;
2556 
2557  /**
2558  * The weight for this event.
2559  */
2560  double XWGTUP;
2561 
2562  /**
2563  * The PDF weights for the two incoming partons. Note that this
2564  * variable is not present in the current LesHouches accord
2565  * (<A HREF="http://arxiv.org/abs/hep-ph/0109068">hep-ph/0109068</A>),
2566  * hopefully it will be present in a future accord.
2567  */
2568  std::pair<double,double> XPDWUP;
2569 
2570  /**
2571  * The scale in GeV used in the calculation of the PDF's in this
2572  * event.
2573  */
2574  double SCALUP;
2575 
2576  /**
2577  * The value of the QED coupling used in this event.
2578  */
2579  double AQEDUP;
2580 
2581  /**
2582  * The value of the QCD coupling used in this event.
2583  */
2584  double AQCDUP;
2585 
2586  /**
2587  * The PDG id's for the particle entries in this event.
2588  */
2589  std::vector<long> IDUP;
2590 
2591  /**
2592  * The status codes for the particle entries in this event.
2593  */
2594  std::vector<int> ISTUP;
2595 
2596  /**
2597  * Indices for the first and last mother for the particle entries in
2598  * this event.
2599  */
2600  std::vector< std::pair<int,int> > MOTHUP;
2601 
2602  /**
2603  * The colour-line indices (first(second) is (anti)colour) for the
2604  * particle entries in this event.
2605  */
2606  std::vector< std::pair<int,int> > ICOLUP;
2607 
2608  /**
2609  * Lab frame momentum (Px, Py, Pz, E and M in GeV) for the particle
2610  * entries in this event.
2611  */
2612  std::vector< std::vector<double> > PUP;
2613 
2614  /**
2615  * Invariant lifetime (c*tau, distance from production to decay in
2616  * mm) for the particle entries in this event.
2617  */
2618  std::vector<double> VTIMUP;
2619 
2620  /**
2621  * Spin info for the particle entries in this event given as the
2622  * cosine of the angle between the spin vector of a particle and the
2623  * 3-momentum of the decaying particle, specified in the lab frame.
2624  */
2625  std::vector<double> SPINUP;
2626 
2627  /**
2628  * A pointer to the current HEPRUP object.
2629  */
2631 
2632  /**
2633  * The current weight info object.
2634  */
2636 
2637  /**
2638  * The weights associated with this event
2639  */
2640  std::vector<Weight> namedweights;
2641 
2642  /**
2643  * The weights for this event and their corresponding WeightInfo object.
2644  */
2645  std::vector< std::pair<double, const WeightInfo *> > weights;
2646 
2647  /**
2648  * Contents of the clustering tag.
2649  */
2650  std::vector<Clus> clustering;
2651 
2652  /**
2653  * Contents of the pdfinfo tag.
2654  */
2656 
2657  /**
2658  * Saved information about pdfs if different in a selected weight.
2659  */
2660  std::pair<int,int> PDFGUPsave;
2661 
2662  /**
2663  * Saved information about pdfs if different in a selected weight.
2664  */
2665  std::pair<int,int> PDFSUPsave;
2666 
2667 
2668  /**
2669  * Contents of the scales tag
2670  */
2672 
2673  /**
2674  * The number of attempts the ME generator did before accepting this
2675  * event.
2676  */
2677  int ntries;
2678 
2679  /**
2680  * Is this an event or an event group?
2681  */
2682  bool isGroup;
2683 
2684  /**
2685  * If this is not a single event, but an event group, the events
2686  * included in the group are in this vector;
2687  */
2689 
2690  /**
2691  * Save junk stuff in events just to be on the safe side
2692  */
2693  std::string junk;
2694 
2695 };
2696 
2697 
2698 // Destructor implemented here.
2699 
2700 inline void EventGroup::clear() {
2701  while ( size() > 0 ) {
2702  delete back();
2703  pop_back();
2704  }
2705 }
2706 
2708  clear();
2709 }
2710 
2712  : std::vector<HEPEUP*>(eg.size()) {
2713  for ( int i = 0, N = eg.size(); i < N; ++i ) at(i) = new HEPEUP(*eg.at(i));
2714 }
2715 
2717  if ( &x == this ) return *this;
2718  clear();
2719  nreal = x.nreal;
2720  ncounter = x.ncounter;
2721  for ( int i = 0, N = x.size(); i < N; ++i ) push_back(new HEPEUP(*x.at(i)));
2722  return *this;
2723 }
2724 
2725 
2726 /**
2727  * The Reader class is initialized with a stream from which to read a
2728  * version 1/2 Les Houches Accord event file. In the constructor of
2729  * the Reader object the optional header information is read and then
2730  * the mandatory init is read. After this the whole header block
2731  * including the enclosing lines with tags are available in the public
2732  * headerBlock member variable. Also the information from the init
2733  * block is available in the heprup member variable and any additional
2734  * comment lines are available in initComments. After each successful
2735  * call to the readEvent() function the standard Les Houches Accord
2736  * information about the event is available in the hepeup member
2737  * variable and any additional comments in the eventComments
2738  * variable. A typical reading sequence would look as follows:
2739  *
2740  *
2741  */
2742 class Reader {
2743 
2744 public:
2745 
2746  /**
2747  * Initialize the Reader with a stream from which to read an event
2748  * file. After the constructor is called the whole header block
2749  * including the enclosing lines with tags are available in the
2750  * public headerBlock member variable. Also the information from the
2751  * init block is available in the heprup member variable and any
2752  * additional comment lines are available in initComments.
2753  *
2754  * @param is the stream to read from.
2755  */
2756  Reader(std::istream & is)
2757  : file(&is), currevent(-1),
2758  curreventfile(-1), currfileevent(-1), dirpath("") {
2759  init();
2760  }
2761 
2762  /**
2763  * Initialize the Reader with a filename from which to read an event
2764  * file. After the constructor is called the whole header block
2765  * including the enclosing lines with tags are available in the
2766  * public headerBlock member variable. Also the information from the
2767  * init block is available in the heprup member variable and any
2768  * additional comment lines are available in initComments.
2769  *
2770  * @param filename the name of the file to read from.
2771  */
2772  Reader(std::string filename)
2773  : intstream(filename.c_str()), file(&intstream), currevent(-1),
2774  curreventfile(-1), currfileevent(-1), dirpath("") {
2775 
2776  size_t slash = filename.find_last_of('/');
2777  if ( slash != std::string::npos ) dirpath = filename.substr(0, slash + 1);
2778  init();
2779  }
2780 
2781 private:
2782 
2783  /**
2784  * Used internally in the constructors to read header and init
2785  * blocks.
2786  */
2787  void init() {
2788 
2789  // initialize reading from multi-file runs.
2790  initfile = file;
2791 
2792  bool readingHeader = false;
2793  bool readingInit = false;
2794 
2795  // Make sure we are reading a LHEF file:
2796  getline();
2797  if ( !currentFind("<LesHouchesEvents") )
2798  throw std::runtime_error
2799  ("Tried to read a file which does not start with the "
2800  "LesHouchesEvents tag.");
2801  version = 1;
2802  if ( currentFind("version=\"3" ) )
2803  version = 3;
2804  else if ( currentFind("version=\"2" ) )
2805  version = 2;
2806  else if ( !currentFind("version=\"1" ) )
2807  throw std::runtime_error
2808  ("Tried to read a LesHouchesEvents file which is above version 3.");
2809 
2810  // Loop over all lines until we hit the </init> tag.
2811  while ( getline() && !currentFind("</init>") ) {
2812  if ( currentFind("<header") ) {
2813  // We have hit the header block, so we should dump this and
2814  // all following lines to headerBlock until we hit the end of
2815  // it.
2816  readingHeader = true;
2817  headerBlock = currentLine + "\n";
2818  }
2819  else if ( currentFind("<init>") ) {
2820  // We have hit the init block
2821  readingInit = true;
2822  initComments = currentLine + "\n";
2823  }
2824  else if ( currentFind("</header>") ) {
2825  // The end of the header block. Dump this line as well to the
2826  // headerBlock and we're done.
2827  readingHeader = false;
2828  headerBlock += currentLine + "\n";
2829  }
2830  else if ( readingHeader ) {
2831  // We are in the process of reading the header block. Dump the
2832  // line to haderBlock.
2833  headerBlock += currentLine + "\n";
2834  }
2835  else if ( readingInit ) {
2836  // Here we found a comment line. Dump it to initComments.
2837  initComments += currentLine + "\n";
2838  }
2839  else {
2840  // We found some other stuff outside the standard tags.
2841  outsideBlock += currentLine + "\n";
2842  }
2843  }
2844  if ( !currentFind("</init>") )
2845  throw std::runtime_error("Found incomplete init tag in "
2846  "Les Houches file.");
2847  initComments += currentLine + "\n";
2848  std::vector<XMLTag*> tags = XMLTag::findXMLTags(initComments);
2849  for ( int i = 0, N = tags.size(); i < N; ++i )
2850  if ( tags[i]->name == "init" ) {
2851  heprup = HEPRUP(*tags[i], version);
2852  break;
2853  }
2854  XMLTag::deleteAll(tags);
2855 
2856  if ( !heprup.eventfiles.empty() ) openeventfile(0);
2857 
2858  }
2859 
2860 public:
2861 
2862  /**
2863  * Read an event from the file and store it in the hepeup
2864  * object. Optional comment lines are stored i the eventComments
2865  * member variable.
2866  * @return true if the read sas successful.
2867  */
2868  bool readEvent() {
2869 
2870  // Check if the initialization was successful. Otherwise we will
2871  // not read any events.
2872  if ( heprup.NPRUP < 0 ) return false;
2873 
2874  std::string eventLines;
2875  int inEvent = 0;
2876 
2877  // Keep reading lines until we hit the end of an event or event group.
2878  while ( getline() ) {
2879  if ( inEvent ) {
2880  eventLines += currentLine + "\n";
2881  if ( inEvent == 1 && currentFind("</event>") ) break;
2882  if ( inEvent == 2 && currentFind("</eventgroup>") ) break;
2883  }
2884  else if ( currentFind("<eventgroup") ) {
2885  eventLines += currentLine + "\n";
2886  inEvent = 2;
2887  }
2888  else if ( currentFind("<event") ) {
2889  eventLines += currentLine + "\n";
2890  inEvent = 1;
2891  }
2892  else {
2893  outsideBlock += currentLine + "\n";
2894  }
2895  }
2896  if ( ( inEvent == 1 && !currentFind("</event>") ) ||
2897  ( inEvent == 2 && !currentFind("</eventgroup>") ) ) {
2898  if ( heprup.eventfiles.empty() ||
2899  ++curreventfile >= int(heprup.eventfiles.size()) ) return false;
2901  return readEvent();
2902  }
2903 
2904  std::vector<XMLTag*> tags = XMLTag::findXMLTags(eventLines);
2905 
2906  for ( int i = 0, N = tags.size(); i < N ; ++i ) {
2907  if ( tags[i]->name == "event" || tags[i]->name == "eventgroup" ) {
2908  hepeup = HEPEUP(*tags[i], heprup);
2909  XMLTag::deleteAll(tags);
2910  ++currevent;
2911  if ( curreventfile >= 0 ) ++currfileevent;
2912  return true;
2913  }
2914  }
2915 
2916  if ( !heprup.eventfiles.empty() &&
2917  ++curreventfile < int(heprup.eventfiles.size()) ) {
2919  return readEvent();
2920  }
2921 
2922  XMLTag::deleteAll(tags);
2923  return false;
2924 
2925  }
2926 
2927  /**
2928  * Open the efentfile with index ifile. If another eventfile is
2929  * being read, its remaining contents is discarded. This is a noop
2930  * if current read session is not a multi-file run.
2931  */
2932  void openeventfile(int ifile) {
2933  std::cerr << "opening file " << ifile << std::endl;
2934  efile.close();
2935  std::string fname = heprup.eventfiles[ifile].filename;
2936  if ( fname[0] != '/' ) fname = dirpath + fname;
2937  efile.open(fname.c_str());
2938  if ( !efile ) throw std::runtime_error("Could not open event file " +
2939  fname);
2940  file = &efile;
2941  curreventfile = ifile;
2942  currfileevent = 0;
2943  }
2944 
2945 protected:
2946 
2947  /**
2948  * Used internally to read a single line from the stream.
2949  */
2950  bool getline() {
2951  return ( (bool)std::getline(*file, currentLine) );
2952  }
2953 
2954  /**
2955  * @return true if the current line contains the given string.
2956  */
2957  bool currentFind(std::string str) const {
2958  return currentLine.find(str) != std::string::npos;
2959  }
2960 
2961 protected:
2962 
2963  /**
2964  * A local stream which is unused if a stream is supplied from the
2965  * outside.
2966  */
2967  std::ifstream intstream;
2968 
2969  /**
2970  * The stream we are reading from. This may be a pointer to an
2971  * external stream or the internal intstream, or a separate event
2972  * file from a multi-file run
2973  */
2974  std::istream * file;
2975 
2976  /**
2977  * The original stream from where we read the init block.
2978  */
2979  std::istream * initfile;
2980 
2981  /**
2982  * A separate stream for reading multi-file runs.
2983  */
2984  std::ifstream efile;
2985 
2986  /**
2987  * The last line read in from the stream in getline().
2988  */
2989  std::string currentLine;
2990 
2991 public:
2992 
2993  /**
2994  * XML file version
2995  */
2996  int version;
2997 
2998  /**
2999  * All lines (since the last readEvent()) outside the header, init
3000  * and event tags.
3001  */
3002  std::string outsideBlock;
3003 
3004  /**
3005  * All lines from the header block.
3006  */
3007  std::string headerBlock;
3008 
3009  /**
3010  * The standard init information.
3011  */
3013 
3014  /**
3015  * Additional comments found in the init block.
3016  */
3017  std::string initComments;
3018 
3019  /**
3020  * The standard information about the last read event.
3021  */
3023 
3024  /**
3025  * Additional comments found with the last read event.
3026  */
3027  std::string eventComments;
3028 
3029  /**
3030  * The number of the current event (starting from 1).
3031  */
3033 
3034  /**
3035  * The current event file being read from (-1 means there are no
3036  * separate event files).
3037  */
3039 
3040  /**
3041  * The number of the current event in the current event file.
3042  */
3044 
3045  /**
3046  * The directory from where we are reading files.
3047  */
3048  std::string dirpath;
3049 
3050 private:
3051 
3052  /**
3053  * The default constructor should never be used.
3054  */
3056 
3057  /**
3058  * The copy constructor should never be used.
3059  */
3060  Reader(const Reader &);
3061 
3062  /**
3063  * The Reader cannot be assigned to.
3064  */
3066 
3067 };
3068 
3069 /**
3070  * The Writer class is initialized with a stream to which to write a
3071  * version 1.0 Les Houches Accord event file. In the constructor of
3072  * the Writer object the main XML tag is written out, with the
3073  * corresponding end tag is written in the destructor. After a Writer
3074  * object has been created, it is possible to assign standard init
3075  * information in the heprup member variable. In addition any XML
3076  * formatted information can be added to the headerBlock member
3077  * variable (directly or via the addHeader() function). Further
3078  * comment line (beginning with a <code>#</code> character) can be
3079  * added to the initComments variable (directly or with the
3080  * addInitComment() function). After this information is set, it
3081  * should be written out to the file with the init() function.
3082  *
3083  * Before each event is written out with the writeEvent() function,
3084  * the standard event information can then be assigned to the hepeup
3085  * variable and optional comment lines (beginning with a
3086  * <code>#</code> character) may be given to the eventComments
3087  * variable (directly or with the addEventComment() function).
3088  *
3089  */
3090 class Writer {
3091 
3092 public:
3093 
3094 #ifndef HEPMC3_PYTHON_BINDINGS
3095  /**
3096  * Create a Writer object giving a stream to write to.
3097  * @param os the stream where the event file is written.
3098  */
3099  Writer(std::ostream & os)
3100  : file(&os), initfile(&os), dirpath("") { }
3101 #endif
3102  /**
3103  * Create a Writer object giving a filename to write to.
3104  * @param filename the name of the event file to be written.
3105  */
3106  Writer(std::string filename)
3107  : intstream(filename.c_str()), file(&intstream), initfile(&intstream),
3108  dirpath("") {
3109  size_t slash = filename.find_last_of('/');
3110  if ( slash != std::string::npos ) dirpath = filename.substr(0, slash + 1);
3111  }
3112 
3113  /**
3114  * The destructor writes out the final XML end-tag.
3115  */
3117  file = initfile;
3118  if ( !heprup.eventfiles.empty() ) {
3119  if ( curreventfile >= 0 &&
3120  curreventfile < int(heprup.eventfiles.size()) &&
3121  heprup.eventfiles[curreventfile].neve < 0 )
3123  writeinit();
3124  }
3125  *file << "</LesHouchesEvents>" << std::endl;
3126  }
3127 #ifndef HEPMC3_PYTHON_BINDINGS
3128  /**
3129  * Add header lines consisting of XML code with this stream.
3130  */
3131  std::ostream & headerBlock() {
3132  return headerStream;
3133  }
3134 
3135  /**
3136  * Add comment lines to the init block with this stream.
3137  */
3138  std::ostream & initComments() {
3139  return initStream;
3140  }
3141 
3142  /**
3143  * Add comment lines to the next event to be written out with this stream.
3144  */
3145  std::ostream & eventComments() {
3146  return eventStream;
3147  }
3148 #endif
3149  /**
3150  * Add header lines consisting of XML code with this stream.
3151  */
3152  void headerBlock(const std::string& a) {
3153  headerStream<<a;;
3154  }
3155 
3156  /**
3157  * Add comment lines to the init block with this stream.
3158  */
3159  void initComments(const std::string& a) {
3160  initStream<<a;
3161  }
3162 
3163  /**
3164  * Add comment lines to the next event to be written out with this stream.
3165  */
3166  void eventComments(const std::string& a) {
3167  eventStream<<a;
3168  }
3169 
3170  /**
3171  * Initialize the writer.
3172  */
3173  void init() {
3174  if ( heprup.eventfiles.empty() ) writeinit();
3175  lastevent = 0;
3177  if ( !heprup.eventfiles.empty() ) openeventfile(0);
3178  }
3179 
3180  /**
3181  * Open a new event file, possibly closing a previous opened one.
3182  */
3183  bool openeventfile(int ifile) {
3184  if ( heprup.eventfiles.empty() ) return false;
3185  if ( ifile < 0 || ifile >= int(heprup.eventfiles.size()) ) return false;
3186  if ( curreventfile >= 0 ) {
3188  if ( ef.neve > 0 && ef.neve != currfileevent )
3189  std::cerr << "LHEF::Writer number of events in event file "
3190  << ef.filename << " does not match the given number."
3191  << std::endl;
3192  ef.neve = currfileevent;
3193  }
3194  efile.close();
3195  std::string fname = heprup.eventfiles[ifile].filename;
3196  if ( fname[0] != '/' ) fname = dirpath + fname;
3197  efile.open(fname.c_str());
3198  if ( !efile ) throw std::runtime_error("Could not open event file " +
3199  fname);
3200  std::cerr << "Opened event file " << fname << std::endl;
3201  file = &efile;
3202  curreventfile = ifile;
3203  currfileevent = 0;
3204  return true;
3205  }
3206 
3207 
3208  /**
3209  * Write out an optional header block followed by the standard init
3210  * block information together with any comment lines.
3211  */
3212  void writeinit() {
3213 
3214  // Write out the standard XML tag for the event file.
3215  if ( heprup.version == 3 )
3216  *file << "<LesHouchesEvents version=\"3.0\">\n";
3217  else if ( heprup.version == 2 )
3218  *file << "<LesHouchesEvents version=\"2.0\">\n";
3219  else
3220  *file << "<LesHouchesEvents version=\"1.0\">\n";
3221 
3222 
3223  *file << std::setprecision(10);
3224 
3225  std::string headBlock = headerStream.str();
3226  if ( headBlock.length() ) {
3227  if ( headBlock.find("<header>") == std::string::npos )
3228  *file << "<header>\n";
3229  if ( headBlock[headBlock.length() - 1] != '\n' )
3230  headBlock += '\n';
3231  *file << headBlock;
3232  if ( headBlock.find("</header>") == std::string::npos )
3233  *file << "</header>\n";
3234  }
3235 
3236  heprup.print(*file);
3237 
3238  }
3239 
3240  /**
3241  * Write the current HEPEUP object to the stream;
3242  */
3243  void writeEvent() {
3244 
3245  if ( !heprup.eventfiles.empty() ) {
3246  if ( currfileevent == heprup.eventfiles[curreventfile].neve &&
3247  curreventfile + 1 < int(heprup.eventfiles.size()) )
3249  }
3250 
3251  hepeup.print(*file);
3252 
3253  ++lastevent;
3254  ++currfileevent;
3255  }
3256 
3257 protected:
3258 
3259  /**
3260  * A local stream which is unused if a stream is supplied from the
3261  * outside.
3262  */
3263  std::ofstream intstream;
3264 
3265  /**
3266  * The stream we are writing to. This may be a reference to an
3267  * external stream or the internal intstream.
3268  */
3269  std::ostream * file;
3270 
3271  /**
3272  * The original stream from where we read the init block.
3273  */
3274  std::ostream * initfile;
3275 
3276  /**
3277  * A separate stream for reading multi-file runs.
3278  */
3279  std::ofstream efile;
3280 
3281  /**
3282  * The number of the last event written (starting from 1).
3283  */
3285 
3286  /**
3287  * The current event file being written to (-1 means there are no
3288  * separate event files).
3289  */
3291 
3292  /**
3293  * The number of the current event in the current event file.
3294  */
3296 
3297  /**
3298  * The directory from where we are reading files.
3299  */
3300  std::string dirpath;
3301 
3302 public:
3303  /**
3304  * The standard init information.
3305  */
3307 
3308 
3309  /**
3310  * The standard information about the event we will write next.
3311  */
3313 
3314 
3315 
3316 private:
3317 
3318  /**
3319  * Stream to add all lines in the header block.
3320  */
3321  std::ostringstream headerStream;
3322 
3323  /**
3324  * Stream to add additional comments to be put in the init block.
3325  */
3326  std::ostringstream initStream;
3327 
3328  /**
3329  * Stream to add additional comments to be written together the next event.
3330  */
3331  std::ostringstream eventStream;
3332 
3333 
3334  /**
3335  * The default constructor should never be used.
3336  */
3338 
3339  /**
3340  * The copy constructor should never be used.
3341  */
3342  Writer(const Writer &);
3343 
3344  /**
3345  * The Writer cannot be assigned to.
3346  */
3348 
3349 };
3350 
3351 }
3352 
3353 /* \example LHEFCat.cc This is a main function which simply reads a
3354  Les Houches Event File from the standard input and writes it again
3355  to the standard output.
3356  This file can be downloaded from
3357  <A HREF="http://www.thep.lu.se/~leif/LHEF/LHEFCat.cc">here</A>.
3358  There are also two sample event files,
3359  <A HREF="http://www.thep.lu.se/~leif/LHEF/ttV_unweighted_events.lhe">ttV_unweighted_events.lhe</A> and <A HREF="http://www.thep.lu.se/~leif/LHEF/testlhef3.lhe">testlhef3.lhe</A>
3360  to try it on.
3361 */
3362 
3363 /* \mainpage Les Houches Event File
3364 
3365 Here are some example classes for reading and writing Les Houches
3366 Event Files according to the
3367 <A HREF="http://www.thep.lu.se/~torbjorn/lhef/lhafile2.pdf">proposal</A>
3368 by Torbj&ouml;rn Sj&ouml;strand discussed at the
3369 <A HREF="http://mc4lhc06.web.cern.ch/mc4lhc06/">MC4LHC</A>
3370 workshop at CERN 2006.
3371 
3372 The classes has now been updated to handle the suggested version 3 of
3373 this file standard as discussed at the <a href="http://phystev.in2p3.fr/wiki/2013:groups:tools:lhef3">Les Houches workshop 2013</a> (The previous suggested version 2.0 was discussed at the <a href="http://www.lpthe.jussieu.fr/LesHouches09Wiki/index.php/LHEF_for_Matching">Les Houches workshop 2009</a>).
3374 
3375 There is a whole set of classes available in a single header file
3376 called <A
3377 HREF="http://www.thep.lu.se/~leif/LHEF/LHEF.h">LHEF.h</A>. The idea is
3378 that matrix element or parton shower generators will implement their
3379 own wrapper using these classes as containers.
3380 
3381 The two classes LHEF::HEPRUP and LHEF::HEPEUP are simple container
3382 classes which contain the same information as the Les Houches standard
3383 Fortran common blocks with the same names. They also contain the extra
3384 information defined in version 3 in the standard. The other two main
3385 classes are called LHEF::Reader and LHEF::Writer and are used to read
3386 and write Les Houches Event Files
3387 
3388 Here are a few <A HREF="examples.html">examples</A> of how to use the
3389 classes:
3390 
3391 \namespace LHEF The LHEF namespace contains some example classes for reading and writing Les Houches
3392 Event Files according to the
3393 <A HREF="http://www.thep.lu.se/~torbjorn/lhef/lhafile2.pdf">proposal</A>
3394 by Torbj&ouml;rn Sj&ouml;strand discussed at the
3395 <A HREF="http://mc4lhc06.web.cern.ch/mc4lhc06/">MC4LHC</A>
3396 workshop at CERN 2006.
3397 
3398 
3399 
3400  */
3401 
3402 
3403 #endif /* HEPMC3_LHEF_H */
LHEF::Reader::getline
bool getline()
Definition: LHEF.h:2950
LHEF::HEPRUP::HEPRUP
HEPRUP()
Definition: LHEF.h:1636
LHEF::XMLTag::AttributeMap
std::map< std::string, std::string > AttributeMap
Definition: LHEF.h:97
LHEF::Reader::intstream
std::ifstream intstream
Definition: LHEF.h:2967
LHEF::PDFInfo::print
void print(std::ostream &file) const
Definition: LHEF.h:1565
LHEF::Cut::type
std::string type
Definition: LHEF.h:879
LHEF::TagBase::getattr
bool getattr(std::string n, bool &v, bool erase=true)
Definition: LHEF.h:382
LHEF::Cut::np1
std::string np1
Definition: LHEF.h:889
LHEF::Writer::currfileevent
int currfileevent
Definition: LHEF.h:3295
LHEF::Writer::eventComments
std::ostream & eventComments()
Definition: LHEF.h:3145
LHEF::EventFile::print
void print(std::ostream &file) const
Definition: LHEF.h:640
LHEF::EventGroup::clear
void clear()
Definition: LHEF.h:2700
LHEF::HEPRUP::resize
void resize(int nrup)
Definition: LHEF.h:1877
LHEF::WeightGroup
Definition: LHEF.h:1128
LHEF::Reader::Reader
Reader(const Reader &)
LHEF::HEPRUP::weightmap
std::map< std::string, int > weightmap
Definition: LHEF.h:2036
LHEF::WeightInfo::pdf
long pdf
Definition: LHEF.h:1115
LHEF::WeightInfo::inGroup
int inGroup
Definition: LHEF.h:1090
LHEF::Cut::deltaR
static double deltaR(const std::vector< double > &p1, const std::vector< double > &p2)
Definition: LHEF.h:860
LHEF::WeightInfo
Definition: LHEF.h:1042
LHEF::Writer::initComments
void initComments(const std::string &a)
Definition: LHEF.h:3159
LHEF::XSecInfos
std::map< std::string, XSecInfo > XSecInfos
Definition: LHEF.h:611
LHEF::HEPEUP::setWeight
void setWeight(int i, double w)
Definition: LHEF.h:2466
LHEF::XSecInfo::meanweight
double meanweight
Definition: LHEF.h:589
LHEF::Reader::operator=
Reader & operator=(const Reader &)
LHEF::EventFile
Definition: LHEF.h:617
LHEF::XSecInfo::negweights
bool negweights
Definition: LHEF.h:594
LHEF::Writer::init
void init()
Definition: LHEF.h:3173
LHEF::HEPEUP::HEPEUP
HEPEUP(const HEPEUP &x)
Definition: LHEF.h:2134
LHEF::Reader::currentLine
std::string currentLine
Definition: LHEF.h:2989
LHEF::HEPEUP::XWGTUP
double XWGTUP
Definition: LHEF.h:2560
LHEF::ProcInfo
Definition: LHEF.h:915
LHEF::Weight::weights
std::vector< double > weights
Definition: LHEF.h:1237
LHEF::Writer::efile
std::ofstream efile
Definition: LHEF.h:3279
LHEF::WeightGroup::type
std::string type
Definition: LHEF.h:1156
LHEF::XMLTag::end
static const pos_t end
Definition: LHEF.h:102
LHEF::Writer
Definition: LHEF.h:3090
LHEF::EventGroup::nreal
int nreal
Definition: LHEF.h:2099
LHEF::TagBase::TagBase
TagBase()
Definition: LHEF.h:355
LHEF::HEPRUP::IDWTUP
int IDWTUP
Definition: LHEF.h:1962
LHEF::WeightGroup::WeightGroup
WeightGroup(const XMLTag &tag, int groupIndex, std::vector< WeightInfo > &wiv)
Definition: LHEF.h:1139
LHEF::HEPRUP::XMAXUP
std::vector< double > XMAXUP
Definition: LHEF.h:1984
LHEF::Writer::curreventfile
int curreventfile
Definition: LHEF.h:3290
LHEF::Reader::outsideBlock
std::string outsideBlock
Definition: LHEF.h:3002
LHEF::PDFInfo::scale
double scale
Definition: LHEF.h:1610
LHEF::ProcInfo::rscheme
std::string rscheme
Definition: LHEF.h:980
LHEF::HEPEUP::operator=
HEPEUP & operator=(const HEPEUP &x)
Definition: LHEF.h:2175
LHEF::HEPEUP::pdfinfo
PDFInfo pdfinfo
Definition: LHEF.h:2655
LHEF::PDFInfo::xf2
double xf2
Definition: LHEF.h:1605
LHEF::Clus::p1
int p1
Definition: LHEF.h:1284
LHEF::WeightInfo::pdf2
long pdf2
Definition: LHEF.h:1121
LHEF::Cut::match
bool match(long id1, long id2=0) const
Definition: LHEF.h:749
LHEF::XMLTag::getattr
bool getattr(std::string n, long &v) const
Definition: LHEF.h:163
LHEF::MergeInfo::mergingscale
double mergingscale
Definition: LHEF.h:1029
LHEF::HEPRUP::nWeights
int nWeights() const
Definition: LHEF.h:1906
LHEF::HEPRUP::generators
std::vector< Generator > generators
Definition: LHEF.h:2026
LHEF::PDFInfo::x2
double x2
Definition: LHEF.h:1595
LHEF::Reader::currevent
int currevent
Definition: LHEF.h:3032
LHEF::Scales::Scales
Scales(double defscale=-1.0, int npart=0)
Definition: LHEF.h:1421
LHEF::Cut::p2
std::set< long > p2
Definition: LHEF.h:894
LHEF::HEPEUP::MOTHUP
std::vector< std::pair< int, int > > MOTHUP
Definition: LHEF.h:2600
LHEF::Writer::initfile
std::ostream * initfile
Definition: LHEF.h:3274
LHEF::Writer::dirpath
std::string dirpath
Definition: LHEF.h:3300
LHEF::HEPRUP::HEPRUP
HEPRUP(const XMLTag &tagin, int versin)
Definition: LHEF.h:1653
LHEF::TagBase
Definition: LHEF.h:345
LHEF::Scale::emitter
int emitter
Definition: LHEF.h:1394
LHEF::XSecInfo::totxsec
double totxsec
Definition: LHEF.h:574
LHEF::HEPEUP::setWeight
bool setWeight(std::string name, double w)
Definition: LHEF.h:2472
LHEF::HEPEUP::SPINUP
std::vector< double > SPINUP
Definition: LHEF.h:2625
LHEF::Writer::Writer
Writer()
LHEF::XMLTag::name
std::string name
Definition: LHEF.h:119
LHEF::Reader::initfile
std::istream * initfile
Definition: LHEF.h:2979
LHEF::Reader::efile
std::ifstream efile
Definition: LHEF.h:2984
LHEF::XMLTag::~XMLTag
~XMLTag()
Definition: LHEF.h:112
LHEF::EventFile::EventFile
EventFile()
Definition: LHEF.h:622
LHEF::Writer::initStream
std::ostringstream initStream
Definition: LHEF.h:3326
M_PI
#define M_PI
Definition of PI. Needed on some platforms.
Definition: FourVector.h:15
LHEF::Writer::Writer
Writer(std::string filename)
Definition: LHEF.h:3106
LHEF::Scale::print
void print(std::ostream &file) const
Definition: LHEF.h:1355
LHEF::operator<<
std::ostream & operator<<(std::ostream &os, const OAttr< T > &oa)
Definition: LHEF.h:76
LHEF::Weight::indices
std::vector< int > indices
Definition: LHEF.h:1242
LHEF::Generator::version
std::string version
Definition: LHEF.h:504
LHEF::EventGroup::ncounter
int ncounter
Definition: LHEF.h:2104
LHEF::MergeInfo::MergeInfo
MergeInfo(const XMLTag &tag)
Definition: LHEF.h:1002
LHEF::TagBase::getattr
bool getattr(std::string n, double &v, bool erase=true)
Definition: LHEF.h:368
LHEF::HEPEUP::resize
void resize(int nup)
Definition: LHEF.h:2424
LHEF::Writer::openeventfile
bool openeventfile(int ifile)
Definition: LHEF.h:3183
LHEF::HEPRUP::procinfo
std::map< long, ProcInfo > procinfo
Definition: LHEF.h:2015
LHEF::TagBase::AttributeMap
XMLTag::AttributeMap AttributeMap
Definition: LHEF.h:350
LHEF::HEPRUP::print
void print(std::ostream &file) const
Definition: LHEF.h:1765
LHEF::HEPRUP::mergeinfo
std::map< long, MergeInfo > mergeinfo
Definition: LHEF.h:2020
LHEF::TagBase::getattr
bool getattr(std::string n, int &v, bool erase=true)
Definition: LHEF.h:410
LHEF::WeightInfo::muf
double muf
Definition: LHEF.h:1105
LHEF::HEPRUP::ptypes
std::map< std::string, std::set< long > > ptypes
Definition: LHEF.h:2010
LHEF::Scales::muf
double muf
Definition: LHEF.h:1512
LHEF::Reader::currfileevent
int currfileevent
Definition: LHEF.h:3043
LHEF::Reader::curreventfile
int curreventfile
Definition: LHEF.h:3038
LHEF::HEPRUP::NPRUP
int NPRUP
Definition: LHEF.h:1967
LHEF::PDFInfo::PDFInfo
PDFInfo(const XMLTag &tag, double defscale=-1.0)
Definition: LHEF.h:1551
LHEF::HEPEUP::reset
void reset()
Definition: LHEF.h:2404
LHEF::Clus::scale
double scale
Definition: LHEF.h:1299
LHEF::OAttr::name
std::string name
Definition: LHEF.h:55
LHEF::HEPRUP::PDFGUP
std::pair< int, int > PDFGUP
Definition: LHEF.h:1949
LHEF::TagBase::closetag
void closetag(std::ostream &file, std::string tag) const
Definition: LHEF.h:445
LHEF::Generator::name
std::string name
Definition: LHEF.h:499
LHEF::XMLTag::getattr
bool getattr(std::string n, int &v) const
Definition: LHEF.h:174
LHEF::Reader::hepeup
HEPEUP hepeup
Definition: LHEF.h:3022
LHEF::Cut::eta
static double eta(const std::vector< double > &p)
Definition: LHEF.h:832
LHEF::Clus
Definition: LHEF.h:1250
LHEF::TagBase::contents
std::string contents
Definition: LHEF.h:462
LHEF::HEPRUP::getXSecInfo
const XSecInfo & getXSecInfo(std::string weightname="") const
Definition: LHEF.h:1925
LHEF::MergeInfo::print
void print(std::ostream &file) const
Definition: LHEF.h:1013
LHEF::Cut::np2
std::string np2
Definition: LHEF.h:899
LHEF::Writer::lastevent
int lastevent
Definition: LHEF.h:3284
LHEF::Scales::SCALUP
double SCALUP
Definition: LHEF.h:1528
LHEF::XSecInfo
Definition: LHEF.h:511
LHEF::HEPRUP::LPRUP
std::vector< int > LPRUP
Definition: LHEF.h:1989
LHEF::HEPRUP::cuts
std::vector< Cut > cuts
Definition: LHEF.h:2005
LHEF::HEPEUP::currentWeight
const WeightInfo * currentWeight
Definition: LHEF.h:2635
LHEF::HEPEUP::setEvent
HEPEUP & setEvent(const HEPEUP &x)
Definition: LHEF.h:2143
LHEF::MergeInfo
Definition: LHEF.h:992
LHEF::hashline
std::string hashline(std::string s)
Definition: LHEF.h:328
LHEF::XMLTag::contents
std::string contents
Definition: LHEF.h:134
LHEF::Writer::headerBlock
void headerBlock(const std::string &a)
Definition: LHEF.h:3152
LHEF::EventGroup::~EventGroup
~EventGroup()
Definition: LHEF.h:2707
LHEF::Reader::heprup
HEPRUP heprup
Definition: LHEF.h:3012
LHEF::Reader::eventComments
std::string eventComments
Definition: LHEF.h:3027
LHEF::TagBase::yes
static std::string yes()
Definition: LHEF.h:467
LHEF::Reader::Reader
Reader(std::istream &is)
Definition: LHEF.h:2756
LHEF::HEPRUP::operator=
HEPRUP & operator=(const HEPRUP &x)=default
LHEF::ProcInfo::eworder
int eworder
Definition: LHEF.h:970
LHEF::Clus::print
void print(std::ostream &file) const
Definition: LHEF.h:1272
LHEF::HEPEUP
Definition: LHEF.h:2117
LHEF::Writer::eventStream
std::ostringstream eventStream
Definition: LHEF.h:3331
LHEF::Weight::iswgt
bool iswgt
Definition: LHEF.h:1222
LHEF::HEPRUP
Definition: LHEF.h:1627
LHEF::ProcInfo::print
void print(std::ostream &file) const
Definition: LHEF.h:940
LHEF::HEPRUP::weightIndex
int weightIndex(std::string name) const
Definition: LHEF.h:1897
LHEF::Reader::initComments
std::string initComments
Definition: LHEF.h:3017
LHEF::XSecInfo::xsecerr
double xsecerr
Definition: LHEF.h:579
LHEF::Scales::scales
std::vector< Scale > scales
Definition: LHEF.h:1533
LHEF::HEPEUP::subevents
EventGroup subevents
Definition: LHEF.h:2688
LHEF::Writer::writeEvent
void writeEvent()
Definition: LHEF.h:3243
LHEF::WeightGroup::combine
std::string combine
Definition: LHEF.h:1161
LHEF::HEPEUP::HEPEUP
HEPEUP(const XMLTag &tagin, HEPRUP &heprupin)
Definition: LHEF.h:2199
LHEF::Cut::Cut
Cut()
Definition: LHEF.h:674
LHEF::Writer::headerStream
std::ostringstream headerStream
Definition: LHEF.h:3321
LHEF::XSecInfo::XSecInfo
XSecInfo()
Definition: LHEF.h:516
LHEF::HEPEUP::clustering
std::vector< Clus > clustering
Definition: LHEF.h:2650
LHEF::Scale::stype
std::string stype
Definition: LHEF.h:1387
LHEF::HEPEUP::totalWeight
double totalWeight(std::string name) const
Definition: LHEF.h:2445
LHEF::Scales::hasInfo
bool hasInfo() const
Definition: LHEF.h:1452
LHEF::PDFInfo::x1
double x1
Definition: LHEF.h:1590
LHEF::ProcInfo::ProcInfo
ProcInfo(const XMLTag &tag)
Definition: LHEF.h:925
LHEF::TagBase::attributes
XMLTag::AttributeMap attributes
Definition: LHEF.h:457
LHEF::Cut::p1
std::set< long > p1
Definition: LHEF.h:884
LHEF::Reader::Reader
Reader(std::string filename)
Definition: LHEF.h:2772
LHEF::XMLTag
Definition: LHEF.h:87
LHEF::MergeInfo::MergeInfo
MergeInfo()
Definition: LHEF.h:997
LHEF::HEPEUP::NUP
int NUP
Definition: LHEF.h:2550
LHEF::Cut::passCuts
bool passCuts(const std::vector< long > &id, const std::vector< std::vector< double > > &p) const
Definition: LHEF.h:765
LHEF::Weight::Weight
Weight()
Definition: LHEF.h:1174
LHEF::ProcInfo::scheme
std::string scheme
Definition: LHEF.h:985
LHEF::Reader::openeventfile
void openeventfile(int ifile)
Definition: LHEF.h:2932
LHEF::PDFInfo::p2
long p2
Definition: LHEF.h:1585
LHEF::OAttr
Definition: LHEF.h:45
LHEF::HEPRUP::resize
void resize()
Definition: LHEF.h:1887
LHEF::XSecInfo::neve
long neve
Definition: LHEF.h:564
LHEF::Scales
Definition: LHEF.h:1416
LHEF::HEPRUP::eventfiles
std::vector< EventFile > eventfiles
Definition: LHEF.h:2000
LHEF::Writer::Writer
Writer(std::ostream &os)
Definition: LHEF.h:3099
LHEF::HEPRUP::junk
std::string junk
Definition: LHEF.h:2046
LHEF::HEPEUP::namedweights
std::vector< Weight > namedweights
Definition: LHEF.h:2640
LHEF::HEPRUP::IDBMUP
std::pair< long, long > IDBMUP
Definition: LHEF.h:1938
LHEF::Writer::operator=
Writer & operator=(const Writer &)
LHEF::Scales::print
void print(std::ostream &file) const
Definition: LHEF.h:1460
LHEF::Cut::min
double min
Definition: LHEF.h:904
LHEF::HEPEUP::IDPRUP
int IDPRUP
Definition: LHEF.h:2555
LHEF::XSecInfo::ntries
long ntries
Definition: LHEF.h:569
LHEF::HEPEUP::weight
double weight(std::string name) const
Definition: LHEF.h:2459
LHEF::Clus::alphas
double alphas
Definition: LHEF.h:1305
LHEF::EventGroup::EventGroup
EventGroup()
Definition: LHEF.h:2074
LHEF::Scale
Definition: LHEF.h:1313
LHEF::Reader::headerBlock
std::string headerBlock
Definition: LHEF.h:3007
LHEF::Writer::eventComments
void eventComments(const std::string &a)
Definition: LHEF.h:3166
LHEF::WeightInfo::print
void print(std::ostream &file) const
Definition: LHEF.h:1070
LHEF::ProcInfo::iproc
int iproc
Definition: LHEF.h:955
LHEF::Clus::Clus
Clus()
Definition: LHEF.h:1255
LHEF::XMLTag::getattr
bool getattr(std::string n, bool &v) const
Definition: LHEF.h:152
LHEF::HEPEUP::isGroup
bool isGroup
Definition: LHEF.h:2682
LHEF::HEPEUP::VTIMUP
std::vector< double > VTIMUP
Definition: LHEF.h:2618
LHEF::Writer::writeinit
void writeinit()
Definition: LHEF.h:3212
LHEF::Scale::scale
double scale
Definition: LHEF.h:1409
LHEF::HEPEUP::ntries
int ntries
Definition: LHEF.h:2677
LHEF::TagBase::TagBase
TagBase(const AttributeMap &attr, std::string conts=std::string())
Definition: LHEF.h:360
LHEF::Weight::Weight
Weight(const XMLTag &tag)
Definition: LHEF.h:1179
LHEF::TagBase::getattr
bool getattr(std::string n, long &v, bool erase=true)
Definition: LHEF.h:396
LHEF::Cut
Definition: LHEF.h:669
LHEF::PDFInfo
Definition: LHEF.h:1540
LHEF::Reader::version
int version
Definition: LHEF.h:2996
LHEF::Weight
Definition: LHEF.h:1169
LHEF::Scale::Scale
Scale(std::string st="veto", int emtr=0, double sc=0.0)
Definition: LHEF.h:1318
LHEF::HEPEUP::IDUP
std::vector< long > IDUP
Definition: LHEF.h:2589
LHEF::EventFile::EventFile
EventFile(const XMLTag &tag)
Definition: LHEF.h:627
LHEF::OAttr::OAttr
OAttr(std::string n, const T &v)
Definition: LHEF.h:50
LHEF::Scales::mur
double mur
Definition: LHEF.h:1517
LHEF::Cut::outside
bool outside(double value) const
Definition: LHEF.h:872
LHEF::XMLTag::getattr
bool getattr(std::string n, std::string &v) const
Definition: LHEF.h:185
LHEF::HEPRUP::dprec
int dprec
Definition: LHEF.h:2056
LHEF::Reader
Definition: LHEF.h:2742
LHEF::Reader::file
std::istream * file
Definition: LHEF.h:2974
LHEF::XSecInfo::weightname
std::string weightname
Definition: LHEF.h:604
LHEF::Writer::Writer
Writer(const Writer &)
LHEF::Weight::born
double born
Definition: LHEF.h:1227
LHEF::Reader::init
void init()
Definition: LHEF.h:2787
LHEF::MergeInfo::iproc
int iproc
Definition: LHEF.h:1024
LHEF::HEPEUP::junk
std::string junk
Definition: LHEF.h:2693
LHEF::HEPRUP::getXSecInfo
XSecInfo & getXSecInfo(std::string weightname="")
Definition: LHEF.h:1914
LHEF::HEPRUP::weightinfo
std::vector< WeightInfo > weightinfo
Definition: LHEF.h:2031
LHEF::XSecInfo::varweights
bool varweights
Definition: LHEF.h:599
LHEF::HEPEUP::ISTUP
std::vector< int > ISTUP
Definition: LHEF.h:2594
LHEF::HEPRUP::clear
void clear()
Definition: LHEF.h:1862
LHEF::HEPEUP::heprup
HEPRUP * heprup
Definition: LHEF.h:2630
LHEF::MergeInfo::maxmult
bool maxmult
Definition: LHEF.h:1034
LHEF::oattr
OAttr< T > oattr(std::string name, const T &value)
Definition: LHEF.h:68
LHEF::Scales::getScale
double getScale(std::string st, int pdgem, int emr, int rec) const
Definition: LHEF.h:1489
LHEF::Writer::initComments
std::ostream & initComments()
Definition: LHEF.h:3138
LHEF::XMLTag::findXMLTags
static std::vector< XMLTag * > findXMLTags(std::string str, std::string *leftover=0)
Definition: LHEF.h:198
LHEF::WeightGroup::WeightGroup
WeightGroup()
Definition: LHEF.h:1133
LHEF::ProcInfo::loops
int loops
Definition: LHEF.h:960
LHEF::EventFile::ntries
long ntries
Definition: LHEF.h:662
LHEF::HEPRUP::PDFSUP
std::pair< int, int > PDFSUP
Definition: LHEF.h:1955
LHEF::PDFInfo::SCALUP
double SCALUP
Definition: LHEF.h:1615
LHEF::HEPEUP::totalWeight
double totalWeight(int i=0) const
Definition: LHEF.h:2433
LHEF::WeightInfo::WeightInfo
WeightInfo()
Definition: LHEF.h:1047
LHEF::EventFile::neve
long neve
Definition: LHEF.h:657
LHEF::XMLTag::deleteAll
static void deleteAll(std::vector< XMLTag * > &tags)
Definition: LHEF.h:293
LHEF::HEPEUP::AQEDUP
double AQEDUP
Definition: LHEF.h:2579
LHEF::XMLTag::print
void print(std::ostream &os) const
Definition: LHEF.h:302
LHEF::Scale::emitted
std::set< int > emitted
Definition: LHEF.h:1404
LHEF::WeightInfo::name
std::string name
Definition: LHEF.h:1100
LHEF::HEPRUP::HEPRUP
HEPRUP(const HEPRUP &)=default
LHEF::WeightInfo::isrwgt
bool isrwgt
Definition: LHEF.h:1095
LHEF::Generator::print
void print(std::ostream &file) const
Definition: LHEF.h:488
LHEF::Writer::intstream
std::ofstream intstream
Definition: LHEF.h:3263
LHEF::ProcInfo::ProcInfo
ProcInfo()
Definition: LHEF.h:920
LHEF::HEPEUP::setSubEvent
bool setSubEvent(unsigned int i)
Definition: LHEF.h:2530
LHEF::EventGroup::operator=
EventGroup & operator=(const EventGroup &)
Definition: LHEF.h:2716
LHEF::Clus::p0
int p0
Definition: LHEF.h:1294
LHEF::XMLTag::attr
AttributeMap attr
Definition: LHEF.h:124
LHEF::HEPEUP::AQCDUP
double AQCDUP
Definition: LHEF.h:2584
LHEF::EventFile::filename
std::string filename
Definition: LHEF.h:652
LHEF::Writer::file
std::ostream * file
Definition: LHEF.h:3269
LHEF::HEPEUP::weight
double weight(int i=0) const
Definition: LHEF.h:2452
LHEF::HEPRUP::version
int version
Definition: LHEF.h:2051
LHEF::HEPEUP::scales
Scales scales
Definition: LHEF.h:2671
LHEF::ProcInfo::qcdorder
int qcdorder
Definition: LHEF.h:965
LHEF::Generator::Generator
Generator(const XMLTag &tag)
Definition: LHEF.h:479
LHEF::XMLTag::tags
std::vector< XMLTag * > tags
Definition: LHEF.h:129
LHEF::Weight::name
std::string name
Definition: LHEF.h:1217
LHEF::HEPEUP::setWeightInfo
bool setWeightInfo(unsigned int i)
Definition: LHEF.h:2499
LHEF::Generator
Definition: LHEF.h:474
LHEF::PDFInfo::xf1
double xf1
Definition: LHEF.h:1600
LHEF::Reader::Reader
Reader()
LHEF::HEPRUP::xsecinfos
XSecInfos xsecinfos
Definition: LHEF.h:1994
LHEF::Reader::currentFind
bool currentFind(std::string str) const
Definition: LHEF.h:2957
LHEF::HEPRUP::XSECUP
std::vector< double > XSECUP
Definition: LHEF.h:1972
LHEF::HEPEUP::PDFGUPsave
std::pair< int, int > PDFGUPsave
Definition: LHEF.h:2660
LHEF::HEPRUP::weightNameHepMC
std::string weightNameHepMC(int i) const
Definition: LHEF.h:1751
LHEF::OAttr::val
T val
Definition: LHEF.h:60
LHEF::Reader::dirpath
std::string dirpath
Definition: LHEF.h:3048
LHEF::HEPRUP::XERRUP
std::vector< double > XERRUP
Definition: LHEF.h:1978
LHEF::PDFInfo::PDFInfo
PDFInfo(double defscale=-1.0)
Definition: LHEF.h:1545
LHEF::Writer::heprup
HEPRUP heprup
Definition: LHEF.h:3306
LHEF::TagBase::getattr
bool getattr(std::string n, std::string &v, bool erase=true)
Definition: LHEF.h:424
LHEF::HEPEUP::weights
std::vector< std::pair< double, const WeightInfo * > > weights
Definition: LHEF.h:2645
LHEF::XMLTag::getattr
bool getattr(std::string n, double &v) const
Definition: LHEF.h:140
LHEF::Writer::~Writer
~Writer()
Definition: LHEF.h:3116
LHEF::HEPEUP::print
void print(std::ostream &file) const
Definition: LHEF.h:2320
LHEF::WeightInfo::mur
double mur
Definition: LHEF.h:1110
LHEF::Weight::print
void print(std::ostream &file) const
Definition: LHEF.h:1197
LHEF::ProcInfo::fscheme
std::string fscheme
Definition: LHEF.h:975
LHEF
Les Houches event file classes.
Definition: LHEF.h:39
LHEF::Clus::Clus
Clus(const XMLTag &tag)
Definition: LHEF.h:1260
LHEF::HEPRUP::EBMUP
std::pair< double, double > EBMUP
Definition: LHEF.h:1943
LHEF::HEPRUP::weightgroup
std::vector< WeightGroup > weightgroup
Definition: LHEF.h:2041
LHEF::WeightInfo::WeightInfo
WeightInfo(const XMLTag &tag)
Definition: LHEF.h:1053
LHEF::Cut::Cut
Cut(const XMLTag &tag, const std::map< std::string, std::set< long > > &ptypes)
Definition: LHEF.h:680
LHEF::EventGroup
Definition: LHEF.h:2069
LHEF::HEPEUP::resize
void resize()
Definition: LHEF.h:2485
LHEF::Scales::Scales
Scales(const XMLTag &tag, double defscale=-1.0, int npart=0)
Definition: LHEF.h:1429
LHEF::Cut::rap
static double rap(const std::vector< double > &p)
Definition: LHEF.h:846
LHEF::Writer::hepeup
HEPEUP hepeup
Definition: LHEF.h:3312
LHEF::Cut::max
double max
Definition: LHEF.h:908
LHEF::HEPEUP::HEPEUP
HEPEUP()
Definition: LHEF.h:2126
LHEF::Weight::sudakov
double sudakov
Definition: LHEF.h:1232
LHEF::XMLTag::XMLTag
XMLTag()
Definition: LHEF.h:107
LHEF::Scale::recoilers
std::set< int > recoilers
Definition: LHEF.h:1399
LHEF::HEPEUP::PDFSUPsave
std::pair< int, int > PDFSUPsave
Definition: LHEF.h:2665
LHEF::Scale::Scale
Scale(const XMLTag &tag)
Definition: LHEF.h:1324
LHEF::HEPEUP::ICOLUP
std::vector< std::pair< int, int > > ICOLUP
Definition: LHEF.h:2606
LHEF::HEPEUP::XPDWUP
std::pair< double, double > XPDWUP
Definition: LHEF.h:2568
LHEF::HEPEUP::clear
void clear()
Definition: LHEF.h:2414
LHEF::TagBase::printattrs
void printattrs(std::ostream &file) const
Definition: LHEF.h:435
LHEF::HEPEUP::~HEPEUP
~HEPEUP()
Definition: LHEF.h:2188
LHEF::Scales::mups
double mups
Definition: LHEF.h:1523
LHEF::Writer::headerBlock
std::ostream & headerBlock()
Definition: LHEF.h:3131
LHEF::XSecInfo::XSecInfo
XSecInfo(const XMLTag &tag)
Definition: LHEF.h:522
LHEF::HEPEUP::SCALUP
double SCALUP
Definition: LHEF.h:2574
LHEF::Clus::p2
int p2
Definition: LHEF.h:1289
LHEF::Cut::print
void print(std::ostream &file) const
Definition: LHEF.h:722
LHEF::PDFInfo::p1
long p1
Definition: LHEF.h:1580
LHEF::XSecInfo::print
void print(std::ostream &file) const
Definition: LHEF.h:546
LHEF::XMLTag::pos_t
std::string::size_type pos_t
Definition: LHEF.h:92
LHEF::Reader::readEvent
bool readEvent()
Definition: LHEF.h:2868
LHEF::HEPEUP::PUP
std::vector< std::vector< double > > PUP
Definition: LHEF.h:2612
LHEF::XSecInfo::maxweight
double maxweight
Definition: LHEF.h:584