PYTHIA  8.312
ExternalMEsMadgraph.h
1 // ExternalMEsMadgraph.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2024 Peter Skands, Stefan Prestel, Philip Ilten, Torbjorn
3 // Sjostrand.
4 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
5 // Please respect the MCnet Guidelines, see GUIDELINES for details.
6 
7 // This file contains the Madgraph parton shower ME plugin class which
8 // interfaces with matrix elements generated with the
9 // PY8Kernels/MG5MES plugin to MadGraph 5.
10 
11 #ifndef Pythia8_ExternalMEsMadgraph_H
12 #define Pythia8_ExternalMEsMadgraph_H
13 
14 // Include Pythia headers.
15 #include "Pythia8/ExternalMEs.h"
16 #include "Pythia8/Plugins.h"
17 
18 // Include Madgraph PY8MEs plugin headers.
19 #include "PY8ME.h"
20 #include "PY8MEs.h"
21 
22 namespace Pythia8 {
23 
24 //==========================================================================
25 
26 // External matrix element class specifically for MadGraph.
27 
29 
30 public:
31 
32  // Constructor.
34  isInitPtr = false; isInit = false; libPtr = nullptr; modelPtr = nullptr;}
35 
36  // Destructor.
37  ~ExternalMEsMadgraph() {if (libPtr != nullptr) delete libPtr;
38  if (modelPtr != nullptr) delete modelPtr;}
39 
40  // Initialisers.
41  bool init() override;
42  bool initVincia(Info* infoPtrIn) override;
43  bool initDire(Info*, string card) override;
44 
45  // Methods to check availability of matrix elements.
46  bool isAvailable(vector<int> idIn, vector<int> idOut) override;
47  bool isAvailable(const Pythia8::Event& event) override;
48  bool isAvailable(const Pythia8::Event& event, int iBeg = 3) override;
49  bool isAvailable(const vector<Particle>& state) override;
50 
51  // Get the matrix element squared for a particle state.
52  double calcME2(const vector<Particle>& state) override;
53  double calcME2(const Event& event, const int iBeg = 3) override;
54 
55 private:
56 
57  // Fill lists of IDs, momenta, colours, and helicities in MG5 format.
58  void fillLists(const vector<Particle>& state, vector<int>& idsIn,
59  vector<int>& idsOut, vector<int>& hels, vector<int>& cols,
60  vector<vector<double>>& pn) const;
61 
62  // Calculate ME2 from pre-set lists.
63  double calcME2(vector<int>& idIn, vector<int>& idOut,
64  vector< vector<double> >& pn, vector<int>& hels, vector<int>& cols,
65  set<int> sChannels = {});
66 
67  PY8MEs_namespace::PY8MEs* libPtr;
68  PARS* modelPtr;
69 
70 };
71 
72 //--------------------------------------------------------------------------
73 
74 // Initialise the Madgraph model, parameters, and couplings for use in Vincia.
75 
76 bool ExternalMEsMadgraph::init() {return true;}
77 
79 
80  // Check if pointers initialized.
81  initPtrs(infoPtrIn);
82  int verbose = settingsPtr->mode("Vincia:verbose");
83  if (verbose > 1)
84  cout << " (ExternalMEsMadgraph::init()) begin -------" << endl;
85  if (!isInitPtr) {
86  loggerPtr->ERROR_MSG("cannot initialize, pointers not set");
87  return false;
88  }
89  isInit = true;
90 
91  // Create new model instance.
92  if (verbose >= 2) cout << " (ExternalMEsMadgraph::init())"
93  << "setting MG5 C++ masses, widths, couplings..." << endl;
94  if (modelPtr != nullptr) delete modelPtr;
95  modelPtr = new PARS();
96  modelPtr->setIndependentParameters(particleDataPtr,coupSMPtr,slhaPtr);
97  modelPtr->setIndependentCouplings();
98  if (verbose >= 3) {
99  modelPtr->printIndependentParameters();
100  modelPtr->printIndependentCouplings();
101  }
102 
103  // In the VINCIA context, alphaS_MGME = 1/4pi (- > gS = 1; we
104  // control the running separately). Thus, even the Madgraph "dependent"
105  // parameters only need to be set once.
106 
107  // Alternatively, we could evaluate the QCD coupling at MZ but should
108  // then use a coupling definition from a Vincia parameter list rather
109  // than PYTHIA's couplings.
110  // double muDefME = 91.2;
111  // double alpS = coupSMPtr->alphaS(pow2(muDefME));
112 
113  // The following is equivalent to
114  // PY8MEs::updateModelDependentCouplingsWithPY8(...)
115  double alpS = 1.0 / ( 4 * M_PI );
116  modelPtr->setDependentParameters(particleDataPtr, coupSMPtr, slhaPtr,
117  alpS);
118  modelPtr->setDependentCouplings();
119 
120  // Construct Madgraph process library.
121  if (verbose >= 3) cout << " (ExternalMEsMadgraph::init())"
122  << " attempting to construct lib" << endl;
123  if (libPtr != nullptr) delete libPtr;
124  libPtr = new PY8MEs_namespace::PY8MEs(modelPtr);
125  // Do not include averaging or symmetry factors in MG5.
126  libPtr->seProcessesIncludeSymmetryFactors(false);
127  libPtr->seProcessesIncludeHelicityAveragingFactors(false);
128  libPtr->seProcessesIncludeColorAveragingFactors(false);
129  // Set whether symmetry and averaging factors are applied in calcME2().
130  inclSymFac = false;
131  inclHelAvgFac = true;
132  inclColAvgFac = true;
133  // Leading-colour colour-ordered amplitude only (can be reset later).
134  colMode = 1;
135  // Implicitly sum over helicities (can be reset later).
136  helMode = 1;
137 
138  return true;
139 
140 }
141 
142 //--------------------------------------------------------------------------
143 
144 // Initialise the Madgraph model, parameters, and couplings for use in Dire.
145 
146 bool ExternalMEsMadgraph::initDire(Info*, string card) {
147 
148  // Redirect output so that Dire can print MG5 initialization.
149  std::streambuf *old = cout.rdbuf();
150  stringstream ss;
151  cout.rdbuf (ss.rdbuf());
152  if (libPtr != nullptr) delete libPtr;
153  libPtr = new PY8MEs_namespace::PY8MEs(card);
154  // Do not include averaging or symmetry factors in MG5.
155  libPtr->seProcessesIncludeSymmetryFactors(false);
156  libPtr->seProcessesIncludeHelicityAveragingFactors(false);
157  libPtr->seProcessesIncludeColorAveragingFactors(false);
158  libPtr->setProcessesExternalMassesMode(1);
159  // Set whether symmetry and averaging factors are applied in calcME2().
160  inclSymFac = false;
161  inclHelAvgFac = true;
162  inclColAvgFac = true;
163  // Leading-colour colour-ordered amplitude only (can be reset later).
164  colMode = 1;
165  // Implicitly sum over helicities (can be reset later).
166  helMode = 1;
167  // Restore print-out.
168  cout.rdbuf (old);
169 
170  return true;
171 
172 }
173 
174 //--------------------------------------------------------------------------
175 
176 // Check if a matrix element is available.
177 
178 bool ExternalMEsMadgraph::isAvailable(vector<int> idIn, vector<int> idOut) {
179  set<int> req_s_channels;
180  PY8MEs_namespace::PY8ME * query
181  = libPtr->getProcess(idIn, idOut, req_s_channels);
182  return (query != nullptr);
183 }
184 
185 bool ExternalMEsMadgraph::isAvailable(const Event& event) {
186 
187  vector <int> in, out;
188  fillIds(event, in, out);
189  set<int> req_s_channels;
190 
191  PY8MEs_namespace::PY8ME* query
192  = libPtr->getProcess(in, out, req_s_channels);
193  return (query != nullptr);
194 }
195 
196 bool ExternalMEsMadgraph::isAvailable(const Event& event, const int iBeg) {
197 
198  vector <int> in, out;
199  fillIds(event, in, out, iBeg);
200  set<int> req_s_channels;
201 
202  PY8MEs_namespace::PY8ME* query
203  = libPtr->getProcess(in, out, req_s_channels);
204  return (query != nullptr);
205 }
206 
207 bool ExternalMEsMadgraph::isAvailable(const vector<Particle>& state) {
208 
209  vector <int> idIn, idOut;
210  for (const Particle& ptcl : state) {
211  if (ptcl.isFinal()) idOut.push_back(ptcl.id());
212  else idIn.push_back(ptcl.id());
213  }
214  set<int> req_s_channels;
215 
216  PY8MEs_namespace::PY8ME* query
217  = libPtr->getProcess(idIn, idOut, req_s_channels);
218  return (query != nullptr);
219 }
220 
221 //--------------------------------------------------------------------------
222 
223 // Calcuate the squared matrix element.
224 
225 double ExternalMEsMadgraph::calcME2(const vector<Particle>& state) {
226 
227  // Prepare lists.
228  vector<int> idIn, idOut, hels, cols;
229  vector<vector<double>> pn;
230  fillLists(state, idIn, idOut, hels, cols, pn);
231  int nIn = idIn.size();
232  if (nIn <= 0 || state.size() - nIn < 1) return -1.;
233 
234  return calcME2(idIn, idOut, pn, hels, cols);
235 
236 }
237 
239  const int iBeg) {
240 
241  // Prepare lists.
242  vector<int> in, out;
243  fillIds(event, in, out, iBeg);
244  vector<int> cols;
245  fillCols(event, cols, iBeg);
246  vector< vector<double> > pvec = fillMoms(event, iBeg);
247  vector<int> helicities;
248 
249  return calcME2(in, out, pvec, helicities, cols);
250 
251 }
252 
253 double ExternalMEsMadgraph::calcME2(vector<int>& idIn, vector<int>& idOut,
254  vector< vector<double> >& pn, vector<int>& hels, vector<int>& cols,
255  set<int> sChannels) {
256 
257  // Access the process.
258  PY8MEs_namespace::process_specifier proc_spec =
259  libPtr->getProcessSpecifier(idIn, idOut, sChannels);
260  PY8MEs_namespace::process_accessor proc_handle =
261  libPtr->getProcess(proc_spec);
262  // Return right away if unavailable.
263  if (proc_handle.second.second < 0) return 0.0;
264  PY8MEs_namespace::PY8ME* proc_ptr = proc_handle.first;
265  proc_ptr->setPermutation(proc_handle.second.first);
266  int procID = proc_handle.second.second;
267  proc_ptr->setProcID(procID);
268  int nIn = idIn.size();
269 
270  // Save current helicity configuration.
271  vector< vector<int> > helConf;
272  helConf.push_back(hels);
273  // Check if we are doing an explicit helicity sum.
274  vector<int> i9;
275  if (helMode == 0) {
276  // Save indices of unpolarised partons.
277  for (int i(0); i<(int)hels.size(); ++i)
278  if (hels[i] == 9) i9.push_back(i);
279  }
280  // Manually calculate helicity average.
281  double helAvgNorm = i9.size()>0 ? 1:proc_ptr->getHelicityAveragingFactor();
282  while (i9.size() > 0) {
283  int i = i9.back();
284  int id = (i < nIn) ? idIn[i] : idOut[i-nIn];
285  // How many spin states.
286  int nS = particleDataPtr->spinType(id);
287  // Massless particles have max 2 (physical) spin states.
288  if (particleDataPtr->m0(id) == 0.0) nS=min(nS,2);
289  if (i < nIn) helAvgNorm *= (double)nS;
290  // Create nS copies of helConf, one for each spin state.
291  int helConfSizeNow = helConf.size();
292  for (int iCopy = 1; iCopy <= nS; ++iCopy) {
293  // Set hel for this particle in this copy.
294  // Start from -1, then 1, then 0 (if 3 states).
295  int h = -1;
296  if (nS == 1) h = 0;
297  else if (iCopy == 2) h = 1;
298  else if (iCopy == 3) h = 0;
299  else if (iCopy == 4) h = -2;
300  else if (iCopy == 5) h = 2;
301  for (int iHelConf=0; iHelConf<helConfSizeNow; ++iHelConf) {
302  vector<int> helNow = helConf[iHelConf];
303  helNow[i] = h;
304  // First copy: use existing.
305  if (iCopy == 1) helConf[iHelConf] = helNow;
306  // Subsequent copies: create new.
307  else helConf.push_back(helNow);
308  }
309  }
310  // Remove the particle whose helicities have been summed over.
311  i9.pop_back();
312  }
313 
314  // Set colour configurations according to requested colour mode.
315  vector<vector<int>> colConf;
316  if (colMode >= 3) {
317  // For FC sum use empty colour vector to communicate it with MG5.
318  cols.clear();
319  colConf.push_back(cols);
320  } else if (colMode == 2) {
321  // For LC sum, fetch all LC colour configurations.
322  vector<vector<int>> colorConfigs = proc_ptr->getColorConfigs(procID);
323  for (const auto& cc : colorConfigs) {
324  int colID = proc_ptr->getColorIDForConfig(cc);
325  if (proc_ptr->getColorFlowRelativeNCPower(colID) == 0)
326  colConf.push_back(cc);
327  }
328  } else colConf.push_back(cols);
329 
330  // Compute sum over colours and helicities.
331  // (Save helicity components if needed later).
332  double me2 = 0.0;
333  me2hel.clear();
334  proc_ptr->setMomenta(pn);
335  for (const auto& cc : colConf) {
336  proc_ptr->setColors(cc);
337  for (int iHC(0); iHC<(int)helConf.size(); ++iHC) {
338  proc_ptr->setHelicities(helConf[iHC]);
339  double me2now = proc_ptr->sigmaKin();
340  // MG5 may produce inf/nan for unphysical hel combinations.
341  if (!isnan(me2now) && !isinf(me2now)) {
342  // Save helicity matrix element for possible later use.
343  me2hel[helConf[iHC]] = me2now;
344  me2 += me2now;
345  }
346  }
347  }
348 
349  // Potentially include symmetry and averaging factors.
350  if (inclSymFac) me2 /= proc_ptr->getSymmetryFactor();
351  if (inclHelAvgFac) me2 /= helAvgNorm;
352  if (inclColAvgFac) me2 /= proc_ptr->getColorAveragingFactor();
353  return me2;
354 
355 }
356 
357 //--------------------------------------------------------------------------
358 
359 // Fill lists.
360 
361 void ExternalMEsMadgraph::fillLists(const vector<Particle>& state,
362  vector<int>& idsIn, vector<int>& idsOut, vector<int>& hels,
363  vector<int>& cols, vector<vector<double>>& pn) const {
364  for (const Particle& ptcl : state) {
365  if (ptcl.isFinal()) idsOut.push_back(ptcl.id());
366  else idsIn.push_back(ptcl.id());
367  vector<double> pNow = {ptcl.e(), ptcl.px(), ptcl.py(), ptcl.pz()};
368  pn.push_back(pNow);
369  hels.push_back(ptcl.pol());
370  cols.push_back(ptcl.col());
371  cols.push_back(ptcl.acol());
372  }
373 }
374 
375 //--------------------------------------------------------------------------
376 
377 // Declare the plugin.
378 
379 PYTHIA8_PLUGIN_CLASS(ExternalMEs, ExternalMEsMadgraph, false, false, false)
380 PYTHIA8_PLUGIN_VERSIONS(PYTHIA_VERSION_INTEGER)
381 
382 //==========================================================================
383 
384 } // end namespace Pythia8
385 
386 #endif // end Pythia8_ExternalMEsMadgraph_H
double calcME2(const vector< Particle > &state) override
Get the matrix element squared for a particle state.
Definition: ExternalMEsMadgraph.h:225
int colMode
Colour mode (0: strict LC, 1: LC, 2: LC sum, 3: FC).
Definition: ExternalMEs.h:95
virtual void initPtrs(Info *infoPtrIn)
Initialisers for pointers.
Definition: ExternalMEs.cc:22
Definition: Info.h:45
int helMode
Helicity mode (0: explicit helicity sum, 1: implicit helicity sum).
Definition: ExternalMEs.h:101
The Event class holds all info on the generated event.
Definition: Event.h:453
void fillCols(const Event &event, vector< int > &colors, int iBeg=3) const
Fill a vector of colors, from an event, starting from entry i = iBeg.
Definition: ExternalMEs.cc:60
void fillMoms(const Event &event, vector< Vec4 > &p, int iBeg=3) const
Fill a vector of momenta, from an event, starting from entry i = iBeg.
Definition: ExternalMEs.cc:50
External matrix element class specifically for MadGraph.
Definition: ExternalMEsMadgraph.h:28
Definition: Logger.h:23
bool isAvailable(vector< int > idIn, vector< int > idOut) override
Methods to check availability of matrix elements.
Definition: ExternalMEsMadgraph.h:178
bool inclSymFac
Symmetry and averaging factors.
Definition: ExternalMEs.h:104
map< vector< int >, double > me2hel
Saved list of helicity components for last ME evaluated.
Definition: ExternalMEs.h:98
~ExternalMEsMadgraph()
Destructor.
Definition: ExternalMEsMadgraph.h:37
Base class for external matrix-element interfaces.
Definition: ExternalMEs.h:30
Definition: Event.h:32
bool init() override
Initialisers.
Definition: ExternalMEsMadgraph.h:76
bool initDire(Info *, string card) override
Initialise the Madgraph model, parameters, and couplings for use in Dire.
Definition: ExternalMEsMadgraph.h:146
bool initVincia(Info *infoPtrIn) override
Definition: ExternalMEsMadgraph.h:78
void fillIds(const Event &event, vector< int > &in, vector< int > &out, int iBeg=3) const
Fill a vector of IDs, from an event, starting from entry i = iBeg.
Definition: ExternalMEs.cc:37
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
bool isInitPtr
Is initialized.
Definition: ExternalMEs.h:116
Definition: Settings.h:195
ExternalMEsMadgraph(Pythia *, Settings *, Logger *)
Constructor.
Definition: ExternalMEsMadgraph.h:33