PYTHIA  8.312
GeneratorInput.h
1 // GeneratorInput.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2024 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // Primary Author: Richard Corke
7 // Secondary Author: Stephen Mrenna
8 // This file provides the following classes:
9 // AlpgenPar: a class for parsing ALPGEN parameter files
10 // and reading back out the values
11 // LHAupAlpgen: an LHAup derived class for reading in ALPGEN
12 // format event files
13 // AlpgenHooks: a UserHooks derived class for providing
14 // 'Alpgen:*' user options
15 // MadgraphPar: a class for parsing MadGraph parameter files
16 // and reading back out the values
17 // Example usage is shown in main32.cc, and further details
18 // can be found in the 'Jet Matching Style' manual page.
19 // Minor changes were made by the secondary author for integration
20 // with Madgraph-style matching, and Madgraph input was added.
21 
22 #ifndef Pythia8_GeneratorInput_H
23 #define Pythia8_GeneratorInput_H
24 
25 // Includes and namespace
26 #include "Pythia8/Pythia.h"
27 
28 namespace Pythia8 {
29 
30 //==========================================================================
31 
32 // AlpgenPar: Class to parse ALPGEN parameter files and make them
33 // available through a simple interface
34 
35 class AlpgenPar {
36 
37 public:
38 
39  // Constructor
40  AlpgenPar() {}
41 
42  // Parse as incoming ALPGEN parameter file (passed as string)
43  bool parse(const string paramStr);
44 
45  // Parse an incoming parameter line
46  void extractRunParam(string line);
47 
48  // Check if a parameter exists
49  bool haveParam(const string &paramIn) {
50  return (params.find(paramIn) == params.end()) ? false : true; }
51 
52  // Get a parameter as a double or integer.
53  // Caller should have already checked existance of the parameter.
54  double getParam(const string &paramIn) {
55  return (haveParam(paramIn)) ? params[paramIn] : 0.; }
56  int getParamAsInt(const string &paramIn) {
57  return (haveParam(paramIn)) ? int(params[paramIn]) : 0.; }
58 
59  // Print parameters read from the '.par' file
60  void printParams();
61 
62 private:
63 
64  // Warn if a parameter is going to be overwriten
65  void warnParamOverwrite(const string &paramIn, double val);
66 
67  // Simple string trimmer
68  static string trim(string s);
69 
70  // Storage for parameters
71  map<string,double> params;
72 
73  // Constants
74  static const double ZEROTHRESHOLD;
75 
76 };
77 
78 //==========================================================================
79 
80 // LHAupAlpgen: LHAup derived class for reading in ALPGEN format
81 // event files.
82 
83 class LHAupAlpgen : public LHAup {
84 
85 public:
86 
87  // Constructor and destructor.
88  LHAupAlpgen(const char *baseFNin);
89  ~LHAupAlpgen() { closeFile(isUnw, ifsUnw); }
90 
91  // Override fileFound routine from LHAup.
92  bool fileFound() { return (isUnw != nullptr); }
93 
94  // Override setInit/setEvent routines from LHAup.
95  bool setInit();
96  bool setEvent(int);
97 
98  // Print list of particles; mainly intended for debugging
99  void printParticles();
100 
101 private:
102 
103  // Add resonances to incoming event.
104  bool addResonances();
105 
106  // Rescale momenta to remove any imbalances.
107  bool rescaleMomenta();
108 
109  // Class variables
110  string baseFN, parFN, unwFN; // Incoming filenames
111  AlpgenPar alpgenPar; // Parameter database
112  int lprup; // Process code
113  double ebmupA, ebmupB; // Beam energies
114  int ihvy1, ihvy2; // Heavy flavours for certain processes
115  double mb; // Bottom mass
116  ifstream ifsUnw; // Input file stream for 'unw' file
117  istream* isUnw; // Input stream for 'unw' file
118  vector<LHAParticle> myParticles; // Local storage for particles
119 
120  // Constants
121  static const bool LHADEBUG, LHADEBUGRESCALE;
122  static const double ZEROTHRESHOLD, EWARNTHRESHOLD, PTWARNTHRESHOLD,
123  INCOMINGMIN;
124 
125 };
126 
127 //==========================================================================
128 
129 // AlpgenHooks: provides Alpgen file reading options.
130 // Note that it is defined with virtual inheritance, so that it can
131 // be combined with other UserHooks classes, see e.g. main32.cc.
132 
133 class AlpgenHooks : virtual public UserHooks {
134 
135 public:
136 
137  // Constructor and destructor
138  AlpgenHooks(Pythia &pythia);
139  ~AlpgenHooks() {}
140 
141  // Override initAfterBeams routine from UserHooks
142  bool initAfterBeams();
143 
144 private:
145 
146  // LHAupAlpgen pointer
147  shared_ptr<LHAupAlpgen> LHAagPtr;
148 
149 };
150 
151 //==========================================================================
152 
153 // MadgraphPar: Class to parse the Madgraph header parameters and
154 // make them available through a simple interface
155 
156 class MadgraphPar {
157 
158 public:
159 
160  // Constructor
162 
163  // Parse an incoming Madgraph parameter file string
164  bool parse(const string paramStr);
165 
166  // Parse an incoming parameter line
167  void extractRunParam(string line);
168 
169  // Check if a parameter exists
170  bool haveParam(const string &paramIn) {
171  return (params.find(paramIn) == params.end()) ? false : true; }
172 
173  // Get a parameter as a double or integer.
174  // Caller should have already checked existance of the parameter.
175  double getParam(const string &paramIn) {
176  return (haveParam(paramIn)) ? params[paramIn] : 0.; }
177  int getParamAsInt(const string &paramIn) {
178  return (haveParam(paramIn)) ? int(params[paramIn]) : 0.; }
179 
180  // Print parameters read from the '.par' file
181  void printParams();
182 
183 private:
184 
185  // Warn if a parameter is going to be overwriten
186  void warnParamOverwrite(const string &paramIn, double val);
187 
188  // Simple string trimmer
189  static string trim(string s);
190 
191  // Storage for parameters
192  map<string,double> params;
193 
194  // Constants
195  static const double ZEROTHRESHOLD;
196 
197 };
198 
199 //==========================================================================
200 
201 // Main implementation of AlpgenPar class.
202 // This may be split out to a separate C++ file if desired,
203 // but currently included here for ease of use.
204 
205 //--------------------------------------------------------------------------
206 
207 // Constants: could be changed here if desired, but normally should not.
208 // These are of technical nature, as described for each.
209 
210 // A zero threshold value for double comparisons.
211 const double AlpgenPar::ZEROTHRESHOLD = 1e-10;
212 
213 //--------------------------------------------------------------------------
214 
215 // Warn if e/pT imbalance greater than these values
216 // Parse an incoming Alpgen parameter file string
217 
218 inline bool AlpgenPar::parse(const string paramStr) {
219 
220  // Read par file in blocks:
221  // 0 - process information
222  // 1 - run parameters
223  // 2 - cross sections
224  int block = 0;
225 
226  // Loop over incoming lines
227  stringstream paramStream(paramStr);
228  string line;
229  while (getline(paramStream, line)) {
230 
231  // Change to 'run parameters' block
232  if (line.find("run parameters") != string::npos) {
233  block = 1;
234 
235  // End of 'run parameters' block
236  } else if (line.find("end parameters") != string::npos) {
237  block = 2;
238 
239  // Do not extract anything from block 0 so far
240  } else if (block == 0) {
241 
242  // Block 1 or 2: extract parameters
243  } else {
244  extractRunParam(line);
245 
246  }
247  } // while (getline(paramStream, line))
248 
249  return true;
250 }
251 
252 //--------------------------------------------------------------------------
253 
254 // Parse an incoming parameter line
255 
256 inline void AlpgenPar::extractRunParam(string line) {
257 
258  // Extract information to the right of the final '!' character
259  size_t idx = line.rfind("!");
260  if (idx == string::npos) return;
261  string paramName = trim(line.substr(idx + 1));
262  string paramVal = trim(line.substr(0, idx));
263  istringstream iss(paramVal);
264 
265  // Special case: 'hard process code' - single integer input
266  double val;
267  if (paramName == "hard process code") {
268  iss >> val;
269  warnParamOverwrite("hpc", val);
270  params["hpc"] = val;
271 
272  // Special case: 'Crosssection +- error (pb)' - two double values
273  } else if (paramName.find("Crosssection") == 0) {
274  double xerrup;
275  iss >> val >> xerrup;
276  warnParamOverwrite("xsecup", val);
277  warnParamOverwrite("xerrup", val);
278  params["xsecup"] = val;
279  params["xerrup"] = xerrup;
280 
281  // Special case: 'unwtd events, lum (pb-1)' - integer and double values
282  } else if (paramName.find("unwtd events") == 0) {
283  int nevent;
284  iss >> nevent >> val;
285  warnParamOverwrite("nevent", val);
286  warnParamOverwrite("lum", val);
287  params["nevent"] = nevent;
288  params["lum"] = val;
289 
290  // Special case: 'mc,mb,...' - split on ',' for name and ' ' for values
291  } else if (paramName.find(",") != string::npos) {
292 
293  // Simple tokeniser
294  string paramNameNow;
295  istringstream issName(paramName);
296  while (getline(issName, paramNameNow, ',')) {
297  iss >> val;
298  warnParamOverwrite(paramNameNow, val);
299  params[paramNameNow] = val;
300  }
301 
302  // Default case: assume integer and double on the left
303  } else {
304  int paramIdx;
305  iss >> paramIdx >> val;
306  warnParamOverwrite(paramName, val);
307  params[paramName] = val;
308  }
309 }
310 
311 //--------------------------------------------------------------------------
312 
313 // Print parameters read from the '.par' file
314 
315 inline void AlpgenPar::printParams() {
316 
317  // Loop over all stored parameters and print
318  cout << fixed << setprecision(3) << endl
319  << " *------- Alpgen parameters -------*" << endl;
320  for (map < string, double >::iterator it = params.begin();
321  it != params.end(); ++it)
322  cout << " | " << left << setw(13) << it->first
323  << " | " << right << setw(13) << it->second
324  << " |" << endl;
325  cout << " *-----------------------------------*" << endl;
326 }
327 
328 //--------------------------------------------------------------------------
329 
330 // Warn if a parameter is going to be overwriten
331 
332 inline void AlpgenPar::warnParamOverwrite(const string &paramIn, double val) {
333 
334  // Check if present and if new value is different
335  if (haveParam(paramIn) &&
336  abs(getParam(paramIn) - val) > ZEROTHRESHOLD) {
337  cout << "Warning in LHAupAlpgen::warnParamOverwrite:"
338  << " overwriting existing parameter" << paramIn << endl;
339  }
340 }
341 
342 //--------------------------------------------------------------------------
343 
344 // Simple string trimmer
345 
346 inline string AlpgenPar::trim(string s) {
347 
348  // Remove whitespace in incoming string
349  size_t i;
350  if ((i = s.find_last_not_of(" \t\r\n")) != string::npos)
351  s = s.substr(0, i + 1);
352  if ((i = s.find_first_not_of(" \t\r\n")) != string::npos)
353  s = s.substr(i);
354  return s;
355 }
356 
357 //==========================================================================
358 
359 // Main implementation of LHAupAlpgen class.
360 // This may be split out to a separate C++ file if desired,
361 // but currently included here for ease of use.
362 
363 // ----------------------------------------------------------------------
364 
365 // Constants: could be changed here if desired, but normally should not.
366 // These are of technical nature, as described for each.
367 
368 // Debug flag to print all particles in each event.
369 const bool LHAupAlpgen::LHADEBUG = false;
370 
371 // Debug flag to print particles when an e/p imbalance is found.
372 const bool LHAupAlpgen::LHADEBUGRESCALE = false;
373 
374 // A zero threshold value for double comparisons.
375 const double LHAupAlpgen::ZEROTHRESHOLD = 1e-10;
376 
377 // Warn if e/pT imbalance greater than these values
378 const double LHAupAlpgen::EWARNTHRESHOLD = 3e-3;
379 const double LHAupAlpgen::PTWARNTHRESHOLD = 1e-3;
380 
381 // If incoming e/pZ is 0, it is reset to this value
382 const double LHAupAlpgen::INCOMINGMIN = 1e-3;
383 
384 // ----------------------------------------------------------------------
385 
386 // Constructor. Opens event file.
387 
388 LHAupAlpgen::LHAupAlpgen(const char* baseFNin)
389  : baseFN(baseFNin), alpgenPar(), isUnw(nullptr) {
390 
391  // Open '.unw' events file (with possible gzip support)
392 #ifdef GZIP
393  unwFN = baseFN + ".unw.gz";
394  isUnw = openFile(unwFN.c_str(), ifsUnw);
395  if (!ifsUnw.is_open()) closeFile(isUnw, ifsUnw);
396 #endif
397  if (isUnw == nullptr) {
398  unwFN = baseFN + ".unw";
399  isUnw = openFile(unwFN.c_str(), ifsUnw);
400  if (!ifsUnw.is_open()) {
401  cout << "Error in LHAupAlpgen::LHAupAlpgen: "
402  << "cannot open event file " << unwFN << endl;
403  closeFile(isUnw, ifsUnw);
404  }
405  }
406 }
407 
408 // ----------------------------------------------------------------------
409 
410 // setInit is a virtual method that must be finalised here.
411 // Opens parameter file and parses it, sets up beams, strategy and processes.
412 
413 inline bool LHAupAlpgen::setInit() {
414 
415  // Read in '_unw.par' file to get parameters
416  ifstream ifsPar;
417  istream* isPar = nullptr;
418 
419  // Try gzip file first then normal file afterwards
420 #ifdef GZIP
421  parFN = baseFN + "_unw.par.gz";
422  isPar = openFile(parFN.c_str(), ifsPar);
423  if (!ifsPar.is_open()) closeFile(isPar, ifsPar);
424 #endif
425  if (isPar == nullptr) {
426  parFN = baseFN + "_unw.par";
427  isPar = openFile(parFN.c_str(), ifsPar);
428  if (!ifsPar.is_open()) {
429  cout << "Error in LHAupAlpgen::LHAupAlpgen: "
430  << "cannot open parameter file " << parFN << endl;
431  closeFile(isPar, ifsPar);
432  return false;
433  }
434  }
435 
436  // Read entire contents into string and close file
437  string paramStr((std::istreambuf_iterator<char>(isPar->rdbuf())),
438  std::istreambuf_iterator<char>());
439 
440  // Make sure we reached EOF and not other error
441  if (ifsPar.bad()) {
442  cout << "Error in LHAupAlpgen::LHAupAlpgen: "
443  << "cannot read parameter file " << parFN << endl;
444  return false;
445  }
446  closeFile(isPar, ifsPar);
447 
448  // Parse file and set LHEF header
449  alpgenPar.parse(paramStr);
450  setInfoHeader("AlpgenPar", paramStr);
451 
452  // Check that all required parameters are present
453  if (!alpgenPar.haveParam("ih2") || !alpgenPar.haveParam("ebeam") ||
454  !alpgenPar.haveParam("hpc") || !alpgenPar.haveParam("xsecup") ||
455  !alpgenPar.haveParam("xerrup")) {
456  cout << "Error in LHAupAlpgen::setInit: "
457  << "missing input parameters" << endl;
458  return false;
459  }
460 
461  // Beam IDs
462  int ih2 = alpgenPar.getParamAsInt("ih2");
463  int idbmupA = 2212;
464  int idbmupB = (ih2 == 1) ? 2212 : -2212;
465 
466  // Beam energies
467  double ebeam = alpgenPar.getParam("ebeam");
468  ebmupA = ebeam;
469  ebmupB = ebmupA;
470 
471  // PDF group and set (at the moment not set)
472  int pdfgupA = 0, pdfsupA = 0;
473  int pdfgupB = 0, pdfsupB = 0;
474 
475  // Strategy is for unweighted events and xmaxup not important
476  int idwtup = 3;
477  double xmaxup = 0.;
478 
479  // Get hard process code
480  lprup = alpgenPar.getParamAsInt("hpc");
481 
482  // Check for unsupported processes
483  if (lprup == 7 || lprup == 8 || lprup == 13) {
484  cout << "Error in LHAupAlpgen::setInit: "
485  << "process not implemented" << endl;
486  return false;
487  }
488 
489  // Depending on the process code, get heavy flavour information:
490  // 6 = QQbar + jets
491  // 7 = QQbar + Q'Qbar' + jets
492  // 8 = QQbar + Higgs + jets
493  // 16 = QQbar + gamma + jets
494  if (lprup == 6 || lprup == 7 || lprup == 8 || lprup == 16) {
495  if (!alpgenPar.haveParam("ihvy")) {
496  cout << "Error in LHAupAlpgen::setInit: "
497  << "heavy flavour information not present" << endl;
498  return false;
499  }
500  ihvy1 = alpgenPar.getParamAsInt("ihvy");
501 
502  } else ihvy1 = -1;
503  if (lprup == 7) {
504  if (!alpgenPar.haveParam("ihvy2")) {
505  cout << "Error in LHAupAlpgen::setInit: "
506  << "heavy flavour information not present" << endl;
507  return false;
508  }
509  ihvy2 = alpgenPar.getParamAsInt("ihvy2");
510  } else ihvy2 = -1;
511  // For single top (process 13), get b mass to set incoming
512  mb = -1.;
513  if (lprup == 13) {
514  if (!alpgenPar.haveParam("mb")) {
515  cout << "Error in LHAupAlpgen::setInit: "
516  << "heavy flavour information not present" << endl;
517  return false;
518  }
519  mb = alpgenPar.getParam("mb");
520  }
521 
522  // Set the beams
523  setBeamA(idbmupA, ebmupA, pdfgupA, pdfsupA);
524  setBeamB(idbmupB, ebmupB, pdfgupB, pdfsupB);
525  setStrategy(idwtup);
526 
527  // Add the process
528  double xsecup = alpgenPar.getParam("xsecup");
529  double xerrup = alpgenPar.getParam("xerrup");
530  addProcess(lprup, xsecup, xerrup, xmaxup);
531  xSecSumSave = xsecup;
532  xErrSumSave = xerrup;
533 
534  // All okay
535  return true;
536 }
537 
538 // ----------------------------------------------------------------------
539 
540 // setEvent is a virtual method that must be finalised here.
541 // Read in an event from the 'unw' file and setup.
542 
543 inline bool LHAupAlpgen::setEvent(int) {
544 
545  // Read in the first line of the event
546  int nEvent, iProc, nParton;
547  double Swgt, Sq;
548  string line;
549  if (!getline(*isUnw, line)) {
550  // Read was bad
551  if (ifsUnw.bad()) {
552  cout << "Error in LHAupAlpgen::setEvent: "
553  << "could not read events from file" << endl;
554  return false;
555  }
556  // End of file reached
557  cout << "Error in LHAupAlpgen::setEvent: "
558  << "end of file reached" << endl;
559  return false;
560  }
561  istringstream iss1(line);
562  iss1 >> nEvent >> iProc >> nParton >> Swgt >> Sq;
563 
564  // Set the process details (ignore alphaQED and alphaQCD parameters)
565  double wgtT = Swgt, scaleT = Sq;
566  setProcess(lprup, wgtT, scaleT);
567 
568  // Incoming flavour and x information for later
569  int id1T, id2T;
570  double x1T, x2T;
571  // Temporary storage for read in parton information
572  int idT, statusT, mother1T, mother2T, col1T, col2T;
573  double pxT, pyT, pzT, eT, mT;
574  // Leave tau and spin as default values
575  double tauT = 0., spinT = 9.;
576 
577  // Store particles locally at first so that resonances can be added
578  myParticles.clear();
579 
580  // Now read in partons
581  for (int i = 0; i < nParton; i++) {
582  // Get the next line
583  if (!getline(*isUnw, line)) {
584  cout << "Error in LHAupAlpgen::setEvent: "
585  << "could not read events from file" << endl;
586  return false;
587  }
588  istringstream iss2(line);
589 
590  // Incoming (flavour, colour, anticolour, pz)
591  if (i < 2) {
592  // Note that mothers will be set automatically by Pythia, and LHA
593  // status -1 is for an incoming parton
594  iss2 >> idT >> col1T >> col2T >> pzT;
595  statusT = -1;
596  mother1T = mother2T = 0;
597  pxT = pyT = mT = 0.;
598  eT = abs(pzT);
599 
600  // Adjust when zero pz/e
601  if (pzT == 0.) {
602  pzT = (i == 0) ? INCOMINGMIN : -INCOMINGMIN;
603  eT = INCOMINGMIN;
604  }
605 
606  // Outgoing (flavour, colour, anticolour, px, py, pz, mass)
607  } else {
608  // Note that mothers 1 and 2 corresport to the incoming partons,
609  // as set above, and LHA status +1 is for outgoing final state
610  iss2 >> idT >> col1T >> col2T >> pxT >> pyT >> pzT >> mT;
611  statusT = 1;
612  mother1T = 1;
613  mother2T = 2;
614  eT = sqrt(max(0., pxT*pxT + pyT*pyT + pzT*pzT + mT*mT));
615  }
616 
617  // Add particle
618  myParticles.push_back(LHAParticle(
619  idT, statusT, mother1T, mother2T, col1T, col2T,
620  pxT, pyT, pzT, eT, mT, tauT, spinT,-1.));
621  }
622 
623  // Add resonances if required
624  if (!addResonances()) return false;
625 
626  // Rescale momenta if required (must be done after full event
627  // reconstruction in addResonances)
628  if (!rescaleMomenta()) return false;
629 
630  // Pass particles on to Pythia
631  for (size_t i = 0; i < myParticles.size(); i++)
632  addParticle(myParticles[i]);
633 
634  // Set incoming flavour/x information and done
635  id1T = myParticles[0].idPart;
636  x1T = myParticles[0].ePart / ebmupA;
637  id2T = myParticles[1].idPart;
638  x2T = myParticles[1].ePart / ebmupA;
639  setIdX(id1T, id2T, x1T, x2T);
640  setPdf(id1T, id2T, x1T, x2T, 0., 0., 0., false);
641  return true;
642 }
643 
644 // ----------------------------------------------------------------------
645 
646 // Print list of particles; mainly intended for debugging
647 
649 
650  cout << endl << "---- LHAupAlpgen particle listing begin ----" << endl;
651  cout << scientific << setprecision(6);
652  for (int i = 0; i < int(myParticles.size()); i++) {
653  cout << setw(5) << i
654  << setw(5) << myParticles[i].idPart
655  << setw(5) << myParticles[i].statusPart
656  << setw(15) << myParticles[i].pxPart
657  << setw(15) << myParticles[i].pyPart
658  << setw(15) << myParticles[i].pzPart
659  << setw(15) << myParticles[i].ePart
660  << setw(15) << myParticles[i].mPart
661  << setw(5) << myParticles[i].mother1Part - 1
662  << setw(5) << myParticles[i].mother2Part - 1
663  << setw(5) << myParticles[i].col1Part
664  << setw(5) << myParticles[i].col2Part
665  << endl;
666  }
667  cout << "---- LHAupAlpgen particle listing end ----" << endl;
668 }
669 
670 // ----------------------------------------------------------------------
671 
672 // Routine to add resonances to an incoming event based on the
673 // hard process code (now stored in lprup).
674 
675 inline bool LHAupAlpgen::addResonances() {
676 
677  // Temporary storage for resonance information
678  int idT, statusT, mother1T, mother2T, col1T, col2T;
679  double pxT, pyT, pzT, eT, mT;
680  // Leave tau and spin as default values
681  double tauT = 0., spinT = 9.;
682 
683  // Alpgen process dependent parts. Processes:
684  // 1 = W + QQbar + jets
685  // 2 = Z/gamma* + QQbar + jets
686  // 3 = W + jets
687  // 4 = Z/gamma* + jets
688  // 10 = W + c + jets
689  // 14 = W + gamma + jets
690  // 15 = W + QQbar + gamma + jets
691  // When QQbar = ttbar, tops are not decayed in these processes.
692  // Explicitly reconstruct W/Z resonances; assumption is that the
693  // decay products are the last two particles.
694  if (lprup <= 4 || lprup == 10 || lprup == 14 || lprup == 15) {
695  // Decay products are the last two entries
696  int i1 = myParticles.size() - 1, i2 = i1 - 1;
697 
698  // Follow 'alplib/alpsho.f' procedure to get ID
699  if (myParticles[i1].idPart + myParticles[i2].idPart == 0)
700  idT = 0;
701  else
702  idT = - (myParticles[i1].idPart % 2) - (myParticles[i2].idPart % 2);
703  idT = (idT > 0) ? 24 : (idT < 0) ? -24 : 23;
704 
705  // Check that we get the expected resonance type; Z/gamma*
706  if (lprup == 2 || lprup == 4) {
707  if (idT != 23) {
708  cout << "Error in "
709  << "LHAupAlpgen::addResonances: wrong resonance type in event"
710  << endl;
711  return false;
712  }
713 
714  // W's
715  } else {
716  if (abs(idT) != 24) {
717  cout << "Error in "
718  << "LHAupAlpgen::addResonances: wrong resonance type in event"
719  << endl;
720  return false;
721  }
722  }
723 
724  // Remaining input
725  statusT = 2;
726  mother1T = 1;
727  mother2T = 2;
728  col1T = col2T = 0;
729  pxT = myParticles[i1].pxPart + myParticles[i2].pxPart;
730  pyT = myParticles[i1].pyPart + myParticles[i2].pyPart;
731  pzT = myParticles[i1].pzPart + myParticles[i2].pzPart;
732  eT = myParticles[i1].ePart + myParticles[i2].ePart;
733  mT = sqrt(eT*eT - pxT*pxT - pyT*pyT - pzT*pzT);
734  myParticles.push_back(LHAParticle(
735  idT, statusT, mother1T, mother2T, col1T, col2T,
736  pxT, pyT, pzT, eT, mT, tauT, spinT, -1.));
737 
738  // Update decay product mothers (note array size as if from 1)
739  myParticles[i1].mother1Part = myParticles[i2].mother1Part =
740  myParticles.size();
741  myParticles[i1].mother2Part = myParticles[i2].mother2Part = 0;
742 
743  // Processes:
744  // 5 = nW + mZ + j gamma + lH + jets
745  // 6 = QQbar + jets (QQbar = ttbar)
746  // 8 = QQbar + Higgs + jets (QQbar = ttbar)
747  // 13 = top + q (topprc = 1)
748  // 13 = top + b (topprc = 2)
749  // 13 = top + W + jets (topprc = 3)
750  // 13 = top + W + b (topprc = 4)
751  // 16 = QQbar + gamma + jets (QQbar = ttbar)
752  //
753  // When tops are present, they are decayed to Wb (both the W and b
754  // are not given), with this W also decaying (decay products given).
755  // The top is marked intermediate, the (intermediate) W is
756  // reconstructed from its decay products, and the decay product mothers
757  // updated. The final-state b is reconstructed from (top - W).
758  //
759  // W/Z resonances are given, as well as their decay products. The
760  // W/Z is marked intermediate, and the decay product mothers updated.
761  //
762  // It is always assumed that decay products are at the end.
763  // For processes 5 and 13, it is also assumed that the decay products
764  // are in the same order as the resonances.
765  // For processes 6, 8 and 16, the possibility of the decay products
766  // being out-of-order is also taken into account.
767  } else if ( ((lprup == 6 || lprup == 8 || lprup == 16) && ihvy1 == 6) ||
768  lprup == 5 || lprup == 13) {
769 
770  // Go backwards through the particles looking for W/Z/top
771  int idx = myParticles.size() - 1;
772  for (int i = myParticles.size() - 1; i > -1; i--) {
773 
774  // W or Z
775  if (myParticles[i].idPart == 23 ||
776  abs(myParticles[i].idPart) == 24) {
777 
778  // Check that decay products and resonance match up
779  int flav;
780  if (myParticles[idx].idPart + myParticles[idx - 1].idPart == 0)
781  flav = 0;
782  else
783  flav = - (myParticles[idx].idPart % 2)
784  - (myParticles[idx - 1].idPart % 2);
785  flav = (flav > 0) ? 24 : (flav < 0) ? -24 : 23;
786  if (flav != myParticles[i].idPart) {
787  if (infoPtr)
788  loggerPtr->ERROR_MSG("resonance does not match decay products");
789  return false;
790  }
791 
792  // Update status/mothers
793  myParticles[i].statusPart = 2;
794  myParticles[idx ].mother1Part = i + 1;
795  myParticles[idx--].mother2Part = 0;
796  myParticles[idx ].mother1Part = i + 1;
797  myParticles[idx--].mother2Part = 0;
798 
799  // Top
800  } else if (abs(myParticles[i].idPart) == 6) {
801 
802  // Check that decay products and resonance match up
803  int flav;
804  if (myParticles[idx].idPart + myParticles[idx - 1].idPart == 0)
805  flav = 0;
806  else
807  flav = - (myParticles[idx].idPart % 2)
808  - (myParticles[idx - 1].idPart % 2);
809  flav = (flav > 0) ? 24 : (flav < 0) ? -24 : 23;
810 
811  bool outOfOrder = false, wrongFlavour = false;;
812  if ( abs(flav) != 24 ||
813  (flav == 24 && myParticles[i].idPart != 6) ||
814  (flav == -24 && myParticles[i].idPart != -6) ) {
815 
816  // Processes 5 and 13, order should always be correct
817  if (lprup == 5 || lprup == 13) {
818  wrongFlavour = true;
819 
820  // Processes 6, 8 and 16, can have out of order decay products
821  } else {
822 
823  // Go back two decay products and retry
824  idx -= 2;
825  if (myParticles[idx].idPart + myParticles[idx - 1].idPart == 0)
826  flav = 0;
827  else
828  flav = - (myParticles[idx].idPart % 2)
829  - (myParticles[idx - 1].idPart % 2);
830  flav = (flav > 0) ? 24 : (flav < 0) ? -24 : 23;
831 
832  // If still the wrong flavour then error
833  if ( abs(flav) != 24 ||
834  (flav == 24 && myParticles[i].idPart != 6) ||
835  (flav == -24 && myParticles[i].idPart != -6) )
836  wrongFlavour = true;
837  else outOfOrder = true;
838  }
839 
840  // Error if wrong flavour
841  if (wrongFlavour) {
842  if (infoPtr)
843  loggerPtr->ERROR_MSG("resonance does not match decay products");
844  return false;
845  }
846  }
847 
848  // Mark t/tbar as now intermediate
849  myParticles[i].statusPart = 2;
850 
851  // New intermediate W+/W-
852  idT = flav;
853  statusT = 2;
854  mother1T = i + 1;
855  mother2T = 0;
856  col1T = col2T = 0;
857  pxT = myParticles[idx].pxPart + myParticles[idx - 1].pxPart;
858  pyT = myParticles[idx].pyPart + myParticles[idx - 1].pyPart;
859  pzT = myParticles[idx].pzPart + myParticles[idx - 1].pzPart;
860  eT = myParticles[idx].ePart + myParticles[idx - 1].ePart;
861  mT = sqrt(eT*eT - pxT*pxT - pyT*pyT - pzT*pzT);
862  myParticles.push_back(LHAParticle(
863  idT, statusT, mother1T, mother2T, col1T, col2T,
864  pxT, pyT, pzT, eT, mT, tauT, spinT, -1.));
865 
866  // Update the decay product mothers
867  myParticles[idx ].mother1Part = myParticles.size();
868  myParticles[idx--].mother2Part = 0;
869  myParticles[idx ].mother1Part = myParticles.size();
870  myParticles[idx--].mother2Part = 0;
871 
872  // New final-state b/bbar
873  idT = (flav == 24) ? 5 : -5;
874  statusT = 1;
875  // Colour from top
876  col1T = myParticles[i].col1Part;
877  col2T = myParticles[i].col2Part;
878  // Momentum from (t/tbar - W+/W-)
879  pxT = myParticles[i].pxPart - myParticles.back().pxPart;
880  pyT = myParticles[i].pyPart - myParticles.back().pyPart;
881  pzT = myParticles[i].pzPart - myParticles.back().pzPart;
882  eT = myParticles[i].ePart - myParticles.back().ePart;
883  mT = sqrt(eT*eT - pxT*pxT - pyT*pyT - pzT*pzT);
884  myParticles.push_back(LHAParticle(
885  idT, statusT, mother1T, mother2T, col1T, col2T,
886  pxT, pyT, pzT, eT, mT, tauT, spinT, -1.));
887 
888  // If decay products were out of order, reset idx to point
889  // at correct decay products
890  if (outOfOrder) idx += 4;
891 
892  } // if (abs(myParticles[i].idPart) == 6)
893  } // for (i)
894 
895 
896  // Processes:
897  // 7 = QQbar + Q'Qbar' + jets (tops are not decayed)
898  // 9 = jets
899  // 11 = gamma + jets
900  // 12 = Higgs + jets
901  } else if (lprup == 7 || lprup == 9 || lprup == 11 || lprup == 12) {
902  // Nothing to do for these processes
903  }
904 
905  // For single top, set incoming b mass if necessary
906  if (lprup == 13) for (int i = 0; i < 2; i++)
907  if (abs(myParticles[i].idPart) == 5) {
908  myParticles[i].mPart = mb;
909  myParticles[i].ePart = sqrt(pow2(myParticles[i].pzPart) + pow2(mb));
910  }
911 
912  // Debug output and done.
913  if (LHADEBUG) printParticles();
914  return true;
915 
916 }
917 
918 // ----------------------------------------------------------------------
919 
920 // Routine to rescale momenta to remove any imbalances. The routine
921 // assumes that any imbalances are due to decimal output/rounding
922 // effects, and are therefore small.
923 //
924 // First any px/py imbalances are fixed by adjusting all outgoing
925 // particles px/py and also updating their energy so mass is fixed.
926 // Because incoming pT is zero, changes should be limited to ~0.001.
927 //
928 // Second, any pz/e imbalances are fixed by scaling the incoming beams
929 // (again, no changes to masses required). Because incoming pz/e is not
930 // zero, effects can be slightly larger ~0.002/0.003.
931 
932 inline bool LHAupAlpgen::rescaleMomenta() {
933 
934  // Total momenta in/out
935  int nOut = 0;
936  Vec4 pIn, pOut;
937  for (int i = 0; i < int(myParticles.size()); i++) {
938  Vec4 pNow = Vec4(myParticles[i].pxPart, myParticles[i].pyPart,
939  myParticles[i].pzPart, myParticles[i].ePart);
940  if (i < 2) pIn += pNow;
941  else if (myParticles[i].statusPart == 1) {
942  nOut++;
943  pOut += pNow;
944  }
945  }
946 
947  // pT out to match pT in. Split any imbalances over all outgoing
948  // particles, and scale energies also to keep m^2 fixed.
949  if (abs(pOut.pT() - pIn.pT()) > ZEROTHRESHOLD) {
950  // Differences in px/py
951  double pxDiff = (pOut.px() - pIn.px()) / nOut,
952  pyDiff = (pOut.py() - pIn.py()) / nOut;
953 
954  // Warn if resulting changes above warning threshold
955  if (pxDiff > PTWARNTHRESHOLD || pyDiff > PTWARNTHRESHOLD) {
956  cout << "Warning in LHAupAlpgen::setEvent: "
957  << "large pT imbalance in incoming event" << endl;
958 
959  // Debug printout
960  if (LHADEBUGRESCALE) {
961  printParticles();
962  cout << "pxDiff = " << pxDiff << ", pyDiff = " << pyDiff << endl;
963  }
964  }
965 
966  // Adjust all final-state outgoing
967  pOut.reset();
968  for (int i = 2; i < int(myParticles.size()); i++) {
969  if (myParticles[i].statusPart != 1) continue;
970  myParticles[i].pxPart -= pxDiff;
971  myParticles[i].pyPart -= pyDiff;
972  myParticles[i].ePart = sqrt(max(0., pow2(myParticles[i].pxPart) +
973  pow2(myParticles[i].pyPart) + pow2(myParticles[i].pzPart) +
974  pow2(myParticles[i].mPart)));
975  pOut += Vec4(myParticles[i].pxPart, myParticles[i].pyPart,
976  myParticles[i].pzPart, myParticles[i].ePart);
977  }
978  }
979 
980  // Differences in E/pZ and scaling factors
981  double de = (pOut.e() - pIn.e());
982  double dp = (pOut.pz() - pIn.pz());
983  double a = 1 + (de + dp) / 2. / myParticles[0].ePart;
984  double b = 1 + (de - dp) / 2. / myParticles[1].ePart;
985 
986  // Warn if resulting energy changes above warning threshold.
987  // Change in pz less than or equal to change in energy (incoming b
988  // quark can have mass at this stage for process 13). Note that for
989  // very small incoming momenta, the relative adjustment may be large,
990  // but still small in absolute terms.
991  if (abs(a - 1.) * myParticles[0].ePart > EWARNTHRESHOLD ||
992  abs(b - 1.) * myParticles[1].ePart > EWARNTHRESHOLD) {
993  cout << "Warning in LHAupAlpgen::setEvent: "
994  << "large rescaling factor" << endl;
995 
996  // Debug printout
997  if (LHADEBUGRESCALE) {
998  printParticles();
999  cout << "de = " << de << ", dp = " << dp
1000  << ", a = " << a << ", b = " << b << endl
1001  << "Absolute energy change for incoming 0 = "
1002  << abs(a - 1.) * myParticles[0].ePart << endl
1003  << "Absolute energy change for incoming 1 = "
1004  << abs(b - 1.) * myParticles[1].ePart << endl;
1005  }
1006  }
1007  myParticles[0].ePart *= a;
1008  myParticles[0].pzPart *= a;
1009  myParticles[1].ePart *= b;
1010  myParticles[1].pzPart *= b;
1011 
1012  // Recalculate resonance four vectors
1013  for (int i = 0; i < int(myParticles.size()); i++) {
1014  if (myParticles[i].statusPart != 2) continue;
1015 
1016  // Only mothers stored in LHA, so go through all
1017  Vec4 resVec;
1018  for (int j = 0; j < int(myParticles.size()); j++) {
1019  if (myParticles[j].mother1Part - 1 != i) continue;
1020  resVec += Vec4(myParticles[j].pxPart, myParticles[j].pyPart,
1021  myParticles[j].pzPart, myParticles[j].ePart);
1022  }
1023 
1024  myParticles[i].pxPart = resVec.px();
1025  myParticles[i].pyPart = resVec.py();
1026  myParticles[i].pzPart = resVec.pz();
1027  myParticles[i].ePart = resVec.e();
1028  }
1029 
1030  return true;
1031 }
1032 
1033 //==========================================================================
1034 
1035 // Main implementation of AlpgenHooks class.
1036 // This may be split out to a separate C++ file if desired,
1037 // but currently included here for ease of use.
1038 
1039 // ----------------------------------------------------------------------
1040 
1041 // Constructor: provides the 'Alpgen:file' option by directly
1042 // changing the Pythia 'Beams' settings
1043 
1045 
1046  // If LHAupAlpgen needed, construct and pass to Pythia
1047  string agFile = pythia.settings.word("Alpgen:file");
1048  if (agFile != "void") {
1049  LHAagPtr = make_shared<LHAupAlpgen>(agFile.c_str());
1050  pythia.settings.mode("Beams:frameType", 5);
1051  pythia.setLHAupPtr(LHAagPtr);
1052  }
1053 }
1054 
1055 // ----------------------------------------------------------------------
1056 
1057 // Initialisation routine which is called by pythia.init().
1058 // This happens after the local pointers have been assigned and after
1059 // Pythia has processed the Beam information (and therefore LHEF header
1060 // information has been read in), but before any other internal
1061 // initialisation. Provides the remaining 'Alpgen:*' options.
1062 
1064 
1065  // Read in ALPGEN specific configuration variables
1066  bool setLightMasses = settingsPtr->flag("Alpgen:setLightMasses");
1067  bool setHeavyMasses = settingsPtr->flag("Alpgen:setHeavyMasses");
1068  bool setNjet = settingsPtr->flag("Alpgen:setNjet");
1069  bool setMLM = settingsPtr->flag("Alpgen:setMLM");
1070 
1071  // If ALPGEN parameters are present, then parse in AlpgenPar object
1072  AlpgenPar par;
1073  string parStr = infoPtr->header("AlpgenPar");
1074  if (!parStr.empty()) {
1075  par.parse(parStr);
1076  par.printParams();
1077  }
1078 
1079  // Set masses if requested
1080  if (setLightMasses) {
1081  if (par.haveParam("mc")) particleDataPtr->m0(4, par.getParam("mc"));
1082  if (par.haveParam("mb")) particleDataPtr->m0(5, par.getParam("mb"));
1083  }
1084  if (setHeavyMasses) {
1085  if (par.haveParam("mt")) particleDataPtr->m0(6, par.getParam("mt"));
1086  if (par.haveParam("mz")) particleDataPtr->m0(23, par.getParam("mz"));
1087  if (par.haveParam("mw")) particleDataPtr->m0(24, par.getParam("mw"));
1088  if (par.haveParam("mh")) particleDataPtr->m0(25, par.getParam("mh"));
1089  }
1090 
1091  // Set MLM:nJets if requested
1092  if (setNjet) {
1093  if (par.haveParam("njets"))
1094  settingsPtr->mode("JetMatching:nJet", par.getParamAsInt("njets"));
1095  else
1096  cout << "Warning in AlpgenHooks:init: "
1097  << "no ALPGEN nJet parameter found" << endl;
1098  }
1099 
1100  // Set MLM merging parameters if requested
1101  if (setMLM) {
1102  if (par.haveParam("ptjmin") && par.haveParam("drjmin") &&
1103  par.haveParam("etajmax")) {
1104  double ptjmin = par.getParam("ptjmin");
1105  ptjmin = max(ptjmin + 5., 1.2 * ptjmin);
1106  settingsPtr->parm("JetMatching:eTjetMin", ptjmin);
1107  settingsPtr->parm("JetMatching:coneRadius", par.getParam("drjmin"));
1108  settingsPtr->parm("JetMatching:etaJetMax", par.getParam("etajmax"));
1109 
1110  // Warn if setMLM requested, but parameters not present
1111  } else {
1112  cout << "Warning in AlpgenHooks:init: "
1113  << "no ALPGEN merging parameters found" << endl;
1114  }
1115  }
1116 
1117  // Initialisation complete.
1118  return true;
1119 }
1120 
1121 //==========================================================================
1122 
1123 // Main implementation of MadgraphPar class.
1124 // This may be split out to a separate C++ file if desired,
1125 // but currently included here for ease of use.
1126 
1127 //--------------------------------------------------------------------------
1128 
1129 // Constants: could be changed here if desired, but normally should not.
1130 // These are of technical nature, as described for each.
1131 
1132 // A zero threshold value for double comparisons.
1133 const double MadgraphPar::ZEROTHRESHOLD = 1e-10;
1134 
1135 //--------------------------------------------------------------------------
1136 
1137 // Parse an incoming Madgraph parameter file string
1138 
1139 inline bool MadgraphPar::parse(const string paramStr) {
1140 
1141  // Loop over incoming lines
1142  stringstream paramStream(paramStr);
1143  string line;
1144  while ( getline(paramStream, line) ) extractRunParam(line);
1145  return true;
1146 
1147 }
1148 
1149 //--------------------------------------------------------------------------
1150 
1151 // Parse an incoming parameter line
1152 
1153 inline void MadgraphPar::extractRunParam(string line) {
1154 
1155  // Extract information to the right of the final '!' character
1156  size_t idz = line.find("#");
1157  if ( !(idz == string::npos) ) return;
1158  size_t idx = line.find("=");
1159  size_t idy = line.find("!");
1160  if (idy == string::npos) idy = line.size();
1161  if (idx == string::npos) return;
1162  string paramName = trim( line.substr( idx + 1, idy - idx - 1) );
1163  string paramVal = trim( line.substr( 0, idx) );
1164  replace( paramVal.begin(), paramVal.end(), 'd', 'e');
1165 
1166  // Simple tokeniser
1167  istringstream iss(paramVal);
1168  double val;
1169  if (paramName.find(",") != string::npos) {
1170  string paramNameNow;
1171  istringstream issName( paramName);
1172  while ( getline(issName, paramNameNow, ',') ) {
1173  iss >> val;
1174  warnParamOverwrite( paramNameNow, val);
1175  params[paramNameNow] = val;
1176  }
1177 
1178  // Default case: assume integer and double on the left
1179  } else {
1180  iss >> val;
1181  warnParamOverwrite( paramName, val);
1182  params[paramName] = val;
1183  }
1184 }
1185 
1186 //--------------------------------------------------------------------------
1187 
1188 // Print parameters read from the '.par' file
1189 
1191 
1192  // Loop over all stored parameters and print
1193  cout << endl
1194  << " *-------- Madgraph parameters --------*" << endl;
1195  for (map<string,double>::iterator it = params.begin();
1196  it != params.end(); ++it)
1197  cout << " | " << left << setw(15) << it->first
1198  << " | " << right << setw(15) << it->second
1199  << " |" << endl;
1200  cout << " *---------------------------------------*" << endl;
1201 }
1202 
1203 //--------------------------------------------------------------------------
1204 
1205 // Warn if a parameter is going to be overwriten
1206 
1207 inline void MadgraphPar::warnParamOverwrite(const string &paramIn,
1208  double val) {
1209 
1210  // Check if present and if new value is different
1211  if (haveParam(paramIn) &&
1212  abs(getParam(paramIn) - val) > ZEROTHRESHOLD) {
1213  cout << "Warning in LHAupAlpgen::"
1214  << "warnParamOverwrite: overwriting existing parameter"
1215  << paramIn << endl;
1216  }
1217 }
1218 
1219 //--------------------------------------------------------------------------
1220 
1221 // Simple string trimmer
1222 
1223 inline string MadgraphPar::trim(string s) {
1224 
1225  // Remove whitespace in incoming string
1226  size_t i;
1227  if ( (i = s.find_last_not_of(" \t\r\n")) != string::npos)
1228  s = s.substr(0, i + 1);
1229  if ( (i = s.find_first_not_of(" \t\r\n")) != string::npos)
1230  s = s.substr(i);
1231  return s;
1232 }
1233 
1234 //==========================================================================
1235 
1236 } // end namespace Pythia8
1237 
1238 #endif // Pythia8_GeneratorInput_H
void setBeamA(int idIn, double eIn, int pdfGroupIn=0, int pdfSetIn=0)
Input beam info.
Definition: LesHouches.h:222
void addProcess(int idProcIn, double xSecIn=1., double xErrIn=0., double xMaxIn=1.)
Input process info.
Definition: LesHouches.h:233
constexpr double pow2(const double &x)
Powers of small integers - for balance speed/code clarity.
Definition: PythiaStdlib.h:182
void printParams()
Print parameters read from the &#39;.par&#39; file.
Definition: GeneratorInput.h:315
void setStrategy(int strategyIn)
Input process weight strategy.
Definition: LesHouches.h:230
istream * openFile(const char *fn, ifstream &ifs)
Definition: LesHouches.cc:648
bool parse(const string paramStr)
Parse as incoming ALPGEN parameter file (passed as string)
Definition: GeneratorInput.h:218
void setIdX(int id1In, int id2In, double x1In, double x2In)
Input info on flavour and x values of hard-process initiators.
Definition: LesHouches.h:263
A class for the particles stored in LHAup.
Definition: LesHouches.h:48
Settings settings
Settings: databases of flags/modes/parms/words to control run.
Definition: Pythia.h:332
bool initAfterBeams()
Override initAfterBeams routine from UserHooks.
Definition: GeneratorInput.h:1063
AlpgenHooks(Pythia &pythia)
Constructor and destructor.
Definition: GeneratorInput.h:1044
bool setLHAupPtr(LHAupPtr lhaUpPtrIn)
Possibility to pass in pointer to external LHA-interfaced generator.
Definition: Pythia.h:133
LHAupAlpgen(const char *baseFNin)
Constructor and destructor.
Definition: GeneratorInput.h:388
string header(const string &key) const
Return an LHEF header.
Definition: Info.h:353
double getParam(const string &paramIn)
Definition: GeneratorInput.h:54
bool setInit()
Override setInit/setEvent routines from LHAup.
Definition: GeneratorInput.h:413
Definition: LesHouches.h:78
void addParticle(LHAParticle particleIn)
Input particle info, one particle at the time.
Definition: LesHouches.h:252
bool haveParam(const string &paramIn)
Check if a parameter exists.
Definition: GeneratorInput.h:49
void setProcess(int idProcIn=0, double weightIn=1., double scaleIn=0., double alphaQEDIn=0.0073, double alphaQCDIn=0.12)
Input info on the selected process.
Definition: LesHouches.h:243
AlpgenPar()
Constructor.
Definition: GeneratorInput.h:40
bool setEvent(int)
Definition: GeneratorInput.h:543
void setPdf(int id1pdfIn, int id2pdfIn, double x1pdfIn, double x2pdfIn, double scalePDFIn, double pdf1In, double pdf2In, bool pdfIsSetIn)
Optionally input info on parton density values of event.
Definition: LesHouches.h:267
Definition: GeneratorInput.h:133
bool parse(const string paramStr)
Parse an incoming Madgraph parameter file string.
Definition: GeneratorInput.h:1139
Definition: GeneratorInput.h:35
void extractRunParam(string line)
Parse an incoming parameter line.
Definition: GeneratorInput.h:1153
void printParams()
Print parameters read from the &#39;.par&#39; file.
Definition: GeneratorInput.h:1190
Definition: GeneratorInput.h:83
void printParticles()
Print list of particles; mainly intended for debugging.
Definition: GeneratorInput.h:648
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
The Pythia class contains the top-level routines to generate an event.
Definition: Pythia.h:71
void setInfoHeader(const string &key, const string &val)
Definition: LesHouches.h:294
Info * infoPtr
Pointer to various information on the generation.
Definition: LesHouches.h:218
UserHooks is base class for user access to program execution.
Definition: UserHooks.h:32
bool haveParam(const string &paramIn)
Check if a parameter exists.
Definition: GeneratorInput.h:170
Definition: GeneratorInput.h:156
MadgraphPar()
Constructor.
Definition: GeneratorInput.h:161
double getParam(const string &paramIn)
Definition: GeneratorInput.h:175
void closeFile(istream *&is, ifstream &ifs)
Definition: LesHouches.cc:659
void reset()
Member functions for input.
Definition: Basics.h:46
bool fileFound()
Override fileFound routine from LHAup.
Definition: GeneratorInput.h:92
void extractRunParam(string line)
Parse an incoming parameter line.
Definition: GeneratorInput.h:256
Definition: Basics.h:32