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