PYTHIA8/pythia8140/include/LesHouches.h

   1 // LesHouches.h 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 // 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 bool fileFound() {return true;}
  86
  87   // A pure virtual method setInit, wherein all initialization information
  88   // is supposed to be set in the derived class. Can do this by reading a
  89   // file or some other way, as desired. Returns false if it did not work.
  90   virtual bool setInit() = 0;
  91
  92   // Give back info on beams.
  93   int    idBeamA()       const {return idBeamASave;}
  94   int    idBeamB()       const {return idBeamBSave;}
  95   double eBeamA()        const {return eBeamASave;}
  96   double eBeamB()        const {return eBeamBSave;}
  97   int    pdfGroupBeamA() const {return pdfGroupBeamASave;}
  98   int    pdfGroupBeamB() const {return pdfGroupBeamBSave;}
  99   int    pdfSetBeamA()   const {return pdfSetBeamASave;}
 100   int    pdfSetBeamB()   const {return pdfSetBeamBSave;}
 101
 102   // Give back weight strategy.
 103   int    strategy()      const {return strategySave;}
 104
 105   // Give back info on processes.
 106   int    sizeProc()      const {return processes.size();}
 107   int    idProcess(int proc) const {return processes[proc].idProc;}
 108   double xSec(int proc)  const {return processes[proc].xSecProc;}
 109   double xErr(int proc)  const {return processes[proc].xErrProc;}
 110   double xMax(int proc)  const {return processes[proc].xMaxProc;}
 111
 112   // Print the initialization info; useful to check that setting it worked.
 113   void   listInit(ostream& os = cout);
 114
 115   // A pure virtual method setEvent, wherein information on the next event
 116   // is supposed to be set in the derived class.
 117   // Strategies +-1 and +-2: idProcIn is the process type, selected by PYTHIA.
 118   // Strategies +-3 and +-4: idProcIn is dummy; process choice is made locally.
 119   // The method can find the next event by a runtime interface to another
 120   // program, or by reading a file, as desired.
 121   // The method should return false if it did not work.
 122   virtual bool setEvent(int idProcIn = 0) = 0;
 123
 124   // Give back process number, weight, scale, alpha_em, alpha_s.
 125   int    idProcess()       const {return idProc;}
 126   double weight()          const {return weightProc;}
 127   double scale()           const {return scaleProc;}
 128   double alphaQED()        const {return alphaQEDProc;}
 129   double alphaQCD()        const {return alphaQCDProc;}
 130
 131   // Give back info on separate particle.
 132   int    sizePart()        const {return particles.size();}
 133   int    id(int part)      const {return particles[part].idPart;}
 134   int    status(int part)  const {return particles[part].statusPart;}
 135   int    mother1(int part) const {return particles[part].mother1Part;}
 136   int    mother2(int part) const {return particles[part].mother2Part;}
 137   int    col1(int part)    const {return particles[part].col1Part;}
 138   int    col2(int part)    const {return particles[part].col2Part;}
 139   double px(int part)      const {return particles[part].pxPart;}
 140   double py(int part)      const {return particles[part].pyPart;}
 141   double pz(int part)      const {return particles[part].pzPart;}
 142   double e(int part)       const {return particles[part].ePart;}
 143   double m(int part)       const {return particles[part].mPart;}
 144   double tau(int part)     const {return particles[part].tauPart;}
 145   double spin(int part)    const {return particles[part].spinPart;}
 146
 147   // Optional: give back info on parton density values of event.
 148   bool   pdfIsSet()        const {return pdfIsSetSave;}
 149   int    id1()             const {return id1Save;}
 150   int    id2()             const {return id2Save;}
 151   double x1()              const {return x1Save;}
 152   double x2()              const {return x2Save;}
 153   double scalePDF()        const {return scalePDFSave;}
 154   double xpdf1()           const {return xpdf1Save;}
 155   double xpdf2()           const {return xpdf2Save;}
 156
 157   // Print the info; useful to check that reading an event worked.
 158   void   listEvent(ostream& os = cout);
 159
 160   // Skip ahead a number of events, which are not considered further.
 161   // Mainly intended for debug when using the LHAupLHEF class.
 162   virtual bool skipEvent(int nSkip) {
 163     for (int iSkip = 0; iSkip < nSkip; ++iSkip)
 164     if (!setEvent()) return false; return true;}
 165
 166   // Four routines to write a Les Houches Event file in steps.
 167   bool   openLHEF(string fileNameIn);
 168   bool   initLHEF();
 169   bool   eventLHEF();
 170   bool   closeLHEF(bool updateInit = false);
 171
 172 protected:
 173
 174   // Constructor. Sets default to be that events come with unit weight.
 175   LHAup(int strategyIn = 3) : strategySave(strategyIn)
 176     { processes.reserve(10); particles.reserve(20); }
 177
 178   // Allow conversion from mb to pb.
 179   static const double CONVERTMB2PB;
 180
 181   // Pointer to various information on the generation.
 182   Info* infoPtr;
 183
 184   // Input beam info.
 185   void setBeamA(int idIn, double eIn, int pdfGroupIn = 0, int pdfSetIn = 0)
 186     { idBeamASave = idIn; eBeamASave = eIn; pdfGroupBeamASave = pdfGroupIn;
 187     pdfSetBeamASave = pdfSetIn;}
 188   void setBeamB(int idIn, double eIn, int pdfGroupIn = 0, int pdfSetIn = 0)
 189     { idBeamBSave = idIn; eBeamBSave = eIn; pdfGroupBeamBSave = pdfGroupIn;
 190     pdfSetBeamBSave = pdfSetIn;}
 191
 192   // Input process weight strategy.
 193   void setStrategy(int strategyIn) {strategySave = strategyIn;}
 194
 195   // Input process info.
 196   void addProcess(int idProcIn, double xSecIn = 1., double xErrIn = 0.,
 197     double xMaxIn = 1.) { processes.push_back( LHAProcess( idProcIn,
 198     xSecIn, xErrIn, xMaxIn)); }
 199
 200   // Possibility to update some cross section info at end of run.
 201   void setXSec(int iP, double xSecIn) {processes[iP].xSecProc = xSecIn;}
 202   void setXErr(int iP, double xErrIn) {processes[iP].xErrProc = xErrIn;}
 203   void setXMax(int iP, double xMaxIn) {processes[iP].xMaxProc = xMaxIn;}
 204
 205   // Input info on the selected process.
 206   void setProcess(int idProcIn = 0, double weightIn = 1., double
 207     scaleIn = 0., double alphaQEDIn = 0.0073, double alphaQCDIn = 0.12) {
 208     idProc = idProcIn; weightProc = weightIn; scaleProc = scaleIn;
 209     alphaQEDProc = alphaQEDIn; alphaQCDProc = alphaQCDIn;
 210     // Clear particle list. Add empty zeroth particle for correct indices.
 211     particles.clear(); addParticle(0); pdfIsSetSave = false;}
 212
 213   // Input particle info, one particle at the time.
 214   void addParticle(LHAParticle particleIn) {
 215     particles.push_back(particleIn);}
 216   void addParticle(int idIn, int statusIn = 0, int mother1In = 0,
 217     int mother2In = 0, int col1In = 0, int col2In = 0, double pxIn = 0.,
 218     double pyIn = 0., double pzIn = 0., double eIn = 0., double mIn = 0.,
 219     double tauIn = 0., double spinIn = 9.) {
 220     particles.push_back( LHAParticle( idIn, statusIn, mother1In, mother2In,
 221     col1In, col2In, pxIn, pyIn, pzIn, eIn, mIn, tauIn, spinIn)); }
 222
 223   // Optionally input info on parton density values of event.
 224   void setPdf(int id1In, int id2In, double x1In, double x2In,
 225     double scalePDFIn, double xpdf1In, double xpdf2In, bool pdfIsSetIn)
 226     { id1Save = id1In; id2Save = id2In; x1Save = x1In; x2Save = x2In;
 227     scalePDFSave = scalePDFIn; xpdf1Save = xpdf1In; xpdf2Save = xpdf2In;
 228     pdfIsSetSave = pdfIsSetIn;}
 229
 230   // Three routines for LHEF files, but put here for flexibility.
 231   bool setInitLHEF(ifstream& is);
 232   bool setNewEventLHEF(ifstream& is);
 233   bool setOldEventLHEF();
 234
 235   // Event properties from LHEF files, for repeated use.
 236   int    nupSave, idprupSave;
 237   double xwgtupSave, scalupSave, aqedupSave, aqcdupSave;
 238   vector<LHAParticle> particlesSave;
 239   bool   getPDFSave;
 240   int    id1InSave, id2InSave;
 241   double x1InSave, x2InSave, scalePDFInSave, xpdf1InSave, xpdf2InSave;
 242
 243 private:
 244
 245   // Event weighting and mixing strategy.
 246   int strategySave;
 247
 248   // Beam particle properties.
 249   int idBeamASave, idBeamBSave;
 250   double eBeamASave, eBeamBSave;
 251   int pdfGroupBeamASave, pdfGroupBeamBSave, pdfSetBeamASave, pdfSetBeamBSave;
 252
 253   // The process list, stored as a vector of processes.
 254   vector<LHAProcess> processes;
 255
 256   // Store info on the selected process.
 257   int idProc;
 258   double weightProc, scaleProc, alphaQEDProc, alphaQCDProc;
 259
 260   // The particle list, stored as a vector of particles.
 261   vector<LHAParticle> particles;
 262
 263   // Optional info on parton density values of event.
 264   bool   pdfIsSetSave;
 265   int    id1Save, id2Save;
 266   double x1Save, x2Save, scalePDFSave, xpdf1Save, xpdf2Save;
 267
 268   // File to which to write Les Houches Event File information.
 269   string fileName;
 270   fstream osLHEF;
 271   char dateNow[12];
 272   char timeNow[9];
 273
 274 };
 275
 276 //==========================================================================
 277
 278 // A derived class with information read from a Les Houches Event File.
 279
 280 class LHAupLHEF : public LHAup {
 281
 282 public:
 283
 284   // Constructor.
 285   LHAupLHEF(const char* fileIn) : is(fileIn) {}
 286
 287   // Destructor.
 288   ~LHAupLHEF() {}
 289
 290   // Confirm that file was found and opened as expected.
 291   bool fileFound() {return is.good();}
 292
 293   // Routine for doing the job of reading and setting initialization info.
 294   bool setInit() {return setInitLHEF(is);}
 295
 296   // Routine for doing the job of reading and setting info on next event.
 297   bool setEvent(int = 0) {if (!setNewEventLHEF(is)) return false;
 298     return setOldEventLHEF();}
 299
 300   // Skip ahead a number of events, which are not considered further.
 301   bool skipEvent(int nSkip) {for (int iSkip = 0; iSkip < nSkip; ++iSkip)
 302     if (!setNewEventLHEF(is)) return false; return true;}
 303
 304 private:
 305
 306   // File from which to read.
 307   ifstream is;
 308
 309 };
 310
 311 //==========================================================================
 312
 313 // A derived class with information read from PYTHIA 8 itself, for output.
 314
 315 class LHAupFromPYTHIA8 : public LHAup {
 316
 317 public:
 318
 319   // Constructor.
 320   LHAupFromPYTHIA8(Event* processPtrIn, Info* infoPtrIn) {
 321     processPtr = processPtrIn; infoPtr = infoPtrIn;}
 322
 323   // Destructor.
 324   ~LHAupFromPYTHIA8() {}
 325
 326   // Routine for doing the job of reading and setting initialization info.
 327   bool setInit();
 328
 329   // Routine for doing the job of reading and setting info on next event.
 330   bool setEvent(int = 0);
 331
 332   // Update cross-section information at the end of the run.
 333   bool updateSigma();
 334
 335 private:
 336
 337   // Pointers to process event record and further information.
 338   Event* processPtr;
 339   Info*  infoPtr;
 340
 341 };
 342
 343 //==========================================================================
 344
 345 } // end namespace Pythia8
 346
 347 #endif // Pythia8_LesHouches_H