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