PYTHIA8/pythia8145/src/LesHouches.cxx

   1 // LesHouches.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2010 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 namespace Pythia8 {
  15
  16 //==========================================================================
  17
  18 // LHAup class.
  19
  20 //--------------------------------------------------------------------------
  21
  22 // Constants: could be changed here if desired, but normally should not.
  23 // These are of technical nature, as described for each.
  24
  25 // LHA convention with cross section in pb may require conversion from mb.
  26 const double LHAup::CONVERTMB2PB = 1e9;
  27
  28 //--------------------------------------------------------------------------
  29
  30 // Print the initialization info; to check it worked.
  31
  32 void LHAup::listInit(ostream& os) {
  33
  34   // Header.
  35   os << "\n --------  LHA initialization information  ------------ \n";
  36
  37   // Beam info.
  38   os << fixed << setprecision(3)
  39      << "\n  beam    kind      energy  pdfgrp  pdfset \n"
  40      << "     A  " << setw(6) << idBeamASave
  41      <<  setw(12) << eBeamASave
  42      << setw(8) << pdfGroupBeamASave
  43      << setw(8) << pdfSetBeamASave << "\n"
  44      << "     B  " << setw(6) << idBeamBSave
  45      <<  setw(12) << eBeamBSave
  46      << setw(8) << pdfGroupBeamBSave
  47      << setw(8) << pdfSetBeamBSave << "\n";
  48
  49   // Event weighting strategy.
  50   os << "\n  Event weighting strategy = " << setw(2)
  51      << strategySave << "\n" ;
  52
  53   // Process list.
  54   os << scientific << setprecision(4)
  55      << "\n  Processes, with strategy-dependent cross section info \n"
  56      << "  number      xsec (pb)      xerr (pb)      xmax (pb) \n" ;
  57   for (int ip = 0; ip < int(processes.size()); ++ip) {
  58     os << setw(8) << processes[ip].idProc
  59        << setw(15) << processes[ip].xSecProc
  60        << setw(15) << processes[ip].xErrProc
  61        << setw(15) << processes[ip].xMaxProc << "\n";
  62   }
  63
  64   // Finished.
  65   os << "\n --------  End LHA initialization information  -------- \n";
  66
  67 }
  68
  69 //--------------------------------------------------------------------------
  70
  71 // Print the event info; to check it worked.
  72
  73 void LHAup::listEvent(ostream& os) {
  74
  75   // Header.
  76   os << "\n --------  LHA event information and listing  -------------"
  77      << "--------------------------------------------------------- \n";
  78
  79   // Basic event info.
  80   os << scientific << setprecision(4)
  81      << "\n    process = " << setw(8) << idProc
  82      << "    weight = " << setw(12) << weightProc
  83      << "     scale = " << setw(12) << scaleProc << " (GeV) \n"
  84      << "                   "
  85      << "     alpha_em = " << setw(12) << alphaQEDProc
  86      << "    alpha_strong = " << setw(12) << alphaQCDProc << "\n";
  87
  88   // Particle list
  89   os << fixed << setprecision(3)
  90      << "\n    Participating Particles \n"
  91      << "    no        id stat     mothers     colours      p_x        "
  92      << "p_y        p_z         e          m        tau    spin \n" ;
  93   for (int ip = 1; ip < int(particles.size()); ++ip) {
  94     os << setw(6) << ip
  95        << setw(10) << particles[ip].idPart
  96        << setw(5) << particles[ip].statusPart
  97        << setw(6) << particles[ip].mother1Part
  98        << setw(6) << particles[ip].mother2Part
  99        << setw(6) << particles[ip].col1Part
 100        << setw(6) << particles[ip].col2Part
 101        << setw(11) << particles[ip].pxPart
 102        << setw(11) << particles[ip].pyPart
 103        << setw(11) << particles[ip].pzPart
 104        << setw(11) << particles[ip].ePart
 105        << setw(11) <<  particles[ip].mPart
 106        << setw(8) <<  particles[ip].tauPart
 107        << setw(8) <<  particles[ip].spinPart << "\n";
 108   }
 109
 110   // PDF info - optional.
 111   if (pdfIsSetSave) os << "\n   pdf: id1 =" << setw(5) << id1Save
 112     << " id2 =" << setw(5) << id2Save
 113     << " x1 ="  << scientific << setw(10) << x1Save
 114     << " x2 =" << setw(10) << x2Save
 115     << " scalePDF =" << setw(10) << scalePDFSave
 116     << " xpdf1 =" << setw(10) << xpdf1Save
 117     << " xpdf2 =" << setw(10) << xpdf2Save << "\n";
 118
 119   // Finished.
 120   os << "\n --------  End LHA event information and listing  ---------"
 121      << "--------------------------------------------------------- \n";
 122
 123 }
 124
 125 //--------------------------------------------------------------------------
 126
 127 // Open and write header to a Les Houches Event File.
 128
 129 bool LHAup::openLHEF(string fileNameIn) {
 130
 131   // Open file for writing. Reset it to be empty.
 132   fileName = fileNameIn;
 133   const char* cstring = fileName.c_str();
 134   osLHEF.open(cstring, ios::out | ios::trunc);
 135   if (!osLHEF) {
 136     infoPtr->errorMsg("Error in LHAup::openLHEF:"
 137       " could not open file", fileName);
 138     return false;
 139   }
 140
 141   // Read out current date and time.
 142   time_t t = time(0);
 143   strftime(dateNow,12,"%d %b %Y",localtime(&t));
 144   strftime(timeNow,9,"%H:%M:%S",localtime(&t));
 145
 146   // Write header.
 147   osLHEF << "<LesHouchesEvents version=\"1.0\">\n"
 148          << "<!--\n"
 149          << "  File written by Pythia8::LHAup on "
 150          << dateNow << " at " << timeNow << "\n"
 151          << "-->" << endl;
 152
 153   // Done.
 154   return true;
 155
 156 }
 157
 158 //--------------------------------------------------------------------------
 159
 160 // Write initialization information to a Les Houches Event File.
 161
 162 bool LHAup::initLHEF() {
 163
 164   // Write information on beams.
 165   osLHEF << "<init>\n" << scientific << setprecision(6)
 166          << "  " << idBeamASave       << "  " << idBeamBSave
 167          << "  " << eBeamASave        << "  " << eBeamBSave
 168          << "  " << pdfGroupBeamASave << "  " << pdfGroupBeamBSave
 169          << "  " << pdfSetBeamASave   << "  " << pdfSetBeamBSave
 170          << "  " << strategySave      << "  " << processes.size() << "\n";
 171
 172   // Write information on all the subprocesses.
 173   for (int ip = 0; ip < int(processes.size()); ++ip)
 174     osLHEF << " " << setw(13) << processes[ip].xSecProc
 175            << " " << setw(13) << processes[ip].xErrProc
 176            << " " << setw(13) << processes[ip].xMaxProc
 177            << " " << setw(6) << processes[ip].idProc << "\n";
 178
 179   // Done.
 180   osLHEF << "</init>" << endl;
 181   return true;
 182
 183 }
 184
 185 //--------------------------------------------------------------------------
 186
 187 // Write event information to a Les Houches Event File.
 188
 189 bool LHAup::eventLHEF() {
 190
 191   // Write information on process as such.
 192   osLHEF << "<event>\n" << scientific << setprecision(6)
 193          << " " << setw(5) << particles.size() - 1
 194          << " " << setw(5) << idProc
 195          << " " << setw(13) << weightProc
 196          << " " << setw(13) << scaleProc
 197          << " " << setw(13) << alphaQEDProc
 198          << " " << setw(13) << alphaQCDProc << "\n";
 199
 200   // Write information on the particles, excluding zeroth.
 201   for (int ip = 1; ip < int(particles.size()); ++ip) {
 202     osLHEF << " " << setw(8) << particles[ip].idPart
 203            << " " << setw(5) << particles[ip].statusPart
 204            << " " << setw(5) << particles[ip].mother1Part
 205            << " " << setw(5) << particles[ip].mother2Part
 206            << " " << setw(5) << particles[ip].col1Part
 207            << " " << setw(5) << particles[ip].col2Part << setprecision(10)
 208            << " " << setw(17) << particles[ip].pxPart
 209            << " " << setw(17) << particles[ip].pyPart
 210            << " " << setw(17) << particles[ip].pzPart
 211            << " " << setw(17) << particles[ip].ePart
 212            << " " << setw(17) <<  particles[ip].mPart << setprecision(6);
 213     if (particles[ip].tauPart == 0.) osLHEF << " 0.";
 214     else osLHEF << " " << setw(13) << particles[ip].tauPart;
 215     if (particles[ip].spinPart == 9.) osLHEF << " 9.";
 216     else osLHEF << " " << setw(13) << particles[ip].spinPart;
 217     osLHEF << "\n";
 218   }
 219
 220   // Optionally write information on PDF values at hard interaction.
 221   if (pdfIsSetSave) osLHEF << "#pdf"
 222            << " " << setw(4) << id1Save
 223            << " " << setw(4) << id2Save
 224            << " " << setw(13) << x1Save
 225            << " " << setw(13) << x2Save
 226            << " " << setw(13) << scalePDFSave
 227            << " " << setw(13) << xpdf1Save
 228            << " " << setw(13) << xpdf2Save << "\n";
 229
 230   // Done.
 231   osLHEF << "</event>" << endl;
 232   return true;
 233
 234 }
 235
 236 //--------------------------------------------------------------------------
 237
 238 // Write end of a Les Houches Event File and close it.
 239
 240 bool LHAup::closeLHEF(bool updateInit) {
 241
 242   // Write an end to the file.
 243   osLHEF << "</LesHouchesEvents>" << endl;
 244   osLHEF.close();
 245
 246   // Optionally update the cross section information.
 247   if (updateInit) {
 248     const char* cstring = fileName.c_str();
 249     osLHEF.open(cstring, ios::in | ios::out);
 250
 251     // Rewrite header; identically with what openLHEF did.
 252     osLHEF << "<LesHouchesEvents version=\"1.0\">\n"
 253            << "<!--\n"
 254            << "  File written by Pythia8::LHAup on "
 255            << dateNow << " at " << timeNow << "\n"
 256            << "-->" << endl;
 257
 258     // Redo initialization information.
 259     initLHEF();
 260     osLHEF.close();
 261   }
 262
 263   // Done.
 264   return true;
 265
 266 }
 267
 268 //--------------------------------------------------------------------------
 269
 270 // Read in initialization information from a Les Houches Event File.
 271
 272 bool LHAup::setInitLHEF(ifstream& is) {
 273
 274   // Check that first line is consistent with proper LHEF file.
 275   string line;
 276   if (!getline(is, line)) return false;
 277   if (line.find("<LesHouchesEvents") == string::npos) return false;
 278   if (line.find("version=\"1.0\"" ) == string::npos ) return false;
 279
 280   // Loop over lines until an <init tag is found first on a line.
 281   string tag = " ";
 282   do {
 283     if (!getline(is, line)) return false;
 284     if (line.find_first_not_of(" \n\t\v\b\r\f\a") != string::npos) {
 285       istringstream getfirst(line);
 286       getfirst >> tag;
 287       if (!getfirst) return false;
 288     }
 289   } while (tag != "<init>" && tag != "<init");
 290
 291   // Read in beam and strategy info, and store it.
 292   int idbmupA, idbmupB;
 293   double ebmupA, ebmupB;
 294   int pdfgupA, pdfgupB, pdfsupA, pdfsupB, idwtup, nprup;
 295   if (!getline(is, line)) return false;
 296   istringstream getbms(line);
 297   getbms >> idbmupA >> idbmupB >> ebmupA >> ebmupB >> pdfgupA
 298      >> pdfgupB >> pdfsupA >> pdfsupB >> idwtup >> nprup;
 299   if (!getbms) return false;
 300   setBeamA(idbmupA, ebmupA, pdfgupA, pdfsupA);
 301   setBeamB(idbmupB, ebmupB, pdfgupB, pdfsupB);
 302   setStrategy(idwtup);
 303
 304   // Read in process info, one process at a time, and store it.
 305   double xsecup, xerrup, xmaxup;
 306   int lprup;
 307   for (int ip = 0; ip < nprup; ++ip) {
 308     if (!getline(is, line)) return false;
 309     istringstream getpro(line);
 310     getpro >> xsecup >> xerrup >> xmaxup >> lprup ;
 311     if (!getpro) return false;
 312     addProcess(lprup, xsecup, xerrup, xmaxup);
 313   }
 314
 315   // Reading worked.
 316   return true;
 317
 318 }
 319
 320 //--------------------------------------------------------------------------
 321
 322 // Read in event information from a Les Houches Event File,
 323 // into a staging area where it can be reused by setOldEventLHEF.
 324
 325 bool LHAup::setNewEventLHEF(ifstream& is) {
 326
 327   // Loop over lines until an <event tag is found first on a line.
 328   string line, tag;
 329   do {
 330     if (!getline(is, line)) return false;
 331     if (line.find_first_not_of(" \n\t\v\b\r\f\a") != string::npos) {
 332       istringstream getfirst(line);
 333       getfirst >> tag;
 334       if (!getfirst) return false;
 335     }
 336   } while (tag != "<event>" && tag != "<event");
 337
 338   // Read in process info and store it.
 339   if (!getline(is, line)) return false;
 340   istringstream getpro(line);
 341   getpro >> nupSave >> idprupSave >> xwgtupSave >> scalupSave
 342     >> aqedupSave >> aqcdupSave;
 343   if (!getpro) return false;
 344
 345   // Reset particlesSave vector, add slot-0 empty particle.
 346   particlesSave.clear();
 347   particlesSave.push_back( LHAParticle() );
 348
 349   // Read in particle info one by one, and store it.
 350   // Note unusual C++ loop range, to better reflect LHA/Fortran standard.
 351   // (Recall that process(...) above added empty particle at index 0.)
 352   int idup, istup, mothup1, mothup2, icolup1, icolup2;
 353   double pup1, pup2, pup3, pup4, pup5, vtimup, spinup;
 354   for (int ip = 1; ip <= nupSave; ++ip) {
 355     if (!getline(is, line)) return false;
 356     istringstream getall(line);
 357     getall >> idup >> istup >> mothup1 >> mothup2 >> icolup1 >> icolup2
 358       >> pup1 >> pup2 >> pup3 >> pup4 >> pup5 >> vtimup >> spinup;
 359     if (!getall) return false;
 360     particlesSave.push_back( LHAParticle( idup, istup, mothup1, mothup2,
 361       icolup1, icolup2, pup1, pup2, pup3, pup4, pup5, vtimup, spinup) );
 362   }
 363
 364   // Continue parsing till </event>. Extract pdf info if present.
 365   getPDFSave = false;
 366   do {
 367     if (!getline(is, line)) return false;
 368     istringstream getpdf(line);
 369     getpdf >> tag;
 370     if (!getpdf) return false;
 371     if (tag == "#pdf") {
 372       getpdf >> id1InSave >> id2InSave >> x1InSave >> x2InSave
 373              >> scalePDFInSave >> xpdf1InSave >> xpdf2InSave;
 374       if (!getpdf) return false;
 375       getPDFSave = true;
 376     }
 377   } while (tag != "</event>" && tag != "</event");
 378
 379   // Need id and x values even when no PDF info. Rest empty.
 380   if (!getPDFSave) {
 381     id1InSave      = particlesSave[1].idPart;
 382     id2InSave      = particlesSave[2].idPart;
 383     x1InSave       = particlesSave[1].ePart / eBeamASave;
 384     x2InSave       = particlesSave[2].ePart / eBeamBSave;
 385     scalePDFInSave = 0.;
 386     xpdf1InSave    = 0.;
 387     xpdf2InSave    = 0.;
 388   }
 389
 390   // Reading worked.
 391   return true;
 392
 393 }
 394
 395 //--------------------------------------------------------------------------
 396
 397 // Make current event information read in by setNewEventLHEF.
 398
 399 bool LHAup::setOldEventLHEF() {
 400
 401   // Store saved event, optionally also parton density information.
 402   setProcess(idprupSave, xwgtupSave, scalupSave, aqedupSave, aqcdupSave);
 403   for (int ip = 1; ip <= nupSave; ++ip) addParticle( particlesSave[ip] );
 404   setPdf(id1InSave, id2InSave, x1InSave, x2InSave, scalePDFInSave,
 405     xpdf1InSave, xpdf2InSave, getPDFSave);
 406
 407   // Done;
 408   return true;
 409
 410 }
 411
 412 //==========================================================================
 413
 414 // LHAupFromPYTHIA8 class.
 415
 416 //--------------------------------------------------------------------------
 417
 418 // Read in initialization information from PYTHIA 8.
 419
 420 bool LHAupFromPYTHIA8::setInit() {
 421
 422   // Read in beam from Info class. Parton density left empty.
 423   int    idbmupA = infoPtr->idA();
 424   int    idbmupB = infoPtr->idB();
 425   double ebmupA  = infoPtr->eA();
 426   double ebmupB  = infoPtr->eB();
 427   int    pdfgupA = 0;
 428   int    pdfgupB = 0;
 429   int    pdfsupA = 0;
 430   int    pdfsupB = 0;
 431   setBeamA(idbmupA, ebmupA, pdfgupA, pdfsupA);
 432   setBeamB(idbmupB, ebmupB, pdfgupB, pdfsupB);
 433
 434   // Currently only one allowed strategy.
 435   int    idwtup = 3;
 436   setStrategy(idwtup);
 437
 438   // Only one process with dummy information. (Can overwrite at the end.)
 439   int    lprup  = 9999;
 440   double xsecup = 1.;
 441   double xerrup = 0.;
 442   double xmaxup = 1.;
 443   addProcess(lprup, xsecup, xerrup, xmaxup);
 444
 445   // Done.
 446   return true;
 447
 448 }
 449
 450 //--------------------------------------------------------------------------
 451
 452 // Read in event information from PYTHIA 8.
 453
 454 bool LHAupFromPYTHIA8::setEvent( int ) {
 455
 456   // Read process information from Info class, and store it.
 457   // Note: renormalization scale here, factorization further down.
 458   int    idprup = infoPtr->code();
 459   // For now always convert to process 9999.
 460   idprup        = 9999;
 461   double xwgtup = infoPtr->weight();
 462   double scalup = infoPtr->QRen();
 463   double aqedup = infoPtr->alphaEM();
 464   double aqcdup = infoPtr->alphaS();
 465   setProcess(idprup, xwgtup, scalup, aqedup, aqcdup);
 466
 467   // Read in particle info one by one, excluding zero and beams, and store it.
 468   // Note unusual C++ loop range, to better reflect LHA/Fortran standard.
 469   int nup   = processPtr->size() - 3;
 470   int    idup, statusup, istup, mothup1, mothup2, icolup1, icolup2;
 471   double pup1, pup2, pup3, pup4, pup5, vtimup, spinup;
 472   for (int ip = 1; ip <= nup; ++ip) {
 473     Particle& particle = (*processPtr)[ip + 2];
 474     idup     = particle.id();
 475     // Convert from PYTHIA8 to LHA status codes.
 476     statusup = particle.status();
 477     if (ip < 3)            istup = -1;
 478     else if (statusup < 0) istup =  2;
 479     else                   istup =  1;
 480     mothup1  = max(0, particle.mother1() - 2);
 481     mothup2  = max(0, particle.mother2() - 2);
 482     icolup1  = particle.col();
 483     icolup2  = particle.acol();
 484     pup1     = particle.px();
 485     pup2     = particle.py();
 486     pup3     = particle.pz();
 487     pup4     = particle.e();
 488     pup5     = particle.m();
 489     vtimup   = particle.tau();
 490     spinup   = 9.;
 491     addParticle(idup, istup, mothup1, mothup2, icolup1, icolup2,
 492       pup1, pup2, pup3, pup4, pup5, vtimup, spinup) ;
 493   }
 494
 495   // Also extract pdf information from Info class, and store it.
 496   int    id1up      = infoPtr->id1();
 497   int    id2up      = infoPtr->id2();
 498   double x1up       = infoPtr->x1();
 499   double x2up       = infoPtr->x2();
 500   double scalePDFup = infoPtr->QFac();
 501   double xpdf1up    = infoPtr->pdf1();
 502   double xpdf2up    = infoPtr->pdf2();
 503   setPdf(id1up, id2up, x1up, x2up, scalePDFup, xpdf1up, xpdf2up, true);
 504
 505   // Done.
 506   return true;
 507
 508 }
 509
 510 //--------------------------------------------------------------------------
 511
 512 //  Update cross-section information at the end of the run.
 513
 514 bool LHAupFromPYTHIA8::updateSigma() {
 515
 516   // Read out information from PYTHIA 8 and send it in to LHA.
 517   double sigGen = CONVERTMB2PB * infoPtr->sigmaGen();
 518   double sigErr = CONVERTMB2PB * infoPtr->sigmaErr();
 519   setXSec(0, sigGen);
 520   setXErr(0, sigErr);
 521
 522   // Done.
 523   return true;
 524
 525 }
 526
 527 //==========================================================================
 528
 529 } // end namespace Pythia8