StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
LesHouches.cc
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) << id1Save
131  << " id2 =" << setw(5) << id2Save
132  << " x1 =" << scientific << setw(10) << x1Save
133  << " x2 =" << setw(10) << x2Save
134  << " scalePDF =" << setw(10) << scalePDFSave
135  << " xpdf1 =" << setw(10) << xpdf1Save
136  << " xpdf2 =" << setw(10) << xpdf2Save << "\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 
208 bool LHAup::eventLHEF() {
209 
210  // Write information on process as such.
211  osLHEF << "<event>\n" << scientific << setprecision(6)
212  << " " << setw(5) << particles.size() - 1
213  << " " << setw(5) << idProc
214  << " " << setw(13) << weightProc
215  << " " << setw(13) << scaleProc
216  << " " << setw(13) << alphaQEDProc
217  << " " << setw(13) << alphaQCDProc << "\n";
218 
219  // Write information on the particles, excluding zeroth.
220  for (int ip = 1; ip < int(particles.size()); ++ip) {
221  osLHEF << " " << setw(8) << particles[ip].idPart
222  << " " << setw(5) << particles[ip].statusPart
223  << " " << setw(5) << particles[ip].mother1Part
224  << " " << setw(5) << particles[ip].mother2Part
225  << " " << setw(5) << particles[ip].col1Part
226  << " " << setw(5) << particles[ip].col2Part << setprecision(10)
227  << " " << setw(17) << particles[ip].pxPart
228  << " " << setw(17) << particles[ip].pyPart
229  << " " << setw(17) << particles[ip].pzPart
230  << " " << setw(17) << particles[ip].ePart
231  << " " << setw(17) << particles[ip].mPart << setprecision(6);
232  if (particles[ip].tauPart == 0.) osLHEF << " 0.";
233  else osLHEF << " " << setw(13) << particles[ip].tauPart;
234  if (particles[ip].spinPart == 9.) osLHEF << " 9.";
235  else osLHEF << " " << setw(13) << particles[ip].spinPart;
236  osLHEF << "\n";
237  }
238 
239  // Optionally write information on PDF values at hard interaction.
240  if (pdfIsSetSave) osLHEF << "#pdf"
241  << " " << setw(4) << id1Save
242  << " " << setw(4) << id2Save
243  << " " << setw(13) << x1Save
244  << " " << setw(13) << x2Save
245  << " " << setw(13) << scalePDFSave
246  << " " << setw(13) << xpdf1Save
247  << " " << setw(13) << xpdf2Save << "\n";
248 
249  // Done.
250  osLHEF << "</event>" << endl;
251  return true;
252 
253 }
254 
255 //--------------------------------------------------------------------------
256 
257 // Write end of a Les Houches Event File and close it.
258 
259 bool LHAup::closeLHEF(bool updateInit) {
260 
261  // Write an end to the file.
262  osLHEF << "</LesHouchesEvents>" << endl;
263  osLHEF.close();
264 
265  // Optionally update the cross section information.
266  if (updateInit) {
267  const char* cstring = fileName.c_str();
268  osLHEF.open(cstring, ios::in | ios::out);
269 
270  // Rewrite header; identically with what openLHEF did.
271  osLHEF << "<LesHouchesEvents version=\"1.0\">\n"
272  << "<!--\n"
273  << " File written by Pythia8::LHAup on "
274  << dateNow << " at " << timeNow << "\n"
275  << "-->" << endl;
276 
277  // Redo initialization information.
278  initLHEF();
279  osLHEF.close();
280  }
281 
282  // Done.
283  return true;
284 
285 }
286 
287 //--------------------------------------------------------------------------
288 
289 // Read in initialization information from a Les Houches Event File.
290 
291 bool LHAup::setInitLHEF(istream& is) {
292 
293  // Check that first line is consistent with proper LHEF file.
294  string line;
295  if (!getline(is, line)) return false;
296  if (line.find("<LesHouchesEvents") == string::npos) return false;
297  if (line.find("version=\"1.0\"" ) == string::npos ) return false;
298 
299  // Loop over lines until an <init tag is found first on a line.
300  string tag = " ";
301  do {
302  if (!getline(is, line)) return false;
303  if (line.find_first_not_of(" \n\t\v\b\r\f\a") != string::npos) {
304  istringstream getfirst(line);
305  getfirst >> tag;
306  if (!getfirst) return false;
307  }
308  } while (tag != "<init>" && tag != "<init");
309 
310  // Read in beam and strategy info, and store it.
311  int idbmupA, idbmupB;
312  double ebmupA, ebmupB;
313  int pdfgupA, pdfgupB, pdfsupA, pdfsupB, idwtup, nprup;
314  if (!getline(is, line)) return false;
315  istringstream getbms(line);
316  getbms >> idbmupA >> idbmupB >> ebmupA >> ebmupB >> pdfgupA
317  >> pdfgupB >> pdfsupA >> pdfsupB >> idwtup >> nprup;
318  if (!getbms) return false;
319  setBeamA(idbmupA, ebmupA, pdfgupA, pdfsupA);
320  setBeamB(idbmupB, ebmupB, pdfgupB, pdfsupB);
321  setStrategy(idwtup);
322 
323  // Read in process info, one process at a time, and store it.
324  double xsecup, xerrup, xmaxup;
325  int lprup;
326  for (int ip = 0; ip < nprup; ++ip) {
327  if (!getline(is, line)) return false;
328  istringstream getpro(line);
329  getpro >> xsecup >> xerrup >> xmaxup >> lprup ;
330  if (!getpro) return false;
331  addProcess(lprup, xsecup, xerrup, xmaxup);
332  }
333 
334  // Reading worked.
335  return true;
336 
337 }
338 
339 //--------------------------------------------------------------------------
340 
341 // Read in event information from a Les Houches Event File,
342 // into a staging area where it can be reused by setOldEventLHEF.
343 
344 bool LHAup::setNewEventLHEF(istream& is) {
345 
346  // Loop over lines until an <event tag is found first on a line.
347  string line, tag;
348  do {
349  if (!getline(is, line)) return false;
350  if (line.find_first_not_of(" \n\t\v\b\r\f\a") != string::npos) {
351  istringstream getfirst(line);
352  getfirst >> tag;
353  if (!getfirst) return false;
354  }
355  } while (tag != "<event>" && tag != "<event");
356 
357  // Read in process info and store it.
358  if (!getline(is, line)) return false;
359  istringstream getpro(line);
360  getpro >> nupSave >> idprupSave >> xwgtupSave >> scalupSave
361  >> aqedupSave >> aqcdupSave;
362  if (!getpro) return false;
363 
364  // Reset particlesSave vector, add slot-0 empty particle.
365  particlesSave.clear();
366  particlesSave.push_back( LHAParticle() );
367 
368  // Read in particle info one by one, and store it.
369  // Note unusual C++ loop range, to better reflect LHA/Fortran standard.
370  // (Recall that process(...) above added empty particle at index 0.)
371  int idup, istup, mothup1, mothup2, icolup1, icolup2;
372  double pup1, pup2, pup3, pup4, pup5, vtimup, spinup;
373  for (int ip = 1; ip <= nupSave; ++ip) {
374  if (!getline(is, line)) return false;
375  istringstream getall(line);
376  getall >> idup >> istup >> mothup1 >> mothup2 >> icolup1 >> icolup2
377  >> pup1 >> pup2 >> pup3 >> pup4 >> pup5 >> vtimup >> spinup;
378  if (!getall) return false;
379  particlesSave.push_back( LHAParticle( idup, istup, mothup1, mothup2,
380  icolup1, icolup2, pup1, pup2, pup3, pup4, pup5, vtimup, spinup) );
381  }
382 
383  // Continue parsing till </event>. Extract pdf info if present.
384  getPDFSave = false;
385  do {
386  if (!getline(is, line)) return false;
387  istringstream getpdf(line);
388  getpdf >> tag;
389  if (!getpdf) return false;
390  if (tag == "#pdf") {
391  getpdf >> id1InSave >> id2InSave >> x1InSave >> x2InSave
392  >> scalePDFInSave >> xpdf1InSave >> xpdf2InSave;
393  if (!getpdf) return false;
394  getPDFSave = true;
395  }
396  } while (tag != "</event>" && tag != "</event");
397 
398  // Need id and x values even when no PDF info. Rest empty.
399  if (!getPDFSave) {
400  id1InSave = particlesSave[1].idPart;
401  id2InSave = particlesSave[2].idPart;
402  x1InSave = particlesSave[1].ePart / eBeamASave;
403  x2InSave = particlesSave[2].ePart / eBeamBSave;
404  scalePDFInSave = 0.;
405  xpdf1InSave = 0.;
406  xpdf2InSave = 0.;
407  }
408 
409  // Reading worked.
410  return true;
411 
412 }
413 
414 //--------------------------------------------------------------------------
415 
416 // Make current event information read in by setNewEventLHEF.
417 
418 bool LHAup::setOldEventLHEF() {
419 
420  // Store saved event, optionally also parton density information.
421  setProcess(idprupSave, xwgtupSave, scalupSave, aqedupSave, aqcdupSave);
422  for (int ip = 1; ip <= nupSave; ++ip) addParticle( particlesSave[ip] );
423  setPdf(id1InSave, id2InSave, x1InSave, x2InSave, scalePDFInSave,
424  xpdf1InSave, xpdf2InSave, getPDFSave);
425 
426  // Done;
427  return true;
428 
429 }
430 
431 
432 //==========================================================================
433 
434 // LHAupLHEF class.
435 
436 //--------------------------------------------------------------------------
437 
438 // Constructor.
439 
440 LHAupLHEF::LHAupLHEF(const char* fileIn) : ifs(fileIn), is(NULL) {
441 
442 // Construct istream without gzip support.
443 #ifndef GZIPSUPPORT
444  is = (istream *) &ifs;
445 
446 // Construct istream with gzip support.
447 #else
448  boost::iostreams::filtering_istream *fis =
449  new boost::iostreams::filtering_istream();
450 
451  // Pass along the 'good()' flag, so code elsewhere works unmodified.
452  if (!ifs.good()) fis->setstate(ios_base::badbit);
453 
454  // Check filename ending to decide which filters to apply.
455  else {
456  const char *last = strrchr(fileIn, '.');
457  if (last && strncmp(last, ".gz", 3) == 0)
458  fis->push(boost::iostreams::gzip_decompressor());
459  fis->push(ifs);
460  }
461  is = (istream *) fis;
462 #endif
463 
464 }
465 
466 //--------------------------------------------------------------------------
467 
468 // Destructor.
469 
470 LHAupLHEF::~LHAupLHEF() {
471 
472 // Delete istream if constructed.
473 #ifdef GZIPSUPPORT
474  if (is) delete is;
475 #endif
476 
477 }
478 
479 //==========================================================================
480 
481 // LHAupFromPYTHIA8 class.
482 
483 //--------------------------------------------------------------------------
484 
485 // Read in initialization information from PYTHIA 8.
486 
487 bool LHAupFromPYTHIA8::setInit() {
488 
489  // Read in beam from Info class. Parton density left empty.
490  int idbmupA = infoPtr->idA();
491  int idbmupB = infoPtr->idB();
492  double ebmupA = infoPtr->eA();
493  double ebmupB = infoPtr->eB();
494  int pdfgupA = 0;
495  int pdfgupB = 0;
496  int pdfsupA = 0;
497  int pdfsupB = 0;
498  setBeamA(idbmupA, ebmupA, pdfgupA, pdfsupA);
499  setBeamB(idbmupB, ebmupB, pdfgupB, pdfsupB);
500 
501  // Currently only one allowed strategy.
502  int idwtup = 3;
503  setStrategy(idwtup);
504 
505  // Only one process with dummy information. (Can overwrite at the end.)
506  int lprup = 9999;
507  double xsecup = 1.;
508  double xerrup = 0.;
509  double xmaxup = 1.;
510  addProcess(lprup, xsecup, xerrup, xmaxup);
511 
512  // Done.
513  return true;
514 
515 }
516 
517 //--------------------------------------------------------------------------
518 
519 // Read in event information from PYTHIA 8.
520 
521 bool LHAupFromPYTHIA8::setEvent( int ) {
522 
523  // Read process information from Info class, and store it.
524  // Note: renormalization scale here, factorization further down.
525  // For now always convert to process 9999, instead of infoPtr->code().
526  int idprup = 9999;
527  double xwgtup = infoPtr->weight();
528  double scalup = infoPtr->QRen();
529  double aqedup = infoPtr->alphaEM();
530  double aqcdup = infoPtr->alphaS();
531  setProcess(idprup, xwgtup, scalup, aqedup, aqcdup);
532 
533  // Read in particle info one by one, excluding zero and beams, and store it.
534  // Note unusual C++ loop range, to better reflect LHA/Fortran standard.
535  int nup = processPtr->size() - 3;
536  int idup, statusup, istup, mothup1, mothup2, icolup1, icolup2;
537  double pup1, pup2, pup3, pup4, pup5, vtimup, spinup;
538  for (int ip = 1; ip <= nup; ++ip) {
539  Particle& particle = (*processPtr)[ip + 2];
540  idup = particle.id();
541  // Convert from PYTHIA8 to LHA status codes.
542  statusup = particle.status();
543  if (ip < 3) istup = -1;
544  else if (statusup < 0) istup = 2;
545  else istup = 1;
546  mothup1 = max(0, particle.mother1() - 2);
547  mothup2 = max(0, particle.mother2() - 2);
548  icolup1 = particle.col();
549  icolup2 = particle.acol();
550  pup1 = particle.px();
551  pup2 = particle.py();
552  pup3 = particle.pz();
553  pup4 = particle.e();
554  pup5 = particle.m();
555  vtimup = particle.tau();
556  spinup = particle.pol();
557  addParticle(idup, istup, mothup1, mothup2, icolup1, icolup2,
558  pup1, pup2, pup3, pup4, pup5, vtimup, spinup) ;
559  }
560 
561  // Also extract pdf information from Info class, and store it.
562  int id1up = infoPtr->id1();
563  int id2up = infoPtr->id2();
564  double x1up = infoPtr->x1();
565  double x2up = infoPtr->x2();
566  double scalePDFup = infoPtr->QFac();
567  double xpdf1up = infoPtr->pdf1();
568  double xpdf2up = infoPtr->pdf2();
569  setPdf(id1up, id2up, x1up, x2up, scalePDFup, xpdf1up, xpdf2up, true);
570 
571  // Done.
572  return true;
573 
574 }
575 
576 //--------------------------------------------------------------------------
577 
578 // Update cross-section information at the end of the run.
579 
580 bool LHAupFromPYTHIA8::updateSigma() {
581 
582  // Read out information from PYTHIA 8 and send it in to LHA.
583  double sigGen = CONVERTMB2PB * infoPtr->sigmaGen();
584  double sigErr = CONVERTMB2PB * infoPtr->sigmaErr();
585  setXSec(0, sigGen);
586  setXErr(0, sigErr);
587 
588  // Done.
589  return true;
590 
591 }
592 
593 //==========================================================================
594 
595 } // end namespace Pythia8
Generated by   doxygen 1.8.5