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