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 |
---|