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