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 name = 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[name] = 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(): 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  using std::setw;
1768  file << std::setprecision(dprec);
1769 
1770  file << "<init>\n"
1771  << " " << setw(8) << IDBMUP.first
1772  << " " << setw(8) << IDBMUP.second
1773  << " " << setw(14) << EBMUP.first
1774  << " " << setw(14) << EBMUP.second
1775  << " " << setw(4) << PDFGUP.first
1776  << " " << setw(4) << PDFGUP.second
1777  << " " << setw(4) << PDFSUP.first
1778  << " " << setw(4) << PDFSUP.second
1779  << " " << setw(4) << IDWTUP
1780  << " " << setw(4) << NPRUP << std::endl;
1781 
1782  for ( int i = 0; i < NPRUP; ++i )
1783  file << " " << setw(14) << XSECUP[i]
1784  << " " << setw(14) << XERRUP[i]
1785  << " " << setw(14) << XMAXUP[i]
1786  << " " << setw(6) << LPRUP[i] << std::endl;
1787 
1788  for ( int i = 0, N = generators.size(); i < N; ++i )
1789  generators[i].print(file);
1790 
1791  if ( !eventfiles.empty() ) {
1792  file << "<eventfiles>\n";
1793  for ( int i = 0, N = eventfiles.size(); i < N; ++i )
1794  eventfiles[i].print(file);
1795  file << "</eventfiles>\n";
1796  }
1797  //AV if ( !xsecinfos.empty() > 0 )
1798  if ( !xsecinfos.empty())
1799  for ( XSecInfos::const_iterator it = xsecinfos.begin();
1800  it != xsecinfos.end(); ++it )
1801  if ( it->second.neve > 0 ) it->second.print(file);
1802 
1803  if ( cuts.size() > 0 ) {
1804  file << "<cutsinfo>" << std::endl;
1805 
1806  for ( std::map<std::string, std::set<long> >::const_iterator ptit =
1807  ptypes.begin(); ptit != ptypes.end(); ++ptit ) {
1808  file << "<ptype" << oattr("name", ptit->first) << ">";
1809  for ( std::set<long>::const_iterator it = ptit->second.begin();
1810  it != ptit->second.end(); ++it )
1811  file << " " << *it;
1812  file << "</ptype>" << std::endl;
1813  }
1814 
1815  for ( int i = 0, N = cuts.size(); i < N; ++i )
1816  cuts[i].print(file);
1817  file << "</cutsinfo>" << std::endl;
1818  }
1819 
1820  for ( std::map<long,ProcInfo>::const_iterator it = procinfo.begin();
1821  it != procinfo.end(); ++it )
1822  it->second.print(file);
1823 
1824  for ( std::map<long,MergeInfo>::const_iterator it = mergeinfo.begin();
1825  it != mergeinfo.end(); ++it )
1826  it->second.print(file);
1827 
1828  bool isrwgt = false;
1829  int ingroup = -1;
1830  for ( int i = 0, N = weightinfo.size(); i < N; ++i ) {
1831  if ( weightinfo[i].isrwgt ) {
1832  if ( !isrwgt ) file << "<initrwgt>\n";
1833  isrwgt = true;
1834  } else {
1835  if ( isrwgt ) file << "</initrwgt>\n";
1836  isrwgt = false;
1837  }
1838  int group = weightinfo[i].inGroup;
1839  if ( group != ingroup ) {
1840  if ( ingroup != -1 ) file << "</weightgroup>\n";
1841  if ( group != -1 ) {
1842  file << "<weightgroup"
1843  << oattr("type", weightgroup[group].type);
1844  if ( !weightgroup[group].combine.empty() )
1845  file << oattr("combine", weightgroup[group].combine);
1846  file << ">\n";
1847  }
1848  ingroup = group;
1849  }
1850  weightinfo[i].print(file);
1851  }
1852  if ( ingroup != -1 ) file << "</weightgroup>\n";
1853  if ( isrwgt ) file << "</initrwgt>\n";
1854 
1855 
1856  file << hashline(junk) << "</init>" << std::endl;
1857 
1858  }
1859 
1860  /**
1861  * Clear all information.
1862  */
1863  void clear() {
1864  procinfo.clear();
1865  mergeinfo.clear();
1866  weightinfo.clear();
1867  weightgroup.clear();
1868  cuts.clear();
1869  ptypes.clear();
1870  junk.clear();
1871  }
1872 
1873  /**
1874  * Set the NPRUP variable, corresponding to the number of
1875  * sub-processes, to \a nrup, and resize all relevant vectors
1876  * accordingly.
1877  */
1878  void resize(int nrup) {
1879  NPRUP = nrup;
1880  resize();
1881  }
1882 
1883  /**
1884  * Assuming the NPRUP variable, corresponding to the number of
1885  * sub-processes, is correctly set, resize the relevant vectors
1886  * accordingly.
1887  */
1888  void resize() {
1889  XSECUP.resize(NPRUP);
1890  XERRUP.resize(NPRUP);
1891  XMAXUP.resize(NPRUP);
1892  LPRUP.resize(NPRUP);
1893  }
1894 
1895  /**
1896  * @return the index of the weight with the given \a name
1897  */
1898  int weightIndex(std::string name) const {
1899  std::map<std::string, int>::const_iterator it = weightmap.find(name);
1900  if ( it != weightmap.end() ) return it->second;
1901  return 0;
1902  }
1903 
1904  /**
1905  * @return the number of weights (including the nominial one).
1906  */
1907  int nWeights() const {
1908  return weightmap.size() + 1;
1909  }
1910 
1911  /**
1912  * @return the XSecInfo object corresponding to the named weight \a
1913  * weithname. If no such object exists, it will be created.
1914  */
1915  XSecInfo & getXSecInfo(std::string weightname = "") {
1916  XSecInfo & xi = xsecinfos[weightname];
1917  xi.weightname = weightname;
1918  return xi;
1919  }
1920 
1921  /**
1922  * @return the XSecInfo object corresponding to the named weight \a
1923  * weithname. If no such object exists, an empty XSecInfo will be
1924  * returned..
1925  */
1926  const XSecInfo & getXSecInfo(std::string weightname = "") const {
1927  static XSecInfo noinfo;
1928  XSecInfos::const_iterator it = xsecinfos.find(weightname);
1929  if ( it != xsecinfos.end() ) return it->second;
1930  else return noinfo;
1931  }
1932 
1933 
1934 public:
1935 
1936  /**
1937  * PDG id's of beam particles. (first/second is in +/-z direction).
1938  */
1939  std::pair<long,long> IDBMUP;
1940 
1941  /**
1942  * Energy of beam particles given in GeV.
1943  */
1944  std::pair<double,double> EBMUP;
1945 
1946  /**
1947  * The author group for the PDF used for the beams according to the
1948  * PDFLib specification.
1949  */
1950  std::pair<int,int> PDFGUP;
1951 
1952  /**
1953  * The id number the PDF used for the beams according to the
1954  * PDFLib specification.
1955  */
1956  std::pair<int,int> PDFSUP;
1957 
1958  /**
1959  * Master switch indicating how the ME generator envisages the
1960  * events weights should be interpreted according to the Les Houches
1961  * accord.
1962  */
1963  int IDWTUP;
1964 
1965  /**
1966  * The number of different subprocesses in this file.
1967  */
1968  int NPRUP;
1969 
1970  /**
1971  * The cross sections for the different subprocesses in pb.
1972  */
1973  std::vector<double> XSECUP;
1974 
1975  /**
1976  * The statistical error in the cross sections for the different
1977  * subprocesses in pb.
1978  */
1979  std::vector<double> XERRUP;
1980 
1981  /**
1982  * The maximum event weights (in HEPEUP::XWGTUP) for different
1983  * subprocesses.
1984  */
1985  std::vector<double> XMAXUP;
1986 
1987  /**
1988  * The subprocess code for the different subprocesses.
1989  */
1990  std::vector<int> LPRUP;
1991 
1992  /**
1993  * Contents of the xsecinfo tags.
1994  */
1996 
1997  /**
1998  * A vector of EventFiles where the events are stored separate fron
1999  * the init block.
2000  */
2001  std::vector<EventFile> eventfiles;
2002 
2003  /**
2004  * Contents of the cuts tag.
2005  */
2006  std::vector<Cut> cuts;
2007 
2008  /**
2009  * A map of codes for different particle types.
2010  */
2011  std::map<std::string, std::set<long> > ptypes;
2012 
2013  /**
2014  * Contents of the procinfo tags
2015  */
2016  std::map<long,ProcInfo> procinfo;
2017 
2018  /**
2019  * Contents of the mergeinfo tags
2020  */
2021  std::map<long,MergeInfo> mergeinfo;
2022 
2023  /**
2024  * The names of the programs and their version information used to
2025  * create this file.
2026  */
2027  std::vector<Generator> generators;
2028 
2029  /**
2030  * The vector of WeightInfo objects for this file.
2031  */
2032  std::vector<WeightInfo> weightinfo;
2033 
2034  /**
2035  * A map relating names of weights to indices of the weightinfo vector.
2036  */
2037  std::map<std::string,int> weightmap;
2038 
2039  /**
2040  * The vector of WeightGroup objects in this file.
2041  */
2042  std::vector<WeightGroup> weightgroup;
2043 
2044  /**
2045  * Just to be on the safe side we save any junk inside the init-tag.
2046  */
2047  std::string junk;
2048 
2049  /**
2050  * The main version of the information stored.
2051  */
2052  int version;
2053 
2054  /**
2055  * The precision used for outputing real numbers.
2056  */
2057  int dprec;
2058 
2059 };
2060 
2061 /**
2062  * Forward declaration of the HEPEUP class.
2063  */
2064 class HEPEUP;
2065 
2066 /**
2067  * The EventGroup represents a set of events which are to be
2068  * considered together.
2069  */
2070 struct EventGroup: public std::vector<HEPEUP*> {
2071 
2072  /**
2073  * Initialize default values.
2074  */
2075  inline EventGroup(): nreal(-1), ncounter(-1) {}
2076 
2077  /**
2078  * The copy constructor also copies the included HEPEUP object.
2079  */
2080  inline EventGroup(const EventGroup &);
2081 
2082  /**
2083  * The assignment also copies the included HEPEUP object.
2084  */
2085  inline EventGroup & operator=(const EventGroup &);
2086 
2087  /**
2088  * Remove all subevents.
2089  */
2090  inline void clear();
2091 
2092  /**
2093  * The destructor deletes the included HEPEUP objects.
2094  */
2095  inline ~EventGroup();
2096 
2097  /**
2098  * The number of real events in this event group.
2099  */
2100  int nreal;
2101 
2102  /**
2103  * The number of counter events in this event group.
2104  */
2106 
2107 };
2108 
2109 
2110 /**
2111  * The HEPEUP class is a simple container corresponding to the Les Houches accord
2112  * (<A HREF="http://arxiv.org/abs/hep-ph/0109068">hep-ph/0109068</A>)
2113  * common block with the same name. The members are named in the same
2114  * way as in the common block. However, fortran arrays are represented
2115  * by vectors, except for the arrays of length two which are
2116  * represented by pair objects.
2117  */
2118 class HEPEUP : public TagBase {
2119 
2120 public:
2121 
2122  /** @name Standard constructors and destructors. */
2123  //@{
2124  /**
2125  * Default constructor.
2126  */
2128  : NUP(0), IDPRUP(0), XWGTUP(0.0), XPDWUP(0.0, 0.0),
2129  SCALUP(0.0), AQEDUP(0.0), AQCDUP(0.0), heprup(0), currentWeight(0),
2130  ntries(1), isGroup(false) {}
2131 
2132  /**
2133  * Copy constructor
2134  */
2135  HEPEUP(const HEPEUP & x)
2136  : TagBase(x), isGroup(false) {
2137  operator=(x);
2138  }
2139 
2140  /**
2141  * Copy information from the given HEPEUP. Sub event information is
2142  * left untouched.
2143  */
2144  HEPEUP & setEvent(const HEPEUP & x) {
2145  NUP = x.NUP;
2146  IDPRUP = x.IDPRUP;
2147  XWGTUP = x.XWGTUP;
2148  XPDWUP = x.XPDWUP;
2149  SCALUP = x.SCALUP;
2150  AQEDUP = x.AQEDUP;
2151  AQCDUP = x.AQCDUP;
2152  IDUP = x.IDUP;
2153  ISTUP = x.ISTUP;
2154  MOTHUP = x.MOTHUP;
2155  ICOLUP = x.ICOLUP;
2156  PUP = x.PUP;
2157  VTIMUP = x.VTIMUP;
2158  SPINUP = x.SPINUP;
2159  heprup = x.heprup;
2161  weights = x.weights;
2162  pdfinfo = x.pdfinfo;
2163  PDFGUPsave = x.PDFGUPsave;
2164  PDFSUPsave = x.PDFSUPsave;
2165  clustering = x.clustering;
2166  scales = x.scales;
2167  junk = x.junk;
2169  ntries = x.ntries;
2170  return *this;
2171  }
2172 
2173  /**
2174  * Assignment operator.
2175  */
2176  HEPEUP & operator=(const HEPEUP & x) {
2177  if ( &x == this ) return *this;
2178  TagBase::operator=(x);
2179  clear();
2180  setEvent(x);
2181  subevents = x.subevents;
2182  isGroup = x.isGroup;
2183  return *this;
2184  }
2185 
2186  /**
2187  * Destructor.
2188  */
2190  clear();
2191  };
2192  //@}
2193 
2194 public:
2195 
2196 
2197  /**
2198  * Constructor from an event or eventgroup tag.
2199  */
2200  HEPEUP(const XMLTag & tagin, HEPRUP & heprupin)
2201  : TagBase(tagin.attr), NUP(0), IDPRUP(0), XWGTUP(0.0), XPDWUP(0.0, 0.0),
2202  SCALUP(0.0), AQEDUP(0.0), AQCDUP(0.0), heprup(&heprupin),
2203  currentWeight(0), ntries(1), isGroup(tagin.name == "eventgroup") {
2204 
2205  if ( heprup->NPRUP < 0 )
2206  throw std::runtime_error("Tried to read events but no processes defined "
2207  "in init block of Les Houches file.");
2208 
2209  std::vector<XMLTag*> tags = tagin.tags;
2210 
2211  if ( isGroup ) {
2212  getattr("nreal", subevents.nreal);
2213  getattr("ncounter", subevents.ncounter);
2214  for ( int i = 0, N = tags.size(); i < N; ++i )
2215  if ( tags[i]->name == "event" )
2216  subevents.push_back(new HEPEUP(*tags[i], heprupin));
2217  return;
2218  } else
2219  getattr("ntries", ntries);
2220 
2221 
2222 
2223  // The event information should be in the first (anonymous) tag
2224  std::istringstream iss(tags[0]->contents);
2225  if ( !( iss >> NUP >> IDPRUP >> XWGTUP >> SCALUP >> AQEDUP >> AQCDUP ) )
2226  throw std::runtime_error("Failed to parse event in Les Houches file.");
2227 
2228  resize();
2229 
2230  // Read all particle lines.
2231  for ( int i = 0; i < NUP; ++i ) {
2232  if ( !( iss >> IDUP[i] >> ISTUP[i] >> MOTHUP[i].first >> MOTHUP[i].second
2233  >> ICOLUP[i].first >> ICOLUP[i].second
2234  >> PUP[i][0] >> PUP[i][1] >> PUP[i][2]
2235  >> PUP[i][3] >> PUP[i][4]
2236  >> VTIMUP[i] >> SPINUP[i] ) )
2237  throw std::runtime_error("Failed to parse event in Les Houches file.");
2238  }
2239 
2240  junk.clear();
2241  std::string ss;
2242  while ( getline(iss, ss) ) junk += ss + '\n';
2243 
2244  scales = Scales(SCALUP, NUP);
2245  pdfinfo = PDFInfo(SCALUP);
2246  namedweights.clear();
2247  weights.clear();
2248  weights.resize(heprup->nWeights(),
2249  std::make_pair(XWGTUP, (WeightInfo*)(0)));
2250  weights.front().first = XWGTUP;
2251  for ( int i = 1, N = weights.size(); i < N; ++i )
2252  weights[i].second = &heprup->weightinfo[i - 1];
2253 
2254  for ( int i = 1, N = tags.size(); i < N; ++i ) {
2255  XMLTag & tag = *tags[i];
2256 
2257  if ( tag.name.empty() ) junk += tag.contents;
2258 
2259  if ( tag.name == "weights" ) {
2260  weights.resize(heprup->nWeights(),
2261  std::make_pair(XWGTUP, (WeightInfo*)(0)));
2262  weights.front().first = XWGTUP;
2263  for ( int ii = 1, NN = weights.size(); ii < NN; ++ii )
2264  weights[ii].second = &heprup->weightinfo[ii - 1];
2265  double w = 0.0;
2266  int iii = 0;
2267  std::istringstream isss(tag.contents);
2268  while ( isss >> w )
2269  if ( ++iii < int(weights.size()) )
2270  weights[iii].first = w;
2271  else
2272  weights.push_back(std::make_pair(w, (WeightInfo*)(0)));
2273  }
2274  if ( tag.name == "weight" ) {
2275  namedweights.push_back(Weight(tag));
2276  }
2277  if ( tag.name == "rwgt" ) {
2278  for ( int j = 0, M = tag.tags.size(); j < M; ++j ) {
2279  if ( tag.tags[j]->name == "wgt" ) {
2280  namedweights.push_back(Weight(*tag.tags[j]));
2281  }
2282  }
2283  }
2284  else if ( tag.name == "clustering" ) {
2285  for ( int j = 0, M= tag.tags.size(); j < M; ++j ) {
2286  if ( tag.tags[j]->name == "clus" )
2287  clustering.push_back(Clus(*tag.tags[j]));
2288  }
2289  }
2290  else if ( tag.name == "pdfinfo" ) {
2291  pdfinfo = PDFInfo(tag, SCALUP);
2292  }
2293  else if ( tag.name == "scales" ) {
2294  scales = Scales(tag, SCALUP, NUP);
2295  }
2296 
2297  }
2298 
2299  for ( int i = 0, N = namedweights.size(); i < N; ++i ) {
2300  int indx = heprup->weightIndex(namedweights[i].name);
2301  if ( indx > 0 ) {
2302  weights[indx].first = namedweights[i].weights[0];
2303  namedweights[i].indices[0] = indx;
2304  } else {
2305  weights.push_back(std::make_pair(namedweights[i].weights[0],
2306  (WeightInfo*)(0)));
2307  namedweights[i].indices[0] = weights.size() - 1;
2308  }
2309  for ( int j = 1, M = namedweights[i].weights.size(); j < M; ++j ) {
2310  weights.push_back(std::make_pair(namedweights[i].weights[j],
2311  (WeightInfo*)(0)));
2312  namedweights[i].indices[j] = weights.size() - 1;
2313  }
2314  }
2315 
2316  }
2317 
2318  /**
2319  * Print out the event (group) as an XML tag.
2320  */
2321  void print(std::ostream & file) const {
2322 
2323  file << std::setprecision(heprup->dprec);
2324 
2325  using std::setw;
2326 
2327  if ( isGroup ) {
2328  file << "<eventgroup";
2329  if ( subevents.nreal > 0 )
2330  file << oattr("nreal", subevents.nreal);
2331  if ( subevents.ncounter > 0 )
2332  file << oattr("ncounter", subevents.ncounter);
2333  printattrs(file);
2334  file << ">\n";
2335  for ( int i = 0, N = subevents.size(); i < N; ++i )
2336  subevents[i]->print(file);
2337  file << "</eventgroup>\n";
2338  return;
2339  }
2340 
2341  file << "<event";
2342  if ( ntries > 1 ) file << oattr("ntries", ntries);
2343  printattrs(file);
2344  file << ">\n";
2345  file << " " << setw(4) << NUP
2346  << " " << setw(6) << IDPRUP
2347  << " " << setw(14) << XWGTUP
2348  << " " << setw(14) << SCALUP
2349  << " " << setw(14) << AQEDUP
2350  << " " << setw(14) << AQCDUP << "\n";
2351 
2352  for ( int i = 0; i < NUP; ++i )
2353  file << " " << setw(8) << IDUP[i]
2354  << " " << setw(2) << ISTUP[i]
2355  << " " << setw(4) << MOTHUP[i].first
2356  << " " << setw(4) << MOTHUP[i].second
2357  << " " << setw(4) << ICOLUP[i].first
2358  << " " << setw(4) << ICOLUP[i].second
2359  << " " << setw(14) << PUP[i][0]
2360  << " " << setw(14) << PUP[i][1]
2361  << " " << setw(14) << PUP[i][2]
2362  << " " << setw(14) << PUP[i][3]
2363  << " " << setw(14) << PUP[i][4]
2364  << " " << setw(1) << VTIMUP[i]
2365  << " " << setw(1) << SPINUP[i] << std::endl;
2366 
2367  if ( weights.size() > 0 ) {
2368  file << "<weights>";
2369  for ( int i = 1, N = weights.size(); i < N; ++i )
2370  file << " " << weights[i].first;
2371  file << "</weights>\n";
2372  }
2373 
2374  bool iswgt = false;
2375  for ( int i = 0, N = namedweights.size(); i < N; ++i ) {
2376  if ( namedweights[i].iswgt ) {
2377  if ( !iswgt ) file << "<rwgt>\n";
2378  iswgt = true;
2379  } else {
2380  if ( iswgt ) file << "</rwgt>\n";
2381  iswgt = false;
2382  }
2383  for ( int j = 0, M = namedweights[i].indices.size(); j < M; ++j )
2384  namedweights[i].weights[j] = weight(namedweights[i].indices[j]);
2385  namedweights[i].print(file);
2386  }
2387  if ( iswgt ) file << "</rwgt>\n";
2388 
2389  if ( !clustering.empty() ) {
2390  file << "<clustering>" << std::endl;
2391  for ( int i = 0, N = clustering.size(); i < N; ++i )
2392  clustering[i].print(file);
2393  file << "</clustering>" << std::endl;
2394  }
2395 
2396  pdfinfo.print(file);
2397  scales.print(file);
2398 
2399  file << hashline(junk) << "</event>\n";
2400 
2401  }
2402 
2403  /**
2404  * Reset the HEPEUP object (does not touch the sub events).
2405  */
2406  void reset() {
2407  setWeightInfo(0);
2408  NUP = 0;
2409  clustering.clear();
2410  weights.clear();
2411  }
2412 
2413  /**
2414  * Clear the HEPEUP object.
2415  */
2416  void clear() {
2417  reset();
2418  subevents.clear();
2419  }
2420 
2421  /**
2422  * Set the NUP variable, corresponding to the number of particles in
2423  * the current event, to \a nup, and resize all relevant vectors
2424  * accordingly.
2425  */
2426  void resize(int nup) {
2427  NUP = nup;
2428  resize();
2429  }
2430 
2431  /**
2432  * Return the total weight for this event (including all sub
2433  * evenets) for the given index.
2434  */
2435  double totalWeight(int i = 0) const {
2436  if ( subevents.empty() ) return weight(i);
2437  double w = 0.0;
2438  for ( int ii = 0, N = subevents.size(); ii < N; ++ii )
2439  w += subevents[ii]->weight(i);
2440  return w;
2441  }
2442 
2443  /**
2444  * Return the total weight for this event (including all sub
2445  * evenets) for the given weight name.
2446  */
2447  double totalWeight(std::string name) const {
2448  return totalWeight(heprup->weightIndex(name));
2449  }
2450 
2451  /**
2452  * Return the weight for the given index.
2453  */
2454  double weight(int i = 0) const {
2455  return weights[i].first;
2456  }
2457 
2458  /**
2459  * Return the weight for the given weight name.
2460  */
2461  double weight(std::string name) const {
2462  return weight(heprup->weightIndex(name));
2463  }
2464 
2465  /**
2466  * Set the weight with the given index.
2467  */
2468  void setWeight(int i, double w) {
2469  weights[i].first = w;
2470  }
2471  /**
2472  * Set the weight with the given name.
2473  */
2474  bool setWeight(std::string name, double w) {
2475  int i = heprup->weightIndex(name);
2476  if ( i >= int(weights.size()) ) return false;
2477  setWeight(i, w);
2478  return true;
2479  }
2480 
2481 
2482  /**
2483  * Assuming the NUP variable, corresponding to the number of
2484  * particles in the current event, is correctly set, resize the
2485  * relevant vectors accordingly.
2486  */
2487  void resize() {
2488  IDUP.resize(NUP);
2489  ISTUP.resize(NUP);
2490  MOTHUP.resize(NUP);
2491  ICOLUP.resize(NUP);
2492  PUP.resize(NUP, std::vector<double>(5));
2493  VTIMUP.resize(NUP);
2494  SPINUP.resize(NUP);
2495  }
2496 
2497  /**
2498  * Setup the current event to use weight i. If zero, the default
2499  * weight will be used.
2500  */
2501  bool setWeightInfo(unsigned int i) {
2502  if ( i >= weights.size() ) return false;
2503  if ( currentWeight ) {
2508  }
2509  XWGTUP = weights[i].first;
2510  currentWeight = weights[i].second;
2511  if ( currentWeight ) {
2516  if ( currentWeight->pdf ) {
2517  heprup->PDFGUP.first = heprup->PDFGUP.second = 0;
2518  heprup->PDFSUP.first = heprup->PDFSUP.second = currentWeight->pdf;
2519  }
2520  if ( currentWeight->pdf2 ) {
2521  heprup->PDFSUP.second = currentWeight->pdf2;
2522  }
2523 
2524  }
2525  return true;
2526  }
2527 
2528  /**
2529  * Setup the current event to use sub event i. If zero, no sub event
2530  * will be chsen.
2531  */
2532  bool setSubEvent(unsigned int i) {
2533  if ( i > subevents.size() || subevents.empty() ) return false;
2534  if ( i == 0 ) {
2535  reset();
2536  weights = subevents[0]->weights;
2537  for ( int ii = 1, N = subevents.size(); ii < N; ++ii )
2538  for ( int j = 0, M = weights.size(); j < M; ++j )
2539  weights[j].first += subevents[ii]->weights[j].first;
2540  currentWeight = 0;
2541  } else {
2542  setEvent(*subevents[i - 1]);
2543  }
2544  return true;
2545  }
2546 
2547 public:
2548 
2549  /**
2550  * The number of particle entries in the current event.
2551  */
2552  int NUP;
2553 
2554  /**
2555  * The subprocess code for this event (as given in LPRUP).
2556  */
2557  int IDPRUP;
2558 
2559  /**
2560  * The weight for this event.
2561  */
2562  double XWGTUP;
2563 
2564  /**
2565  * The PDF weights for the two incoming partons. Note that this
2566  * variable is not present in the current LesHouches accord
2567  * (<A HREF="http://arxiv.org/abs/hep-ph/0109068">hep-ph/0109068</A>),
2568  * hopefully it will be present in a future accord.
2569  */
2570  std::pair<double,double> XPDWUP;
2571 
2572  /**
2573  * The scale in GeV used in the calculation of the PDF's in this
2574  * event.
2575  */
2576  double SCALUP;
2577 
2578  /**
2579  * The value of the QED coupling used in this event.
2580  */
2581  double AQEDUP;
2582 
2583  /**
2584  * The value of the QCD coupling used in this event.
2585  */
2586  double AQCDUP;
2587 
2588  /**
2589  * The PDG id's for the particle entries in this event.
2590  */
2591  std::vector<long> IDUP;
2592 
2593  /**
2594  * The status codes for the particle entries in this event.
2595  */
2596  std::vector<int> ISTUP;
2597 
2598  /**
2599  * Indices for the first and last mother for the particle entries in
2600  * this event.
2601  */
2602  std::vector< std::pair<int,int> > MOTHUP;
2603 
2604  /**
2605  * The colour-line indices (first(second) is (anti)colour) for the
2606  * particle entries in this event.
2607  */
2608  std::vector< std::pair<int,int> > ICOLUP;
2609 
2610  /**
2611  * Lab frame momentum (Px, Py, Pz, E and M in GeV) for the particle
2612  * entries in this event.
2613  */
2614  std::vector< std::vector<double> > PUP;
2615 
2616  /**
2617  * Invariant lifetime (c*tau, distance from production to decay in
2618  * mm) for the particle entries in this event.
2619  */
2620  std::vector<double> VTIMUP;
2621 
2622  /**
2623  * Spin info for the particle entries in this event given as the
2624  * cosine of the angle between the spin vector of a particle and the
2625  * 3-momentum of the decaying particle, specified in the lab frame.
2626  */
2627  std::vector<double> SPINUP;
2628 
2629  /**
2630  * A pointer to the current HEPRUP object.
2631  */
2633 
2634  /**
2635  * The current weight info object.
2636  */
2638 
2639  /**
2640  * The weights associated with this event
2641  */
2642  std::vector<Weight> namedweights;
2643 
2644  /**
2645  * The weights for this event and their corresponding WeightInfo object.
2646  */
2647  std::vector< std::pair<double, const WeightInfo *> > weights;
2648 
2649  /**
2650  * Contents of the clustering tag.
2651  */
2652  std::vector<Clus> clustering;
2653 
2654  /**
2655  * Contents of the pdfinfo tag.
2656  */
2658 
2659  /**
2660  * Saved information about pdfs if different in a selected weight.
2661  */
2662  std::pair<int,int> PDFGUPsave;
2663 
2664  /**
2665  * Saved information about pdfs if different in a selected weight.
2666  */
2667  std::pair<int,int> PDFSUPsave;
2668 
2669 
2670  /**
2671  * Contents of the scales tag
2672  */
2674 
2675  /**
2676  * The number of attempts the ME generator did before accepting this
2677  * event.
2678  */
2679  int ntries;
2680 
2681  /**
2682  * Is this an event or an event group?
2683  */
2684  bool isGroup;
2685 
2686  /**
2687  * If this is not a single event, but an event group, the events
2688  * included in the group are in this vector;
2689  */
2691 
2692  /**
2693  * Save junk stuff in events just to be on the safe side
2694  */
2695  std::string junk;
2696 
2697 };
2698 
2699 
2700 // Destructor implemented here.
2701 
2702 inline void EventGroup::clear() {
2703  while ( size() > 0 ) {
2704  delete back();
2705  pop_back();
2706  }
2707 }
2708 
2710  clear();
2711 }
2712 
2714  : std::vector<HEPEUP*>(eg.size()) {
2715  for ( int i = 0, N = eg.size(); i < N; ++i ) at(i) = new HEPEUP(*eg.at(i));
2716 }
2717 
2719  if ( &x == this ) return *this;
2720  clear();
2721  nreal = x.nreal;
2722  ncounter = x.ncounter;
2723  for ( int i = 0, N = x.size(); i < N; ++i ) push_back(new HEPEUP(*x.at(i)));
2724  return *this;
2725 }
2726 
2727 
2728 /**
2729  * The Reader class is initialized with a stream from which to read a
2730  * version 1/2 Les Houches Accord event file. In the constructor of
2731  * the Reader object the optional header information is read and then
2732  * the mandatory init is read. After this the whole header block
2733  * including the enclosing lines with tags are available in the public
2734  * headerBlock member variable. Also the information from the init
2735  * block is available in the heprup member variable and any additional
2736  * comment lines are available in initComments. After each successful
2737  * call to the readEvent() function the standard Les Houches Accord
2738  * information about the event is available in the hepeup member
2739  * variable and any additional comments in the eventComments
2740  * variable. A typical reading sequence would look as follows:
2741  *
2742  *
2743  */
2744 class Reader {
2745 
2746 public:
2747 
2748  /**
2749  * Initialize the Reader with a stream from which to read an event
2750  * file. After the constructor is called the whole header block
2751  * including the enclosing lines with tags are available in the
2752  * public headerBlock member variable. Also the information from the
2753  * init block is available in the heprup member variable and any
2754  * additional comment lines are available in initComments.
2755  *
2756  * @param is the stream to read from.
2757  */
2758  Reader(std::istream & is)
2759  : file(&is), currevent(-1),
2760  curreventfile(-1), currfileevent(-1), dirpath("") {
2761  init();
2762  }
2763 
2764  /**
2765  * Initialize the Reader with a filename from which to read an event
2766  * file. After the constructor is called the whole header block
2767  * including the enclosing lines with tags are available in the
2768  * public headerBlock member variable. Also the information from the
2769  * init block is available in the heprup member variable and any
2770  * additional comment lines are available in initComments.
2771  *
2772  * @param filename the name of the file to read from.
2773  */
2774  Reader(std::string filename)
2775  : intstream(filename.c_str()), file(&intstream), currevent(-1),
2776  curreventfile(-1), currfileevent(-1), dirpath("") {
2777 
2778  size_t slash = filename.find_last_of('/');
2779  if ( slash != std::string::npos ) dirpath = filename.substr(0, slash + 1);
2780  init();
2781  }
2782 
2783 private:
2784 
2785  /**
2786  * Used internally in the constructors to read header and init
2787  * blocks.
2788  */
2789  void init() {
2790 
2791  // initialize reading from multi-file runs.
2792  initfile = file;
2793 
2794  bool readingHeader = false;
2795  bool readingInit = false;
2796 
2797  // Make sure we are reading a LHEF file:
2798  getline();
2799  if ( !currentFind("<LesHouchesEvents") )
2800  throw std::runtime_error
2801  ("Tried to read a file which does not start with the "
2802  "LesHouchesEvents tag.");
2803  version = 1;
2804  if ( currentFind("version=\"3" ) )
2805  version = 3;
2806  else if ( currentFind("version=\"2" ) )
2807  version = 2;
2808  else if ( !currentFind("version=\"1" ) )
2809  throw std::runtime_error
2810  ("Tried to read a LesHouchesEvents file which is above version 3.");
2811 
2812  // Loop over all lines until we hit the </init> tag.
2813  while ( getline() && !currentFind("</init>") ) {
2814  if ( currentFind("<header") ) {
2815  // We have hit the header block, so we should dump this and
2816  // all following lines to headerBlock until we hit the end of
2817  // it.
2818  readingHeader = true;
2819  headerBlock = currentLine + "\n";
2820  }
2821  else if ( currentFind("<init>") ) {
2822  // We have hit the init block
2823  readingInit = true;
2824  initComments = currentLine + "\n";
2825  }
2826  else if ( currentFind("</header>") ) {
2827  // The end of the header block. Dump this line as well to the
2828  // headerBlock and we're done.
2829  readingHeader = false;
2830  headerBlock += currentLine + "\n";
2831  }
2832  else if ( readingHeader ) {
2833  // We are in the process of reading the header block. Dump the
2834  // line to haderBlock.
2835  headerBlock += currentLine + "\n";
2836  }
2837  else if ( readingInit ) {
2838  // Here we found a comment line. Dump it to initComments.
2839  initComments += currentLine + "\n";
2840  }
2841  else {
2842  // We found some other stuff outside the standard tags.
2843  outsideBlock += currentLine + "\n";
2844  }
2845  }
2846  if ( !currentFind("</init>") )
2847  throw std::runtime_error("Found incomplete init tag in "
2848  "Les Houches file.");
2849  initComments += currentLine + "\n";
2850  std::vector<XMLTag*> tags = XMLTag::findXMLTags(initComments);
2851  for ( int i = 0, N = tags.size(); i < N; ++i )
2852  if ( tags[i]->name == "init" ) {
2853  heprup = HEPRUP(*tags[i], version);
2854  break;
2855  }
2856  XMLTag::deleteAll(tags);
2857 
2858  if ( !heprup.eventfiles.empty() ) openeventfile(0);
2859 
2860  }
2861 
2862 public:
2863 
2864  /**
2865  * Read an event from the file and store it in the hepeup
2866  * object. Optional comment lines are stored i the eventComments
2867  * member variable.
2868  * @return true if the read sas successful.
2869  */
2870  bool readEvent() {
2871 
2872  // Check if the initialization was successful. Otherwise we will
2873  // not read any events.
2874  if ( heprup.NPRUP < 0 ) return false;
2875 
2876  std::string eventLines;
2877  int inEvent = 0;
2878 
2879  // Keep reading lines until we hit the end of an event or event group.
2880  while ( getline() ) {
2881  if ( inEvent ) {
2882  eventLines += currentLine + "\n";
2883  if ( inEvent == 1 && currentFind("</event>") ) break;
2884  if ( inEvent == 2 && currentFind("</eventgroup>") ) break;
2885  }
2886  else if ( currentFind("<eventgroup") ) {
2887  eventLines += currentLine + "\n";
2888  inEvent = 2;
2889  }
2890  else if ( currentFind("<event") ) {
2891  eventLines += currentLine + "\n";
2892  inEvent = 1;
2893  }
2894  else {
2895  outsideBlock += currentLine + "\n";
2896  }
2897  }
2898  if ( ( inEvent == 1 && !currentFind("</event>") ) ||
2899  ( inEvent == 2 && !currentFind("</eventgroup>") ) ) {
2900  if ( heprup.eventfiles.empty() ||
2901  ++curreventfile >= int(heprup.eventfiles.size()) ) return false;
2903  return readEvent();
2904  }
2905 
2906  std::vector<XMLTag*> tags = XMLTag::findXMLTags(eventLines);
2907 
2908  for ( int i = 0, N = tags.size(); i < N ; ++i ) {
2909  if ( tags[i]->name == "event" || tags[i]->name == "eventgroup" ) {
2910  hepeup = HEPEUP(*tags[i], heprup);
2911  XMLTag::deleteAll(tags);
2912  ++currevent;
2913  if ( curreventfile >= 0 ) ++currfileevent;
2914  return true;
2915  }
2916  }
2917 
2918  if ( !heprup.eventfiles.empty() &&
2919  ++curreventfile < int(heprup.eventfiles.size()) ) {
2921  return readEvent();
2922  }
2923 
2924  XMLTag::deleteAll(tags);
2925  return false;
2926 
2927  }
2928 
2929  /**
2930  * Open the efentfile with index ifile. If another eventfile is
2931  * being read, its remaining contents is discarded. This is a noop
2932  * if current read session is not a multi-file run.
2933  */
2934  void openeventfile(int ifile) {
2935  std::cerr << "opening file " << ifile << std::endl;
2936  efile.close();
2937  std::string fname = heprup.eventfiles[ifile].filename;
2938  if ( fname[0] != '/' ) fname = dirpath + fname;
2939  efile.open(fname.c_str());
2940  if ( !efile ) throw std::runtime_error("Could not open event file " +
2941  fname);
2942  file = &efile;
2943  curreventfile = ifile;
2944  currfileevent = 0;
2945  }
2946 
2947 protected:
2948 
2949  /**
2950  * Used internally to read a single line from the stream.
2951  */
2952  bool getline() {
2953  return ( (bool)std::getline(*file, currentLine) );
2954  }
2955 
2956  /**
2957  * @return true if the current line contains the given string.
2958  */
2959  bool currentFind(std::string str) const {
2960  return currentLine.find(str) != std::string::npos;
2961  }
2962 
2963 protected:
2964 
2965  /**
2966  * A local stream which is unused if a stream is supplied from the
2967  * outside.
2968  */
2969  std::ifstream intstream;
2970 
2971  /**
2972  * The stream we are reading from. This may be a pointer to an
2973  * external stream or the internal intstream, or a separate event
2974  * file from a multi-file run
2975  */
2976  std::istream * file;
2977 
2978  /**
2979  * The original stream from where we read the init block.
2980  */
2981  std::istream * initfile;
2982 
2983  /**
2984  * A separate stream for reading multi-file runs.
2985  */
2986  std::ifstream efile;
2987 
2988  /**
2989  * The last line read in from the stream in getline().
2990  */
2991  std::string currentLine;
2992 
2993 public:
2994 
2995  /**
2996  * XML file version
2997  */
2998  int version;
2999 
3000  /**
3001  * All lines (since the last readEvent()) outside the header, init
3002  * and event tags.
3003  */
3004  std::string outsideBlock;
3005 
3006  /**
3007  * All lines from the header block.
3008  */
3009  std::string headerBlock;
3010 
3011  /**
3012  * The standard init information.
3013  */
3015 
3016  /**
3017  * Additional comments found in the init block.
3018  */
3019  std::string initComments;
3020 
3021  /**
3022  * The standard information about the last read event.
3023  */
3025 
3026  /**
3027  * Additional comments found with the last read event.
3028  */
3029  std::string eventComments;
3030 
3031  /**
3032  * The number of the current event (starting from 1).
3033  */
3035 
3036  /**
3037  * The current event file being read from (-1 means there are no
3038  * separate event files).
3039  */
3041 
3042  /**
3043  * The number of the current event in the current event file.
3044  */
3046 
3047  /**
3048  * The directory from where we are reading files.
3049  */
3050  std::string dirpath;
3051 
3052 private:
3053 
3054  /**
3055  * The default constructor should never be used.
3056  */
3057  Reader();
3058 
3059  /**
3060  * The copy constructor should never be used.
3061  */
3062  Reader(const Reader &);
3063 
3064  /**
3065  * The Reader cannot be assigned to.
3066  */
3067  Reader & operator=(const Reader &);
3068 
3069 };
3070 
3071 /**
3072  * The Writer class is initialized with a stream to which to write a
3073  * version 1.0 Les Houches Accord event file. In the constructor of
3074  * the Writer object the main XML tag is written out, with the
3075  * corresponding end tag is written in the destructor. After a Writer
3076  * object has been created, it is possible to assign standard init
3077  * information in the heprup member variable. In addition any XML
3078  * formatted information can be added to the headerBlock member
3079  * variable (directly or via the addHeader() function). Further
3080  * comment line (beginning with a <code>#</code> character) can be
3081  * added to the initComments variable (directly or with the
3082  * addInitComment() function). After this information is set, it
3083  * should be written out to the file with the init() function.
3084  *
3085  * Before each event is written out with the writeEvent() function,
3086  * the standard event information can then be assigned to the hepeup
3087  * variable and optional comment lines (beginning with a
3088  * <code>#</code> character) may be given to the eventComments
3089  * variable (directly or with the addEventComment() function).
3090  *
3091  */
3092 class Writer {
3093 
3094 public:
3095 
3096 #ifndef HEPMC3_PYTHON_BINDINGS
3097  /**
3098  * Create a Writer object giving a stream to write to.
3099  * @param os the stream where the event file is written.
3100  */
3101  Writer(std::ostream & os)
3102  : file(&os), initfile(&os), dirpath("") { }
3103 #endif
3104  /**
3105  * Create a Writer object giving a filename to write to.
3106  * @param filename the name of the event file to be written.
3107  */
3108  Writer(std::string filename)
3109  : intstream(filename.c_str()), file(&intstream), initfile(&intstream),
3110  dirpath("") {
3111  size_t slash = filename.find_last_of('/');
3112  if ( slash != std::string::npos ) dirpath = filename.substr(0, slash + 1);
3113  }
3114 
3115  /**
3116  * The destructor writes out the final XML end-tag.
3117  */
3119  file = initfile;
3120  if ( !heprup.eventfiles.empty() ) {
3121  if ( curreventfile >= 0 &&
3122  curreventfile < int(heprup.eventfiles.size()) &&
3123  heprup.eventfiles[curreventfile].neve < 0 )
3125  writeinit();
3126  }
3127  *file << "</LesHouchesEvents>" << std::endl;
3128  }
3129 #ifndef HEPMC3_PYTHON_BINDINGS
3130  /**
3131  * Add header lines consisting of XML code with this stream.
3132  */
3133  std::ostream & headerBlock() {
3134  return headerStream;
3135  }
3136 
3137  /**
3138  * Add comment lines to the init block with this stream.
3139  */
3140  std::ostream & initComments() {
3141  return initStream;
3142  }
3143 
3144  /**
3145  * Add comment lines to the next event to be written out with this stream.
3146  */
3147  std::ostream & eventComments() {
3148  return eventStream;
3149  }
3150 #endif
3151  /**
3152  * Add header lines consisting of XML code with this stream.
3153  */
3154  void headerBlock(const std::string& a) {
3155  headerStream<<a;;
3156  }
3157 
3158  /**
3159  * Add comment lines to the init block with this stream.
3160  */
3161  void initComments(const std::string& a) {
3162  initStream<<a;
3163  }
3164 
3165  /**
3166  * Add comment lines to the next event to be written out with this stream.
3167  */
3168  void eventComments(const std::string& a) {
3169  eventStream<<a;
3170  }
3171 
3172  /**
3173  * Initialize the writer.
3174  */
3175  void init() {
3176  if ( heprup.eventfiles.empty() ) writeinit();
3177  lastevent = 0;
3179  if ( !heprup.eventfiles.empty() ) openeventfile(0);
3180  }
3181 
3182  /**
3183  * Open a new event file, possibly closing a previous opened one.
3184  */
3185  bool openeventfile(int ifile) {
3186  if ( heprup.eventfiles.empty() ) return false;
3187  if ( ifile < 0 || ifile >= int(heprup.eventfiles.size()) ) return false;
3188  if ( curreventfile >= 0 ) {
3190  if ( ef.neve > 0 && ef.neve != currfileevent )
3191  std::cerr << "LHEF::Writer number of events in event file "
3192  << ef.filename << " does not match the given number."
3193  << std::endl;
3194  ef.neve = currfileevent;
3195  }
3196  efile.close();
3197  std::string fname = heprup.eventfiles[ifile].filename;
3198  if ( fname[0] != '/' ) fname = dirpath + fname;
3199  efile.open(fname.c_str());
3200  if ( !efile ) throw std::runtime_error("Could not open event file " +
3201  fname);
3202  std::cerr << "Opened event file " << fname << std::endl;
3203  file = &efile;
3204  curreventfile = ifile;
3205  currfileevent = 0;
3206  return true;
3207  }
3208 
3209 
3210  /**
3211  * Write out an optional header block followed by the standard init
3212  * block information together with any comment lines.
3213  */
3214  void writeinit() {
3215 
3216  // Write out the standard XML tag for the event file.
3217  if ( heprup.version == 3 )
3218  *file << "<LesHouchesEvents version=\"3.0\">\n";
3219  else if ( heprup.version == 2 )
3220  *file << "<LesHouchesEvents version=\"2.0\">\n";
3221  else
3222  *file << "<LesHouchesEvents version=\"1.0\">\n";
3223 
3224 
3225  *file << std::setprecision(10);
3226 
3227  using std::setw;
3228 
3229  std::string headBlock = headerStream.str();
3230  if ( headBlock.length() ) {
3231  if ( headBlock.find("<header>") == std::string::npos )
3232  *file << "<header>\n";
3233  if ( headBlock[headBlock.length() - 1] != '\n' )
3234  headBlock += '\n';
3235  *file << headBlock;
3236  if ( headBlock.find("</header>") == std::string::npos )
3237  *file << "</header>\n";
3238  }
3239 
3240  heprup.print(*file);
3241 
3242  }
3243 
3244  /**
3245  * Write the current HEPEUP object to the stream;
3246  */
3247  void writeEvent() {
3248 
3249  if ( !heprup.eventfiles.empty() ) {
3250  if ( currfileevent == heprup.eventfiles[curreventfile].neve &&
3251  curreventfile + 1 < int(heprup.eventfiles.size()) )
3253  }
3254 
3255  hepeup.print(*file);
3256 
3257  ++lastevent;
3258  ++currfileevent;
3259  }
3260 
3261 protected:
3262 
3263  /**
3264  * A local stream which is unused if a stream is supplied from the
3265  * outside.
3266  */
3267  std::ofstream intstream;
3268 
3269  /**
3270  * The stream we are writing to. This may be a reference to an
3271  * external stream or the internal intstream.
3272  */
3273  std::ostream * file;
3274 
3275  /**
3276  * The original stream from where we read the init block.
3277  */
3278  std::ostream * initfile;
3279 
3280  /**
3281  * A separate stream for reading multi-file runs.
3282  */
3283  std::ofstream efile;
3284 
3285  /**
3286  * The number of the last event written (starting from 1).
3287  */
3289 
3290  /**
3291  * The current event file being written to (-1 means there are no
3292  * separate event files).
3293  */
3295 
3296  /**
3297  * The number of the current event in the current event file.
3298  */
3300 
3301  /**
3302  * The directory from where we are reading files.
3303  */
3304  std::string dirpath;
3305 
3306 public:
3307  /**
3308  * The standard init information.
3309  */
3311 
3312 
3313  /**
3314  * The standard information about the event we will write next.
3315  */
3317 
3318 
3319 
3320 private:
3321 
3322  /**
3323  * Stream to add all lines in the header block.
3324  */
3325  std::ostringstream headerStream;
3326 
3327  /**
3328  * Stream to add additional comments to be put in the init block.
3329  */
3330  std::ostringstream initStream;
3331 
3332  /**
3333  * Stream to add additional comments to be written together the next event.
3334  */
3335  std::ostringstream eventStream;
3336 
3337 
3338  /**
3339  * The default constructor should never be used.
3340  */
3341  Writer();
3342 
3343  /**
3344  * The copy constructor should never be used.
3345  */
3346  Writer(const Writer &);
3347 
3348  /**
3349  * The Writer cannot be assigned to.
3350  */
3351  Writer & operator=(const Writer &);
3352 
3353 };
3354 
3355 }
3356 
3357 /* \example LHEFCat.cc This is a main function which simply reads a
3358  Les Houches Event File from the standard input and writes it again
3359  to the standard output.
3360  This file can be downloaded from
3361  <A HREF="http://www.thep.lu.se/~leif/LHEF/LHEFCat.cc">here</A>.
3362  There are also two sample event files,
3363  <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>
3364  to try it on.
3365 */
3366 
3367 /* \mainpage Les Houches Event File
3368 
3369 Here are some example classes for reading and writing Les Houches
3370 Event Files according to the
3371 <A HREF="http://www.thep.lu.se/~torbjorn/lhef/lhafile2.pdf">proposal</A>
3372 by Torbj&ouml;rn Sj&ouml;strand discussed at the
3373 <A HREF="http://mc4lhc06.web.cern.ch/mc4lhc06/">MC4LHC</A>
3374 workshop at CERN 2006.
3375 
3376 The classes has now been updated to handle the suggested version 3 of
3377 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>).
3378 
3379 There is a whole set of classes available in a single header file
3380 called <A
3381 HREF="http://www.thep.lu.se/~leif/LHEF/LHEF.h">LHEF.h</A>. The idea is
3382 that matrix element or parton shower generators will implement their
3383 own wrapper using these classes as containers.
3384 
3385 The two classes LHEF::HEPRUP and LHEF::HEPEUP are simple container
3386 classes which contain the same information as the Les Houches standard
3387 Fortran common blocks with the same names. They also contain the extra
3388 information defined in version 3 in the standard. The other two main
3389 classes are called LHEF::Reader and LHEF::Writer and are used to read
3390 and write Les Houches Event Files
3391 
3392 Here are a few <A HREF="examples.html">examples</A> of how to use the
3393 classes:
3394 
3395 \namespace LHEF The LHEF namespace contains some example classes for reading and writing Les Houches
3396 Event Files according to the
3397 <A HREF="http://www.thep.lu.se/~torbjorn/lhef/lhafile2.pdf">proposal</A>
3398 by Torbj&ouml;rn Sj&ouml;strand discussed at the
3399 <A HREF="http://mc4lhc06.web.cern.ch/mc4lhc06/">MC4LHC</A>
3400 workshop at CERN 2006.
3401 
3402 
3403 
3404  */
3405 
3406 
3407 #endif /* HEPMC3_LHEF_H */
XSecInfo & getXSecInfo(std::string weightname="")
Definition: LHEF.h:1915
MergeInfo(const XMLTag &tag)
Definition: LHEF.h:1002
void print(std::ostream &file) const
Definition: LHEF.h:1197
std::vector< XMLTag * > tags
Definition: LHEF.h:129
int weightIndex(std::string name) const
Definition: LHEF.h:1898
std::ofstream intstream
Definition: LHEF.h:3267
std::string weightNameHepMC(int i) const
Definition: LHEF.h:1751
void setWeight(int i, double w)
Definition: LHEF.h:2468
bool getattr(std::string n, int &v) const
Definition: LHEF.h:174
double x1
Definition: LHEF.h:1590
std::vector< double > VTIMUP
Definition: LHEF.h:2620
bool getattr(std::string n, double &v, bool erase=true)
Definition: LHEF.h:368
double totalWeight(std::string name) const
Definition: LHEF.h:2447
Scales(const XMLTag &tag, double defscale=-1.0, int npart=0)
Definition: LHEF.h:1429
void print(std::ostream &file) const
Definition: LHEF.h:1013
int nWeights() const
Definition: LHEF.h:1907
void print(std::ostream &file) const
Definition: LHEF.h:722
bool getattr(std::string n, int &v, bool erase=true)
Definition: LHEF.h:410
std::ofstream efile
Definition: LHEF.h:3283
std::pair< int, int > PDFGUPsave
Definition: LHEF.h:2662
std::string np1
Definition: LHEF.h:889
static void deleteAll(std::vector< XMLTag *> &tags)
Definition: LHEF.h:293
bool setWeight(std::string name, double w)
Definition: LHEF.h:2474
Reader(std::string filename)
Definition: LHEF.h:2774
static std::vector< XMLTag * > findXMLTags(std::string str, std::string *leftover=0)
Definition: LHEF.h:198
double mups
Definition: LHEF.h:1523
std::vector< int > ISTUP
Definition: LHEF.h:2596
bool hasInfo() const
Definition: LHEF.h:1452
HEPEUP hepeup
Definition: LHEF.h:3024
double XWGTUP
Definition: LHEF.h:2562
std::string name
Definition: LHEF.h:55
double mergingscale
Definition: LHEF.h:1029
std::string contents
Definition: LHEF.h:462
double totxsec
Definition: LHEF.h:574
XSecInfo(const XMLTag &tag)
Definition: LHEF.h:522
XMLTag()
Definition: LHEF.h:107
double scale
Definition: LHEF.h:1299
Clus(const XMLTag &tag)
Definition: LHEF.h:1260
void print(std::ostream &os) const
Definition: LHEF.h:302
HEPEUP & setEvent(const HEPEUP &x)
Definition: LHEF.h:2144
bool getattr(std::string n, std::string &v, bool erase=true)
Definition: LHEF.h:424
std::string rscheme
Definition: LHEF.h:980
double maxweight
Definition: LHEF.h:584
double muf
Definition: LHEF.h:1512
double min
Definition: LHEF.h:904
std::istream * initfile
Definition: LHEF.h:2981
bool negweights
Definition: LHEF.h:594
std::string hashline(std::string s)
Definition: LHEF.h:328
double getScale(std::string st, int pdgem, int emr, int rec) const
Definition: LHEF.h:1489
double AQCDUP
Definition: LHEF.h:2586
std::ostream & eventComments()
Definition: LHEF.h:3147
std::string headerBlock
Definition: LHEF.h:3009
std::ostringstream headerStream
Definition: LHEF.h:3325
double mur
Definition: LHEF.h:1517
std::vector< EventFile > eventfiles
Definition: LHEF.h:2001
std::string stype
Definition: LHEF.h:1387
std::set< long > p2
Definition: LHEF.h:894
double xsecerr
Definition: LHEF.h:579
std::ostream & headerBlock()
Definition: LHEF.h:3133
void resize()
Definition: LHEF.h:1888
double mur
Definition: LHEF.h:1110
std::vector< Scale > scales
Definition: LHEF.h:1533
Reader & operator=(const Reader &)
STL namespace.
std::vector< std::pair< int, int > > MOTHUP
Definition: LHEF.h:2602
bool readEvent()
Definition: LHEF.h:2870
void reset()
Definition: LHEF.h:2406
std::pair< double, double > EBMUP
Definition: LHEF.h:1944
#define M_PI
Definition of PI. Needed on some platforms.
Definition: FourVector.h:15
double SCALUP
Definition: LHEF.h:1615
std::string initComments
Definition: LHEF.h:3019
int dprec
Definition: LHEF.h:2057
void init()
Definition: LHEF.h:3175
std::ostringstream initStream
Definition: LHEF.h:3330
~XMLTag()
Definition: LHEF.h:112
bool maxmult
Definition: LHEF.h:1034
std::map< long, ProcInfo > procinfo
Definition: LHEF.h:2016
int NUP
Definition: LHEF.h:2552
double born
Definition: LHEF.h:1227
const XSecInfo & getXSecInfo(std::string weightname="") const
Definition: LHEF.h:1926
bool match(long id1, long id2=0) const
Definition: LHEF.h:749
std::vector< double > XERRUP
Definition: LHEF.h:1979
int IDWTUP
Definition: LHEF.h:1963
double sudakov
Definition: LHEF.h:1232
OAttr< T > oattr(std::string name, const T &value)
Definition: LHEF.h:68
std::ostream & initComments()
Definition: LHEF.h:3140
AttributeMap attr
Definition: LHEF.h:124
double muf
Definition: LHEF.h:1105
std::pair< double, double > XPDWUP
Definition: LHEF.h:2570
Writer(std::ostream &os)
Definition: LHEF.h:3101
bool passCuts(const std::vector< long > &id, const std::vector< std::vector< double > > &p) const
Definition: LHEF.h:765
OAttr(std::string n, const T &v)
Definition: LHEF.h:50
WeightGroup(const XMLTag &tag, int groupIndex, std::vector< WeightInfo > &wiv)
Definition: LHEF.h:1139
int eworder
Definition: LHEF.h:970
double max
Definition: LHEF.h:908
std::ifstream intstream
Definition: LHEF.h:2969
std::map< std::string, XSecInfo > XSecInfos
Definition: LHEF.h:611
std::string name
Definition: LHEF.h:499
Clus()
Definition: LHEF.h:1255
HEPRUP heprup
Definition: LHEF.h:3310
std::pair< int, int > PDFSUPsave
Definition: LHEF.h:2667
bool isGroup
Definition: LHEF.h:2684
std::string dirpath
Definition: LHEF.h:3050
int version
Definition: LHEF.h:2998
double xf1
Definition: LHEF.h:1600
bool getattr(std::string n, std::string &v) const
Definition: LHEF.h:185
Writer(std::string filename)
Definition: LHEF.h:3108
std::string contents
Definition: LHEF.h:134
Generator(const XMLTag &tag)
Definition: LHEF.h:479
XMLTag::AttributeMap AttributeMap
Definition: LHEF.h:350
HEPEUP(const HEPEUP &x)
Definition: LHEF.h:2135
int currevent
Definition: LHEF.h:3034
bool getattr(std::string n, bool &v) const
Definition: LHEF.h:152
int NPRUP
Definition: LHEF.h:1968
static std::string yes()
Definition: LHEF.h:467
PDFInfo(double defscale=-1.0)
Definition: LHEF.h:1545
ProcInfo(const XMLTag &tag)
Definition: LHEF.h:925
Les Houches event file classes.
Definition: LHEF.h:39
void print(std::ostream &file) const
Definition: LHEF.h:1565
void print(std::ostream &file) const
Definition: LHEF.h:940
EventFile(const XMLTag &tag)
Definition: LHEF.h:627
std::string name
Definition: LHEF.h:1217
std::string junk
Definition: LHEF.h:2047
HEPRUP(const XMLTag &tagin, int versin)
Definition: LHEF.h:1653
XMLTag::AttributeMap attributes
Definition: LHEF.h:457
int curreventfile
Definition: LHEF.h:3294
bool setSubEvent(unsigned int i)
Definition: LHEF.h:2532
double weight(std::string name) const
Definition: LHEF.h:2461
Scale(std::string st="veto", int emtr=0, double sc=0.0)
Definition: LHEF.h:1318
XSecInfos xsecinfos
Definition: LHEF.h:1995
bool getattr(std::string n, long &v, bool erase=true)
Definition: LHEF.h:396
bool iswgt
Definition: LHEF.h:1222
std::pair< int, int > PDFGUP
Definition: LHEF.h:1950
const WeightInfo * currentWeight
Definition: LHEF.h:2637
long p1
Definition: LHEF.h:1580
std::string type
Definition: LHEF.h:879
int currfileevent
Definition: LHEF.h:3299
int lastevent
Definition: LHEF.h:3288
PDFInfo pdfinfo
Definition: LHEF.h:2657
std::ostringstream eventStream
Definition: LHEF.h:3335
double alphas
Definition: LHEF.h:1305
void print(std::ostream &file) const
Definition: LHEF.h:488
void print(std::ostream &file) const
Definition: LHEF.h:640
std::ostream * initfile
Definition: LHEF.h:3278
double scale
Definition: LHEF.h:1610
HEPEUP(const XMLTag &tagin, HEPRUP &heprupin)
Definition: LHEF.h:2200
int p0
Definition: LHEF.h:1294
void resize(int nrup)
Definition: LHEF.h:1878
WeightInfo(const XMLTag &tag)
Definition: LHEF.h:1053
double xf2
Definition: LHEF.h:1605
std::ifstream efile
Definition: LHEF.h:2986
std::string outsideBlock
Definition: LHEF.h:3004
std::vector< std::pair< double, const WeightInfo * > > weights
Definition: LHEF.h:2647
std::string::size_type pos_t
Definition: LHEF.h:92
std::string name
Definition: LHEF.h:1100
double x2
Definition: LHEF.h:1595
void clear()
Definition: LHEF.h:1863
Scales scales
Definition: LHEF.h:2673
int loops
Definition: LHEF.h:960
bool getline()
Definition: LHEF.h:2952
TagBase(const AttributeMap &attr, std::string conts=std::string())
Definition: LHEF.h:360
bool currentFind(std::string str) const
Definition: LHEF.h:2959
std::set< int > recoilers
Definition: LHEF.h:1399
std::string scheme
Definition: LHEF.h:985
bool outside(double value) const
Definition: LHEF.h:872
int version
Definition: LHEF.h:2052
std::vector< Generator > generators
Definition: LHEF.h:2027
void openeventfile(int ifile)
Definition: LHEF.h:2934
void initComments(const std::string &a)
Definition: LHEF.h:3161
std::vector< std::vector< double > > PUP
Definition: LHEF.h:2614
std::string np2
Definition: LHEF.h:899
HEPRUP & operator=(const HEPRUP &x)=default
void print(std::ostream &file) const
Definition: LHEF.h:1765
void writeEvent()
Definition: LHEF.h:3247
void headerBlock(const std::string &a)
Definition: LHEF.h:3154
HEPEUP hepeup
Definition: LHEF.h:3316
std::map< long, MergeInfo > mergeinfo
Definition: LHEF.h:2021
std::vector< WeightGroup > weightgroup
Definition: LHEF.h:2042
std::vector< long > IDUP
Definition: LHEF.h:2591
int curreventfile
Definition: LHEF.h:3040
void print(std::ostream &file) const
Definition: LHEF.h:1355
long neve
Definition: LHEF.h:564
std::string dirpath
Definition: LHEF.h:3304
std::vector< Weight > namedweights
Definition: LHEF.h:2642
EventGroup & operator=(const EventGroup &)
Definition: LHEF.h:2718
Weight(const XMLTag &tag)
Definition: LHEF.h:1179
std::vector< Cut > cuts
Definition: LHEF.h:2006
std::vector< double > SPINUP
Definition: LHEF.h:2627
void closetag(std::ostream &file, std::string tag) const
Definition: LHEF.h:445
std::istream * file
Definition: LHEF.h:2976
long ntries
Definition: LHEF.h:569
T val
Definition: LHEF.h:60
int iproc
Definition: LHEF.h:955
void print(std::ostream &file) const
Definition: LHEF.h:1460
long ntries
Definition: LHEF.h:662
int emitter
Definition: LHEF.h:1394
std::string filename
Definition: LHEF.h:652
double AQEDUP
Definition: LHEF.h:2581
std::string eventComments
Definition: LHEF.h:3029
std::string junk
Definition: LHEF.h:2695
Reader(std::istream &is)
Definition: LHEF.h:2758
bool getattr(std::string n, bool &v, bool erase=true)
Definition: LHEF.h:382
void clear()
Definition: LHEF.h:2416
std::map< std::string, std::set< long > > ptypes
Definition: LHEF.h:2011
std::string combine
Definition: LHEF.h:1161
bool setWeightInfo(unsigned int i)
Definition: LHEF.h:2501
PDFInfo(const XMLTag &tag, double defscale=-1.0)
Definition: LHEF.h:1551
Cut()
Definition: LHEF.h:674
std::string weightname
Definition: LHEF.h:604
void clear()
Definition: LHEF.h:2702
void print(std::ostream &file) const
Definition: LHEF.h:1272
bool openeventfile(int ifile)
Definition: LHEF.h:3185
std::string currentLine
Definition: LHEF.h:2991
void writeinit()
Definition: LHEF.h:3214
static double eta(const std::vector< double > &p)
Definition: LHEF.h:832
int qcdorder
Definition: LHEF.h:965
std::string name
Definition: LHEF.h:119
std::set< int > emitted
Definition: LHEF.h:1404
int ntries
Definition: LHEF.h:2679
void eventComments(const std::string &a)
Definition: LHEF.h:3168
std::ostream * file
Definition: LHEF.h:3273
Writer & operator=(const Writer &)
int currfileevent
Definition: LHEF.h:3045
void printattrs(std::ostream &file) const
Definition: LHEF.h:435
double scale
Definition: LHEF.h:1409
bool varweights
Definition: LHEF.h:599
std::vector< double > XMAXUP
Definition: LHEF.h:1985
void resize(int nup)
Definition: LHEF.h:2426
std::vector< WeightInfo > weightinfo
Definition: LHEF.h:2032
std::vector< int > LPRUP
Definition: LHEF.h:1990
std::string fscheme
Definition: LHEF.h:975
std::string version
Definition: LHEF.h:504
void init()
Definition: LHEF.h:2789
std::set< long > p1
Definition: LHEF.h:884
std::vector< int > indices
Definition: LHEF.h:1242
double SCALUP
Definition: LHEF.h:2576
std::string type
Definition: LHEF.h:1156
Cut(const XMLTag &tag, const std::map< std::string, std::set< long > > &ptypes)
Definition: LHEF.h:680
std::vector< Clus > clustering
Definition: LHEF.h:2652
HEPRUP heprup
Definition: LHEF.h:3014
bool getattr(std::string n, long &v) const
Definition: LHEF.h:163
Scale(const XMLTag &tag)
Definition: LHEF.h:1324
std::vector< std::pair< int, int > > ICOLUP
Definition: LHEF.h:2608
std::map< std::string, int > weightmap
Definition: LHEF.h:2037
static double deltaR(const std::vector< double > &p1, const std::vector< double > &p2)
Definition: LHEF.h:860
double SCALUP
Definition: LHEF.h:1528
HEPEUP & operator=(const HEPEUP &x)
Definition: LHEF.h:2176
std::vector< double > weights
Definition: LHEF.h:1237
static double rap(const std::vector< double > &p)
Definition: LHEF.h:846
void print(std::ostream &file) const
Definition: LHEF.h:1070
std::vector< double > XSECUP
Definition: LHEF.h:1973
void print(std::ostream &file) const
Definition: LHEF.h:546
Scales(double defscale=-1.0, int npart=0)
Definition: LHEF.h:1421
static const pos_t end
Definition: LHEF.h:102
double meanweight
Definition: LHEF.h:589
std::pair< int, int > PDFSUP
Definition: LHEF.h:1956
long neve
Definition: LHEF.h:657
int p1
Definition: LHEF.h:1284
double totalWeight(int i=0) const
Definition: LHEF.h:2435
void resize()
Definition: LHEF.h:2487
EventGroup subevents
Definition: LHEF.h:2690
void print(std::ostream &file) const
Definition: LHEF.h:2321
long p2
Definition: LHEF.h:1585
bool getattr(std::string n, double &v) const
Definition: LHEF.h:140
std::pair< long, long > IDBMUP
Definition: LHEF.h:1939
HEPRUP * heprup
Definition: LHEF.h:2632
int p2
Definition: LHEF.h:1289
double weight(int i=0) const
Definition: LHEF.h:2454
std::map< std::string, std::string > AttributeMap
Definition: LHEF.h:97
int IDPRUP
Definition: LHEF.h:2557