PYTHIA8/pythia8130/src/LesHouches.cxx

   1 // LesHouches.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2008 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 // LHAupLHEF class.
 271
 272 //*********
 273
 274 // Read in initialization information from a Les Houches Event File.
 275
 276 bool LHAupLHEF::setInit() {
 277
 278   // Check that first line is consistent with proper LHEF file.
 279   string line;
 280   if (!getline(is, line)) return false;
 281   if (line.find("<LesHouchesEvents") == string::npos) return false;
 282   if (line.find("version=\"1.0\"" ) == string::npos ) return false;
 283
 284   // Loop over lines until an <init tag is found first on a line.
 285   string tag = " ";
 286   do {
 287     if (!getline(is, line)) return false;
 288     if (line.find_first_not_of(" ") != string::npos) {
 289       istringstream getfirst(line);
 290       getfirst >> tag;
 291       if (!getfirst) return false;
 292     }
 293   } while (tag != "<init>" && tag != "<init");
 294
 295   // Read in beam and strategy info, and store it.
 296   int idbmupA, idbmupB;
 297   double ebmupA, ebmupB;
 298   int pdfgupA, pdfgupB, pdfsupA, pdfsupB, idwtup, nprup;
 299   if (!getline(is, line)) return false;
 300   istringstream getbms(line);
 301   getbms >> idbmupA >> idbmupB >> ebmupA >> ebmupB >> pdfgupA
 302      >> pdfgupB >> pdfsupA >> pdfsupB >> idwtup >> nprup;
 303   if (!getbms) return false;
 304   setBeamA(idbmupA, ebmupA, pdfgupA, pdfsupA);
 305   setBeamB(idbmupB, ebmupB, pdfgupB, pdfsupB);
 306   setStrategy(idwtup);
 307
 308   // Read in process info, one process at a time, and store it.
 309   double xsecup, xerrup, xmaxup;
 310   int lprup;
 311   for (int ip = 0; ip < nprup; ++ip) {
 312     if (!getline(is, line)) return false;
 313     istringstream getpro(line);
 314     getpro >> xsecup >> xerrup >> xmaxup >> lprup ;
 315     if (!getpro) return false;
 316     addProcess(lprup, xsecup, xerrup, xmaxup);
 317   }
 318
 319   // Reading worked.
 320   return true;
 321
 322 }
 323
 324 //*********
 325
 326 // Read in event information from a Les Houches Event File.
 327
 328 bool LHAupLHEF::setEvent( int ) {
 329
 330   // Loop over lines until an <event tag is found first on a line.
 331   string line, tag;
 332   do {
 333     if (!getline(is, line)) return false;
 334     if (line.find_first_not_of(" ") != string::npos) {
 335       istringstream getfirst(line);
 336       getfirst >> tag;
 337       if (!getfirst) return false;
 338     }
 339   } while (tag != "<event>" && tag != "<event");
 340
 341   // Read in process info and store it.
 342   int nup, idprup;
 343   double xwgtup, scalup, aqedup, aqcdup;
 344   if (!getline(is, line)) return false;
 345   istringstream getpro(line);
 346   getpro >> nup >> idprup >> xwgtup >> scalup >> aqedup >> aqcdup;
 347   if (!getpro) return false;
 348   setProcess(idprup, xwgtup, scalup, aqedup, aqcdup);
 349
 350   // Read in particle info one by one, and store it.
 351   // Note unusual C++ loop range, to better reflect LHA/Fortran standard.
 352   // (Recall that process(...) above added empty particle at index 0.)
 353   int idup, istup, mothup1, mothup2, icolup1, icolup2;
 354   double pup1, pup2, pup3, pup4, pup5, vtimup, spinup;
 355   for (int ip = 1; ip <= nup; ++ip) {
 356     if (!getline(is, line)) return false;
 357     istringstream getall(line);
 358     getall >> idup >> istup >> mothup1 >> mothup2 >> icolup1 >> icolup2
 359       >> pup1 >> pup2 >> pup3 >> pup4 >> pup5 >> vtimup >> spinup;
 360     if (!getall) return false;
 361     addParticle(idup, istup, mothup1, mothup2, icolup1, icolup2,
 362       pup1, pup2, pup3, pup4, pup5, vtimup, spinup) ;
 363   }
 364
 365   // Continue parsing till </event>. Extract pdf info if present.
 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       int id1In, id2In;
 373       double x1In, x2In, scalePDFIn, xpdf1In, xpdf2In;
 374       getpdf >> id1In >> id2In >>  x1In >> x2In >> scalePDFIn
 375              >> xpdf1In >> xpdf2In;
 376       if (!getpdf) return false;
 377       setPdf(id1In, id2In, x1In, x2In, scalePDFIn, xpdf1In, xpdf2In);
 378     }
 379   } while (tag != "</event>" && tag != "</event");
 380
 381   // Reading worked.
 382   return true;
 383
 384 }
 385
 386 //**************************************************************************
 387
 388 // LHAupFromPYTHIA8 class.
 389
 390 //*********
 391
 392 // Read in initialization information from PYTHIA 8.
 393
 394 bool LHAupFromPYTHIA8::setInit() {
 395
 396   // Read in beam from Info class. Parton density left empty.
 397   int    idbmupA = infoPtr->idA();
 398   int    idbmupB = infoPtr->idB();
 399   double ebmupA  = infoPtr->eA();
 400   double ebmupB  = infoPtr->eB();
 401   int    pdfgupA = 0;
 402   int    pdfgupB = 0;
 403   int    pdfsupA = 0;
 404   int    pdfsupB = 0;
 405   setBeamA(idbmupA, ebmupA, pdfgupA, pdfsupA);
 406   setBeamB(idbmupB, ebmupB, pdfgupB, pdfsupB);
 407
 408   // Currently only one allowed strategy.
 409   int    idwtup = 3;
 410   setStrategy(idwtup);
 411
 412   // Only one process with dummy information. (Can overwrite at the end.)
 413   int    lprup  = 9999;
 414   double xsecup = 1.;
 415   double xerrup = 0.;
 416   double xmaxup = 1.;
 417   addProcess(lprup, xsecup, xerrup, xmaxup);
 418
 419   // Done.
 420   return true;
 421
 422 }
 423
 424 //*********
 425
 426 // Read in event information from PYTHIA 8.
 427
 428 bool LHAupFromPYTHIA8::setEvent( int ) {
 429
 430   // Read process information from Info class, and store it.
 431   // Note: renormalization scale here, factorization further down.
 432   int    idprup = infoPtr->code();
 433   // For now always convert to process 9999.
 434   idprup        = 9999;
 435   double xwgtup = infoPtr->weight();
 436   double scalup = infoPtr->QRen();
 437   double aqedup = infoPtr->alphaEM();
 438   double aqcdup = infoPtr->alphaS();
 439   setProcess(idprup, xwgtup, scalup, aqedup, aqcdup);
 440
 441   // Read in particle info one by one, excluding zero and beams, and store it.
 442   // Note unusual C++ loop range, to better reflect LHA/Fortran standard.
 443   int nup   = processPtr->size() - 3;
 444   int    idup, statusup, istup, mothup1, mothup2, icolup1, icolup2;
 445   double pup1, pup2, pup3, pup4, pup5, vtimup, spinup;
 446   for (int ip = 1; ip <= nup; ++ip) {
 447     Particle& particle = (*processPtr)[ip + 2];
 448     idup     = particle.id();
 449     // Convert from PYTHIA8 to LHA status codes.
 450     statusup = particle.status();
 451     if (ip < 3)            istup = -1;
 452     else if (statusup < 0) istup =  2;
 453     else                   istup =  1;
 454     mothup1  = max(0, particle.mother1() - 2);
 455     mothup2  = max(0, particle.mother2() - 2);
 456     icolup1  = particle.col();
 457     icolup2  = particle.acol();
 458     pup1     = particle.px();
 459     pup2     = particle.py();
 460     pup3     = particle.pz();
 461     pup4     = particle.e();
 462     pup5     = particle.m();
 463     vtimup   = particle.tau();
 464     spinup   = 9.;
 465     addParticle(idup, istup, mothup1, mothup2, icolup1, icolup2,
 466       pup1, pup2, pup3, pup4, pup5, vtimup, spinup) ;
 467   }
 468
 469   // Also extract pdf information from Info class, and store it.
 470   int    id1up      = infoPtr->id1();
 471   int    id2up      = infoPtr->id2();
 472   double x1up       = infoPtr->x1();
 473   double x2up       = infoPtr->x2();
 474   double scalePDFup = infoPtr->QFac();
 475   double xpdf1up    = infoPtr->pdf1();
 476   double xpdf2up    = infoPtr->pdf2();
 477   setPdf(id1up, id2up, x1up, x2up, scalePDFup, xpdf1up, xpdf2up);
 478
 479   // Done.
 480   return true;
 481
 482 }
 483
 484 //*********
 485
 486 //  Update cross-section information at the end of the run.
 487
 488 bool LHAupFromPYTHIA8::updateSigma() {
 489
 490   // Read out information from PYTHIA 8 and send it in to LHA.
 491   double sigGen = CONVERTMB2PB * infoPtr->sigmaGen();
 492   double sigErr = CONVERTMB2PB * infoPtr->sigmaErr();
 493   setXSec(0, sigGen);
 494   setXErr(0, sigErr);
 495
 496   // Done.
 497   return true;
 498
 499 }
 500
 501 //**************************************************************************
 502
 503 } // end namespace Pythia8