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  */
124  AttributeMap attr;
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  */
1995  XSecInfos xsecinfos;
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;
2160  namedweights = x.namedweights;
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;
2168  currentWeight = x.currentWeight;
2169  ntries = x.ntries;
2170  return *this;
2171  }
2172 
2173  /**
2174  * Assignment operator.
2175  */
2176  HEPEUP & operator=(const HEPEUP & x) {
2177  clear();
2178  setEvent(x);
2179  subevents = x.subevents;
2180  isGroup = x.isGroup;
2181  return *this;
2182  }
2183 
2184  /**
2185  * Destructor.
2186  */
2188  clear();
2189  };
2190  //@}
2191 
2192 public:
2193 
2194 
2195  /**
2196  * Constructor from an event or eventgroup tag.
2197  */
2198  HEPEUP(const XMLTag & tagin, HEPRUP & heprupin)
2199  : TagBase(tagin.attr), NUP(0), IDPRUP(0), XWGTUP(0.0), XPDWUP(0.0, 0.0),
2200  SCALUP(0.0), AQEDUP(0.0), AQCDUP(0.0), heprup(&heprupin),
2201  currentWeight(0), ntries(1), isGroup(tagin.name == "eventgroup") {
2202 
2203  if ( heprup->NPRUP < 0 )
2204  throw std::runtime_error("Tried to read events but no processes defined "
2205  "in init block of Les Houches file.");
2206 
2207  std::vector<XMLTag*> tags = tagin.tags;
2208 
2209  if ( isGroup ) {
2210  getattr("nreal", subevents.nreal);
2211  getattr("ncounter", subevents.ncounter);
2212  for ( int i = 0, N = tags.size(); i < N; ++i )
2213  if ( tags[i]->name == "event" )
2214  subevents.push_back(new HEPEUP(*tags[i], heprupin));
2215  return;
2216  } else
2217  getattr("ntries", ntries);
2218 
2219 
2220 
2221  // The event information should be in the first (anonymous) tag
2222  std::istringstream iss(tags[0]->contents);
2223  if ( !( iss >> NUP >> IDPRUP >> XWGTUP >> SCALUP >> AQEDUP >> AQCDUP ) )
2224  throw std::runtime_error("Failed to parse event in Les Houches file.");
2225 
2226  resize();
2227 
2228  // Read all particle lines.
2229  for ( int i = 0; i < NUP; ++i ) {
2230  if ( !( iss >> IDUP[i] >> ISTUP[i] >> MOTHUP[i].first >> MOTHUP[i].second
2231  >> ICOLUP[i].first >> ICOLUP[i].second
2232  >> PUP[i][0] >> PUP[i][1] >> PUP[i][2]
2233  >> PUP[i][3] >> PUP[i][4]
2234  >> VTIMUP[i] >> SPINUP[i] ) )
2235  throw std::runtime_error("Failed to parse event in Les Houches file.");
2236  }
2237 
2238  junk.clear();
2239  std::string ss;
2240  while ( getline(iss, ss) ) junk += ss + '\n';
2241 
2242  scales = Scales(SCALUP, NUP);
2243  pdfinfo = PDFInfo(SCALUP);
2244  namedweights.clear();
2245  weights.clear();
2246  weights.resize(heprup->nWeights(),
2247  std::make_pair(XWGTUP, (WeightInfo*)(0)));
2248  weights.front().first = XWGTUP;
2249  for ( int i = 1, N = weights.size(); i < N; ++i )
2250  weights[i].second = &heprup->weightinfo[i - 1];
2251 
2252  for ( int i = 1, N = tags.size(); i < N; ++i ) {
2253  XMLTag & tag = *tags[i];
2254 
2255  if ( tag.name.empty() ) junk += tag.contents;
2256 
2257  if ( tag.name == "weights" ) {
2258  weights.resize(heprup->nWeights(),
2259  std::make_pair(XWGTUP, (WeightInfo*)(0)));
2260  weights.front().first = XWGTUP;
2261  for ( int ii = 1, NN = weights.size(); ii < NN; ++ii )
2262  weights[ii].second = &heprup->weightinfo[ii - 1];
2263  double w = 0.0;
2264  int iii = 0;
2265  std::istringstream isss(tag.contents);
2266  while ( isss >> w )
2267  if ( ++iii < int(weights.size()) )
2268  weights[iii].first = w;
2269  else
2270  weights.push_back(std::make_pair(w, (WeightInfo*)(0)));
2271  }
2272  if ( tag.name == "weight" ) {
2273  namedweights.push_back(Weight(tag));
2274  }
2275  if ( tag.name == "rwgt" ) {
2276  for ( int j = 0, M = tag.tags.size(); j < M; ++j ) {
2277  if ( tag.tags[j]->name == "wgt" ) {
2278  namedweights.push_back(Weight(*tag.tags[j]));
2279  }
2280  }
2281  }
2282  else if ( tag.name == "clustering" ) {
2283  for ( int j = 0, M= tag.tags.size(); j < M; ++j ) {
2284  if ( tag.tags[j]->name == "clus" )
2285  clustering.push_back(Clus(*tag.tags[j]));
2286  }
2287  }
2288  else if ( tag.name == "pdfinfo" ) {
2289  pdfinfo = PDFInfo(tag, SCALUP);
2290  }
2291  else if ( tag.name == "scales" ) {
2292  scales = Scales(tag, SCALUP, NUP);
2293  }
2294 
2295  }
2296 
2297  for ( int i = 0, N = namedweights.size(); i < N; ++i ) {
2298  int indx = heprup->weightIndex(namedweights[i].name);
2299  if ( indx > 0 ) {
2300  weights[indx].first = namedweights[i].weights[0];
2301  namedweights[i].indices[0] = indx;
2302  } else {
2303  weights.push_back(std::make_pair(namedweights[i].weights[0],
2304  (WeightInfo*)(0)));
2305  namedweights[i].indices[0] = weights.size() - 1;
2306  }
2307  for ( int j = 1, M = namedweights[i].weights.size(); j < M; ++j ) {
2308  weights.push_back(std::make_pair(namedweights[i].weights[j],
2309  (WeightInfo*)(0)));
2310  namedweights[i].indices[j] = weights.size() - 1;
2311  }
2312  }
2313 
2314  }
2315 
2316  /**
2317  * Print out the event (group) as an XML tag.
2318  */
2319  void print(std::ostream & file) const {
2320 
2321  file << std::setprecision(heprup->dprec);
2322 
2323  using std::setw;
2324 
2325  if ( isGroup ) {
2326  file << "<eventgroup";
2327  if ( subevents.nreal > 0 )
2328  file << oattr("nreal", subevents.nreal);
2329  if ( subevents.ncounter > 0 )
2330  file << oattr("ncounter", subevents.ncounter);
2331  printattrs(file);
2332  file << ">\n";
2333  for ( int i = 0, N = subevents.size(); i < N; ++i )
2334  subevents[i]->print(file);
2335  file << "</eventgroup>\n";
2336  return;
2337  }
2338 
2339  file << "<event";
2340  if ( ntries > 1 ) file << oattr("ntries", ntries);
2341  printattrs(file);
2342  file << ">\n";
2343  file << " " << setw(4) << NUP
2344  << " " << setw(6) << IDPRUP
2345  << " " << setw(14) << XWGTUP
2346  << " " << setw(14) << SCALUP
2347  << " " << setw(14) << AQEDUP
2348  << " " << setw(14) << AQCDUP << "\n";
2349 
2350  for ( int i = 0; i < NUP; ++i )
2351  file << " " << setw(8) << IDUP[i]
2352  << " " << setw(2) << ISTUP[i]
2353  << " " << setw(4) << MOTHUP[i].first
2354  << " " << setw(4) << MOTHUP[i].second
2355  << " " << setw(4) << ICOLUP[i].first
2356  << " " << setw(4) << ICOLUP[i].second
2357  << " " << setw(14) << PUP[i][0]
2358  << " " << setw(14) << PUP[i][1]
2359  << " " << setw(14) << PUP[i][2]
2360  << " " << setw(14) << PUP[i][3]
2361  << " " << setw(14) << PUP[i][4]
2362  << " " << setw(1) << VTIMUP[i]
2363  << " " << setw(1) << SPINUP[i] << std::endl;
2364 
2365  if ( weights.size() > 0 ) {
2366  file << "<weights>";
2367  for ( int i = 1, N = weights.size(); i < N; ++i )
2368  file << " " << weights[i].first;
2369  file << "</weights>\n";
2370  }
2371 
2372  bool iswgt = false;
2373  for ( int i = 0, N = namedweights.size(); i < N; ++i ) {
2374  if ( namedweights[i].iswgt ) {
2375  if ( !iswgt ) file << "<rwgt>\n";
2376  iswgt = true;
2377  } else {
2378  if ( iswgt ) file << "</rwgt>\n";
2379  iswgt = false;
2380  }
2381  for ( int j = 0, M = namedweights[i].indices.size(); j < M; ++j )
2382  namedweights[i].weights[j] = weight(namedweights[i].indices[j]);
2383  namedweights[i].print(file);
2384  }
2385  if ( iswgt ) file << "</rwgt>\n";
2386 
2387  if ( !clustering.empty() ) {
2388  file << "<clustering>" << std::endl;
2389  for ( int i = 0, N = clustering.size(); i < N; ++i )
2390  clustering[i].print(file);
2391  file << "</clustering>" << std::endl;
2392  }
2393 
2394  pdfinfo.print(file);
2395  scales.print(file);
2396 
2397  file << hashline(junk) << "</event>\n";
2398 
2399  }
2400 
2401  /**
2402  * Reset the HEPEUP object (does not touch the sub events).
2403  */
2404  void reset() {
2405  setWeightInfo(0);
2406  NUP = 0;
2407  clustering.clear();
2408  weights.clear();
2409  }
2410 
2411  /**
2412  * Clear the HEPEUP object.
2413  */
2414  void clear() {
2415  reset();
2416  subevents.clear();
2417  }
2418 
2419  /**
2420  * Set the NUP variable, corresponding to the number of particles in
2421  * the current event, to \a nup, and resize all relevant vectors
2422  * accordingly.
2423  */
2424  void resize(int nup) {
2425  NUP = nup;
2426  resize();
2427  }
2428 
2429  /**
2430  * Return the total weight for this event (including all sub
2431  * evenets) for the given index.
2432  */
2433  double totalWeight(int i = 0) const {
2434  if ( subevents.empty() ) return weight(i);
2435  double w = 0.0;
2436  for ( int ii = 0, N = subevents.size(); ii < N; ++ii )
2437  w += subevents[ii]->weight(i);
2438  return w;
2439  }
2440 
2441  /**
2442  * Return the total weight for this event (including all sub
2443  * evenets) for the given weight name.
2444  */
2445  double totalWeight(std::string name) const {
2446  return totalWeight(heprup->weightIndex(name));
2447  }
2448 
2449  /**
2450  * Return the weight for the given index.
2451  */
2452  double weight(int i = 0) const {
2453  return weights[i].first;
2454  }
2455 
2456  /**
2457  * Return the weight for the given weight name.
2458  */
2459  double weight(std::string name) const {
2460  return weight(heprup->weightIndex(name));
2461  }
2462 
2463  /**
2464  * Set the weight with the given index.
2465  */
2466  void setWeight(int i, double w) {
2467  weights[i].first = w;
2468  }
2469  /**
2470  * Set the weight with the given name.
2471  */
2472  bool setWeight(std::string name, double w) {
2473  int i = heprup->weightIndex(name);
2474  if ( i >= int(weights.size()) ) return false;
2475  setWeight(i, w);
2476  return true;
2477  }
2478 
2479 
2480  /**
2481  * Assuming the NUP variable, corresponding to the number of
2482  * particles in the current event, is correctly set, resize the
2483  * relevant vectors accordingly.
2484  */
2485  void resize() {
2486  IDUP.resize(NUP);
2487  ISTUP.resize(NUP);
2488  MOTHUP.resize(NUP);
2489  ICOLUP.resize(NUP);
2490  PUP.resize(NUP, std::vector<double>(5));
2491  VTIMUP.resize(NUP);
2492  SPINUP.resize(NUP);
2493  }
2494 
2495  /**
2496  * Setup the current event to use weight i. If zero, the default
2497  * weight will be used.
2498  */
2499  bool setWeightInfo(unsigned int i) {
2500  if ( i >= weights.size() ) return false;
2501  if ( currentWeight ) {
2502  scales.mur /= currentWeight->mur;
2503  scales.muf /= currentWeight->muf;
2504  heprup->PDFGUP = PDFGUPsave;
2505  heprup->PDFSUP = PDFSUPsave;
2506  }
2507  XWGTUP = weights[i].first;
2508  currentWeight = weights[i].second;
2509  if ( currentWeight ) {
2510  scales.mur *= currentWeight->mur;
2511  scales.muf *= currentWeight->muf;
2512  PDFGUPsave = heprup->PDFGUP;
2513  PDFSUPsave = heprup->PDFSUP;
2514  if ( currentWeight->pdf ) {
2515  heprup->PDFGUP.first = heprup->PDFGUP.second = 0;
2516  heprup->PDFSUP.first = heprup->PDFSUP.second = currentWeight->pdf;
2517  }
2518  if ( currentWeight->pdf2 ) {
2519  heprup->PDFSUP.second = currentWeight->pdf2;
2520  }
2521 
2522  }
2523  return true;
2524  }
2525 
2526  /**
2527  * Setup the current event to use sub event i. If zero, no sub event
2528  * will be chsen.
2529  */
2530  bool setSubEvent(unsigned int i) {
2531  if ( i > subevents.size() || subevents.empty() ) return false;
2532  if ( i == 0 ) {
2533  reset();
2534  weights = subevents[0]->weights;
2535  for ( int ii = 1, N = subevents.size(); ii < N; ++ii )
2536  for ( int j = 0, M = weights.size(); j < M; ++j )
2537  weights[j].first += subevents[ii]->weights[j].first;
2538  currentWeight = 0;
2539  } else {
2540  setEvent(*subevents[i - 1]);
2541  }
2542  return true;
2543  }
2544 
2545 public:
2546 
2547  /**
2548  * The number of particle entries in the current event.
2549  */
2550  int NUP;
2551 
2552  /**
2553  * The subprocess code for this event (as given in LPRUP).
2554  */
2555  int IDPRUP;
2556 
2557  /**
2558  * The weight for this event.
2559  */
2560  double XWGTUP;
2561 
2562  /**
2563  * The PDF weights for the two incoming partons. Note that this
2564  * variable is not present in the current LesHouches accord
2565  * (<A HREF="http://arxiv.org/abs/hep-ph/0109068">hep-ph/0109068</A>),
2566  * hopefully it will be present in a future accord.
2567  */
2568  std::pair<double,double> XPDWUP;
2569 
2570  /**
2571  * The scale in GeV used in the calculation of the PDF's in this
2572  * event.
2573  */
2574  double SCALUP;
2575 
2576  /**
2577  * The value of the QED coupling used in this event.
2578  */
2579  double AQEDUP;
2580 
2581  /**
2582  * The value of the QCD coupling used in this event.
2583  */
2584  double AQCDUP;
2585 
2586  /**
2587  * The PDG id's for the particle entries in this event.
2588  */
2589  std::vector<long> IDUP;
2590 
2591  /**
2592  * The status codes for the particle entries in this event.
2593  */
2594  std::vector<int> ISTUP;
2595 
2596  /**
2597  * Indices for the first and last mother for the particle entries in
2598  * this event.
2599  */
2600  std::vector< std::pair<int,int> > MOTHUP;
2601 
2602  /**
2603  * The colour-line indices (first(second) is (anti)colour) for the
2604  * particle entries in this event.
2605  */
2606  std::vector< std::pair<int,int> > ICOLUP;
2607 
2608  /**
2609  * Lab frame momentum (Px, Py, Pz, E and M in GeV) for the particle
2610  * entries in this event.
2611  */
2612  std::vector< std::vector<double> > PUP;
2613 
2614  /**
2615  * Invariant lifetime (c*tau, distance from production to decay in
2616  * mm) for the particle entries in this event.
2617  */
2618  std::vector<double> VTIMUP;
2619 
2620  /**
2621  * Spin info for the particle entries in this event given as the
2622  * cosine of the angle between the spin vector of a particle and the
2623  * 3-momentum of the decaying particle, specified in the lab frame.
2624  */
2625  std::vector<double> SPINUP;
2626 
2627  /**
2628  * A pointer to the current HEPRUP object.
2629  */
2631 
2632  /**
2633  * The current weight info object.
2634  */
2636 
2637  /**
2638  * The weights associated with this event
2639  */
2640  std::vector<Weight> namedweights;
2641 
2642  /**
2643  * The weights for this event and their corresponding WeightInfo object.
2644  */
2645  std::vector< std::pair<double, const WeightInfo *> > weights;
2646 
2647  /**
2648  * Contents of the clustering tag.
2649  */
2650  std::vector<Clus> clustering;
2651 
2652  /**
2653  * Contents of the pdfinfo tag.
2654  */
2656 
2657  /**
2658  * Saved information about pdfs if different in a selected weight.
2659  */
2660  std::pair<int,int> PDFGUPsave;
2661 
2662  /**
2663  * Saved information about pdfs if different in a selected weight.
2664  */
2665  std::pair<int,int> PDFSUPsave;
2666 
2667 
2668  /**
2669  * Contents of the scales tag
2670  */
2672 
2673  /**
2674  * The number of attempts the ME generator did before accepting this
2675  * event.
2676  */
2677  int ntries;
2678 
2679  /**
2680  * Is this an event or an event group?
2681  */
2682  bool isGroup;
2683 
2684  /**
2685  * If this is not a single event, but an event group, the events
2686  * included in the group are in this vector;
2687  */
2689 
2690  /**
2691  * Save junk stuff in events just to be on the safe side
2692  */
2693  std::string junk;
2694 
2695 };
2696 
2697 
2698 // Destructor implemented here.
2699 
2700 inline void EventGroup::clear() {
2701  while ( size() > 0 ) {
2702  delete back();
2703  pop_back();
2704  }
2705 }
2706 
2708  clear();
2709 }
2710 
2712  : std::vector<HEPEUP*>(eg.size()) {
2713  for ( int i = 0, N = eg.size(); i < N; ++i ) at(i) = new HEPEUP(*eg.at(i));
2714 }
2715 
2717  if ( &x == this ) return *this;
2718  clear();
2719  nreal = x.nreal;
2720  ncounter = x.ncounter;
2721  for ( int i = 0, N = x.size(); i < N; ++i ) push_back(new HEPEUP(*x.at(i)));
2722  return *this;
2723 }
2724 
2725 
2726 /**
2727  * The Reader class is initialized with a stream from which to read a
2728  * version 1/2 Les Houches Accord event file. In the constructor of
2729  * the Reader object the optional header information is read and then
2730  * the mandatory init is read. After this the whole header block
2731  * including the enclosing lines with tags are available in the public
2732  * headerBlock member variable. Also the information from the init
2733  * block is available in the heprup member variable and any additional
2734  * comment lines are available in initComments. After each successful
2735  * call to the readEvent() function the standard Les Houches Accord
2736  * information about the event is available in the hepeup member
2737  * variable and any additional comments in the eventComments
2738  * variable. A typical reading sequence would look as follows:
2739  *
2740  *
2741  */
2742 class Reader {
2743 
2744 public:
2745 
2746  /**
2747  * Initialize the Reader with a stream from which to read an event
2748  * file. After the constructor is called the whole header block
2749  * including the enclosing lines with tags are available in the
2750  * public headerBlock member variable. Also the information from the
2751  * init block is available in the heprup member variable and any
2752  * additional comment lines are available in initComments.
2753  *
2754  * @param is the stream to read from.
2755  */
2756  Reader(std::istream & is)
2757  : file(&is), currevent(-1),
2758  curreventfile(-1), currfileevent(-1), dirpath("") {
2759  init();
2760  }
2761 
2762  /**
2763  * Initialize the Reader with a filename from which to read an event
2764  * file. After the constructor is called the whole header block
2765  * including the enclosing lines with tags are available in the
2766  * public headerBlock member variable. Also the information from the
2767  * init block is available in the heprup member variable and any
2768  * additional comment lines are available in initComments.
2769  *
2770  * @param filename the name of the file to read from.
2771  */
2772  Reader(std::string filename)
2773  : intstream(filename.c_str()), file(&intstream), currevent(-1),
2774  curreventfile(-1), currfileevent(-1), dirpath("") {
2775 
2776  size_t slash = filename.find_last_of('/');
2777  if ( slash != std::string::npos ) dirpath = filename.substr(0, slash + 1);
2778  init();
2779  }
2780 
2781 private:
2782 
2783  /**
2784  * Used internally in the constructors to read header and init
2785  * blocks.
2786  */
2787  void init() {
2788 
2789  // initialize reading from multi-file runs.
2790  initfile = file;
2791 
2792  bool readingHeader = false;
2793  bool readingInit = false;
2794 
2795  // Make sure we are reading a LHEF file:
2796  getline();
2797  if ( !currentFind("<LesHouchesEvents") )
2798  throw std::runtime_error
2799  ("Tried to read a file which does not start with the "
2800  "LesHouchesEvents tag.");
2801  version = 1;
2802  if ( currentFind("version=\"3" ) )
2803  version = 3;
2804  else if ( currentFind("version=\"2" ) )
2805  version = 2;
2806  else if ( !currentFind("version=\"1" ) )
2807  throw std::runtime_error
2808  ("Tried to read a LesHouchesEvents file which is above version 3.");
2809 
2810  // Loop over all lines until we hit the </init> tag.
2811  while ( getline() && !currentFind("</init>") ) {
2812  if ( currentFind("<header") ) {
2813  // We have hit the header block, so we should dump this and
2814  // all following lines to headerBlock until we hit the end of
2815  // it.
2816  readingHeader = true;
2817  headerBlock = currentLine + "\n";
2818  }
2819  else if ( currentFind("<init>") ) {
2820  // We have hit the init block
2821  readingInit = true;
2822  initComments = currentLine + "\n";
2823  }
2824  else if ( currentFind("</header>") ) {
2825  // The end of the header block. Dump this line as well to the
2826  // headerBlock and we're done.
2827  readingHeader = false;
2828  headerBlock += currentLine + "\n";
2829  }
2830  else if ( readingHeader ) {
2831  // We are in the process of reading the header block. Dump the
2832  // line to haderBlock.
2833  headerBlock += currentLine + "\n";
2834  }
2835  else if ( readingInit ) {
2836  // Here we found a comment line. Dump it to initComments.
2837  initComments += currentLine + "\n";
2838  }
2839  else {
2840  // We found some other stuff outside the standard tags.
2841  outsideBlock += currentLine + "\n";
2842  }
2843  }
2844  if ( !currentFind("</init>") )
2845  throw std::runtime_error("Found incomplete init tag in "
2846  "Les Houches file.");
2847  initComments += currentLine + "\n";
2848  std::vector<XMLTag*> tags = XMLTag::findXMLTags(initComments);
2849  for ( int i = 0, N = tags.size(); i < N; ++i )
2850  if ( tags[i]->name == "init" ) {
2851  heprup = HEPRUP(*tags[i], version);
2852  break;
2853  }
2854  XMLTag::deleteAll(tags);
2855 
2856  if ( !heprup.eventfiles.empty() ) openeventfile(0);
2857 
2858  }
2859 
2860 public:
2861 
2862  /**
2863  * Read an event from the file and store it in the hepeup
2864  * object. Optional comment lines are stored i the eventComments
2865  * member variable.
2866  * @return true if the read sas successful.
2867  */
2868  bool readEvent() {
2869 
2870  // Check if the initialization was successful. Otherwise we will
2871  // not read any events.
2872  if ( heprup.NPRUP < 0 ) return false;
2873 
2874  std::string eventLines;
2875  int inEvent = 0;
2876 
2877  // Keep reading lines until we hit the end of an event or event group.
2878  while ( getline() ) {
2879  if ( inEvent ) {
2880  eventLines += currentLine + "\n";
2881  if ( inEvent == 1 && currentFind("</event>") ) break;
2882  if ( inEvent == 2 && currentFind("</eventgroup>") ) break;
2883  }
2884  else if ( currentFind("<eventgroup") ) {
2885  eventLines += currentLine + "\n";
2886  inEvent = 2;
2887  }
2888  else if ( currentFind("<event") ) {
2889  eventLines += currentLine + "\n";
2890  inEvent = 1;
2891  }
2892  else {
2893  outsideBlock += currentLine + "\n";
2894  }
2895  }
2896  if ( ( inEvent == 1 && !currentFind("</event>") ) ||
2897  ( inEvent == 2 && !currentFind("</eventgroup>") ) ) {
2898  if ( heprup.eventfiles.empty() ||
2899  ++curreventfile >= int(heprup.eventfiles.size()) ) return false;
2900  openeventfile(curreventfile);
2901  return readEvent();
2902  }
2903 
2904  std::vector<XMLTag*> tags = XMLTag::findXMLTags(eventLines);
2905 
2906  for ( int i = 0, N = tags.size(); i < N ; ++i ) {
2907  if ( tags[i]->name == "event" || tags[i]->name == "eventgroup" ) {
2908  hepeup = HEPEUP(*tags[i], heprup);
2909  XMLTag::deleteAll(tags);
2910  ++currevent;
2911  if ( curreventfile >= 0 ) ++currfileevent;
2912  return true;
2913  }
2914  }
2915 
2916  if ( !heprup.eventfiles.empty() &&
2917  ++curreventfile < int(heprup.eventfiles.size()) ) {
2918  openeventfile(curreventfile);
2919  return readEvent();
2920  }
2921 
2922  XMLTag::deleteAll(tags);
2923  return false;
2924 
2925  }
2926 
2927  /**
2928  * Open the efentfile with index ifile. If another eventfile is
2929  * being read, its remaining contents is discarded. This is a noop
2930  * if current read session is not a multi-file run.
2931  */
2932  void openeventfile(int ifile) {
2933  std::cerr << "opening file " << ifile << std::endl;
2934  efile.close();
2935  std::string fname = heprup.eventfiles[ifile].filename;
2936  if ( fname[0] != '/' ) fname = dirpath + fname;
2937  efile.open(fname.c_str());
2938  if ( !efile ) throw std::runtime_error("Could not open event file " +
2939  fname);
2940  file = &efile;
2941  curreventfile = ifile;
2942  currfileevent = 0;
2943  }
2944 
2945 protected:
2946 
2947  /**
2948  * Used internally to read a single line from the stream.
2949  */
2950  bool getline() {
2951  return ( (bool)std::getline(*file, currentLine) );
2952  }
2953 
2954  /**
2955  * @return true if the current line contains the given string.
2956  */
2957  bool currentFind(std::string str) const {
2958  return currentLine.find(str) != std::string::npos;
2959  }
2960 
2961 protected:
2962 
2963  /**
2964  * A local stream which is unused if a stream is supplied from the
2965  * outside.
2966  */
2967  std::ifstream intstream;
2968 
2969  /**
2970  * The stream we are reading from. This may be a pointer to an
2971  * external stream or the internal intstream, or a separate event
2972  * file from a multi-file run
2973  */
2974  std::istream * file;
2975 
2976  /**
2977  * The original stream from where we read the init block.
2978  */
2979  std::istream * initfile;
2980 
2981  /**
2982  * A separate stream for reading multi-file runs.
2983  */
2984  std::ifstream efile;
2985 
2986  /**
2987  * The last line read in from the stream in getline().
2988  */
2989  std::string currentLine;
2990 
2991 public:
2992 
2993  /**
2994  * XML file version
2995  */
2996  int version;
2997 
2998  /**
2999  * All lines (since the last readEvent()) outside the header, init
3000  * and event tags.
3001  */
3002  std::string outsideBlock;
3003 
3004  /**
3005  * All lines from the header block.
3006  */
3007  std::string headerBlock;
3008 
3009  /**
3010  * The standard init information.
3011  */
3013 
3014  /**
3015  * Additional comments found in the init block.
3016  */
3017  std::string initComments;
3018 
3019  /**
3020  * The standard information about the last read event.
3021  */
3023 
3024  /**
3025  * Additional comments found with the last read event.
3026  */
3027  std::string eventComments;
3028 
3029  /**
3030  * The number of the current event (starting from 1).
3031  */
3033 
3034  /**
3035  * The current event file being read from (-1 means there are no
3036  * separate event files).
3037  */
3039 
3040  /**
3041  * The number of the current event in the current event file.
3042  */
3044 
3045  /**
3046  * The directory from where we are reading files.
3047  */
3048  std::string dirpath;
3049 
3050 private:
3051 
3052  /**
3053  * The default constructor should never be used.
3054  */
3055  Reader();
3056 
3057  /**
3058  * The copy constructor should never be used.
3059  */
3060  Reader(const Reader &);
3061 
3062  /**
3063  * The Reader cannot be assigned to.
3064  */
3065  Reader & operator=(const Reader &);
3066 
3067 };
3068 
3069 /**
3070  * The Writer class is initialized with a stream to which to write a
3071  * version 1.0 Les Houches Accord event file. In the constructor of
3072  * the Writer object the main XML tag is written out, with the
3073  * corresponding end tag is written in the destructor. After a Writer
3074  * object has been created, it is possible to assign standard init
3075  * information in the heprup member variable. In addition any XML
3076  * formatted information can be added to the headerBlock member
3077  * variable (directly or via the addHeader() function). Further
3078  * comment line (beginning with a <code>#</code> character) can be
3079  * added to the initComments variable (directly or with the
3080  * addInitComment() function). After this information is set, it
3081  * should be written out to the file with the init() function.
3082  *
3083  * Before each event is written out with the writeEvent() function,
3084  * the standard event information can then be assigned to the hepeup
3085  * variable and optional comment lines (beginning with a
3086  * <code>#</code> character) may be given to the eventComments
3087  * variable (directly or with the addEventComment() function).
3088  *
3089  */
3090 class Writer {
3091 
3092 public:
3093 
3094  /**
3095  * Create a Writer object giving a stream to write to.
3096  * @param os the stream where the event file is written.
3097  */
3098  Writer(std::ostream & os)
3099  : file(&os), initfile(&os), dirpath("") { }
3100 
3101  /**
3102  * Create a Writer object giving a filename to write to.
3103  * @param filename the name of the event file to be written.
3104  */
3105  Writer(std::string filename)
3106  : intstream(filename.c_str()), file(&intstream), initfile(&intstream),
3107  dirpath("") {
3108  size_t slash = filename.find_last_of('/');
3109  if ( slash != std::string::npos ) dirpath = filename.substr(0, slash + 1);
3110  }
3111 
3112  /**
3113  * The destructor writes out the final XML end-tag.
3114  */
3116  file = initfile;
3117  if ( !heprup.eventfiles.empty() ) {
3118  if ( curreventfile >= 0 &&
3119  curreventfile < int(heprup.eventfiles.size()) &&
3120  heprup.eventfiles[curreventfile].neve < 0 )
3121  heprup.eventfiles[curreventfile].neve = currfileevent;
3122  writeinit();
3123  }
3124  *file << "</LesHouchesEvents>" << std::endl;
3125  }
3126 
3127  /**
3128  * Add header lines consisting of XML code with this stream.
3129  */
3130  std::ostream & headerBlock() {
3131  return headerStream;
3132  }
3133 
3134  /**
3135  * Add comment lines to the init block with this stream.
3136  */
3137  std::ostream & initComments() {
3138  return initStream;
3139  }
3140 
3141  /**
3142  * Add comment lines to the next event to be written out with this stream.
3143  */
3144  std::ostream & eventComments() {
3145  return eventStream;
3146  }
3147 
3148  /**
3149  * Initialize the writer.
3150  */
3151  void init() {
3152  if ( heprup.eventfiles.empty() ) writeinit();
3153  lastevent = 0;
3154  curreventfile = currfileevent = -1;
3155  if ( !heprup.eventfiles.empty() ) openeventfile(0);
3156  }
3157 
3158  /**
3159  * Open a new event file, possibly closing a previous opened one.
3160  */
3161  bool openeventfile(int ifile) {
3162  if ( heprup.eventfiles.empty() ) return false;
3163  if ( ifile < 0 || ifile >= int(heprup.eventfiles.size()) ) return false;
3164  if ( curreventfile >= 0 ) {
3165  EventFile & ef = heprup.eventfiles[curreventfile];
3166  if ( ef.neve > 0 && ef.neve != currfileevent )
3167  std::cerr << "LHEF::Writer number of events in event file "
3168  << ef.filename << " does not match the given number."
3169  << std::endl;
3170  ef.neve = currfileevent;
3171  }
3172  efile.close();
3173  std::string fname = heprup.eventfiles[ifile].filename;
3174  if ( fname[0] != '/' ) fname = dirpath + fname;
3175  efile.open(fname.c_str());
3176  if ( !efile ) throw std::runtime_error("Could not open event file " +
3177  fname);
3178  std::cerr << "Opened event file " << fname << std::endl;
3179  file = &efile;
3180  curreventfile = ifile;
3181  currfileevent = 0;
3182  return true;
3183  }
3184 
3185 
3186  /**
3187  * Write out an optional header block followed by the standard init
3188  * block information together with any comment lines.
3189  */
3190  void writeinit() {
3191 
3192  // Write out the standard XML tag for the event file.
3193  if ( heprup.version == 3 )
3194  *file << "<LesHouchesEvents version=\"3.0\">\n";
3195  else if ( heprup.version == 2 )
3196  *file << "<LesHouchesEvents version=\"2.0\">\n";
3197  else
3198  *file << "<LesHouchesEvents version=\"1.0\">\n";
3199 
3200 
3201  *file << std::setprecision(10);
3202 
3203  using std::setw;
3204 
3205  std::string headBlock = headerStream.str();
3206  if ( headBlock.length() ) {
3207  if ( headBlock.find("<header>") == std::string::npos )
3208  *file << "<header>\n";
3209  if ( headBlock[headBlock.length() - 1] != '\n' )
3210  headBlock += '\n';
3211  *file << headBlock;
3212  if ( headBlock.find("</header>") == std::string::npos )
3213  *file << "</header>\n";
3214  }
3215 
3216  heprup.print(*file);
3217 
3218  }
3219 
3220  /**
3221  * Write the current HEPEUP object to the stream;
3222  */
3223  void writeEvent() {
3224 
3225  if ( !heprup.eventfiles.empty() ) {
3226  if ( currfileevent == heprup.eventfiles[curreventfile].neve &&
3227  curreventfile + 1 < int(heprup.eventfiles.size()) )
3228  openeventfile(curreventfile + 1);
3229  }
3230 
3231  hepeup.print(*file);
3232 
3233  ++lastevent;
3234  ++currfileevent;
3235  }
3236 
3237 protected:
3238 
3239  /**
3240  * A local stream which is unused if a stream is supplied from the
3241  * outside.
3242  */
3243  std::ofstream intstream;
3244 
3245  /**
3246  * The stream we are writing to. This may be a reference to an
3247  * external stream or the internal intstream.
3248  */
3249  std::ostream * file;
3250 
3251  /**
3252  * The original stream from where we read the init block.
3253  */
3254  std::ostream * initfile;
3255 
3256  /**
3257  * A separate stream for reading multi-file runs.
3258  */
3259  std::ofstream efile;
3260 
3261  /**
3262  * The number of the last event written (starting from 1).
3263  */
3265 
3266  /**
3267  * The current event file being written to (-1 means there are no
3268  * separate event files).
3269  */
3271 
3272  /**
3273  * The number of the current event in the current event file.
3274  */
3276 
3277  /**
3278  * The directory from where we are reading files.
3279  */
3280  std::string dirpath;
3281 
3282 public:
3283 
3284  /**
3285  * Stream to add all lines in the header block.
3286  */
3287  std::ostringstream headerStream;
3288 
3289  /**
3290  * The standard init information.
3291  */
3293 
3294  /**
3295  * Stream to add additional comments to be put in the init block.
3296  */
3297  std::ostringstream initStream;
3298 
3299  /**
3300  * The standard information about the event we will write next.
3301  */
3303 
3304  /**
3305  * Stream to add additional comments to be written together the next event.
3306  */
3307  std::ostringstream eventStream;
3308 
3309 private:
3310 
3311  /**
3312  * The default constructor should never be used.
3313  */
3314  Writer();
3315 
3316  /**
3317  * The copy constructor should never be used.
3318  */
3319  Writer(const Writer &);
3320 
3321  /**
3322  * The Writer cannot be assigned to.
3323  */
3324  Writer & operator=(const Writer &);
3325 
3326 };
3327 
3328 }
3329 
3330 /* \example LHEFCat.cc This is a main function which simply reads a
3331  Les Houches Event File from the standard input and writes it again
3332  to the standard output.
3333  This file can be downloaded from
3334  <A HREF="http://www.thep.lu.se/~leif/LHEF/LHEFCat.cc">here</A>.
3335  There are also two sample event files,
3336  <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>
3337  to try it on.
3338 */
3339 
3340 /* \mainpage Les Houches Event File
3341 
3342 Here are some example classes for reading and writing Les Houches
3343 Event Files according to the
3344 <A HREF="http://www.thep.lu.se/~torbjorn/lhef/lhafile2.pdf">proposal</A>
3345 by Torbj&ouml;rn Sj&ouml;strand discussed at the
3346 <A HREF="http://mc4lhc06.web.cern.ch/mc4lhc06/">MC4LHC</A>
3347 workshop at CERN 2006.
3348 
3349 The classes has now been updated to handle the suggested version 3 of
3350 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>).
3351 
3352 There is a whole set of classes available in a single header file
3353 called <A
3354 HREF="http://www.thep.lu.se/~leif/LHEF/LHEF.h">LHEF.h</A>. The idea is
3355 that matrix element or parton shower generators will implement their
3356 own wrapper using these classes as containers.
3357 
3358 The two classes LHEF::HEPRUP and LHEF::HEPEUP are simple container
3359 classes which contain the same information as the Les Houches standard
3360 Fortran common blocks with the same names. They also contain the extra
3361 information defined in version 3 in the standard. The other two main
3362 classes are called LHEF::Reader and LHEF::Writer and are used to read
3363 and write Les Houches Event Files
3364 
3365 Here are a few <A HREF="examples.html">examples</A> of how to use the
3366 classes:
3367 
3368 \namespace LHEF The LHEF namespace contains some example classes for reading and writing Les Houches
3369 Event Files according to the
3370 <A HREF="http://www.thep.lu.se/~torbjorn/lhef/lhafile2.pdf">proposal</A>
3371 by Torbj&ouml;rn Sj&ouml;strand discussed at the
3372 <A HREF="http://mc4lhc06.web.cern.ch/mc4lhc06/">MC4LHC</A>
3373 workshop at CERN 2006.
3374 
3375 
3376 
3377  */
3378 
3379 
3380 #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:3243
std::string weightNameHepMC(int i) const
Definition: LHEF.h:1751
void setWeight(int i, double w)
Definition: LHEF.h:2466
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:2618
bool getattr(std::string n, double &v, bool erase=true)
Definition: LHEF.h:368
double totalWeight(std::string name) const
Definition: LHEF.h:2445
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:3259
std::pair< int, int > PDFGUPsave
Definition: LHEF.h:2660
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:2472
Reader(std::string filename)
Definition: LHEF.h:2772
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:2594
bool hasInfo() const
Definition: LHEF.h:1452
HEPEUP hepeup
Definition: LHEF.h:3022
double XWGTUP
Definition: LHEF.h:2560
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:2979
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:2584
std::ostream & eventComments()
Definition: LHEF.h:3144
std::string headerBlock
Definition: LHEF.h:3007
std::ostringstream headerStream
Definition: LHEF.h:3287
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:3130
void resize()
Definition: LHEF.h:1888
double mur
Definition: LHEF.h:1110
std::vector< Scale > scales
Definition: LHEF.h:1533
STL namespace.
std::vector< std::pair< int, int > > MOTHUP
Definition: LHEF.h:2600
bool readEvent()
Definition: LHEF.h:2868
void reset()
Definition: LHEF.h:2404
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:3017
int dprec
Definition: LHEF.h:2057
void init()
Definition: LHEF.h:3151
std::ostringstream initStream
Definition: LHEF.h:3297
~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:2550
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:3137
AttributeMap attr
Definition: LHEF.h:124
double muf
Definition: LHEF.h:1105
std::pair< double, double > XPDWUP
Definition: LHEF.h:2568
Writer(std::ostream &os)
Definition: LHEF.h:3098
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:2967
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:3292
std::pair< int, int > PDFSUPsave
Definition: LHEF.h:2665
bool isGroup
Definition: LHEF.h:2682
std::string dirpath
Definition: LHEF.h:3048
int version
Definition: LHEF.h:2996
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:3105
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:3032
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:3270
bool setSubEvent(unsigned int i)
Definition: LHEF.h:2530
double weight(std::string name) const
Definition: LHEF.h:2459
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:2635
long p1
Definition: LHEF.h:1580
std::string type
Definition: LHEF.h:879
int currfileevent
Definition: LHEF.h:3275
int lastevent
Definition: LHEF.h:3264
PDFInfo pdfinfo
Definition: LHEF.h:2655
std::ostringstream eventStream
Definition: LHEF.h:3307
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:3254
double scale
Definition: LHEF.h:1610
HEPEUP(const XMLTag &tagin, HEPRUP &heprupin)
Definition: LHEF.h:2198
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:2984
std::string outsideBlock
Definition: LHEF.h:3002
std::vector< std::pair< double, const WeightInfo * > > weights
Definition: LHEF.h:2645
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:2671
int loops
Definition: LHEF.h:960
bool getline()
Definition: LHEF.h:2950
TagBase(const AttributeMap &attr, std::string conts=std::string())
Definition: LHEF.h:360
bool currentFind(std::string str) const
Definition: LHEF.h:2957
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:2932
std::vector< std::vector< double > > PUP
Definition: LHEF.h:2612
std::string np2
Definition: LHEF.h:899
void print(std::ostream &file) const
Definition: LHEF.h:1765
void writeEvent()
Definition: LHEF.h:3223
HEPEUP hepeup
Definition: LHEF.h:3302
std::map< long, MergeInfo > mergeinfo
Definition: LHEF.h:2021
std::vector< WeightGroup > weightgroup
Definition: LHEF.h:2042
std::vector< long > IDUP
Definition: LHEF.h:2589
int curreventfile
Definition: LHEF.h:3038
void print(std::ostream &file) const
Definition: LHEF.h:1355
long neve
Definition: LHEF.h:564
std::string dirpath
Definition: LHEF.h:3280
std::vector< Weight > namedweights
Definition: LHEF.h:2640
EventGroup & operator=(const EventGroup &)
Definition: LHEF.h:2716
Weight(const XMLTag &tag)
Definition: LHEF.h:1179
std::vector< Cut > cuts
Definition: LHEF.h:2006
std::vector< double > SPINUP
Definition: LHEF.h:2625
void closetag(std::ostream &file, std::string tag) const
Definition: LHEF.h:445
std::istream * file
Definition: LHEF.h:2974
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:2579
std::string eventComments
Definition: LHEF.h:3027
std::string junk
Definition: LHEF.h:2693
Reader(std::istream &is)
Definition: LHEF.h:2756
bool getattr(std::string n, bool &v, bool erase=true)
Definition: LHEF.h:382
void clear()
Definition: LHEF.h:2414
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:2499
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:2700
void print(std::ostream &file) const
Definition: LHEF.h:1272
bool openeventfile(int ifile)
Definition: LHEF.h:3161
std::string currentLine
Definition: LHEF.h:2989
void writeinit()
Definition: LHEF.h:3190
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:2677
std::ostream * file
Definition: LHEF.h:3249
int currfileevent
Definition: LHEF.h:3043
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:2424
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:2787
std::set< long > p1
Definition: LHEF.h:884
std::vector< int > indices
Definition: LHEF.h:1242
double SCALUP
Definition: LHEF.h:2574
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:2650
HEPRUP heprup
Definition: LHEF.h:3012
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:2606
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
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:2433
void resize()
Definition: LHEF.h:2485
EventGroup subevents
Definition: LHEF.h:2688
void print(std::ostream &file) const
Definition: LHEF.h:2319
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:2630
int p2
Definition: LHEF.h:1289
double weight(int i=0) const
Definition: LHEF.h:2452
std::map< std::string, std::string > AttributeMap
Definition: LHEF.h:97
int IDPRUP
Definition: LHEF.h:2555