Update to pythi8.170
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8170 / src / LesHouches.cxx
1 // LesHouches.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2012 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5
6 // Function definitions (not found in the header) for the LHAup and
7 // LHAupLHEF classes.
8
9 #include "LesHouches.h"
10
11 // Access time information.
12 #include <ctime>
13
14 // GZIP support.
15 #ifdef GZIPSUPPORT
16
17 // For GCC versions >= 4.6, can switch off shadow warnings.
18 #if (defined GZIPSUPPORT && ((__GNUC__ * 100) + __GNUC_MINOR__) >= 406)
19 #pragma GCC diagnostic ignored "-Wshadow"
20 #endif
21
22 // Boost includes.
23 #include "boost/iostreams/filtering_stream.hpp"
24 #include "boost/iostreams/filter/gzip.hpp"
25
26 // Switch shadow warnings back on.
27 #if (defined GZIPSUPPORT && ((__GNUC__ * 100) + __GNUC_MINOR__) >= 406)
28 #pragma GCC diagnostic warning "-Wshadow"
29 #endif
30
31 #endif // GZIPSUPPORT
32
33 namespace Pythia8 {
34
35 //==========================================================================
36
37 // LHAup class.
38
39 //--------------------------------------------------------------------------
40
41 // Constants: could be changed here if desired, but normally should not.
42 // These are of technical nature, as described for each.
43
44 // LHA convention with cross section in pb may require conversion from mb.
45 const double LHAup::CONVERTMB2PB = 1e9;
46
47 //--------------------------------------------------------------------------
48
49 // Print the initialization info; to check it worked.
50
51 void LHAup::listInit(ostream& os) {
52
53   // Header.
54   os << "\n --------  LHA initialization information  ------------ \n"; 
55
56   // Beam info.
57   os << fixed << setprecision(3) 
58      << "\n  beam    kind      energy  pdfgrp  pdfset \n" 
59      << "     A  " << setw(6) << idBeamASave 
60      <<  setw(12) << eBeamASave 
61      << setw(8) << pdfGroupBeamASave 
62      << setw(8) << pdfSetBeamASave << "\n" 
63      << "     B  " << setw(6) << idBeamBSave 
64      <<  setw(12) << eBeamBSave 
65      << setw(8) << pdfGroupBeamBSave 
66      << setw(8) << pdfSetBeamBSave << "\n"; 
67
68   // Event weighting strategy.
69   os << "\n  Event weighting strategy = " << setw(2) 
70      << strategySave << "\n" ;
71
72   // Process list.
73   os << scientific << setprecision(4) 
74      << "\n  Processes, with strategy-dependent cross section info \n" 
75      << "  number      xsec (pb)      xerr (pb)      xmax (pb) \n" ;
76   for (int ip = 0; ip < int(processes.size()); ++ip) {
77     os << setw(8) << processes[ip].idProc 
78        << setw(15) << processes[ip].xSecProc 
79        << setw(15) << processes[ip].xErrProc 
80        << setw(15) << processes[ip].xMaxProc << "\n";
81   }
82
83   // Finished.
84   os << "\n --------  End LHA initialization information  -------- \n"; 
85
86 }
87
88 //--------------------------------------------------------------------------
89
90 // Print the event info; to check it worked.
91
92 void LHAup::listEvent(ostream& os) {
93
94   // Header.
95   os << "\n --------  LHA event information and listing  -------------"
96      << "--------------------------------------------------------- \n"; 
97
98   // Basic event info.
99   os << scientific << setprecision(4) 
100      << "\n    process = " << setw(8) << idProc 
101      << "    weight = " << setw(12) << weightProc 
102      << "     scale = " << setw(12) << scaleProc << " (GeV) \n"
103      << "                   "
104      << "     alpha_em = " << setw(12) << alphaQEDProc 
105      << "    alpha_strong = " << setw(12) << alphaQCDProc << "\n";
106
107   // Particle list
108   os << fixed << setprecision(3) 
109      << "\n    Participating Particles \n" 
110      << "    no        id stat     mothers     colours      p_x        "
111      << "p_y        p_z         e          m        tau    spin \n" ;
112   for (int ip = 1; ip < int(particles.size()); ++ip) {
113     os << setw(6) << ip 
114        << setw(10) << particles[ip].idPart 
115        << setw(5) << particles[ip].statusPart 
116        << setw(6) << particles[ip].mother1Part
117        << setw(6) << particles[ip].mother2Part 
118        << setw(6) << particles[ip].col1Part
119        << setw(6) << particles[ip].col2Part 
120        << setw(11) << particles[ip].pxPart
121        << setw(11) << particles[ip].pyPart
122        << setw(11) << particles[ip].pzPart 
123        << setw(11) << particles[ip].ePart 
124        << setw(11) <<  particles[ip].mPart 
125        << setw(8) <<  particles[ip].tauPart 
126        << setw(8) <<  particles[ip].spinPart << "\n";
127   }
128
129   // PDF info - optional.
130   if (pdfIsSetSave) os << "\n     pdf: id1 =" << setw(5) << id1pdfSave  
131     << " id2 =" << setw(5) << id2pdfSave 
132     << " x1 ="  << scientific << setw(10) << x1pdfSave    
133     << " x2 =" << setw(10) << x2pdfSave 
134     << " scalePDF =" << setw(10) << scalePDFSave 
135     << " pdf1 =" << setw(10) << pdf1Save    
136     << " pdf2 =" << setw(10) << pdf2Save << "\n";    
137
138   // Finished.
139   os << "\n --------  End LHA event information and listing  ---------"
140      << "--------------------------------------------------------- \n"; 
141
142 }
143
144 //--------------------------------------------------------------------------
145
146 // Open and write header to a Les Houches Event File.
147
148 bool LHAup::openLHEF(string fileNameIn) {
149
150   // Open file for writing. Reset it to be empty.
151   fileName = fileNameIn;
152   const char* cstring = fileName.c_str();
153   osLHEF.open(cstring, ios::out | ios::trunc);  
154   if (!osLHEF) {
155     infoPtr->errorMsg("Error in LHAup::openLHEF:"
156       " could not open file", fileName);
157     return false;
158   }
159
160   // Read out current date and time.
161   time_t t = time(0);
162   strftime(dateNow,12,"%d %b %Y",localtime(&t));
163   strftime(timeNow,9,"%H:%M:%S",localtime(&t));
164
165   // Write header.
166   osLHEF << "<LesHouchesEvents version=\"1.0\">\n" 
167          << "<!--\n"
168          << "  File written by Pythia8::LHAup on " 
169          << dateNow << " at " << timeNow << "\n" 
170          << "-->" << endl;
171
172   // Done.
173   return true;
174
175 }
176
177 //--------------------------------------------------------------------------
178
179 // Write initialization information to a Les Houches Event File.
180
181 bool LHAup::initLHEF() {
182
183   // Write information on beams. 
184   osLHEF << "<init>\n" << scientific << setprecision(6)
185          << "  " << idBeamASave       << "  " << idBeamBSave 
186          << "  " << eBeamASave        << "  " << eBeamBSave 
187          << "  " << pdfGroupBeamASave << "  " << pdfGroupBeamBSave 
188          << "  " << pdfSetBeamASave   << "  " << pdfSetBeamBSave 
189          << "  " << strategySave      << "  " << processes.size() << "\n";
190
191   // Write information on all the subprocesses.
192   for (int ip = 0; ip < int(processes.size()); ++ip) 
193     osLHEF << " " << setw(13) << processes[ip].xSecProc 
194            << " " << setw(13) << processes[ip].xErrProc 
195            << " " << setw(13) << processes[ip].xMaxProc 
196            << " " << setw(6) << processes[ip].idProc << "\n";
197
198   // Done.
199   osLHEF << "</init>" << endl;
200   return true;
201
202 }
203
204 //--------------------------------------------------------------------------
205
206 // Write event information to a Les Houches Event File.
207 // Normal mode is to line up event info in columns, but the non-verbose
208 // altnernative saves space at the expense of human readability.
209
210 bool LHAup::eventLHEF(bool verbose) {
211
212   // Default verbose option.
213   if (verbose) { 
214
215     // Write information on process as such. 
216     osLHEF << "<event>\n" << scientific << setprecision(6)
217            << " " << setw(5) << particles.size() - 1 
218            << " " << setw(5) << idProc
219            << " " << setw(13) << weightProc 
220            << " " << setw(13) << scaleProc
221            << " " << setw(13) << alphaQEDProc
222            << " " << setw(13) << alphaQCDProc << "\n";
223
224     // Write information on the particles, excluding zeroth. 
225     for (int ip = 1; ip < int(particles.size()); ++ip) {
226       LHAParticle& ptNow = particles[ip];
227       osLHEF << " " << setw(8) << ptNow.idPart 
228              << " " << setw(5) << ptNow.statusPart 
229              << " " << setw(5) << ptNow.mother1Part
230              << " " << setw(5) << ptNow.mother2Part 
231              << " " << setw(5) << ptNow.col1Part
232              << " " << setw(5) << ptNow.col2Part << setprecision(10)
233              << " " << setw(17) << ptNow.pxPart
234              << " " << setw(17) << ptNow.pyPart
235              << " " << setw(17) << ptNow.pzPart 
236              << " " << setw(17) << ptNow.ePart 
237              << " " << setw(17) <<  ptNow.mPart << setprecision(6);
238       if (ptNow.tauPart == 0.) osLHEF << " 0.";
239       else osLHEF << " " << setw(13) << ptNow.tauPart;
240       if (ptNow.spinPart == 9.) osLHEF << " 9.";
241       else osLHEF << " " << setw(13) << ptNow.spinPart;
242       osLHEF << "\n";
243     }
244
245     // Optionally write information on PDF values at hard interaction.
246     if (pdfIsSetSave) osLHEF << "#pdf" 
247              << " " << setw(4) << id1pdfSave
248              << " " << setw(4) << id2pdfSave
249              << " " << setw(13) << x1pdfSave 
250              << " " << setw(13) << x2pdfSave 
251              << " " << setw(13) << scalePDFSave 
252              << " " << setw(13) << pdf1Save 
253              << " " << setw(13) << pdf2Save << "\n"; 
254
255   // Alternative non-verbose option. 
256   } else {
257
258     // Write information on process as such. 
259     osLHEF << "<event>\n" << scientific << setprecision(6)
260            << particles.size() - 1 << " " << idProc       << " "
261            << weightProc           << " " << scaleProc    << " "
262            << alphaQEDProc         << " " << alphaQCDProc << "\n";
263
264     // Write information on the particles, excluding zeroth. 
265     for (int ip = 1; ip < int(particles.size()); ++ip) {
266       LHAParticle& ptNow = particles[ip];
267       osLHEF        << ptNow.idPart      << " " << ptNow.statusPart 
268              << " " << ptNow.mother1Part << " " << ptNow.mother2Part 
269              << " " << ptNow.col1Part    << " " << ptNow.col2Part 
270              << setprecision(10)         << " " << ptNow.pxPart
271              << " " << ptNow.pyPart      << " " << ptNow.pzPart 
272              << " " << ptNow.ePart       << " " << ptNow.mPart 
273              << setprecision(6);
274       if (ptNow.tauPart == 0.) osLHEF << " 0.";
275       else osLHEF << " " << setw(13) << ptNow.tauPart;
276       if (ptNow.spinPart == 9.) osLHEF << " 9.";
277       else osLHEF << " " << setw(13) << ptNow.spinPart;
278       osLHEF << "\n";
279     }
280
281     // Optionally write information on PDF values at hard interaction.
282     if (pdfIsSetSave) osLHEF << "#pdf" << " " << id1pdfSave
283              << " " << id2pdfSave << " " << x1pdfSave << " " << x2pdfSave 
284              << " " << scalePDFSave << " " << pdf1Save << " " << pdf2Save 
285              << "\n"; 
286   }
287
288   // Done.
289   osLHEF << "</event>" << endl;
290   return true;
291
292 }
293
294 //--------------------------------------------------------------------------
295
296 // Write end of a Les Houches Event File and close it.
297
298 bool LHAup::closeLHEF(bool updateInit) {
299
300   // Write an end to the file.
301   osLHEF << "</LesHouchesEvents>" << endl;
302   osLHEF.close();
303
304   // Optionally update the cross section information.
305   if (updateInit) {
306     const char* cstring = fileName.c_str();
307     osLHEF.open(cstring, ios::in | ios::out); 
308   
309     // Rewrite header; identically with what openLHEF did.
310     osLHEF << "<LesHouchesEvents version=\"1.0\">\n" 
311            << "<!--\n"
312            << "  File written by Pythia8::LHAup on " 
313            << dateNow << " at " << timeNow << "\n" 
314            << "-->" << endl;
315
316     // Redo initialization information.
317     initLHEF();
318     osLHEF.close();
319   }  
320
321   // Done.
322   return true;
323
324 }
325
326 //--------------------------------------------------------------------------
327
328 // Read in initialization information from a Les Houches Event File.
329
330 bool LHAup::setInitLHEF(istream& is, bool readHeaders) {
331
332   // Check that first line is consistent with proper LHEF file.
333   string line;
334   if (!getline(is, line)) return false;
335   if (line.find("<LesHouchesEvents") == string::npos) return false;  
336   if (line.find("version=\"1.0\"" ) == string::npos ) return false;
337
338   // What to search for if reading headers; if not reading
339   // headers then return to default behaviour
340   string headerTag = (readHeaders) ? "<header>" : "<init";
341
342   // Loop over lines until an <init (or optionally <header>) tag
343   // is found first on a line.
344   string tag = " ";
345   do { 
346     if (!getline(is, line)) return false;
347     if (line.find_first_not_of(" \n\t\v\b\r\f\a") != string::npos) {
348       istringstream getfirst(line);
349       getfirst >> tag;
350       if (!getfirst) return false;
351     }
352   } while (tag != "<init>" && tag != "<init" && tag != headerTag);
353
354   // If header tag found, process if required
355   if (readHeaders == true && tag == headerTag) {
356     // Temporary local storage
357     map < string, string > headerMap;
358
359     // Loop over lines until an <init> tag is found.
360     bool read = true, newKey = false;
361     string key = "base";
362     vector < string > keyVec;
363     while (true) { 
364       if (!getline(is, line)) return false;
365
366       // Check if this line is a tag; '<' as first character,
367       // '>' as last character, exclusing whitespace
368       size_t pos1 = line.find_first_not_of(" \n\t\v\b\r\f\a");
369       size_t pos2 = line.find_last_not_of(" \n\t\v\b\r\f\a");
370       if (pos1 != string::npos && line[pos1] == '<' &&
371           pos2 != string::npos && line[pos2] == '>' &&
372           pos1 < pos2) {
373
374         // Only take the first word of the tag
375         tag = line.substr(pos1 + 1, pos2 - pos1 - 1);
376         istringstream getfirst(tag);
377         getfirst >> tag;
378
379         // Tag present, so handle here
380         if (getfirst) {
381
382           // Exit condition 
383           if (tag == "init") break;
384
385           // End of header block; keep reading until <init> tag,
386           // but do not store any further information
387           else if (tag == "/header") {
388             read = false;
389             continue;
390
391           // Opening tag
392           } else if (tag[0] != '/') {
393             keyVec.push_back(tag);
394             newKey = true;
395             continue;
396
397           // Closing tag that matches current key
398           } else if (tag == "/" + keyVec.back()) {
399             keyVec.pop_back();
400             newKey = true;
401             continue;
402           }
403
404         } // if (getfirst)
405       }
406
407       // At this point we have a line that is not a tag; if no longer
408       // reading headers then keep going
409       if (!read) continue;
410       
411       // Check for key change
412       if (newKey) {
413         if (keyVec.empty()) key = "base";
414         else                key = keyVec[0];
415         for (size_t i = 1; i < keyVec.size(); i++)
416           key += "." + keyVec[i];
417         newKey = false;
418       }
419
420       // Append information to local storage
421       headerMap[key] += line + "\n";
422
423     } // while (true)
424
425     // Copy information to info using LHAup::setInfoHeader
426     for (map < string, string >::iterator it = headerMap.begin();
427         it != headerMap.end(); it++)
428       setInfoHeader(it->first, it->second);
429
430   } // if (readHeaders == true && tag == headerTag)
431   
432   // Read in beam and strategy info, and store it. 
433   int idbmupA, idbmupB;
434   double ebmupA, ebmupB;
435   int pdfgupA, pdfgupB, pdfsupA, pdfsupB, idwtup, nprup;
436   if (!getline(is, line)) return false;
437   istringstream getbms(line);
438   getbms >> idbmupA >> idbmupB >> ebmupA >> ebmupB >> pdfgupA 
439      >> pdfgupB >> pdfsupA >> pdfsupB >> idwtup >> nprup;
440   if (!getbms) return false;
441   setBeamA(idbmupA, ebmupA, pdfgupA, pdfsupA);
442   setBeamB(idbmupB, ebmupB, pdfgupB, pdfsupB);
443   setStrategy(idwtup);
444
445   // Read in process info, one process at a time, and store it.
446   double xsecup, xerrup, xmaxup;
447   xSecSumSave = 0.;
448   xErrSumSave = 0.;
449   int lprup; 
450   for (int ip = 0; ip < nprup; ++ip) { 
451     if (!getline(is, line)) return false;
452     istringstream getpro(line);
453     getpro >> xsecup >> xerrup >> xmaxup >> lprup ;
454     if (!getpro) return false;
455     addProcess(lprup, xsecup, xerrup, xmaxup);
456     xSecSumSave += xsecup;
457     xErrSumSave += pow2(xerrup);
458   }
459   xErrSumSave = sqrt(xErrSumSave);
460
461   // Reading worked.
462   return true;
463
464 }
465
466 //--------------------------------------------------------------------------
467
468 // Read in event information from a Les Houches Event File,
469 // into a staging area where it can be reused by setOldEventLHEF.
470
471 bool LHAup::setNewEventLHEF(istream& is) {
472   
473   // Loop over lines until an <event tag is found first on a line.
474   string line, tag;
475   do { 
476     if (!getline(is, line)) return false;
477     if (line.find_first_not_of(" \n\t\v\b\r\f\a") != string::npos) {
478       istringstream getfirst(line);
479       getfirst >> tag;
480       if (!getfirst) return false;
481     }
482   } while (tag != "<event>" && tag != "<event"); 
483
484   // Read in process info and store it.
485   if (!getline(is, line)) return false;
486   istringstream getpro(line);
487   getpro >> nupSave >> idprupSave >> xwgtupSave >> scalupSave
488     >> aqedupSave >> aqcdupSave;
489   if (!getpro) return false;
490
491   // Reset particlesSave vector, add slot-0 empty particle.
492   particlesSave.clear(); 
493   particlesSave.push_back( LHAParticle() );
494
495   // Read in particle info one by one, and store it.
496   // Note unusual C++ loop range, to better reflect LHA/Fortran standard.
497   // (Recall that process(...) above added empty particle at index 0.) 
498   int idup, istup, mothup1, mothup2, icolup1, icolup2; 
499   double pup1, pup2, pup3, pup4, pup5, vtimup, spinup;
500   for (int ip = 1; ip <= nupSave; ++ip) { 
501     if (!getline(is, line)) return false;
502     istringstream getall(line);
503     getall >> idup >> istup >> mothup1 >> mothup2 >> icolup1 >> icolup2 
504       >> pup1 >> pup2 >> pup3 >> pup4 >> pup5 >> vtimup >> spinup;
505     if (!getall) return false;   
506     particlesSave.push_back( LHAParticle( idup, istup, mothup1, mothup2, 
507       icolup1, icolup2, pup1, pup2, pup3, pup4, pup5, vtimup, spinup) );
508   }
509
510   // Flavour and x values of hard-process initiators.
511   id1InSave = particlesSave[1].idPart;
512   id2InSave = particlesSave[2].idPart;
513   x1InSave  = (eBeamASave > 0.) ? particlesSave[1].ePart / eBeamASave : 0.; 
514   x2InSave  = (eBeamBSave > 0.) ? particlesSave[2].ePart / eBeamBSave : 0.; 
515
516   // Continue parsing till </event>. Extract pdf info if present.
517   getPDFSave = false;
518   do { 
519     if (!getline(is, line)) return false;
520     istringstream getpdf(line);
521     getpdf >> tag;
522     if (!getpdf) return false;
523     if (tag == "#pdf") {
524       getpdf >> id1pdfInSave >> id2pdfInSave >> x1pdfInSave >> x2pdfInSave 
525              >> scalePDFInSave >> pdf1InSave >> pdf2InSave;
526       if (!getpdf) return false;
527       getPDFSave = true;
528     }
529   } while (tag != "</event>" && tag != "</event"); 
530
531   // Need id and x values even when no PDF info. Rest empty.
532   if (!getPDFSave) {
533     id1pdfInSave   = id1InSave;
534     id2pdfInSave   = id2InSave;
535     x1pdfInSave    = x1InSave; 
536     x2pdfInSave    = x2InSave; 
537     scalePDFInSave = 0.;
538     pdf1InSave     = 0.;
539     pdf2InSave     = 0.;
540   }
541   
542   // Reading worked.
543   return true;
544
545 }
546
547 //--------------------------------------------------------------------------
548
549 // Make current event information read in by setNewEventLHEF.
550
551 bool LHAup::setOldEventLHEF() {
552
553   // Store saved event, optionally also parton density information.
554   setProcess( idprupSave, xwgtupSave, scalupSave, aqedupSave, aqcdupSave);
555   for (int ip = 1; ip <= nupSave; ++ip) addParticle( particlesSave[ip] );
556   setIdX( id1InSave, id2InSave, x1InSave, x2InSave);  
557   setPdf( id1pdfInSave, id2pdfInSave, x1pdfInSave, x2pdfInSave, 
558     scalePDFInSave, pdf1InSave, pdf2InSave, getPDFSave);  
559
560   // Done;
561   return true;
562
563 }
564
565 //--------------------------------------------------------------------------
566
567 // Open a file using provided ifstream and return a pointer to an istream
568 // that can be used to process the file. This is designed to handle
569 // GZIPSUPPORT in a transparent manner:
570 //  - no GZIPSUPPORT, the istream pointer is exactly the ifstream object.
571 //  - with GZIPSUPPORT, the istream pointer is a Boost filtering istream
572 //    (memory allocated), which reads from the ifstream and decompresses.
573 // The companion 'closeFile' can be used to correctly close a file and
574 // deallocate memory if needed. Note that a gzip filter is only applied
575 // if the final three characters of the filename are '.gz'.
576
577 istream* LHAup::openFile(const char *fn, ifstream &ifs) {
578   // Open the file
579   ifs.open(fn);
580
581 // No gzip support, so just return pointer to the istream
582 #ifndef GZIPSUPPORT
583   return (istream *) &ifs;
584
585 // Gzip support, so construct istream with gzip support
586 #else
587   // Boost filtering istream
588   boost::iostreams::filtering_istream *fis =
589     new boost::iostreams::filtering_istream();
590
591   // Pass along the 'good()' flag, so code elsewhere works unmodified.
592   if (!ifs.good()) fis->setstate(ios_base::badbit);
593
594   // Check filename ending to decide which filters to apply.
595   else {
596     const char *last = strrchr(fn, '.');
597     if (last && strncmp(last, ".gz", 3) == 0)
598       fis->push(boost::iostreams::gzip_decompressor());
599     fis->push(ifs);
600   }
601   return (istream *) fis;
602
603 #endif // GZIPSUPPORT
604 }
605
606 //--------------------------------------------------------------------------
607
608 // Companion method to 'openFile', above.
609 // Correctly deallocates memory if required before closing the file.
610
611 void LHAup::closeFile(istream *&is, ifstream &ifs) {
612   // If the istream pointer is not NULL and is not the
613   // same as the ifstream, then delete pointer.
614   if (is && is != &ifs) delete is;
615   is = NULL;
616
617   // Close the file
618   if (ifs.is_open()) ifs.close();
619 }
620
621
622 //==========================================================================
623
624 // LHAupLHEF class.
625
626 //--------------------------------------------------------------------------
627
628 // Constructor.
629
630 LHAupLHEF::LHAupLHEF(const char* fileIn, const char* headerIn,
631                      bool readHeadersIn)
632     : is(NULL), isHead(NULL), readHeaders(readHeadersIn) {
633
634   // Open LHEF and optionally header file as well. Note that both
635   // are opened here so that initialisation can be aborted if
636   // either of the files is missing, see fileFound().
637   is     = openFile(fileIn, ifs);
638   isHead = (headerIn == NULL) ? is : openFile(headerIn, ifsHead);
639 }
640
641 //--------------------------------------------------------------------------
642
643 // Destructor.
644
645 LHAupLHEF::~LHAupLHEF() {
646   // Close files
647   closeAllFiles();
648 }
649
650 //==========================================================================
651
652 // LHAupFromPYTHIA8 class.
653
654 //--------------------------------------------------------------------------
655
656 // Read in initialization information from PYTHIA 8.
657
658 bool LHAupFromPYTHIA8::setInit() {
659   
660   // Read in beam from Info class. Parton density left empty. 
661   int    idbmupA = infoPtr->idA();
662   int    idbmupB = infoPtr->idB();
663   double ebmupA  = infoPtr->eA();
664   double ebmupB  = infoPtr->eB();
665   int    pdfgupA = 0; 
666   int    pdfgupB = 0; 
667   int    pdfsupA = 0; 
668   int    pdfsupB = 0; 
669   setBeamA(idbmupA, ebmupA, pdfgupA, pdfsupA);
670   setBeamB(idbmupB, ebmupB, pdfgupB, pdfsupB);
671
672   // Currently only one allowed strategy.
673   int    idwtup = 3;
674   setStrategy(idwtup);
675
676   // Only one process with dummy information. (Can overwrite at the end.)
677   int    lprup  = 9999; 
678   double xsecup = 1.;
679   double xerrup = 0.;
680   double xmaxup = 1.;
681   addProcess(lprup, xsecup, xerrup, xmaxup);
682
683   // Done.
684   return true;
685
686 }
687
688 //--------------------------------------------------------------------------
689
690 // Read in event information from PYTHIA 8.
691
692 bool LHAupFromPYTHIA8::setEvent( int ) {
693
694   // Read process information from Info class, and store it.
695   // Note: renormalization scale here, factorization further down.
696   // For now always convert to process 9999, instead of infoPtr->code(). 
697   int    idprup = 9999;
698   double xwgtup = infoPtr->weight();
699   double scalup = infoPtr->QRen();
700   double aqedup = infoPtr->alphaEM();
701   double aqcdup = infoPtr->alphaS();
702   setProcess(idprup, xwgtup, scalup, aqedup, aqcdup);
703
704   // Read in particle info one by one, excluding zero and beams, and store it.
705   // Note unusual C++ loop range, to better reflect LHA/Fortran standard.
706   int nup   = processPtr->size() - 3;
707   int    idup, statusup, istup, mothup1, mothup2, icolup1, icolup2; 
708   double pup1, pup2, pup3, pup4, pup5, vtimup, spinup;
709   for (int ip = 1; ip <= nup; ++ip) {
710     Particle& particle = (*processPtr)[ip + 2];
711     idup     = particle.id(); 
712     // Convert from PYTHIA8 to LHA status codes.
713     statusup = particle.status();
714     if (ip < 3)            istup = -1;
715     else if (statusup < 0) istup =  2;
716     else                   istup =  1;
717     mothup1  = max(0, particle.mother1() - 2); 
718     mothup2  = max(0, particle.mother2() - 2); 
719     icolup1  = particle.col();
720     icolup2  = particle.acol();
721     pup1     = particle.px();
722     pup2     = particle.py();
723     pup3     = particle.pz();
724     pup4     = particle.e();
725     pup5     = particle.m();
726     vtimup   = particle.tau(); 
727     spinup   = particle.pol();
728     addParticle(idup, istup, mothup1, mothup2, icolup1, icolup2,
729       pup1, pup2, pup3, pup4, pup5, vtimup, spinup) ;
730   }
731
732   // Extract hard-process initiator information from Info class, and store it.
733   int    id1up      = infoPtr->id1();
734   int    id2up      = infoPtr->id2();
735   double x1up       = infoPtr->x1();
736   double x2up       = infoPtr->x2();
737   setIdX( id1up, id2up, x1up, x2up);
738
739   // Also extract pdf information from Info class, and store it.
740   int    id1pdfup   = infoPtr->id1pdf();
741   int    id2pdfup   = infoPtr->id2pdf();
742   double x1pdfup    = infoPtr->x1pdf();
743   double x2pdfup    = infoPtr->x2pdf();
744   double scalePDFup = infoPtr->QFac();
745   double pdf1up     = infoPtr->pdf1();
746   double pdf2up     = infoPtr->pdf2();
747   setPdf(id1pdfup, id2pdfup, x1pdfup, x2pdfup, scalePDFup, pdf1up, 
748     pdf2up, true);  
749
750   // Done.
751   return true;
752
753 }
754
755 //--------------------------------------------------------------------------
756
757 //  Update cross-section information at the end of the run.
758
759 bool LHAupFromPYTHIA8::updateSigma() {
760
761   // Read out information from PYTHIA 8 and send it in to LHA.
762   double sigGen = CONVERTMB2PB * infoPtr->sigmaGen();
763   double sigErr = CONVERTMB2PB * infoPtr->sigmaErr(); 
764   setXSec(0, sigGen);
765   setXErr(0, sigErr);
766
767   // Done.
768   return true;
769
770 }
771  
772 //==========================================================================
773
774 } // end namespace Pythia8