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