PYTHIA8/pythia8130/include/LesHouches.h

   1 // LesHouches.h 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 // Header file for Les Houches Accord user process information.
   7 // LHAup: base class for initialization and event information.
   8 // LHAupLHEF: derived class for reading from an Les Houches Event File.
   9 // Code for interfacing with Fortran commonblocks is found in LHAFortran.h.
  10
  11 #ifndef Pythia8_LesHouches_H
  12 #define Pythia8_LesHouches_H
  13
  14 #include "Event.h"
  15 #include "Info.h"
  16 #include "PythiaStdlib.h"
  17 #include "Settings.h"
  18
  19 namespace Pythia8 {
  20
  21 //**************************************************************************
  22
  23 // LHAup is base class for initialization and event information
  24 // from an external parton-level generator.
  25
  26 class LHAup {
  27
  28 public:
  29
  30   // Destructor.
  31   virtual ~LHAup() {}
  32
  33   // Set info pointer.
  34   void setPtr(Info* infoPtrIn) {infoPtr = infoPtrIn;}
  35
  36   // A pure virtual method setInit, wherein all initialization information
  37   // is supposed to be set in the derived class. Can do this by reading a
  38   // file or some other way, as desired. Returns false if it did not work.
  39   virtual bool setInit() = 0;
  40
  41   // Give back info on beams.
  42   int    idBeamA()       const {return idBeamASave;}
  43   int    idBeamB()       const {return idBeamBSave;}
  44   double eBeamA()        const {return eBeamASave;}
  45   double eBeamB()        const {return eBeamBSave;}
  46   int    pdfGroupBeamA() const {return pdfGroupBeamASave;}
  47   int    pdfGroupBeamB() const {return pdfGroupBeamBSave;}
  48   int    pdfSetBeamA()   const {return pdfSetBeamASave;}
  49   int    pdfSetBeamB()   const {return pdfSetBeamBSave;}
  50
  51   // Give back weight strategy.
  52   int    strategy()      const {return strategySave;}
  53
  54   // Give back info on processes.
  55   int    sizeProc()      const {return processes.size();}
  56   int    idProcess(int proc) const {return processes[proc].idProc;}
  57   double xSec(int proc)  const {return processes[proc].xSecProc;}
  58   double xErr(int proc)  const {return processes[proc].xErrProc;}
  59   double xMax(int proc)  const {return processes[proc].xMaxProc;}
  60
  61   // Print the initialization info; useful to check that setting it worked.
  62   void   listInit(ostream& os = cout);
  63
  64   // A pure virtual method setEvent, wherein information on the next event
  65   // is supposed to be set in the derived class.
  66   // Strategies +-1 and +-2: idProcIn is the process type, selected by PYTHIA.
  67   // Strategies +-3 and +-4: idProcIn is dummy; process choice is made locally.
  68   // The method can find the next event by a runtime interface to another
  69   // program, or by reading a file, as desired.
  70   // The method should return false if it did not work.
  71   virtual bool setEvent(int idProcIn = 0) = 0;
  72
  73   // Give back process number, weight, scale, alpha_em, alpha_s.
  74   int    idProcess()       const {return idProc;}
  75   double weight()          const {return weightProc;}
  76   double scale()           const {return scaleProc;}
  77   double alphaQED()        const {return alphaQEDProc;}
  78   double alphaQCD()        const {return alphaQCDProc;}
  79
  80   // Give back info on separate particle.
  81   int    sizePart()        const {return particles.size();}
  82   int    id(int part)      const {return particles[part].idPart;}
  83   int    status(int part)  const {return particles[part].statusPart;}
  84   int    mother1(int part) const {return particles[part].mother1Part;}
  85   int    mother2(int part) const {return particles[part].mother2Part;}
  86   int    col1(int part)    const {return particles[part].col1Part;}
  87   int    col2(int part)    const {return particles[part].col2Part;}
  88   double px(int part)      const {return particles[part].pxPart;}
  89   double py(int part)      const {return particles[part].pyPart;}
  90   double pz(int part)      const {return particles[part].pzPart;}
  91   double e(int part)       const {return particles[part].ePart;}
  92   double m(int part)       const {return particles[part].mPart;}
  93   double tau(int part)     const {return particles[part].tauPart;}
  94   double spin(int part)    const {return particles[part].spinPart;}
  95
  96   // Optional: give back info on parton density values of event.
  97   bool   pdfIsSet()        const {return pdfIsSetSave;}
  98   int    id1()             const {return id1Save;}
  99   int    id2()             const {return id2Save;}
 100   double x1()              const {return x1Save;}
 101   double x2()              const {return x2Save;}
 102   double scalePDF()        const {return scalePDFSave;}
 103   double xpdf1()           const {return xpdf1Save;}
 104   double xpdf2()           const {return xpdf2Save;}
 105
 106   // Print the info; useful to check that reading an event worked.
 107   void   listEvent(ostream& os = cout);
 108
 109   // Four routines to write a Les Houches Event file in steps.
 110   bool   openLHEF(string fileNameIn);
 111   bool   initLHEF();
 112   bool   eventLHEF();
 113   bool   closeLHEF(bool updateInit = false);
 114
 115 protected:
 116
 117   // Constructor. Sets default to be that events come with unit weight.
 118   LHAup(int strategyIn = 3) : strategySave(strategyIn)
 119     { processes.reserve(10); particles.reserve(20); }
 120
 121   // Allow conversion from mb to pb.
 122   static const double CONVERTMB2PB;
 123
 124   // Pointer to various information on the generation.
 125   Info* infoPtr;
 126
 127   // Input beam info.
 128   void setBeamA(int idIn, double eIn, int pdfGroupIn = 0, int pdfSetIn = 0)
 129     { idBeamASave = idIn; eBeamASave = eIn; pdfGroupBeamASave = pdfGroupIn;
 130     pdfSetBeamASave = pdfSetIn;}
 131   void setBeamB(int idIn, double eIn, int pdfGroupIn = 0, int pdfSetIn = 0)
 132     { idBeamBSave = idIn; eBeamBSave = eIn; pdfGroupBeamBSave = pdfGroupIn;
 133     pdfSetBeamBSave = pdfSetIn;}
 134
 135   // Input process weight strategy.
 136   void setStrategy(int strategyIn) {strategySave = strategyIn;}
 137
 138   // Input process info.
 139   void addProcess(int idProcIn, double xSecIn = 1., double xErrIn = 0.,
 140     double xMaxIn = 1.) { processes.push_back( LHAProcess( idProcIn,
 141     xSecIn, xErrIn, xMaxIn)); }
 142
 143   // Possibility to update some cross section info at end of run.
 144   void setXSec(int iP, double xSecIn) {processes[iP].xSecProc = xSecIn;}
 145   void setXErr(int iP, double xErrIn) {processes[iP].xErrProc = xErrIn;}
 146   void setXMax(int iP, double xMaxIn) {processes[iP].xMaxProc = xMaxIn;}
 147
 148   // Input info on the selected process.
 149   void setProcess(int idProcIn = 0, double weightIn = 1., double
 150     scaleIn = 0., double alphaQEDIn = 0.0073, double alphaQCDIn = 0.12) {
 151     idProc = idProcIn; weightProc = weightIn; scaleProc = scaleIn;
 152     alphaQEDProc = alphaQEDIn; alphaQCDProc = alphaQCDIn;
 153     // Clear particle list. Add empty zeroth particle for correct indices.
 154     particles.clear(); addParticle(0); pdfIsSetSave = false;}
 155
 156   // Input particle info, one particle at the time.
 157   void addParticle(int idIn, int statusIn = 0, int mother1In = 0,
 158     int mother2In = 0, int col1In = 0, int col2In = 0, double pxIn = 0.,
 159     double pyIn = 0., double pzIn = 0., double eIn = 0., double mIn = 0.,
 160     double tauIn = 0., double spinIn = 9.) {
 161     particles.push_back( LHAParticle( idIn, statusIn, mother1In, mother2In,
 162     col1In, col2In, pxIn, pyIn, pzIn, eIn, mIn, tauIn, spinIn)); }
 163
 164   // Optionally input info on parton density values of event.
 165   void setPdf(int id1In, int id2In, double x1In, double x2In,
 166     double scalePDFIn, double xpdf1In, double xpdf2In)
 167     { id1Save = id1In; id2Save = id2In; x1Save = x1In; x2Save = x2In;
 168     scalePDFSave = scalePDFIn; xpdf1Save = xpdf1In; xpdf2Save = xpdf2In;
 169     pdfIsSetSave = true;}
 170
 171 private:
 172
 173   // Event weighting and mixing strategy.
 174   int strategySave;
 175
 176   // Beam particle properties.
 177   int idBeamASave, idBeamBSave;
 178   double eBeamASave, eBeamBSave;
 179   int pdfGroupBeamASave, pdfGroupBeamBSave, pdfSetBeamASave, pdfSetBeamBSave;
 180
 181   // A nested class for processes...
 182   class LHAProcess {
 183   public:
 184     LHAProcess() : idProc(0), xSecProc(0.), xErrProc(0.), xMaxProc(0.) { }
 185     LHAProcess(int idProcIn, double xSecIn, double xErrIn, double xMaxIn) :
 186       idProc(idProcIn), xSecProc(xSecIn), xErrProc(xErrIn),
 187       xMaxProc(xMaxIn) { }
 188     int idProc;
 189     double xSecProc, xErrProc, xMaxProc;
 190   } ;
 191
 192   // ...so that the process list can be kept as a vector.
 193   vector<LHAProcess> processes;
 194
 195   // Store info on the selected process.
 196   int idProc;
 197   double weightProc, scaleProc, alphaQEDProc, alphaQCDProc;
 198
 199   // A nested class for particles...
 200   class LHAParticle {
 201   public:
 202     LHAParticle() : idPart(0), statusPart(0), mother1Part(0),
 203       mother2Part(0), col1Part(0), col2Part(0), pxPart(0.), pyPart(0.),
 204       pzPart(0.), ePart(0.), mPart(0.), tauPart(0.), spinPart(9.) { }
 205     LHAParticle(int idIn, int statusIn, int mother1In, int mother2In,
 206       int col1In, int col2In, double pxIn, double pyIn, double pzIn,
 207       double eIn, double mIn, double tauIn, double spinIn) :
 208       idPart(idIn), statusPart(statusIn), mother1Part(mother1In),
 209       mother2Part(mother2In), col1Part(col1In), col2Part(col2In),
 210       pxPart(pxIn), pyPart(pyIn), pzPart(pzIn), ePart(eIn), mPart(mIn),
 211       tauPart(tauIn), spinPart(spinIn) { }
 212     int idPart, statusPart, mother1Part, mother2Part, col1Part, col2Part;
 213     double pxPart, pyPart, pzPart, ePart, mPart, tauPart, spinPart;
 214   } ;
 215
 216   // ...so that the particle list can be kept as a vector.
 217   vector<LHAParticle> particles;
 218
 219   // Optional info on parton density values of event.
 220   bool   pdfIsSetSave;
 221   int    id1Save, id2Save;
 222   double x1Save, x2Save, scalePDFSave, xpdf1Save, xpdf2Save;
 223
 224   // File to which to write Les Houches Event File information.
 225   string fileName;
 226   fstream osLHEF;
 227   char dateNow[12];
 228   char timeNow[9];
 229
 230 };
 231
 232 //**************************************************************************
 233
 234 // A derived class with information read from a Les Houches Event File.
 235
 236 class LHAupLHEF : public LHAup {
 237
 238 public:
 239
 240   // Constructor.
 241   LHAupLHEF(const char* fileIn) : is(fileIn) {}
 242
 243   // Destructor.
 244   ~LHAupLHEF() {}
 245
 246   // Routine for doing the job of reading and setting initialization info.
 247   bool setInit();
 248
 249   // Routine for doing the job of reading and setting info on next event.
 250   bool setEvent(int = 0);
 251
 252 private:
 253
 254   // File from which to read.
 255   ifstream is;
 256
 257 };
 258
 259 //**************************************************************************
 260
 261 // A derived class with information read from PYTHIA 8 itself, for output.
 262
 263 class LHAupFromPYTHIA8 : public LHAup {
 264
 265 public:
 266
 267   // Constructor.
 268   LHAupFromPYTHIA8(Event* processPtrIn, Info* infoPtrIn) {
 269     processPtr = processPtrIn; infoPtr = infoPtrIn;}
 270
 271   // Destructor.
 272   ~LHAupFromPYTHIA8() {}
 273
 274   // Routine for doing the job of reading and setting initialization info.
 275   bool setInit();
 276
 277   // Routine for doing the job of reading and setting info on next event.
 278   bool setEvent(int = 0);
 279
 280   // Update cross-section information at the end of the run.
 281   bool updateSigma();
 282
 283 private:
 284
 285   // Pointers to process event record and further information.
 286   Event* processPtr;
 287   Info*  infoPtr;
 288
 289 };
 290
 291 //**************************************************************************
 292
 293 } // end namespace Pythia8
 294
 295 #endif // Pythia8_LesHouches_H