File indexing completed on 2025-01-18 10:06:26
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013 #ifndef Pythia8_LesHouches_H
0014 #define Pythia8_LesHouches_H
0015
0016 #include "Pythia8/Event.h"
0017 #include "Pythia8/Info.h"
0018 #include "Pythia8/PythiaStdlib.h"
0019 #include "Pythia8/Settings.h"
0020 #include "Pythia8/Streams.h"
0021
0022 namespace Pythia8 {
0023
0024
0025
0026
0027
0028 class LHAProcess {
0029
0030 public:
0031
0032
0033 LHAProcess() : idProc(0), xSecProc(0.), xErrProc(0.), xMaxProc(0.) { }
0034 LHAProcess(int idProcIn, double xSecIn, double xErrIn, double xMaxIn) :
0035 idProc(idProcIn), xSecProc(xSecIn), xErrProc(xErrIn),
0036 xMaxProc(xMaxIn) { }
0037
0038
0039 int idProc;
0040 double xSecProc, xErrProc, xMaxProc;
0041
0042 } ;
0043
0044
0045
0046
0047
0048 class LHAParticle {
0049
0050 public:
0051
0052
0053 LHAParticle() : idPart(0), statusPart(0), mother1Part(0),
0054 mother2Part(0), col1Part(0), col2Part(0), pxPart(0.), pyPart(0.),
0055 pzPart(0.), ePart(0.), mPart(0.), tauPart(0.), spinPart(9.),
0056 scalePart(-1.) { }
0057 LHAParticle(int idIn, int statusIn, int mother1In, int mother2In,
0058 int col1In, int col2In, double pxIn, double pyIn, double pzIn,
0059 double eIn, double mIn, double tauIn, double spinIn,
0060 double scaleIn) :
0061 idPart(idIn), statusPart(statusIn), mother1Part(mother1In),
0062 mother2Part(mother2In), col1Part(col1In), col2Part(col2In),
0063 pxPart(pxIn), pyPart(pyIn), pzPart(pzIn), ePart(eIn), mPart(mIn),
0064 tauPart(tauIn), spinPart(spinIn), scalePart(scaleIn) { }
0065
0066
0067 int idPart, statusPart, mother1Part, mother2Part, col1Part, col2Part;
0068 double pxPart, pyPart, pzPart, ePart, mPart, tauPart, spinPart,
0069 scalePart;
0070
0071 } ;
0072
0073
0074
0075
0076
0077
0078 class LHAup {
0079
0080 public:
0081
0082
0083 virtual ~LHAup() {}
0084
0085
0086 void setPtr(Info* infoPtrIn) {infoPtr = infoPtrIn;
0087 loggerPtr = infoPtrIn->loggerPtr;}
0088
0089
0090 virtual void newEventFile(const char*) {}
0091 virtual bool fileFound() {return true;}
0092 virtual bool useExternal() {return false;}
0093
0094
0095
0096
0097 virtual bool setInit() = 0;
0098
0099
0100 int idBeamA() const {return idBeamASave;}
0101 int idBeamB() const {return idBeamBSave;}
0102 double eBeamA() const {return eBeamASave;}
0103 double eBeamB() const {return eBeamBSave;}
0104 int pdfGroupBeamA() const {return pdfGroupBeamASave;}
0105 int pdfGroupBeamB() const {return pdfGroupBeamBSave;}
0106 int pdfSetBeamA() const {return pdfSetBeamASave;}
0107 int pdfSetBeamB() const {return pdfSetBeamBSave;}
0108
0109
0110 int strategy() const {return strategySave;}
0111
0112
0113 int sizeProc() const {return processes.size();}
0114 int idProcess(int proc) const {return processes[proc].idProc;}
0115 double xSec(int proc) const {return processes[proc].xSecProc;}
0116 double xErr(int proc) const {return processes[proc].xErrProc;}
0117 double xMax(int proc) const {return processes[proc].xMaxProc;}
0118 double xSecSum() const {return xSecSumSave;}
0119 double xErrSum() const {return xErrSumSave;}
0120
0121
0122 void listInit();
0123
0124
0125
0126
0127
0128
0129
0130
0131 virtual bool setEvent(int idProcIn = 0) = 0;
0132
0133
0134 int idProcess() const {return idProc;}
0135 double weight() const {return weightProc;}
0136 double scale() const {return scaleProc;}
0137 double alphaQED() const {return alphaQEDProc;}
0138 double alphaQCD() const {return alphaQCDProc;}
0139
0140
0141 int sizePart() const {return particles.size();}
0142 int id(int part) const {return particles[part].idPart;}
0143 int status(int part) const {return particles[part].statusPart;}
0144 int mother1(int part) const {return particles[part].mother1Part;}
0145 int mother2(int part) const {return particles[part].mother2Part;}
0146 int col1(int part) const {return particles[part].col1Part;}
0147 int col2(int part) const {return particles[part].col2Part;}
0148 double px(int part) const {return particles[part].pxPart;}
0149 double py(int part) const {return particles[part].pyPart;}
0150 double pz(int part) const {return particles[part].pzPart;}
0151 double e(int part) const {return particles[part].ePart;}
0152 double m(int part) const {return particles[part].mPart;}
0153 double tau(int part) const {return particles[part].tauPart;}
0154 double spin(int part) const {return particles[part].spinPart;}
0155 double scale(int part) const {return particles[part].scalePart;}
0156
0157
0158 int id1() const {return id1Save;}
0159 int id2() const {return id2Save;}
0160 double x1() const {return x1Save;}
0161 double x2() const {return x2Save;}
0162
0163
0164 bool pdfIsSet() const {return pdfIsSetSave;}
0165 int id1pdf() const {return id1pdfSave;}
0166 int id2pdf() const {return id2pdfSave;}
0167 double x1pdf() const {return x1pdfSave;}
0168 double x2pdf() const {return x2pdfSave;}
0169 double scalePDF() const {return scalePDFSave;}
0170 double pdf1() const {return pdf1Save;}
0171 double pdf2() const {return pdf2Save;}
0172
0173
0174 bool scaleShowersIsSet() const {return scaleShowersIsSetSave;}
0175 double scaleShowers(int i) const {return scaleShowersSave[i];}
0176
0177
0178 void listEvent();
0179
0180
0181
0182 virtual bool skipEvent(int nSkip) {
0183 for (int iSkip = 0; iSkip < nSkip; ++iSkip) if (!setEvent()) return false;
0184 return true;}
0185
0186
0187 virtual bool openLHEF(string fileNameIn);
0188 virtual bool closeLHEF(bool updateInit = false);
0189 bool initLHEF();
0190 bool eventLHEF(bool verbose = true);
0191
0192
0193 string getFileName() const {return fileName;}
0194
0195 protected:
0196
0197
0198 LHAup(int strategyIn = 3) : infoPtr(), nupSave(), idprupSave(),
0199 xwgtupSave(), scalupSave(), aqedupSave(), aqcdupSave(), xSecSumSave(),
0200 xErrSumSave(), getPDFSave(), getScale(), getScaleShowers(),
0201 id1InSave(), id2InSave(), id1pdfInSave(), id2pdfInSave(), x1InSave(),
0202 x2InSave(), x1pdfInSave(), x2pdfInSave(), scalePDFInSave(),
0203 pdf1InSave(), pdf2InSave(), scaleShowersInSave(), fileName("void"),
0204 dateNow(), timeNow(), strategySave(strategyIn), idBeamASave(),
0205 idBeamBSave(), eBeamASave(), eBeamBSave(), pdfGroupBeamASave(),
0206 pdfGroupBeamBSave(), pdfSetBeamASave(), pdfSetBeamBSave(), idProc(),
0207 weightProc(), scaleProc(), alphaQEDProc(), alphaQCDProc(),
0208 pdfIsSetSave(), scaleShowersIsSetSave(false), id1Save(),
0209 id2Save(), id1pdfSave(), id2pdfSave(), x1Save(), x2Save(), x1pdfSave(),
0210 x2pdfSave(), scalePDFSave(), pdf1Save(), pdf2Save(), scaleShowersSave() {
0211 processes.reserve(10); particles.reserve(20);
0212 setBeamA( 0, 0., 0, 0); setBeamB( 0, 0., 0, 0); }
0213
0214
0215 static const double CONVERTMB2PB;
0216
0217
0218 Info* infoPtr;
0219 Logger* loggerPtr;
0220
0221
0222 void setBeamA(int idIn, double eIn, int pdfGroupIn = 0, int pdfSetIn = 0)
0223 { idBeamASave = idIn; eBeamASave = eIn; pdfGroupBeamASave = pdfGroupIn;
0224 pdfSetBeamASave = pdfSetIn;}
0225 void setBeamB(int idIn, double eIn, int pdfGroupIn = 0, int pdfSetIn = 0)
0226 { idBeamBSave = idIn; eBeamBSave = eIn; pdfGroupBeamBSave = pdfGroupIn;
0227 pdfSetBeamBSave = pdfSetIn;}
0228
0229
0230 void setStrategy(int strategyIn) {strategySave = strategyIn;}
0231
0232
0233 void addProcess(int idProcIn, double xSecIn = 1., double xErrIn = 0.,
0234 double xMaxIn = 1.) { processes.push_back( LHAProcess( idProcIn,
0235 xSecIn, xErrIn, xMaxIn)); }
0236
0237
0238 void setXSec(int iP, double xSecIn) {processes[iP].xSecProc = xSecIn;}
0239 void setXErr(int iP, double xErrIn) {processes[iP].xErrProc = xErrIn;}
0240 void setXMax(int iP, double xMaxIn) {processes[iP].xMaxProc = xMaxIn;}
0241
0242
0243 void setProcess(int idProcIn = 0, double weightIn = 1., double
0244 scaleIn = 0., double alphaQEDIn = 0.0073, double alphaQCDIn = 0.12) {
0245 idProc = idProcIn; weightProc = weightIn; scaleProc = scaleIn;
0246 alphaQEDProc = alphaQEDIn; alphaQCDProc = alphaQCDIn;
0247
0248 particles.clear(); addParticle(0); pdfIsSetSave = false;
0249 scaleShowersIsSetSave = false;}
0250
0251
0252 void addParticle(LHAParticle particleIn) {
0253 particles.push_back(particleIn);}
0254 void addParticle(int idIn, int statusIn = 0, int mother1In = 0,
0255 int mother2In = 0, int col1In = 0, int col2In = 0, double pxIn = 0.,
0256 double pyIn = 0., double pzIn = 0., double eIn = 0., double mIn = 0.,
0257 double tauIn = 0., double spinIn = 9., double scaleIn = -1.) {
0258 particles.push_back( LHAParticle( idIn, statusIn, mother1In, mother2In,
0259 col1In, col2In, pxIn, pyIn, pzIn, eIn, mIn, tauIn, spinIn,
0260 scaleIn) ); }
0261
0262
0263 void setIdX(int id1In, int id2In, double x1In, double x2In)
0264 { id1Save = id1In; id2Save = id2In; x1Save = x1In; x2Save = x2In;}
0265
0266
0267 void setPdf(int id1pdfIn, int id2pdfIn, double x1pdfIn, double x2pdfIn,
0268 double scalePDFIn, double pdf1In, double pdf2In, bool pdfIsSetIn)
0269 { id1pdfSave = id1pdfIn; id2pdfSave = id2pdfIn; x1pdfSave = x1pdfIn;
0270 x2pdfSave = x2pdfIn; scalePDFSave = scalePDFIn; pdf1Save = pdf1In;
0271 pdf2Save = pdf2In; pdfIsSetSave = pdfIsSetIn;}
0272
0273
0274 void setScaleShowers( double scaleIn1, double scaleIn2 = 0.)
0275 { scaleShowersSave[0] = scaleIn1; scaleShowersSave[1] = scaleIn2;
0276 scaleShowersIsSetSave = true;}
0277
0278
0279 bool setInitLHEF(istream& is, bool readHeaders = false);
0280 bool setNewEventLHEF(istream& is);
0281 bool setOldEventLHEF();
0282
0283
0284
0285
0286
0287
0288 istream* openFile(const char *fn, ifstream &ifs);
0289 void closeFile(istream *&is, ifstream &ifs);
0290
0291
0292
0293
0294 void setInfoHeader(const string &key, const string &val) {
0295 infoPtr->setHeader(key, val); }
0296
0297
0298 int nupSave, idprupSave;
0299 double xwgtupSave, scalupSave, aqedupSave, aqcdupSave, xSecSumSave,
0300 xErrSumSave;
0301 vector<LHAParticle> particlesSave;
0302 bool getPDFSave, getScale, getScaleShowers;
0303 int id1InSave, id2InSave, id1pdfInSave, id2pdfInSave;
0304 double x1InSave, x2InSave, x1pdfInSave, x2pdfInSave, scalePDFInSave,
0305 pdf1InSave, pdf2InSave, scaleShowersInSave[2];
0306
0307
0308 string fileName;
0309 fstream osLHEF;
0310 char dateNow[12];
0311 char timeNow[9];
0312
0313 private:
0314
0315
0316 int strategySave;
0317
0318
0319 int idBeamASave, idBeamBSave;
0320 double eBeamASave, eBeamBSave;
0321 int pdfGroupBeamASave, pdfGroupBeamBSave, pdfSetBeamASave,
0322 pdfSetBeamBSave;
0323
0324
0325 vector<LHAProcess> processes;
0326
0327
0328 int idProc;
0329 double weightProc, scaleProc, alphaQEDProc, alphaQCDProc;
0330
0331
0332 vector<LHAParticle> particles;
0333
0334
0335 bool pdfIsSetSave, scaleShowersIsSetSave;
0336 int id1Save, id2Save, id1pdfSave, id2pdfSave;
0337 double x1Save, x2Save, x1pdfSave, x2pdfSave, scalePDFSave, pdf1Save,
0338 pdf2Save, scaleShowersSave[2];
0339
0340 };
0341
0342
0343
0344
0345
0346 class LHAupLHEF : public LHAup {
0347
0348 public:
0349
0350
0351 LHAupLHEF(Pythia8::Info* infoPtrIn, istream* isIn, istream* isHeadIn,
0352 bool readHeadersIn = false, bool setScalesFromLHEFIn = false ) :
0353 filename(""), headerfile(""),
0354 is(isIn), is_gz(nullptr), isHead(isHeadIn), isHead_gz(nullptr),
0355 readHeaders(readHeadersIn), reader(is),
0356 setScalesFromLHEF(setScalesFromLHEFIn), hasExtFileStream(true),
0357 hasExtHeaderStream(true) {setPtr(infoPtrIn);}
0358
0359 LHAupLHEF(Pythia8::Info* infoPtrIn, const char* filenameIn,
0360 const char* headerIn = nullptr, bool readHeadersIn = false,
0361 bool setScalesFromLHEFIn = false ) :
0362 filename(filenameIn), headerfile(headerIn),
0363 is(nullptr), is_gz(nullptr), isHead(nullptr), isHead_gz(nullptr),
0364 readHeaders(readHeadersIn), reader(filenameIn),
0365 setScalesFromLHEF(setScalesFromLHEFIn), hasExtFileStream(false),
0366 hasExtHeaderStream(false) {
0367 setPtr(infoPtrIn);
0368 is = (openFile(filenameIn, ifs));
0369 isHead = (headerfile == nullptr) ? is : openFile(headerfile, ifsHead);
0370 is_gz = new igzstream(filename);
0371 isHead_gz = (headerfile == nullptr) ? is_gz : new igzstream(headerfile);
0372 }
0373
0374
0375 ~LHAupLHEF() {
0376
0377 closeAllFiles();
0378 }
0379
0380
0381 void closeAllFiles() {
0382
0383 if (!hasExtHeaderStream && isHead_gz != is_gz) isHead_gz->close();
0384 if (isHead_gz != is_gz) delete isHead_gz;
0385 if (is_gz) is_gz->close();
0386 if (is_gz) delete is_gz;
0387
0388
0389 if (!hasExtHeaderStream && isHead != is) closeFile(isHead, ifsHead);
0390 if (!hasExtFileStream) closeFile(is, ifs);
0391 }
0392
0393
0394 void newEventFile(const char* filenameIn) {
0395
0396 closeAllFiles();
0397 is = (openFile(filenameIn, ifs));
0398 is_gz = new igzstream(filenameIn);
0399
0400 reader.setup(filenameIn);
0401
0402
0403 isHead = is;
0404 isHead_gz = is_gz;
0405 }
0406
0407
0408 bool fileFound() {return (useExternal() || (isHead->good() && is->good()));}
0409 bool useExternal() {return (hasExtHeaderStream && hasExtFileStream);}
0410
0411
0412 bool setInit() {
0413 return setInitLHEF(*isHead, readHeaders);
0414 }
0415
0416
0417 bool setInitLHEF( istream & isIn, bool readHead);
0418
0419
0420 bool setEvent(int = 0) {
0421 if (!setNewEventLHEF()) return false;
0422 return setOldEventLHEF();
0423 }
0424
0425
0426 bool skipEvent(int nSkip) {
0427 for (int iSkip = 0; iSkip < nSkip; ++iSkip)
0428 if (!setNewEventLHEF()) return false;
0429 return true;
0430 }
0431
0432
0433 bool setNewEventLHEF();
0434
0435
0436 bool updateSigma() {return true;}
0437
0438 protected:
0439
0440
0441 bool getLine(string & line, bool header = true) {
0442 #ifdef GZIP
0443 if ( isHead_gz && header && !getline(*isHead_gz, line)) return false;
0444 else if ( is_gz && !header && !getline(*is_gz, line)) return false;
0445 if (header && !getline(*isHead, line)) return false;
0446 else if (!header && !getline(*is, line)) return false;
0447 #else
0448 if (header && !getline(*isHead, line)) return false;
0449 else if (!header && !getline(*is, line)) return false;
0450 #endif
0451
0452 replace(line.begin(),line.end(),'\'','\"');
0453 return true;
0454 }
0455
0456 private:
0457
0458 const char* filename;
0459 const char* headerfile;
0460
0461
0462
0463 istream *is;
0464 igzstream *is_gz;
0465 ifstream ifs;
0466 istream *isHead;
0467 igzstream *isHead_gz;
0468 ifstream ifsHead;
0469
0470
0471 bool readHeaders;
0472
0473 Reader reader;
0474
0475
0476 bool setScalesFromLHEF, hasExtFileStream, hasExtHeaderStream;
0477
0478 };
0479
0480
0481
0482
0483
0484 class LHAupFromPYTHIA8 : public LHAup {
0485
0486 public:
0487
0488
0489 LHAupFromPYTHIA8(Event* processPtrIn, const Info* infoPtrIn) {
0490 processPtr = processPtrIn; infoPtr = infoPtrIn;}
0491
0492
0493 ~LHAupFromPYTHIA8() {}
0494
0495
0496 bool setInit();
0497
0498
0499 bool setEvent(int = 0);
0500
0501
0502 bool updateSigma();
0503
0504 private:
0505
0506
0507 Event* processPtr;
0508
0509
0510 const Info* infoPtr;
0511
0512 };
0513
0514
0515
0516
0517
0518
0519 class LHEF3FromPythia8 : public LHAup {
0520
0521 public:
0522
0523
0524 LHEF3FromPythia8(Event* eventPtrIn, const Info* infoPtrIn,
0525 int pDigitsIn = 15, bool writeToFileIn = true) :
0526 eventPtr(eventPtrIn),infoPtr(infoPtrIn),
0527 particleDataPtr(infoPtrIn->particleDataPtr),
0528 settingsPtr(infoPtrIn->settingsPtr), writer(osLHEF),
0529 pDigits(pDigitsIn), writeToFile(writeToFileIn) {}
0530
0531
0532 bool setInit();
0533
0534
0535 void setEventPtr(Event* evPtr) { eventPtr = evPtr; }
0536 bool setEvent(int = 0);
0537 string getEventString() { return writer.getEventString(&hepeup); }
0538
0539
0540 bool openLHEF(string fileNameIn);
0541
0542
0543 bool closeLHEF(bool updateInit = false);
0544
0545
0546 HEPRUP heprup;
0547 HEPEUP hepeup;
0548
0549
0550 Event* eventPtr;
0551
0552
0553 const Info* infoPtr;
0554
0555 ParticleData* particleDataPtr;
0556
0557 private:
0558
0559
0560 Settings* settingsPtr;
0561
0562
0563 Writer writer;
0564
0565
0566 int pDigits;
0567 bool writeToFile;
0568
0569 };
0570
0571
0572
0573 }
0574
0575 #endif