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