PYTHIA8/pythia8170/include/LesHouches.h

   1 // LesHouches.h 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 // Header file for Les Houches Accord user process information.
   7 // LHAProcess: stores a single process; used by the other classes.
   8 // LHAParticle: stores a single particle; used by the other classes.
   9 // LHAup: base class for initialization and event information.
  10 // LHAupLHEF: derived class for reading from an Les Houches Event File.
  11 // Code for interfacing with Fortran commonblocks is found in LHAFortran.h.
  12
  13 #ifndef Pythia8_LesHouches_H
  14 #define Pythia8_LesHouches_H
  15
  16 #include "Event.h"
  17 #include "Info.h"
  18 #include "PythiaStdlib.h"
  19 #include "Settings.h"
  20
  21 namespace Pythia8 {
  22
  23 //==========================================================================
  24
  25 // A class for the processes stored in LHAup.
  26
  27 class LHAProcess {
  28
  29 public:
  30
  31   // Constructors.
  32   LHAProcess() : idProc(0), xSecProc(0.), xErrProc(0.), xMaxProc(0.) { }
  33   LHAProcess(int idProcIn, double xSecIn, double xErrIn, double xMaxIn) :
  34     idProc(idProcIn), xSecProc(xSecIn), xErrProc(xErrIn),
  35     xMaxProc(xMaxIn) { }
  36
  37   // Process properties.
  38   int idProc;
  39   double xSecProc, xErrProc, xMaxProc;
  40
  41 } ;
  42
  43 //==========================================================================
  44
  45 // A class for the particles stored in LHAup.
  46
  47 class LHAParticle {
  48
  49 public:
  50
  51   // Constructors.
  52   LHAParticle() : idPart(0), statusPart(0), mother1Part(0),
  53     mother2Part(0), col1Part(0), col2Part(0), pxPart(0.), pyPart(0.),
  54     pzPart(0.), ePart(0.), mPart(0.), tauPart(0.), spinPart(9.) { }
  55   LHAParticle(int idIn, int statusIn, int mother1In, int mother2In,
  56     int col1In, int col2In, double pxIn, double pyIn, double pzIn,
  57     double eIn, double mIn, double tauIn, double spinIn) :
  58     idPart(idIn), statusPart(statusIn), mother1Part(mother1In),
  59     mother2Part(mother2In), col1Part(col1In), col2Part(col2In),
  60     pxPart(pxIn), pyPart(pyIn), pzPart(pzIn), ePart(eIn), mPart(mIn),
  61     tauPart(tauIn), spinPart(spinIn) { }
  62
  63   // Particle properties.
  64   int idPart, statusPart, mother1Part, mother2Part, col1Part, col2Part;
  65   double pxPart, pyPart, pzPart, ePart, mPart, tauPart, spinPart;
  66
  67 } ;
  68
  69 //==========================================================================
  70
  71 // LHAup is base class for initialization and event information
  72 // from an external parton-level generator.
  73
  74 class LHAup {
  75
  76 public:
  77
  78   // Destructor.
  79   virtual ~LHAup() {}
  80
  81   // Set info pointer.
  82   void setPtr(Info* infoPtrIn) {infoPtr = infoPtrIn;}
  83
  84   // Method to be used for LHAupLHEF derived class.
  85   virtual void newEventFile(const char*) {}
  86   virtual bool fileFound() {return true;}
  87
  88   // A pure virtual method setInit, wherein all initialization information
  89   // is supposed to be set in the derived class. Can do this by reading a
  90   // file or some other way, as desired. Returns false if it did not work.
  91   virtual bool setInit() = 0;
  92
  93   // Give back info on beams.
  94   int    idBeamA()       const {return idBeamASave;}
  95   int    idBeamB()       const {return idBeamBSave;}
  96   double eBeamA()        const {return eBeamASave;}
  97   double eBeamB()        const {return eBeamBSave;}
  98   int    pdfGroupBeamA() const {return pdfGroupBeamASave;}
  99   int    pdfGroupBeamB() const {return pdfGroupBeamBSave;}
 100   int    pdfSetBeamA()   const {return pdfSetBeamASave;}
 101   int    pdfSetBeamB()   const {return pdfSetBeamBSave;}
 102
 103   // Give back weight strategy.
 104   int    strategy()      const {return strategySave;}
 105
 106   // Give back info on processes.
 107   int    sizeProc()      const {return processes.size();}
 108   int    idProcess(int proc) const {return processes[proc].idProc;}
 109   double xSec(int proc)  const {return processes[proc].xSecProc;}
 110   double xErr(int proc)  const {return processes[proc].xErrProc;}
 111   double xMax(int proc)  const {return processes[proc].xMaxProc;}
 112   double xSecSum()       const {return xSecSumSave;}
 113   double xErrSum()       const {return xErrSumSave;}
 114
 115   // Print the initialization info; useful to check that setting it worked.
 116   void   listInit(ostream& os = cout);
 117
 118   // A pure virtual method setEvent, wherein information on the next event
 119   // is supposed to be set in the derived class.
 120   // Strategies +-1 and +-2: idProcIn is the process type, selected by PYTHIA.
 121   // Strategies +-3 and +-4: idProcIn is dummy; process choice is made locally.
 122   // The method can find the next event by a runtime interface to another
 123   // program, or by reading a file, as desired.
 124   // The method should return false if it did not work.
 125   virtual bool setEvent(int idProcIn = 0) = 0;
 126
 127   // Give back process number, weight, scale, alpha_em, alpha_s.
 128   int    idProcess()       const {return idProc;}
 129   double weight()          const {return weightProc;}
 130   double scale()           const {return scaleProc;}
 131   double alphaQED()        const {return alphaQEDProc;}
 132   double alphaQCD()        const {return alphaQCDProc;}
 133
 134   // Give back info on separate particle.
 135   int    sizePart()        const {return particles.size();}
 136   int    id(int part)      const {return particles[part].idPart;}
 137   int    status(int part)  const {return particles[part].statusPart;}
 138   int    mother1(int part) const {return particles[part].mother1Part;}
 139   int    mother2(int part) const {return particles[part].mother2Part;}
 140   int    col1(int part)    const {return particles[part].col1Part;}
 141   int    col2(int part)    const {return particles[part].col2Part;}
 142   double px(int part)      const {return particles[part].pxPart;}
 143   double py(int part)      const {return particles[part].pyPart;}
 144   double pz(int part)      const {return particles[part].pzPart;}
 145   double e(int part)       const {return particles[part].ePart;}
 146   double m(int part)       const {return particles[part].mPart;}
 147   double tau(int part)     const {return particles[part].tauPart;}
 148   double spin(int part)    const {return particles[part].spinPart;}
 149
 150   // Give back info on flavour and x values of hard-process initiators.
 151   int    id1()             const {return id1Save;}
 152   int    id2()             const {return id2Save;}
 153   double x1()              const {return x1Save;}
 154   double x2()              const {return x2Save;}
 155
 156   // Optional: give back info on parton density values of event.
 157   bool   pdfIsSet()        const {return pdfIsSetSave;}
 158   int    id1pdf()          const {return id1pdfSave;}
 159   int    id2pdf()          const {return id2pdfSave;}
 160   double x1pdf()           const {return x1pdfSave;}
 161   double x2pdf()           const {return x2pdfSave;}
 162   double scalePDF()        const {return scalePDFSave;}
 163   double pdf1()            const {return pdf1Save;}
 164   double pdf2()            const {return pdf2Save;}
 165
 166   // Print the info; useful to check that reading an event worked.
 167   void   listEvent(ostream& os = cout);
 168
 169   // Skip ahead a number of events, which are not considered further.
 170   // Mainly intended for debug when using the LHAupLHEF class.
 171   virtual bool skipEvent(int nSkip) {
 172     for (int iSkip = 0; iSkip < nSkip; ++iSkip)
 173     if (!setEvent()) return false; return true;}
 174
 175   // Four routines to write a Les Houches Event file in steps.
 176   bool   openLHEF(string fileNameIn);
 177   bool   initLHEF();
 178   bool   eventLHEF(bool verbose = true);
 179   bool   closeLHEF(bool updateInit = false);
 180
 181 protected:
 182
 183   // Constructor. Sets default to be that events come with unit weight.
 184   LHAup(int strategyIn = 3) : strategySave(strategyIn)
 185     { processes.reserve(10); particles.reserve(20); }
 186
 187   // Allow conversion from mb to pb.
 188   static const double CONVERTMB2PB;
 189
 190   // Pointer to various information on the generation.
 191   Info* infoPtr;
 192
 193   // Input beam info.
 194   void setBeamA(int idIn, double eIn, int pdfGroupIn = 0, int pdfSetIn = 0)
 195     { idBeamASave = idIn; eBeamASave = eIn; pdfGroupBeamASave = pdfGroupIn;
 196     pdfSetBeamASave = pdfSetIn;}
 197   void setBeamB(int idIn, double eIn, int pdfGroupIn = 0, int pdfSetIn = 0)
 198     { idBeamBSave = idIn; eBeamBSave = eIn; pdfGroupBeamBSave = pdfGroupIn;
 199     pdfSetBeamBSave = pdfSetIn;}
 200
 201   // Input process weight strategy.
 202   void setStrategy(int strategyIn) {strategySave = strategyIn;}
 203
 204   // Input process info.
 205   void addProcess(int idProcIn, double xSecIn = 1., double xErrIn = 0.,
 206     double xMaxIn = 1.) { processes.push_back( LHAProcess( idProcIn,
 207     xSecIn, xErrIn, xMaxIn)); }
 208
 209   // Possibility to update some cross section info at end of run.
 210   void setXSec(int iP, double xSecIn) {processes[iP].xSecProc = xSecIn;}
 211   void setXErr(int iP, double xErrIn) {processes[iP].xErrProc = xErrIn;}
 212   void setXMax(int iP, double xMaxIn) {processes[iP].xMaxProc = xMaxIn;}
 213
 214   // Input info on the selected process.
 215   void setProcess(int idProcIn = 0, double weightIn = 1., double
 216     scaleIn = 0., double alphaQEDIn = 0.0073, double alphaQCDIn = 0.12) {
 217     idProc = idProcIn; weightProc = weightIn; scaleProc = scaleIn;
 218     alphaQEDProc = alphaQEDIn; alphaQCDProc = alphaQCDIn;
 219     // Clear particle list. Add empty zeroth particle for correct indices.
 220     particles.clear(); addParticle(0); pdfIsSetSave = false;}
 221
 222   // Input particle info, one particle at the time.
 223   void addParticle(LHAParticle particleIn) {
 224     particles.push_back(particleIn);}
 225   void addParticle(int idIn, int statusIn = 0, int mother1In = 0,
 226     int mother2In = 0, int col1In = 0, int col2In = 0, double pxIn = 0.,
 227     double pyIn = 0., double pzIn = 0., double eIn = 0., double mIn = 0.,
 228     double tauIn = 0., double spinIn = 9.) {
 229     particles.push_back( LHAParticle( idIn, statusIn, mother1In, mother2In,
 230     col1In, col2In, pxIn, pyIn, pzIn, eIn, mIn, tauIn, spinIn)); }
 231
 232   // Input info on flavour and x values of hard-process initiators.
 233   void setIdX(int id1In, int id2In, double x1In, double x2In)
 234     { id1Save = id1In; id2Save = id2In; x1Save = x1In; x2Save = x2In;}
 235
 236   // Optionally input info on parton density values of event.
 237   void setPdf(int id1pdfIn, int id2pdfIn, double x1pdfIn, double x2pdfIn,
 238     double scalePDFIn, double pdf1In, double pdf2In, bool pdfIsSetIn)
 239     { id1pdfSave = id1pdfIn; id2pdfSave = id2pdfIn; x1pdfSave = x1pdfIn;
 240     x2pdfSave = x2pdfIn; scalePDFSave = scalePDFIn; pdf1Save = pdf1In;
 241     pdf2Save = pdf2In; pdfIsSetSave = pdfIsSetIn;}
 242
 243   // Three routines for LHEF files, but put here for flexibility.
 244   bool setInitLHEF(istream& is, bool readHeaders = false);
 245   bool setNewEventLHEF(istream& is);
 246   bool setOldEventLHEF();
 247
 248   // Helper routines to open and close a file handling GZIPSUPPORT:
 249   //   ifstream ifs;
 250   //   istream *is = openFile("myFile.txt", ifs);
 251   //   -- Process file using is --
 252   //   closeFile(is, ifs);
 253   istream* openFile(const char *fn, ifstream &ifs);
 254   void     closeFile(istream *&is, ifstream &ifs);
 255
 256   // LHAup is a friend class to infoPtr, but derived classes
 257   // are not. This wrapper function can be used by derived classes
 258   // to set headers in the Info class.
 259   void setInfoHeader(const string &key, const string &val) {
 260     infoPtr->setHeader(key, val); }
 261
 262   // Event properties from LHEF files, for repeated use.
 263   int    nupSave, idprupSave;
 264   double xwgtupSave, scalupSave, aqedupSave, aqcdupSave, xSecSumSave,
 265          xErrSumSave;
 266   vector<LHAParticle> particlesSave;
 267   bool   getPDFSave;
 268   int    id1InSave, id2InSave, id1pdfInSave, id2pdfInSave;
 269   double x1InSave, x2InSave, x1pdfInSave, x2pdfInSave, scalePDFInSave,
 270          pdf1InSave, pdf2InSave;
 271
 272 private:
 273
 274   // Event weighting and mixing strategy.
 275   int strategySave;
 276
 277   // Beam particle properties.
 278   int    idBeamASave, idBeamBSave;
 279   double eBeamASave, eBeamBSave;
 280   int    pdfGroupBeamASave, pdfGroupBeamBSave, pdfSetBeamASave, pdfSetBeamBSave;
 281
 282   // The process list, stored as a vector of processes.
 283   vector<LHAProcess> processes;
 284
 285   // Store info on the selected process.
 286   int    idProc;
 287   double weightProc, scaleProc, alphaQEDProc, alphaQCDProc;
 288
 289   // The particle list, stored as a vector of particles.
 290   vector<LHAParticle> particles;
 291
 292   // Info on initiators and optionally on parton density values of event.
 293   bool   pdfIsSetSave;
 294   int    id1Save, id2Save, id1pdfSave, id2pdfSave;
 295   double x1Save, x2Save, x1pdfSave, x2pdfSave, scalePDFSave, pdf1Save,
 296          pdf2Save;
 297
 298   // File to which to write Les Houches Event File information.
 299   string fileName;
 300   fstream osLHEF;
 301   char dateNow[12];
 302   char timeNow[9];
 303
 304 };
 305
 306 //==========================================================================
 307
 308 // A derived class with information read from a Les Houches Event File.
 309
 310 class LHAupLHEF : public LHAup {
 311
 312 public:
 313
 314   // Constructor.
 315   LHAupLHEF(const char* fileIn, const char* headerIn = NULL,
 316             bool readHeaders = false);
 317
 318   // Destructor.
 319   ~LHAupLHEF();
 320
 321   // Helper routine to correctly close files
 322   void closeAllFiles() {
 323     // Close header file if separate, and close main file
 324     if (isHead != is) closeFile(isHead, ifsHead);
 325     closeFile(is, ifs);
 326   }
 327
 328   // Want to use new file with events, but without reinitialization.
 329   void newEventFile(const char* fileIn) {
 330     // Close files and then open new file
 331     closeAllFiles();
 332     is = openFile(fileIn, ifs);
 333
 334     // Set isHead to is to keep expected behaviour in
 335     // fileFound() and closeAllFiles()
 336     isHead = is;
 337   }
 338
 339   // Confirm that file was found and opened as expected.
 340   bool fileFound() { return (isHead->good() && is->good()); }
 341
 342   // Routine for doing the job of reading and setting initialization info.
 343   bool setInit() {return setInitLHEF(*isHead, readHeaders);}
 344
 345   // Routine for doing the job of reading and setting info on next event.
 346   bool setEvent(int = 0) {if (!setNewEventLHEF(*is)) return false;
 347     return setOldEventLHEF();}
 348
 349   // Skip ahead a number of events, which are not considered further.
 350   bool skipEvent(int nSkip) {for (int iSkip = 0; iSkip < nSkip; ++iSkip)
 351     if (!setNewEventLHEF(*is)) return false; return true;}
 352
 353 private:
 354
 355   // File from which to read (or a stringstream).
 356   // Optionally also a file from which to read the LHEF header.
 357   ifstream  ifs,  ifsHead;
 358   istream  *is,  *isHead;
 359
 360   // Flag to read headers or not
 361   bool readHeaders;
 362 };
 363
 364 //==========================================================================
 365
 366 // A derived class with information read from PYTHIA 8 itself, for output.
 367
 368 class LHAupFromPYTHIA8 : public LHAup {
 369
 370 public:
 371
 372   // Constructor.
 373   LHAupFromPYTHIA8(Event* processPtrIn, Info* infoPtrIn) {
 374     processPtr = processPtrIn; infoPtr = infoPtrIn;}
 375
 376   // Destructor.
 377   ~LHAupFromPYTHIA8() {}
 378
 379   // Routine for doing the job of reading and setting initialization info.
 380   bool setInit();
 381
 382   // Routine for doing the job of reading and setting info on next event.
 383   bool setEvent(int = 0);
 384
 385   // Update cross-section information at the end of the run.
 386   bool updateSigma();
 387
 388 private:
 389
 390   // Pointers to process event record and further information.
 391   Event* processPtr;
 392   Info*  infoPtr;
 393
 394 };
 395
 396 //==========================================================================
 397
 398 } // end namespace Pythia8
 399
 400 #endif // Pythia8_LesHouches_H