StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
LHEF3.cc
1 // LHEF3.cc 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.
7 // It contains the main class for LHEF 3.0 functionalities.
8 // Function definitions.
9 
10 #include "Pythia8/LHEF3.h"
11 
12 namespace Pythia8 {
13 
14 //==========================================================================
15 
16 // The XMLTag struct is used to represent all information within an XML tag.
17 // It contains the attributes as a map, any sub-tags as a vector of pointers
18 // to other XMLTag objects, and any other information as a single string.
19 
20 //--------------------------------------------------------------------------
21 
22 // Constants.
23 const XMLTag::pos_t XMLTag::end = string::npos;
24 
25 //==========================================================================
26 
27 // The LHAweights struct.
28 
29 //--------------------------------------------------------------------------
30 
31 // Construct from XML tag.
32 
33 LHAweights::LHAweights(const XMLTag & tag) {
34  for ( map<string,string>::const_iterator it = tag.attr.begin();
35  it != tag.attr.end(); ++it ) {
36  string v = it->second.c_str();
37  attributes[it->first] = v;
38  }
39 
40  contents = tag.contents;
41 
42  istringstream iss(tag.contents);
43  double w;
44  while ( iss >> w ) weights.push_back(w);
45 }
46 
47 //--------------------------------------------------------------------------
48 
49 // Print out.
50 
51 void LHAweights::list(ostream & file) const {
52  file << "<weights";
53  for ( map<string,string>::const_iterator it = attributes.begin();
54  it != attributes.end(); ++it )
55  file << " " << it->first << "=\"" << it->second << "\"";
56  file << ">";
57  for ( int j = 0, M = weights.size(); j < M; ++j ) file << " " << weights[j];
58  file << "</weights>" << endl;
59 }
60 
61 //==========================================================================
62 
63 // The LHAscales struct: Collect different scales relevant for an event.
64 
65 //--------------------------------------------------------------------------
66 
67 // Construct from an XML-tag.
68 
69 LHAscales::LHAscales(const XMLTag & tag, double defscale)
70  : muf(defscale), mur(defscale), mups(defscale), SCALUP(defscale) {
71  for ( map<string,string>::const_iterator it = tag.attr.begin();
72  it != tag.attr.end(); ++it ) {
73  double v = atof(it->second.c_str());
74  if ( it->first == "muf" ) muf = v;
75  else if ( it->first == "mur" ) mur = v;
76  else if ( it->first == "mups" ) mups = v;
77  else attributes.insert(make_pair(it->first, v));
78  }
79  contents = tag.contents;
80 }
81 
82 //--------------------------------------------------------------------------
83 
84 // Print out the corresponding XML-tag.
85 
86 void LHAscales::list(ostream & file) const {
87  file << "<scales";
88  file << " muf=\"" << muf << "\"";
89  file << " mur=\"" << mur << "\"";
90  file << " mups=\"" << mups << "\"";
91  for ( map<string,double>::const_iterator it = attributes.begin();
92  it != attributes.end(); ++it )
93  file << " " << it->first << "=\"" << it->second << "\"";
94  file << ">" << contents;
95  file << "</scales>" << endl;
96 }
97 
98 //==========================================================================
99 
100 // The LHAgenerator struct: Collect generator information for an event file.
101 
102 //--------------------------------------------------------------------------
103 
104 // Construct from an XML-tag
105 
106 LHAgenerator::LHAgenerator(const XMLTag & tag, string defname)
107  : name(defname), version(defname), contents(defname) {
108  for ( map<string,string>::const_iterator it = tag.attr.begin();
109  it != tag.attr.end(); ++it ) {
110  if ( it->first == "name" ) name = it->second;
111  else if ( it->first == "version" ) version = it->second;
112  else attributes.insert(make_pair(it->first, it->second));
113  }
114  contents = tag.contents;
115 }
116 
117 //--------------------------------------------------------------------------
118 
119 // Print out the corresponding XML-tag.
120 
121 void LHAgenerator::list(ostream & file) const {
122  file << "<generator";
123  if ( name != "" ) file << " name=\"" << name << "\"";
124  if ( version != "" ) file << " version=\"" << version << "\"";
125  for ( map<string,string>::const_iterator it = attributes.begin();
126  it != attributes.end(); ++it )
127  file << " " << it->first << "=\"" << it->second << "\"";
128  file << " >";
129  file << contents;
130  file << "</generator>" << endl;
131 }
132 
133 //==========================================================================
134 
135 // The LHAwgt struct: Collect the wgt information.
136 
137 //--------------------------------------------------------------------------
138 
139 // Construct from an XML-tag
140 
141 LHAwgt::LHAwgt(const XMLTag & tag, double defwgt)
142  : id(""), contents(defwgt) {
143  for ( map<string,string>::const_iterator it = tag.attr.begin();
144  it != tag.attr.end(); ++it ) {
145  if ( it->first == "id" ) id = it->second;
146  else attributes.insert(make_pair(it->first, it->second));
147  }
148  contents = atof(tag.contents.c_str());
149 }
150 
151 //--------------------------------------------------------------------------
152 
153 // Print out the corresponding XML-tag.
154 
155 void LHAwgt::list(ostream & file) const {
156  file << "<wgt";
157  if ( id != "" ) file << " id=\"" << id << "\"";
158  for ( map<string,string>::const_iterator it = attributes.begin();
159  it != attributes.end(); ++it )
160  file << " " << it->first << "=\"" << it->second << "\"";
161  file << " >";
162  file << contents;
163  file << "</wgt>" << endl;
164 }
165 
166 //==========================================================================
167 
168 // The LHAweight struct: Collect the weight information.
169 
170 //--------------------------------------------------------------------------
171 
172 // Construct from an XML-tag.
173 
174 LHAweight::LHAweight(const XMLTag & tag, string defname)
175  : id(defname), contents(defname) {
176  for ( map<string,string>::const_iterator it = tag.attr.begin();
177  it != tag.attr.end(); ++it ) {
178  if ( it->first == "id" ) id = it->second;
179  else attributes.insert(make_pair(it->first, it->second));
180  }
181  contents = tag.contents;
182 }
183 
184 //--------------------------------------------------------------------------
185 
186 // Print out the corresponding XML-tag.
187 
188 void LHAweight::list(ostream & file) const {
189  file << "<weight";
190  if ( id != "" ) file << " id=\"" << id << "\"";
191  for ( map<string,string>::const_iterator it = attributes.begin();
192  it != attributes.end(); ++it )
193  file << " " << it->first << "=\"" << it->second << "\"";
194  file << " >";
195  file << contents;
196  file << "</weight>" << endl;
197 }
198 
199 //==========================================================================
200 
201 // The LHAweightgroup struct: The LHAweightgroup assigns a group-name to a set
202 // of LHAweight objects.
203 
204 //--------------------------------------------------------------------------
205 
206 // Construct a group of LHAweight objects from an XML tag and
207 // insert them in the given vector.
208 
209 LHAweightgroup::LHAweightgroup(const XMLTag & tag) {
210 
211  for ( map<string,string>::const_iterator it = tag.attr.begin();
212  it != tag.attr.end(); ++it ) {
213  if ( it->first == "name" ) name = it->second;
214  else attributes.insert(make_pair(it->first,it->second));
215  }
216  if ( name=="" ) {
217  string key("type");
218  if( attributes.find(key) != attributes.end() ) {
219  name = attributes[key];
220  }
221  }
222 
223  contents = tag.contents;
224 
225  // Now add the weight's step by step.
226  string s;
227  vector<XMLTag*> tags = XMLTag::findXMLTags(tag.contents, &s);
228  for ( int i = 0, N = tags.size(); i < N; ++i ) {
229  const XMLTag & tagnow = *tags[i];
230  LHAweight wt(tagnow);
231  weights.insert(make_pair(wt.id, wt));
232  weightsKeys.push_back(wt.id);
233  }
234  for ( int i = 0, N = tag.tags.size(); i < N; ++i ) {
235  const XMLTag & tagnow = *tag.tags[i];
236  const LHAweight & wt(tagnow);
237  weights.insert(make_pair(wt.id, wt));
238  weightsKeys.push_back(wt.id);
239  }
240 
241  for ( int i = 0, N = tags.size(); i < N; ++i ) if (tags[i]) delete tags[i];
242 
243 }
244 
245 //--------------------------------------------------------------------------
246 
247 // Print out the corresponding XML-tag.
248 
249 void LHAweightgroup::list(ostream & file) const {
250  file << "<weightgroup";
251  if ( name != "" ) file << " name=\"" << name << "\"";
252  for ( map<string,string>::const_iterator it = attributes.begin();
253  it != attributes.end(); ++it )
254  file << " " << it->first << "=\"" << it->second << "\"";
255  file << " >\n";
256  for ( map<string,LHAweight>::const_iterator it = weights.begin();
257  it != weights.end(); ++it ) it->second.list(file);
258  file << "</weightgroup>" << endl;
259 }
260 
261 //==========================================================================
262 
263 // The LHArwgt struct: Assigns a group-name to a set of LHAwgt objects.
264 
265 //--------------------------------------------------------------------------
266 
267 // Construct a group of LHAwgt objects from an XML tag and
268 // insert them in the given vector.
269 
270 LHArwgt::LHArwgt(const XMLTag & tag) {
271 
272  for ( map<string,string>::const_iterator it = tag.attr.begin();
273  it != tag.attr.end(); ++it ) {
274  string v = it->second.c_str();
275  attributes[it->first] = v;
276  }
277  contents = tag.contents;
278 
279  // Now add the wgt's step by step.
280  string s;
281  vector<XMLTag*> tags = XMLTag::findXMLTags(tag.contents, &s);
282  for ( int i = 0, N = tags.size(); i < N; ++i ) {
283  const XMLTag & tagnow = *tags[i];
284  LHAwgt wt(tagnow);
285  wgts.insert(make_pair(wt.id, wt));
286  wgtsKeys.push_back(wt.id);
287  }
288  for ( int i = 0, N = tag.tags.size(); i < N; ++i ) {
289  const XMLTag & tagnow = *tag.tags[i];
290  LHAwgt wt(tagnow);
291  wgts.insert(make_pair(wt.id, wt));
292  wgtsKeys.push_back(wt.id);
293  }
294 
295  for ( int i = 0, N = tags.size(); i < N; ++i ) if (tags[i]) delete tags[i];
296 
297 }
298 
299 //--------------------------------------------------------------------------
300 
301 // Print out the corresponding XML-tag.
302 
303 void LHArwgt::list(ostream & file) const {
304  file << "<rwgt";
305  for ( map<string,string>::const_iterator it = attributes.begin();
306  it != attributes.end(); ++it )
307  file << " " << it->first << "=\"" << it->second << "\"";
308  file << " >\n";
309  for ( map<string,LHAwgt>::const_iterator it = wgts.begin();
310  it != wgts.end(); ++it ) it->second.list(file);
311  file << "</rwgt>" << endl;
312 }
313 
314 //==========================================================================
315 
316 // The LHAinitrwgt assigns a group-name to a set of LHAweightgroup objects.
317 
318 //--------------------------------------------------------------------------
319 
320 // Construct a group of LHAweightgroup objects from an XML tag and
321 // insert them in the given vector.
322 
323 LHAinitrwgt::LHAinitrwgt(const XMLTag & tag) {
324  for ( map<string,string>::const_iterator it = tag.attr.begin();
325  it != tag.attr.end(); ++it ) {
326  string v = it->second.c_str();
327  attributes[it->first] = v;
328  }
329  contents = tag.contents;
330 
331  // Now add the wgt's step by step.
332  string s;
333  vector<XMLTag*> tags = XMLTag::findXMLTags(tag.contents, &s);
334  for ( int i = 0, N = tags.size(); i < N; ++i ) {
335  const XMLTag & tagnow = *tags[i];
336  if ( tagnow.name == "weightgroup" ) {
337  LHAweightgroup wgroup(tagnow);
338  string wgname = wgroup.name;
339  // if still no name, use integer as a key
340  if (wgname=="") {
341  stringstream iss;
342  iss << i;
343  wgname=iss.str();
344  }
345  weightgroups.insert(make_pair(wgname, wgroup));
346  weightgroupsKeys.push_back(wgname);
347  string ss;
348  vector<XMLTag*> tags2 = XMLTag::findXMLTags(tagnow.contents, &ss);
349  for ( int k = 0, M = tags2.size(); k < M; ++k ) {
350  const XMLTag & tagnow2 = *tags2[k];
351  if ( tagnow2.name == "weight" ) {
352  LHAweight wt(tagnow2);
353  string wtname = wt.id;
354  weights.insert(make_pair(wtname, wt));
355  weightsKeys.push_back(wtname);
356  }
357  }
358  for ( int j = 0, M = tags2.size(); j < M; ++j )
359  if (tags2[j]) delete tags2[j];
360  } else if ( tagnow.name == "weight" ) {
361  LHAweight wt(tagnow);
362  string wtname = wt.id;
363  weights.insert(make_pair(wtname, wt));
364  weightsKeys.push_back(wtname);
365  }
366  }
367 
368  // Now add the wgt's step by step.
369  for ( int i = 0, N = tag.tags.size(); i < N; ++i ) {
370  const XMLTag & tagnow = *tag.tags[i];
371  if ( tagnow.name == "weightgroup" ) {
372  LHAweightgroup wgroup(tagnow);
373  string wgname = wgroup.name;
374  weightgroups.insert(make_pair(wgname, wgroup));
375  weightgroupsKeys.push_back(wgname);
376  string ss;
377  vector<XMLTag*> tags2 = XMLTag::findXMLTags(tagnow.contents, &ss);
378  for ( int k = 0, M = tags2.size(); k < M; ++k ) {
379  const XMLTag & tagnow2 = *tags2[k];
380  if ( tagnow2.name == "weight" ) {
381  LHAweight wt(tagnow2);
382  string wtname = wt.id;
383  weights.insert(make_pair(wtname, wt));
384  weightsKeys.push_back(wtname);
385  }
386  }
387  for ( int k = 0, M = tagnow.tags.size(); k < M; ++k ) {
388  const XMLTag & tagnow2 = *tagnow.tags[k];
389  if ( tagnow2.name == "weight" ) {
390  LHAweight wt(tagnow2);
391  string wtname = wt.id;
392  weights.insert(make_pair(wtname, wt));
393  weightsKeys.push_back(wtname);
394  }
395  }
396  for ( int j = 0, M = tags2.size(); j < M; ++j )
397  if (tags2[j]) delete tags2[j];
398  } else if ( tagnow.name == "weight" ) {
399  LHAweight wt(tagnow);
400  string wtname = wt.id;
401  weights.insert(make_pair(wtname, wt));
402  weightsKeys.push_back(wtname);
403  }
404  }
405 
406  for ( int i = 0, N = tags.size(); i < N; ++i ) if (tags[i]) delete tags[i];
407 
408 }
409 
410 //--------------------------------------------------------------------------
411 
412 // Print out the corresponding XML-tag.
413 
414 void LHAinitrwgt::list(ostream & file) const {
415  file << "<initrwgt";
416  for ( map<string,string>::const_iterator it = attributes.begin();
417  it != attributes.end(); ++it )
418  file << " " << it->first << "=\"" << it->second << "\"";
419  file << " >\n";
420  for ( map<string,LHAweightgroup>::const_iterator it = weightgroups.begin();
421  it != weightgroups.end(); ++it ) it->second.list(file);
422  for ( map<string,LHAweight>::const_iterator it = weights.begin();
423  it != weights.end(); ++it ) it->second.list(file);
424  file << "</initrwgt>" << endl;
425 }
426 
427 //==========================================================================
428 
429 // The HEPRUP class is a simple container for the Les Houches file init block.
430 
431 void HEPRUP::clear() {
432  IDBMUP = make_pair(0,0);
433  EBMUP = make_pair(0,0);
434  PDFGUP = make_pair(0,0);
435  PDFSUP = make_pair(0,0);
436  IDWTUP = -1;
437  NPRUP = 0;
438  XSECUP.resize(0);
439  XERRUP.resize(0);
440  XMAXUP.resize(0);
441  LPRUP.resize(0);
442  initrwgt.clear();
443  generators.resize(0);
444  weightgroups.clear();
445  weights.clear();
446 
447 }
448 
449 //==========================================================================
450 
451 // The HEPEUP class is a simple container corresponding to the Les Houches
452 // accord (<A HREF="http://arxiv.org/abs/hep-ph/0109068">hep-ph/0109068</A>)
453 // common block with the same name. The members are named in the same
454 // way as in the common block. However, fortran arrays are represented
455 // by vectors, except for the arrays of length two which are
456 // represented by pair objects.
457 
458 //--------------------------------------------------------------------------
459 
460 // Copy information from the given HEPEUP.
461 
462 HEPEUP & HEPEUP::setEvent(const HEPEUP & x) {
463 
464  NUP = x.NUP;
465  IDPRUP = x.IDPRUP;
466  XWGTUP = x.XWGTUP;
467  XPDWUP = x.XPDWUP;
468  SCALUP = x.SCALUP;
469  AQEDUP = x.AQEDUP;
470  AQCDUP = x.AQCDUP;
471  IDUP = x.IDUP;
472  ISTUP = x.ISTUP;
473  MOTHUP = x.MOTHUP;
474  ICOLUP = x.ICOLUP;
475  PUP = x.PUP;
476  VTIMUP = x.VTIMUP;
477  SPINUP = x.SPINUP;
478  heprup = x.heprup;
479  scalesSave = x.scalesSave;
480  weightsSave = x.weightsSave;
481  weights_detailed = x.weights_detailed;
482  weights_compressed = x.weights_compressed;
483  rwgtSave = x.rwgtSave;
484  attributes = x.attributes;
485  return *this;
486 
487 }
488 
489 //--------------------------------------------------------------------------
490 
491 // Reset the HEPEUP object.
492 
493 void HEPEUP::reset() {
494  NUP = 0;
495  weights_detailed.clear();
496  weights_compressed.clear();
497  weightsSave.clear();
498  rwgtSave.clear();
499  scalesSave.clear();
500  attributes.clear();
501 }
502 
503 //--------------------------------------------------------------------------
504 
505 // Assuming the NUP variable, corresponding to the number of
506 // particles in the current event, is correctly set, resize the
507 // relevant vectors accordingly.
508 
509 void HEPEUP::resize() {
510  IDUP.resize(NUP);
511  ISTUP.resize(NUP);
512  MOTHUP.resize(NUP);
513  ICOLUP.resize(NUP);
514  PUP.resize(NUP, vector<double>(5));
515  VTIMUP.resize(NUP);
516  SPINUP.resize(NUP);
517 }
518 
519 //==========================================================================
520 
521 // The Reader class is initialized with a stream from which to read a
522 // version 1/2 Les Houches Accord event file. In the constructor of
523 // the Reader object the optional header information is read and then
524 // the mandatory init is read. After this the whole header block
525 // including the enclosing lines with tags are available in the public
526 // headerBlock member variable. Also the information from the init
527 // block is available in the heprup member variable and any additional
528 // comment lines are available in initComments. After each successful
529 // call to the readEvent() function the standard Les Houches Accord
530 // information about the event is available in the hepeup member
531 // variable and any additional comments in the eventComments
532 // variable. A typical reading sequence would look as follows:
533 
534 //--------------------------------------------------------------------------
535 
536 // Used internally in the constructors to read header and init blocks.
537 bool Reader::init() {
538 
539  bool readingHeader = false;
540  bool readingInit = false;
541 
542  // Make sure we are reading a LHEF file:
543  getLine();
544 
545  if ( currentLine.find("<LesHouchesEvents" ) == string::npos )
546  return false;
547  version = 0;
548  if ( currentLine.find("version=\"1" ) != string::npos )
549  version = 1;
550  else if ( currentLine.find("version=\"2" ) != string::npos )
551  version = 2;
552  else if ( currentLine.find("version=\"3" ) != string::npos )
553  version = 3;
554  else
555  return false;
556 
557  // Clear all members.
558  outsideBlock="";
559  headerBlock="";
560  headerComments="";
561  heprup.clear();
562  initComments="";
563  hepeup.clear();
564  eventComments="";
565 
566  // Loop over all lines until we hit the </init> tag.
567  while ( getLine() && currentLine.find("</init>") == string::npos ) {
568  if ( currentLine.find("<header") != string::npos
569  && currentLine.find("#") == string::npos) {
570  // We have hit the header block, so we should dump this and
571  // all following lines to headerBlock until we hit the end of
572  // it.
573  readingHeader = true;
574  headerBlock = currentLine + "\n";
575  }
576  else if ( ( currentLine.find("<init>") != string::npos
577  || currentLine.find("<init ") != string::npos )
578  && currentLine.find("#") == string::npos) {
579  // We have hit the init block, so we should expect to find the
580  // standard information in the following.
581  readingInit = true;
582 
583  // The first line tells us how many lines to read next.
584  getLine();
585  istringstream iss(currentLine);
586  if ( !( iss >> heprup.IDBMUP.first >> heprup.IDBMUP.second
587  >> heprup.EBMUP.first >> heprup.EBMUP.second
588  >> heprup.PDFGUP.first >> heprup.PDFGUP.second
589  >> heprup.PDFSUP.first >> heprup.PDFSUP.second
590  >> heprup.IDWTUP >> heprup.NPRUP ) ) {
591  heprup.NPRUP = -42;
592  return false;
593  }
594  heprup.resize();
595 
596  for ( int i = 0; i < heprup.NPRUP; ++i ) {
597  getLine();
598  istringstream isss(currentLine);
599  if ( !( isss >> heprup.XSECUP[i] >> heprup.XERRUP[i]
600  >> heprup.XMAXUP[i] >> heprup.LPRUP[i] ) ) {
601  heprup.NPRUP = -42;
602  return false;
603  }
604  }
605  }
606  else if ( currentLine.find("</header>") != string::npos
607  && currentLine.find("#") == string::npos) {
608  // The end of the header block. Dump this line as well to the
609  // headerBlock and we're done.
610  readingHeader = false;
611  headerBlock += currentLine + "\n";
612  }
613  else if ( readingHeader ) {
614  // We are in the process of reading the header block. Dump the
615  // line to headerBlock.
616  headerBlock += currentLine + "\n";
617  headerComments += currentLine + "\n";
618  }
619  else if ( readingInit ) {
620  // Here we found a comment line. Dump it to initComments.
621  initComments += currentLine + "\n";
622  }
623  else {
624  // We found some other stuff outside the standard tags.
625  outsideBlock += currentLine + "\n";
626  }
627  }
628 
629  if ( file == NULL ) heprup.NPRUP = -42;
630 
631  // Scan the header block for XML tags
632  string leftovers;
633  vector<XMLTag*> tags1 = XMLTag::findXMLTags(headerComments, &leftovers);
634  if ( leftovers.find_first_not_of(" \t\n") == string::npos )
635  leftovers="";
636 
637  for ( int i = 0, N = tags1.size(); i < N; ++i ) {
638  const XMLTag & tag = *tags1[i];
639 
640  if ( tag.name == "initrwgt" ) {
641  LHAinitrwgt irwgt(tag);
642  heprup.initrwgt = irwgt;
643  for ( int j = 0, M = tag.tags.size(); j < M; ++j ) {
644  XMLTag & ctag = *tag.tags[j];
645  if ( ctag.name == "weightgroup" ) {
646  LHAweightgroup wgroup(ctag);
647  string wgname = wgroup.name;
648  heprup.weightgroups.insert(make_pair(wgname, wgroup));
649 
650  string ss;
651  vector<XMLTag*> tags2 = XMLTag::findXMLTags(ctag.contents, &ss);
652  for ( int k = 0, O = tags2.size(); k < O; ++k ) {
653  const XMLTag & tagnow2 = *tags2[k];
654  if ( tagnow2.name == "weight" ) {
655  LHAweight wt(tagnow2);
656  string wtname = wt.id;
657  heprup.weights.insert(make_pair(wtname, wt));
658  }
659  }
660  for ( int k = 0, O = ctag.tags.size(); k < O; ++k ) {
661  const XMLTag & tagnow2 = *ctag.tags[k];
662  if ( tagnow2.name == "weight" ) {
663  LHAweight wt(tagnow2);
664  string wtname = wt.id;
665  heprup.weights.insert(make_pair(wtname, wt));
666  }
667  }
668  } else if ( ctag.name == "weight" ) {
669  string tname = ctag.attr["id"];
670  heprup.weights.insert(make_pair(tname, LHAweight(ctag)));
671  }
672  }
673  }
674  }
675 
676  heprup.generators.clear();
677  // Scan the init block for XML tags
678  leftovers="";
679  vector<XMLTag*> tags2 = XMLTag::findXMLTags(initComments, &leftovers);
680  if ( leftovers.find_first_not_of(" \t\n") == string::npos )
681  leftovers="";
682 
683  for ( int i = 0, N = tags2.size(); i < N; ++i ) {
684  const XMLTag & tag = *tags2[i];
685  if ( tag.name == "generator" ) {
686  heprup.generators.push_back(LHAgenerator(tag));
687  }
688  }
689 
690  for ( int i = 0, N = tags1.size(); i < N; ++i )
691  if (tags1[i]) delete tags1[i];
692  for ( int i = 0, N = tags2.size(); i < N; ++i )
693  if (tags2[i]) delete tags2[i];
694 
695  // Done
696  return true;
697 
698 }
699 
700 //--------------------------------------------------------------------------
701 
702 // Read an event from the file and store it in the hepeup
703 // object. Optional comment lines are stored in the eventComments
704 // member variable. return true if the read was successful.
705 
706 bool Reader::readEvent(HEPEUP * peup) {
707 
708  HEPEUP & eup = (peup? *peup: hepeup);
709  eup.clear();
710  eup.heprup = &heprup;
711  weights_detailed_vec.clear();
712  weightnames_detailed_vec.clear();
713 
714  // Check if the initialization was successful. Otherwise we will
715  // not read any events.
716  if ( heprup.NPRUP < 0 ) return false;
717  eventComments = "";
718  outsideBlock = "";
719  eup.NUP = 0;
720 
721  // Keep reading lines until we hit the next event or the end of
722  // the event block. Save any inbetween lines. Exit if we didn't
723  // find an event.
724  while ( getLine() && currentLine.find("<event") == string::npos )
725  outsideBlock += currentLine + "\n";
726 
727  // Get event attributes.
728  if (currentLine != "") {
729  string eventLine(currentLine);
730  eventLine += "</event>";
731  vector<XMLTag*> evtags = XMLTag::findXMLTags(eventLine);
732  XMLTag & evtag = *evtags[0];
733  for ( map<string,string>::const_iterator it = evtag.attr.begin();
734  it != evtag.attr.end(); ++it ) {
735  eup.attributes.insert(make_pair(it->first,it->second));
736  }
737  for ( int i = 0, N = evtags.size(); i < N; ++i )
738  if (evtags[i]) delete evtags[i];
739  }
740 
741  if ( !getLine() ) return false;
742 
743  // We found an event. The first line determines how many
744  // subsequent particle lines we have.
745  istringstream iss(currentLine);
746  if ( !( iss >> eup.NUP >> eup.IDPRUP >> eup.XWGTUP
747  >> eup.SCALUP >> eup.AQEDUP >> eup.AQCDUP ) )
748  return false;
749  eup.resize();
750 
751  // Read all particle lines.
752  for ( int i = 0; i < eup.NUP; ++i ) {
753  if ( !getLine() ) return false;
754  istringstream isss(currentLine);
755  if ( !( isss >> eup.IDUP[i] >> eup.ISTUP[i]
756  >> eup.MOTHUP[i].first >> eup.MOTHUP[i].second
757  >> eup.ICOLUP[i].first >> eup.ICOLUP[i].second
758  >> eup.PUP[i][0] >> eup.PUP[i][1] >> eup.PUP[i][2]
759  >> eup.PUP[i][3] >> eup.PUP[i][4]
760  >> eup.VTIMUP[i] >> eup.SPINUP[i] ) )
761  return false;
762  }
763 
764  // Now read any additional comments.
765  while ( getLine() && currentLine.find("</event>") == string::npos )
766  eventComments += currentLine + "\n";
767 
768  if ( file == NULL ) return false;
769 
770  eup.scalesSave = LHAscales(eup.SCALUP);
771 
772  // Scan the init block for XML tags
773  string leftovers;
774  vector<XMLTag*> tags = XMLTag::findXMLTags(eventComments, &leftovers);
775  if ( leftovers.find_first_not_of(" \t\n") == string::npos )
776  leftovers="";
777 
778  eventComments = "";
779  istringstream f(leftovers);
780  string l;
781  while (getline(f, l)) {
782  size_t p = l.find_first_not_of(" \t");
783  l.erase(0, p);
784  p = l.find_last_not_of(" \t");
785  if (string::npos != p) l.erase(p+1);
786  if (l.find_last_not_of("\n") != string::npos)
787  eventComments += l + "\n";
788  }
789 
790  for ( int i = 0, N = tags.size(); i < N; ++i ) {
791  XMLTag & tag = *tags[i];
792 
793  if ( tag.name == "weights" ) {
794  LHAweights wts(tag);
795  eup.weightsSave = wts;
796 
797  for ( int k = 0, M = int(wts.weights.size()); k < M; ++k ) {
798  eup.weights_compressed.push_back(wts.weights[k]);
799  }
800 
801  }
802  else if ( tag.name == "scales" ) {
803  eup.scalesSave = LHAscales(tag, eup.SCALUP);
804  }
805  else if ( tag.name == "rwgt" ) {
806  LHArwgt rwgt0(tag);
807  eup.rwgtSave = rwgt0;
808  string s;
809  vector<XMLTag*> tags2 = XMLTag::findXMLTags(rwgt0.contents, &s);
810  for ( int k = 0, M = tags2.size(); k < M; ++k ) {
811  const XMLTag & tagnow = *tags2[k];
812  if ( tagnow.name == "wgt" ) {
813  LHAwgt wt(tagnow);
814  eup.weights_detailed.insert(make_pair(wt.id, wt.contents));
815  weights_detailed_vec.push_back(wt.contents);
816  weightnames_detailed_vec.push_back(wt.id);
817  }
818  }
819  for ( int k = 0, M = tag.tags.size(); k < M; ++k ) {
820  const XMLTag & tagnow = *tag.tags[k];
821  if ( tagnow.name == "wgt" ) {
822  LHAwgt wt(tagnow);
823  eup.weights_detailed.insert(make_pair(wt.id, wt.contents));
824  weights_detailed_vec.push_back(wt.contents);
825  weightnames_detailed_vec.push_back(wt.id);
826  }
827  }
828  }
829  }
830 
831  for ( int i = 0, N = tags.size(); i < N; ++i ) if (tags[i]) delete tags[i];
832 
833  return true;
834 
835 }
836 
837 //==========================================================================
838 
839 // The Writer class is initialized with a stream to which to write a
840 // version 3.0 Les Houches Accord event file.
841 
842 //--------------------------------------------------------------------------
843 
844 // Write out an optional header block followed by the standard init
845 // block information together with any comment lines.
846 
847 void Writer::init() {
848 
849  // Write out the standard XML tag for the event file.
850  if ( version == 1 )
851  file << "<LesHouchesEvents version=\"1.0\">" << endl;
852  else
853  file << "<LesHouchesEvents version=\"3.0\">" << endl;
854 
855  file << setprecision(8);
856 
857  // Print headercomments and header init information.
858  file << "<header>" << endl;
859  file << hashline(headerStream.str(),true) << std::flush;
860  if ( version != 1 ) heprup.initrwgt.list(file);
861  file << "</header>" << endl;
862 
863  file << "<init>"<< endl
864  << " " << setw(8) << heprup.IDBMUP.first
865  << " " << setw(8) << heprup.IDBMUP.second
866  << " " << setw(14) << heprup.EBMUP.first
867  << " " << setw(14) << heprup.EBMUP.second
868  << " " << setw(4) << heprup.PDFGUP.first
869  << " " << setw(4) << heprup.PDFGUP.second
870  << " " << setw(4) << heprup.PDFSUP.first
871  << " " << setw(4) << heprup.PDFSUP.second
872  << " " << setw(4) << heprup.IDWTUP
873  << " " << setw(4) << heprup.NPRUP << endl;
874  heprup.resize();
875  for ( int i = 0; i < heprup.NPRUP; ++i )
876  file << " " << setw(14) << heprup.XSECUP[i]
877  << " " << setw(14) << heprup.XERRUP[i]
878  << " " << setw(14) << heprup.XMAXUP[i]
879  << " " << setw(6) << heprup.LPRUP[i] << endl;
880 
881  if ( version == 1 ) {
882  file << hashline(initStream.str(),true) << std::flush
883  << "</init>" << endl;
884  initStream.str("");
885  return;
886  }
887 
888  for ( int i = 0, N = heprup.generators.size(); i < N; ++i ) {
889  heprup.generators[i].list(file);
890  }
891 
892  file << hashline(initStream.str(),true) << std::flush
893  << "</init>" << endl;
894  initStream.str("");
895 }
896 
897 //--------------------------------------------------------------------------
898 
899 // Write out the event stored in hepeup, followed by optional
900 // comment lines.
901 
902 bool Writer::writeEvent(HEPEUP * peup, int pDigits) {
903 
904  HEPEUP & eup = (peup? *peup: hepeup);
905 
906  file << "<event";
907  for ( map<string,string>::const_iterator it = eup.attributes.begin();
908  it != eup.attributes.end(); ++it )
909  file << " " << it->first << "=\"" << it->second << "\"";
910  file << ">" << std::flush << endl;
911  file << " " << setw(4) << eup.NUP
912  << " " << setw(6) << eup.IDPRUP
913  << " " << setw(14) << eup.XWGTUP
914  << " " << setw(14) << eup.SCALUP
915  << " " << setw(14) << eup.AQEDUP
916  << " " << setw(14) << eup.AQCDUP << endl;
917  eup.resize();
918 
919  for ( int i = 0; i < eup.NUP; ++i )
920  file << " " << setw(8) << eup.IDUP[i]
921  << " " << setw(2) << eup.ISTUP[i]
922  << " " << setw(4) << eup.MOTHUP[i].first
923  << " " << setw(4) << eup.MOTHUP[i].second
924  << " " << setw(4) << eup.ICOLUP[i].first
925  << " " << setw(4) << eup.ICOLUP[i].second
926  << " " << setw(pDigits) << eup.PUP[i][0]
927  << " " << setw(pDigits) << eup.PUP[i][1]
928  << " " << setw(pDigits) << eup.PUP[i][2]
929  << " " << setw(pDigits) << eup.PUP[i][3]
930  << " " << setw(pDigits) << eup.PUP[i][4]
931  << " " << setw(1) << eup.VTIMUP[i]
932  << " " << setw(1) << eup.SPINUP[i] << endl;
933 
934  // Write event comments.
935  file << hashline(eventStream.str()) << std::flush;
936  eventStream.str("");
937 
938  if ( version != 1 ) {
939  eup.rwgtSave.list(file);
940  eup.weightsSave.list(file);
941  eup.scalesSave.list(file);
942  }
943 
944  file << "</event>" << endl;
945 
946  if ( !file ) return false;
947 
948  return true;
949 
950 }
951 
952 //--------------------------------------------------------------------------
953 
954 // Write out an event as a string.
955 
956 string Writer::getEventString(HEPEUP * peup) {
957 
958  HEPEUP & eup = (peup? *peup: hepeup);
959 
960  stringstream helper;
961 
962  helper << "<event";
963  for ( map<string,string>::const_iterator it = eup.attributes.begin();
964  it != eup.attributes.end(); ++it )
965  helper << " " << it->first << "=\"" << it->second << "\"";
966  helper << ">" << std::flush << endl;
967  helper << " " << setw(4) << eup.NUP
968  << " " << setw(6) << eup.IDPRUP
969  << " " << setw(14) << eup.XWGTUP
970  << " " << setw(14) << eup.SCALUP
971  << " " << setw(14) << eup.AQEDUP
972  << " " << setw(14) << eup.AQCDUP << endl;
973  eup.resize();
974 
975  for ( int i = 0; i < eup.NUP; ++i ) {
976  helper << " " << setw(8) << eup.IDUP[i]
977  << " " << setw(2) << eup.ISTUP[i]
978  << " " << setw(4) << eup.MOTHUP[i].first
979  << " " << setw(4) << eup.MOTHUP[i].second
980  << " " << setw(6) << eup.ICOLUP[i].first
981  << " " << setw(6) << eup.ICOLUP[i].second
982  << fixed
983  << setprecision(15)
984  << " " << setw(22) << eup.PUP[i][0]
985  << " " << setw(22) << eup.PUP[i][1]
986  << " " << setw(22) << eup.PUP[i][2]
987  << " " << setw(22) << eup.PUP[i][3]
988  << " " << setw(22) << eup.PUP[i][4]
989  << " " << setw(6) << eup.VTIMUP[i]
990  << " " << setw(6) << eup.SPINUP[i] << endl;
991  }
992 
993  // Write event comments.
994  helper << hashline(eventStream.str()) << std::flush;
995  eventStream.str("");
996 
997  if ( version != 1 ) {
998  eup.rwgtSave.list(helper);
999  eup.weightsSave.list(helper);
1000  eup.scalesSave.list(helper);
1001  }
1002 
1003  helper << "</event>" << endl;
1004 
1005  string helperString = helper.str();
1006  //event = helperString.c_str();
1007 
1008  return helperString;
1009 
1010 }
1011 
1012 //--------------------------------------------------------------------------
1013 
1014 // Make sure that each line in the string s starts with a
1015 // #-character and that the string ends with a new-line.
1016 
1017 string Writer::hashline(string s, bool comment) {
1018  string ret;
1019  istringstream is(s);
1020  string ss;
1021  while ( getline(is, ss) ) {
1022  if ( comment )
1023  ss = "# " + ss;
1024  ret += ss + '\n';
1025  }
1026  return ret;
1027 }
1028 
1029 //==========================================================================
1030 
1031 }
void version(std::ostream &os=std::cout)
print HepMC version
Definition: Version.h:27