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