StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
LHEF3.h
1 // LHEF3.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // This file is written by Stefan Prestel. The code evolved from
7 // a LHEF 2.0 reader supplied by Leif Lonnblad.
8 // LHEF3.h contains the main class for LHEF 3.0 functionalities.
9 // Header file.
10 
11 #ifndef Pythia8_LHEF3_H
12 #define Pythia8_LHEF3_H
13 
14 #include "Pythia8/PythiaStdlib.h"
15 #include "Pythia8/Streams.h"
16 #include <stdexcept>
17 
18 namespace Pythia8 {
19 
20 //==========================================================================
21 
22 // The XMLTag struct is used to represent all information within an XML tag.
23 // It contains the attributes as a map, any sub-tags as a vector of pointers
24 // to other XMLTag objects, and any other information as a single string.
25 // The XMLTag struct written by Leif Lonnblad.
26 
27 struct XMLTag {
28 
29  // Convenient typdef.
30  typedef string::size_type pos_t;
31 
32  // Convenient alias for npos.
33  static const pos_t end;
34 
35  // The destructor also destroys any sub-tags.
36  ~XMLTag() {
37  for ( int i = 0, N = tags.size(); i < N; ++i )
38  if (tags[i]) delete tags[i];
39  }
40 
41  // The name of this tag.
42  string name;
43 
44  // The attributes of this tag.
45  map<string,string> attr;
46 
47  // A vector of sub-tags.
48  vector<XMLTag*> tags;
49 
50  // The contents of this tag.
51  string contents;
52 
53  // Find an attribute named n and set the double variable v to
54  // the corresponding value. Return false if no attribute was found.
55  bool getattr(string n, double & v) const {
56  map<string,string>::const_iterator it = attr.find(n);
57  if ( it == attr.end() ) return false;
58  v = atof(it->second.c_str());
59  return true;
60  }
61 
62  // Find an attribute named n and set the bool variable v to true if the
63  // corresponding value is "yes". Return false if no attribute was found.
64  bool getattr(string n, bool & v) const {
65  map<string,string>::const_iterator it = attr.find(n);
66  if ( it == attr.end() ) return false;
67  if ( it->second == "yes" ) v = true;
68  return true;
69  }
70 
71  // Find an attribute named n and set the long variable v to the
72  // corresponding value. Return false if no attribute was found.
73  bool getattr(string n, long & v) const {
74  map<string,string>::const_iterator it = attr.find(n);
75  if ( it == attr.end() ) return false;
76  v = atoi(it->second.c_str());
77  return true;
78  }
79 
80  // Find an attribute named n and set the long variable v to the
81  // corresponding value. Return false if no attribute was found.
82  bool getattr(string n, int & v) const {
83  map<string,string>::const_iterator it = attr.find(n);
84  if ( it == attr.end() ) return false;
85  v = int(atoi(it->second.c_str()));
86  return true;
87  }
88 
89  // Find an attribute named n and set the string variable v to the
90  // corresponding value. Return false if no attribute was found.
91  bool getattr(string n, string & v) const {
92  map<string,string>::const_iterator it = attr.find(n);
93  if ( it == attr.end() ) return false;
94  v = it->second;
95  return true;
96  }
97 
98  // Scan the given string and return all XML tags found as a vector
99  // of pointers to XMLTag objects.
100  static vector<XMLTag*> findXMLTags(string str,
101  string * leftover = 0) {
102  vector<XMLTag*> tags;
103  pos_t curr = 0;
104 
105  while ( curr != end ) {
106 
107  // Find the first tag.
108  pos_t begin = str.find("<", curr);
109 
110  // Skip tags in lines beginning with #.
111  pos_t lastbreak_before_begin = str.find_last_of("\n",begin);
112  pos_t lastpound_before_begin = str.find_last_of("#",begin);
113  // Logic: Last newline before begin was before last pound sign (or there
114  // was no last newline at all, i.e. this is the special case of the
115  // first line), and hence the pound sign was before the tag was opened
116  // (at begin) with '<'. Thus, skip forward to next new line.
117  if ( (lastbreak_before_begin < lastpound_before_begin
118  || lastbreak_before_begin == end)
119  && begin > lastpound_before_begin) {
120  pos_t endcom = str.find_first_of("\n",begin);
121  if ( endcom == end ) {
122  if ( leftover ) *leftover += str.substr(curr);
123  return tags;
124  }
125  if ( leftover ) *leftover += str.substr(curr, endcom - curr);
126  curr = endcom;
127  continue;
128  }
129 
130  // Skip xml-style comments.
131  if ( begin != end && str.find("<!--", curr) == begin ) {
132  pos_t endcom = str.find("-->", begin);
133  if ( endcom == end ) {
134  if ( leftover ) *leftover += str.substr(curr);
135  return tags;
136  }
137  if ( leftover ) *leftover += str.substr(curr, endcom - curr);
138  curr = endcom;
139  continue;
140  }
141 
142  // Also skip CDATA statements.
143  // Used for text data that should not be parsed by the XML parser.
144  // (e.g., JavaScript code contains a lot of "<" or "&" characters
145  // which XML would erroneously interpret as the start of a new
146  // element or the start of a character entity, respectively.)
147  // See eg http://www.w3schools.com/xml/xml_cdata.asp
148  if ( str.find("<![CDATA[", curr) == begin ) {
149  pos_t endcom = str.find("]]>", begin);
150  if ( endcom == end ) {
151  if ( leftover ) *leftover += str.substr(curr);
152  return tags;
153  }
154  if ( leftover ) *leftover += str.substr(curr, endcom - curr);
155  curr = endcom;
156  continue;
157  }
158 
159  if ( leftover ) *leftover += str.substr(curr, begin - curr);
160  if ( begin == end || begin > str.length() - 3 || str[begin + 1] == '/' )
161  return tags;
162 
163  pos_t close = str.find(">", curr);
164  if ( close == end ) return tags;
165 
166  // Find the tag name.
167  curr = str.find_first_of(" \t\n/>", begin);
168  tags.push_back(new XMLTag());
169  tags.back()->name = str.substr(begin + 1, curr - begin - 1);
170 
171  while ( true ) {
172 
173  // Now skip some white space to see if we can find an attribute.
174  curr = str.find_first_not_of(" \t\n", curr);
175  if ( curr == end || curr >= close ) break;
176 
177  pos_t tend = str.find_first_of("= \t\n", curr);
178  if ( tend == end || tend >= close ) break;
179 
180  string name = str.substr(curr, tend - curr);
181  curr = str.find("=", curr) + 1;
182 
183  // OK now find the beginning and end of the atribute.
184  curr = str.find("\"", curr);
185  if ( curr == end || curr >= close ) break;
186  pos_t bega = ++curr;
187  curr = str.find("\"", curr);
188  while ( curr != end && str[curr - 1] == '\\' )
189  curr = str.find("\"", curr + 1);
190 
191  string value = str.substr(bega, curr == end? end: curr - bega);
192 
193  tags.back()->attr[name] = value;
194 
195  ++curr;
196 
197  }
198 
199  curr = close + 1;
200  if ( str[close - 1] == '/' ) continue;
201 
202  pos_t endtag = str.find("</" + tags.back()->name + ">", curr);
203  if ( endtag == end ) {
204  tags.back()->contents = str.substr(curr);
205  curr = endtag;
206  } else {
207  tags.back()->contents = str.substr(curr, endtag - curr);
208  curr = endtag + tags.back()->name.length() + 3;
209  }
210 
211  string leftovers;
212  tags.back()->tags = findXMLTags(tags.back()->contents, &leftovers);
213  if ( leftovers.find_first_not_of(" \t\n") == end ) leftovers="";
214  tags.back()->contents = leftovers;
215 
216  }
217 
218  return tags;
219 
220  }
221 
222  // Print out this tag to a stream.
223  void list(ostream & os) const {
224  os << "<" << name;
225  for ( map<string,string>::const_iterator it = attr.begin();
226  it != attr.end(); ++it )
227  os << " " << it->first << "=\"" << it->second << "\"";
228  if ( contents.empty() && tags.empty() ) {
229  os << "/>" << endl;
230  return;
231  }
232  os << ">" << endl;
233  for ( int i = 0, N = tags.size(); i < N; ++i )
234  tags[i]->list(os);
235 
236  os << "````" << contents << "''''</" << name << ">" << endl;
237  }
238 
239 };
240 
241 //==========================================================================
242 
243 // The LHAweights struct represents the information in a weights tag.
244 
245 struct LHAweights {
246 
247  // Initialize default values.
248  LHAweights() {}
249 
250  // Construct from XML tag
251  LHAweights(const XMLTag & tag);
252 
253  // Print out an XML tag.
254  void list(ostream & file) const;
255 
256  // Function to reset this object.
257  void clear() {
258  contents="";
259  weights.clear();
260  attributes.clear();
261  }
262 
263  // The weights of this event.
264  vector<double> weights;
265 
266  // Any other attributes.
267  map<string,string> attributes;
268 
269  // The contents of the tag.
270  string contents;
271 
272  // Return number of weights.
273  int size() const { return int(weights.size()); }
274 
275 };
276 
277 //==========================================================================
278 
279 // Collect different scales relevant for an event.
280 
281 struct LHAscales {
282 
283  // Empty constructor.
284  LHAscales(double defscale = -1.0)
285  : muf(defscale), mur(defscale), mups(defscale), SCALUP(defscale) {}
286 
287  // Construct from an XML-tag
288  LHAscales(const XMLTag & tag, double defscale = -1.0);
289 
290  // Print out the corresponding XML-tag.
291  void list(ostream & file) const;
292 
293  // Function to reset this object.
294  void clear() {
295  contents="";
296  muf=mur=mups=SCALUP;
297  attributes.clear();
298  }
299 
300  // The factorization scale used for this event.
301  double muf;
302 
303  // The renormalization scale used for this event.
304  double mur;
305 
306  // The starting scale for the parton shower as suggested by the
307  // matrix element generator.
308  double mups;
309 
310  // Any other scales reported by the matrix element generator.
311  map<string,double> attributes;
312 
313  // The default scale in this event.
314  double SCALUP;
315 
316  // The contents of the tag.
317  string contents;
318 
319 };
320 
321 //==========================================================================
322 
323 // Collect generator information for an event file.
324 
325 struct LHAgenerator {
326 
327  // Empty constructor.
328  LHAgenerator()
329  : name(""), version(""), contents("") {}
330 
331  // Construct from an XML-tag
332  LHAgenerator(const XMLTag & tag, string defname = "");
333 
334  // Print out the corresponding XML-tag.
335  void list(ostream & file) const;
336 
337  // Function to reset this object.
338  void clear() {
339  contents="";
340  name="";
341  version="";
342  attributes.clear();
343  }
344 
345  // The generator name used for this file.
346  string name;
347 
348  // The generator version used for this file.
349  string version;
350 
351  // Any other attributes.
352  map<string,string> attributes;
353 
354  // The contents of the tag.
355  string contents;
356 
357 };
358 
359 //==========================================================================
360 
361 // Collect the wgt information.
362 
363 struct LHAwgt {
364 
365  // Empty constructor.
366  LHAwgt(double defwgt = 1.0)
367  : id(""), contents(defwgt) {}
368 
369  // Construct from an XML-tag
370  LHAwgt(const XMLTag & tag, double defwgt = 1.0);
371 
372  // Print out the corresponding XML-tag.
373  void list(ostream & file) const;
374 
375  // Function to reset this object.
376  void clear() {
377  contents=0.0;
378  id="";
379  attributes.clear();
380  }
381 
382  // The identification number of this wgt tag.
383  string id;
384 
385  // Any other attributes.
386  map<string,string> attributes;
387 
388  // The weight associated to this tag.
389  double contents;
390 
391 };
392 
393 //==========================================================================
394 
395 // Collect the wgt information.
396 
397 struct LHAweight {
398 
399  // Empty constructor.
400  LHAweight(string defname = "")
401  : id(defname), contents(defname) {}
402 
403  // Construct from an XML-tag
404  LHAweight(const XMLTag & tag, string defname = "");
405 
406  // Print out the corresponding XML-tag.
407  void list(ostream & file) const;
408 
409  // Function to reset this object.
410  void clear() {
411  contents="";
412  id="";
413  attributes.clear();
414  }
415 
416  // The identification number of this weight tag.
417  string id;
418 
419  // Any other attributes.
420  map<string,string> attributes;
421 
422  // The weight description associated to this tag.
423  string contents;
424 
425 };
426 
427 //==========================================================================
428 
429 // The LHAweightgroup assigns a group-name to a set of LHAweight objects.
430 
431 struct LHAweightgroup {
432 
433  // Default constructor;
434  LHAweightgroup() {}
435 
436  // Construct a group of LHAweight objects from an XML tag and
437  // insert them in the given vector.
438  LHAweightgroup(const XMLTag & tag);
439 
440  // Print out the corresponding XML-tag.
441  void list(ostream & file) const;
442 
443  // Function to reset this object.
444  void clear() {
445  contents="";
446  name="";
447  weights.clear();
448  attributes.clear();
449  }
450 
451  // The contents of the tag.
452  string contents;
453 
454  // The name.
455  string name;
456 
457  // The vector of weights.
458  map<string, LHAweight> weights;
459  vector<string> weightsKeys;
460 
461  // Any other attributes.
462  map<string,string> attributes;
463 
464  // Return number of weights.
465  int size() const { return int(weights.size()); }
466 
467 };
468 
469 //==========================================================================
470 
471 // The LHArwgt assigns a group-name to a set of LHAwgt objects.
472 
473 struct LHArwgt {
474 
475  // Default constructor;
476  LHArwgt() {}
477 
478  // Construct a group of LHAwgt objects from an XML tag and
479  // insert them in the given vector.
480  LHArwgt(const XMLTag & tag);
481 
482  // Print out the corresponding XML-tag.
483  void list(ostream & file) const;
484 
485  // Function to reset this object.
486  void clear() {
487  contents="";
488  wgts.clear();
489  attributes.clear();
490  }
491 
492  // The contents of the tag.
493  string contents;
494 
495  // The map of weights.
496  map<string, LHAwgt> wgts;
497  vector<string> wgtsKeys;
498 
499  // Any other attributes.
500  map<string,string> attributes;
501 
502  // Return number of weights.
503  int size() const { return int(wgts.size()); }
504 
505 };
506 
507 //==========================================================================
508 
509 // The LHAinitrwgt assigns a group-name to a set of LHAweightgroup objects.
510 
511 struct LHAinitrwgt {
512 
513  // Default constructor;
514  LHAinitrwgt() {}
515 
516  // Construct a group of LHAweightgroup objects from an XML tag and
517  // insert them in the given vector.
518  LHAinitrwgt(const XMLTag & tag);
519 
520  // Print out the corresponding XML-tag.
521  void list(ostream & file) const;
522 
523  // Function to reset this object.
524  void clear() {
525  contents="";
526  weights.clear();
527  weightgroups.clear();
528  attributes.clear();
529  }
530 
531  // The contents of the tag.
532  string contents;
533 
534  // The vector of weight's.
535  map<string, LHAweight> weights;
536  vector<string> weightsKeys;
537 
538  // The vector of weightgroup's.
539  map<string, LHAweightgroup> weightgroups;
540  vector<string> weightgroupsKeys;
541 
542  // Any other attributes.
543  map<string,string> attributes;
544 
545  // Return number of weights.
546  int size() const { return int(weights.size());}
547 
548  // Return number of weights.
549  int sizeWeightGroups() const { return int(weightgroups.size()); }
550 
551 };
552 
553 //==========================================================================
554 
555 // The HEPRUP class is a simple container corresponding to the Les Houches
556 // accord (<A HREF="http://arxiv.org/abs/hep-ph/0109068">hep-ph/0109068</A>)
557 // common block with the same name. The members are named in the same
558 // way as in the common block. However, fortran arrays are represented
559 // by vectors, except for the arrays of length two which are
560 // represented by pair objects.
561 
562 class HEPRUP {
563 
564 public:
565 
566  // Default constructor.
567  HEPRUP() : IDWTUP(0), NPRUP(0) {}
568 
569  // Assignment operator.
570  HEPRUP & operator=(const HEPRUP & x) {
571  IDBMUP = x.IDBMUP;
572  EBMUP = x.EBMUP;
573  PDFGUP = x.PDFGUP;
574  PDFSUP = x.PDFSUP;
575  IDWTUP = x.IDWTUP;
576  NPRUP = x.NPRUP;
577  XSECUP = x.XSECUP;
578  XERRUP = x.XERRUP;
579  XMAXUP = x.XMAXUP;
580  LPRUP = x.LPRUP;
581  initrwgt = x.initrwgt;
582  generators = x.generators;
583  weightgroups = x.weightgroups;
584  weights = x.weights;
585  return *this;
586  }
587 
588  // Destructor.
589  ~HEPRUP() {}
590 
591  // Set the NPRUP variable, corresponding to the number of
592  // sub-processes, to \a nrup, and resize all relevant vectors
593  // accordingly.
594  void resize(int nrup) {
595  NPRUP = nrup;
596  resize();
597  }
598 
599  // Assuming the NPRUP variable, corresponding to the number of
600  // sub-processes, is correctly set, resize the relevant vectors
601  // accordingly.
602  void resize() {
603  XSECUP.resize(NPRUP);
604  XERRUP.resize(NPRUP);
605  XMAXUP.resize(NPRUP);
606  LPRUP.resize(NPRUP);
607  }
608 
609  // Clear all members.
610  void clear();
611 
612  // PDG id's of beam particles. (first/second is in +/-z direction).
613  pair<long,long> IDBMUP;
614 
615  // Energy of beam particles given in GeV.
616  pair<double,double> EBMUP;
617 
618  // The author group for the PDF used for the beams according to the
619  // PDFLib specification.
620  pair<int,int> PDFGUP;
621 
622  // The id number the PDF used for the beams according to the
623  // PDFLib specification.
624  pair<int,int> PDFSUP;
625 
626  // Master switch indicating how the ME generator envisages the
627  // events weights should be interpreted according to the Les Houches
628  // accord.
629  int IDWTUP;
630 
631  // The number of different subprocesses in this file.
632  int NPRUP;
633 
634  // The cross sections for the different subprocesses in pb.
635  vector<double> XSECUP;
636 
637  // The statistical error in the cross sections for the different
638  // subprocesses in pb.
639  vector<double> XERRUP;
640 
641  // The maximum event weights (in HEPEUP::XWGTUP) for different
642  // subprocesses.
643  vector<double> XMAXUP;
644 
645  // The subprocess code for the different subprocesses.
646  vector<int> LPRUP;
647 
648  // Contents of the LHAinitrwgt tag
649  LHAinitrwgt initrwgt;
650 
651  // Contents of the LHAgenerator tags.
652  vector<LHAgenerator> generators;
653 
654  // A map of the LHAweightgroup tags, indexed by name.
655  map<string,LHAweightgroup> weightgroups;
656 
657  // A map of the LHAweight tags, indexed by name.
658  map<string,LHAweight> weights;
659 
660 };
661 
662 //==========================================================================
663 
664 // The HEPEUP class is a simple container corresponding to the Les Houches
665 // accord (<A HREF="http://arxiv.org/abs/hep-ph/0109068">hep-ph/0109068</A>)
666 // common block with the same name. The members are named in the same
667 // way as in the common block. However, fortran arrays are represented
668 // by vectors, except for the arrays of length two which are
669 // represented by pair objects.
670 
671 class HEPEUP {
672 
673 public:
674 
675  // Default constructor.
676  HEPEUP()
677  : NUP(0), IDPRUP(0), XWGTUP(0.0), XPDWUP(0.0, 0.0),
678  SCALUP(0.0), AQEDUP(0.0), AQCDUP(0.0), heprup(0) {}
679 
680  // Copy constructor
681  HEPEUP(const HEPEUP & x) { operator=(x); }
682 
683  // Copy information from the given HEPEUP.
684  HEPEUP & setEvent(const HEPEUP & x);
685 
686  // Assignment operator.
687  HEPEUP & operator=(const HEPEUP & x) {
688  clear();
689  setEvent(x);
690  return *this;
691  }
692 
693  // Destructor.
694  ~HEPEUP() {
695  clear();
696  };
697 
698  // Reset the HEPEUP object.
699  void reset();
700 
701  // Clear the HEPEUP object.
702  void clear() {
703  reset();
704  }
705 
706  // Set the NUP variable, corresponding to the number of particles in
707  // the current event, to \a nup, and resize all relevant vectors
708  // accordingly.
709  void resize(int nup) {
710  NUP = nup;
711  resize();
712  }
713 
714  // Return the main weight for this event.
715  double weight() const {
716  return XWGTUP;
717  }
718 
719  // Assuming the NUP variable, corresponding to the number of
720  // particles in the current event, is correctly set, resize the
721  // relevant vectors accordingly.
722  void resize();
723 
724  // The number of particle entries in the current event.
725  int NUP;
726 
727  // The subprocess code for this event (as given in LPRUP).
728  int IDPRUP;
729 
730  // The weight for this event.
731  double XWGTUP;
732 
733  // The PDF weights for the two incoming partons. Note that this
734  // variable is not present in the current LesHouches accord
735  // (<A HREF="http://arxiv.org/abs/hep-ph/0109068">hep-ph/0109068</A>),
736  // hopefully it will be present in a future accord.
737  pair<double,double> XPDWUP;
738 
739  // The scale in GeV used in the calculation of the PDF's in this
740  // event.
741  double SCALUP;
742 
743  // The value of the QED coupling used in this event.
744  double AQEDUP;
745 
746  // The value of the QCD coupling used in this event.
747  double AQCDUP;
748 
749  // The PDG id's for the particle entries in this event.
750  vector<long> IDUP;
751 
752  // The status codes for the particle entries in this event.
753  vector<int> ISTUP;
754 
755  // Indices for the first and last mother for the particle entries in
756  // this event.
757  vector< pair<int,int> > MOTHUP;
758 
759  // The colour-line indices (first(second) is (anti)colour) for the
760  // particle entries in this event.
761  vector< pair<int,int> > ICOLUP;
762 
763  // Lab frame momentum (Px, Py, Pz, E and M in GeV) for the particle
764  // entries in this event.
765  vector< vector<double> > PUP;
766 
767  // Invariant lifetime (c*tau, distance from production to decay in
768  // mm) for the particle entries in this event.
769  vector<double> VTIMUP;
770 
771  // Spin info for the particle entries in this event given as the
772  // cosine of the angle between the spin vector of a particle and the
773  // 3-momentum of the decaying particle, specified in the lab frame.
774  vector<double> SPINUP;
775 
776  // A pointer to the current HEPRUP object.
777  HEPRUP * heprup;
778 
779  // The weights associated with this event, as given by the LHAwgt tags.
780  map<string,double> weights_detailed;
781 
782  // The weights associated with this event, as given by the LHAweights tags.
783  vector<double> weights_compressed;
784 
785  // Contents of the LHAscales tag
786  LHAscales scalesSave;
787 
788  // Contents of the LHAweights tag (compressed format)
789  LHAweights weightsSave;
790 
791  // Contents of the LHArwgt tag (detailed format)
792  LHArwgt rwgtSave;
793 
794  // Any other attributes.
795  map<string,string> attributes;
796 
797 };
798 
799 //==========================================================================
800 
801 // The Reader class is initialized with a stream from which to read a
802 // version 1/3 Les Houches Accord event file. In the constructor of
803 // the Reader object the optional header information is read and then
804 // the mandatory init is read. After this the whole header block
805 // including the enclosing lines with tags are available in the public
806 // headerBlock member variable. Also the information from the init
807 // block is available in the heprup member variable and any additional
808 // comment lines are available in initComments. After each successful
809 // call to the readEvent() function the standard Les Houches Accord
810 // information about the event is available in the hepeup member
811 // variable and any additional comments in the eventComments
812 // variable. A typical reading sequence would look as follows:
813 
814 class Reader {
815 
816 public:
817 
818  // Initialize the Reader with a filename from which to read an event
819  // file. After the constructor is called the whole header block
820  // including the enclosing lines with tags are available in the
821  // public headerBlock member variable. Also the information from the
822  // init block is available in the heprup member variable and any
823  // additional comment lines are available in initComments.
824  //
825  // filename: the name of the file to read from.
826  //
827  Reader(string filenameIn)
828  : filename(filenameIn), intstream(NULL), file(NULL), version() {
829  intstream = new igzstream(filename.c_str());
830  file = intstream;
831  isGood = init();
832  }
833 
834  Reader(istream* is)
835  : filename(""), intstream(NULL), file(is), version() {
836  isGood = init();
837  }
838 
839  // Clean up
840  ~Reader() {
841  if (intstream) delete intstream;
842  }
843 
844 
845 
846  // (Re)initialize the Reader with a filename from which to read an event
847  // file. After this, all information from the header and init block is
848  // available.
849  //
850  // filename: File name (not used as the input file stream is given)
851  // isIn : Name of the input file stream.
852  bool setup(string filenameIn) {
853  filename = filenameIn;
854  if (intstream) delete intstream;
855  intstream = new igzstream(filename.c_str());
856  file = intstream;
857  isGood = init();
858  return isGood;
859  }
860 
861 private:
862 
863  // Used internally in the constructors to read header and init blocks.
864  bool init();
865 
866 public:
867 
868  // Read an event from the file and store it in the hepeup
869  // object. Optional comment lines are stored in the eventComments
870  // member variable.
871  bool readEvent(HEPEUP * peup = 0);
872 
873  // Reset values of all event-related members to their defaults.
874  void clearEvent() {
875  currentLine = "";
876  hepeup.clear();
877  eventComments = "";
878  weights_detailed_vec.resize(0);
879  weightnames_detailed_vec.resize(0);
880  }
881 
882 protected:
883 
884  // Used internally to read a single line from the stream.
885  bool getLine() {
886  currentLine = "";
887  if(!getline(*file, currentLine)) return false;
888  // Replace single by double quotes
889  replace(currentLine.begin(),currentLine.end(),'\'','\"');
890  return true;
891  }
892 
893 protected:
894 
895  // Name of file-to-be-read.
896  string filename;
897 
898  // A local stream which is unused if a stream is supplied from the
899  // outside.
900  igzstream* intstream;
901 
902  // The stream we are reading from. This may be a pointer to an
903  // external stream or the internal intstream.
904  istream * file;
905 
906  // The last line read in from the stream in getline().
907  string currentLine;
908 
909 public:
910 
911  // Save if the initialisation worked.
912  bool isGood;
913 
914  // XML file version
915  int version;
916 
917  // All lines (since the last readEvent()) outside the header, init
918  // and event tags.
919  string outsideBlock;
920 
921  // All lines from the header block.
922  string headerBlock;
923  string headerComments;
924 
925  // The standard init information.
926  HEPRUP heprup;
927 
928  // Additional comments found in the init block.
929  string initComments;
930 
931  // The standard information about the last read event.
932  HEPEUP hepeup;
933 
934  // Additional comments found with the last read event.
935  string eventComments;
936 
937  // The detailed weights associated with this event, linearized to a vector.
938  vector<double> weights_detailed_vec;
939  vector<double> weights_detailed_vector()
940  { return weights_detailed_vec; }
941  vector<string> weightnames_detailed_vec;
942  vector<string> weightnames_detailed_vector()
943  { return weightnames_detailed_vec; }
944 
945 private:
946 
947  // The default constructor should never be used.
948  Reader();
949 
950  // The copy constructor should never be used.
951  Reader(const Reader &);
952 
953  // The Reader cannot be assigned to.
954  Reader & operator=(const Reader &);
955 
956 };
957 
958 //==========================================================================
959 
960 // The Writer class is initialized with a stream to which to write a
961 // version 1.0 or 3.0 Les Houches Accord event file. In the init() function of
962 // the Writer object the main XML tag, header and init blocks are written,
963 // with the corresponding end tag is written by list_end_tag().
964 // After a Writer object (in the following called "writer") has been created,
965 // it is possible to assign version (3 by default) information by
966 //
967 // writer.version = <value>;
968 //
969 // The header block (called "someHeaderString" below) is assigned by
970 //
971 // writer.headerBlock() << someHeaderString;
972 //
973 // and the init block comments (called "someInitString" below) are assigned via
974 //
975 // writer.initComments() << someInitString;
976 //
977 // The standard init information (including amendments for LHEF 3.0) can
978 // be assigned by the heprup member variable:
979 //
980 // writer.heprup = heprup;
981 //
982 // where heprup is an object of type HEPRUP. All of the above information
983 // will be writen by calling the init() function.
984 //
985 // Before each event is written out with the writeEvent() function,
986 // the standard event information can be assigned to the hepeup
987 // variable by
988 //
989 // writer.hepeup = hepeup;
990 //
991 // where hepeup is of type HEPEUP. Event comments (called
992 // "someCommentString" below) can be assigned through
993 //
994 // writer.eventComments() << someCommentString;
995 //
996 // All of this event information is written by the writeEvent() function.
997 
998 class Writer {
999 
1000 public:
1001 
1002  // Create a Writer object giving a stream to write to.
1003  // @param os the stream where the event file is written.
1004  Writer(ostream & os)
1005  : file(os), version(3) {}
1006 
1007  // Create a Writer object giving a filename to write to.
1008  // @param filename the name of the event file to be written.
1009  Writer(string filename)
1010  : intstream(filename.c_str()), file(intstream), version(3) {}
1011 
1012  // The destructor.
1013  ~Writer() {}
1014 
1015  // Add header lines consisting of XML code with this stream.
1016  ostream & headerBlock() {
1017  return headerStream;
1018  }
1019 
1020  // Add comment lines to the init block with this stream.
1021  ostream & initComments() {
1022  return initStream;
1023  }
1024 
1025  // Add comment lines to the next event to be written out with this stream.
1026  ostream & eventComments() {
1027  return eventStream;
1028  }
1029 
1030  // Write out the final XML end-tag.
1031  void list_end_tag() {
1032  file << "</LesHouchesEvents>" << endl;
1033  }
1034 
1035  // Write out an optional header block followed by the standard init
1036  // block information together with any comment lines.
1037  void init();
1038 
1039  // Write out the event stored in hepeup, followed by optional
1040  // comment lines.
1041  bool writeEvent(HEPEUP * peup = 0, int pDigits = 15);
1042  // Write out an event as a string.
1043  string getEventString(HEPEUP * peup = 0);
1044 
1045 protected:
1046 
1047  // Make sure that each line in the string \a s starts with a
1048  // #-character and that the string ends with a new-line.
1049  string hashline(string s, bool comment = false);
1050 
1051 protected:
1052 
1053  // A local stream which is unused if a stream is supplied from the
1054  // outside.
1055  ofstream intstream;
1056 
1057  // The stream we are writing to. This may be a reference to an
1058  // external stream or the internal intstream.
1059  ostream & file;
1060 
1061 public:
1062 
1063  // Stream to add all lines in the header block.
1064  ostringstream headerStream;
1065 
1066  // The standard init information.
1067  HEPRUP heprup;
1068 
1069  // Stream to add additional comments to be put in the init block.
1070  ostringstream initStream;
1071 
1072  // The standard information about the event we will write next.
1073  HEPEUP hepeup;
1074 
1075  // Stream to add additional comments to be written together the next event.
1076  ostringstream eventStream;
1077 
1078  // XML file version
1079  int version;
1080 
1081 private:
1082 
1083  // The default constructor should never be used.
1084  Writer();
1085 
1086  // The copy constructor should never be used.
1087  Writer(const Writer &);
1088 
1089  // The Writer cannot be assigned to.
1090  Writer & operator=(const Writer &);
1091 
1092 };
1093 
1094 //==========================================================================
1095 
1096 } // end namespace Pythia8
1097 
1098 #endif // end Pythia8_LHEF3_H