include/Pythia8Plugins/ExternalMEsMadgraph.h

0001 // ExternalMEsMadgraph.h is a part of the PYTHIA event generator.
0002 // Copyright (C) 2025 Peter Skands, Stefan Prestel, Philip Ilten, Torbjorn
0003 // Sjostrand.
0004 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
0005 // Please respect the MCnet Guidelines, see GUIDELINES for details.
0006
0007 // This file contains the Madgraph parton shower ME plugin class which
0008 // interfaces with matrix elements generated with the
0009 // PY8Kernels/MG5MES plugin to MadGraph 5.
0010
0011 #ifndef Pythia8_ExternalMEsMadgraph_H
0012 #define Pythia8_ExternalMEsMadgraph_H
0013
0014 // Include Pythia headers.
0015 #include "Pythia8/ExternalMEs.h"
0016 #include "Pythia8/Plugins.h"
0017
0018 // Include Madgraph PY8MEs plugin headers.
0019 #include "PY8ME.h"
0020 #include "PY8MEs.h"
0021
0022 namespace Pythia8 {
0023
0024 //==========================================================================
0025
0026 // External matrix element class specifically for MadGraph.
0027
0028 class ExternalMEsMadgraph : public ExternalMEs {
0029
0030 public:
0031
0032   // Constructor.
0033   ExternalMEsMadgraph(Pythia*, Settings*, Logger*) {
0034     isInitPtr = false; isInit = false; libPtr = nullptr; modelPtr = nullptr;}
0035
0036   // Destructor.
0037   ~ExternalMEsMadgraph() {if (libPtr != nullptr) delete libPtr;
0038     if (modelPtr != nullptr) delete modelPtr;}
0039
0040   // Initialisers.
0041   bool init() override;
0042   bool initVincia(Info* infoPtrIn) override;
0043   bool initDire(Info*, string card) override;
0044
0045   // Methods to check availability of matrix elements.
0046   bool isAvailable(vector<int> idIn, vector<int> idOut) override;
0047   bool isAvailable(const Pythia8::Event& event) override;
0048   bool isAvailable(const Pythia8::Event& event, int iBeg = 3) override;
0049   bool isAvailable(const vector<Particle>& state) override;
0050
0051   // Get the matrix element squared for a particle state.
0052   double calcME2(const vector<Particle>& state) override;
0053   double calcME2(const Event& event, const int iBeg = 3) override;
0054
0055 private:
0056
0057   // Fill lists of IDs, momenta, colours, and helicities in MG5 format.
0058   void fillLists(const vector<Particle>& state, vector<int>& idsIn,
0059     vector<int>& idsOut, vector<int>& hels, vector<int>& cols,
0060     vector<vector<double>>& pn) const;
0061
0062   // Calculate ME2 from pre-set lists.
0063   double calcME2(vector<int>& idIn, vector<int>& idOut,
0064     vector< vector<double> >& pn, vector<int>& hels, vector<int>& cols,
0065     set<int> sChannels = {});
0066
0067   PY8MEs_namespace::PY8MEs* libPtr;
0068   PARS* modelPtr;
0069
0070 };
0071
0072 //--------------------------------------------------------------------------
0073
0074 // Initialise the Madgraph model, parameters, and couplings for use in Vincia.
0075
0076 bool ExternalMEsMadgraph::init() {return true;}
0077
0078 bool ExternalMEsMadgraph::initVincia(Info* infoPtrIn) {
0079
0080   // Check if pointers initialized.
0081   initPtrs(infoPtrIn);
0082   int verbose = settingsPtr->mode("Vincia:verbose");
0083   if (verbose > 1)
0084     cout << " (ExternalMEsMadgraph::init()) begin -------" << endl;
0085   if (!isInitPtr) {
0086     loggerPtr->ERROR_MSG("cannot initialize, pointers not set");
0087     return false;
0088   }
0089   isInit = true;
0090
0091   // Create new model instance.
0092   if (verbose >= 2) cout << " (ExternalMEsMadgraph::init())"
0093     << "setting MG5 C++ masses, widths, couplings..." << endl;
0094   if (modelPtr != nullptr) delete modelPtr;
0095   modelPtr = new PARS();
0096   modelPtr->setIndependentParameters(particleDataPtr,coupSMPtr,slhaPtr);
0097   modelPtr->setIndependentCouplings();
0098   if (verbose >= 3) {
0099     modelPtr->printIndependentParameters();
0100     modelPtr->printIndependentCouplings();
0101   }
0102
0103   // In the VINCIA context, alphaS_MGME = 1/4pi (- > gS = 1; we
0104   // control the running separately). Thus, even the Madgraph "dependent"
0105   // parameters only need to be set once.
0106
0107   // Alternatively, we could evaluate the QCD coupling at MZ but should
0108   // then use a coupling definition from a Vincia parameter list rather
0109   // than PYTHIA's couplings.
0110   //    double muDefME  = 91.2;
0111   //    double alpS = coupSMPtr->alphaS(pow2(muDefME));
0112
0113   // The following is equivalent to
0114   // PY8MEs::updateModelDependentCouplingsWithPY8(...)
0115   double alpS = 1.0 / ( 4 * M_PI );
0116   modelPtr->setDependentParameters(particleDataPtr, coupSMPtr, slhaPtr,
0117     alpS);
0118   modelPtr->setDependentCouplings();
0119
0120   // Construct Madgraph process library.
0121   if (verbose >= 3) cout << " (ExternalMEsMadgraph::init())"
0122     << " attempting to construct lib" << endl;
0123   if (libPtr != nullptr) delete libPtr;
0124   libPtr = new PY8MEs_namespace::PY8MEs(modelPtr);
0125   // Do not include averaging or symmetry factors in MG5.
0126   libPtr->seProcessesIncludeSymmetryFactors(false);
0127   libPtr->seProcessesIncludeHelicityAveragingFactors(false);
0128   libPtr->seProcessesIncludeColorAveragingFactors(false);
0129   // Set whether symmetry and averaging factors are applied in calcME2().
0130   inclSymFac    = false;
0131   inclHelAvgFac = true;
0132   inclColAvgFac = true;
0133   // Leading-colour colour-ordered amplitude only (can be reset later).
0134   colMode = 1;
0135   // Implicitly sum over helicities (can be reset later).
0136   helMode = 1;
0137
0138   return true;
0139
0140 }
0141
0142 //--------------------------------------------------------------------------
0143
0144 // Initialise the Madgraph model, parameters, and couplings for use in Dire.
0145
0146 bool ExternalMEsMadgraph::initDire(Info*, string card) {
0147
0148   // Redirect output so that Dire can print MG5 initialization.
0149   std::streambuf *old = cout.rdbuf();
0150   stringstream ss;
0151   cout.rdbuf (ss.rdbuf());
0152   if (libPtr != nullptr) delete libPtr;
0153   libPtr = new PY8MEs_namespace::PY8MEs(card);
0154   // Do not include averaging or symmetry factors in MG5.
0155   libPtr->seProcessesIncludeSymmetryFactors(false);
0156   libPtr->seProcessesIncludeHelicityAveragingFactors(false);
0157   libPtr->seProcessesIncludeColorAveragingFactors(false);
0158   libPtr->setProcessesExternalMassesMode(1);
0159   // Set whether symmetry and averaging factors are applied in calcME2().
0160   inclSymFac    = false;
0161   inclHelAvgFac = true;
0162   inclColAvgFac = true;
0163   // Leading-colour colour-ordered amplitude only (can be reset later).
0164   colMode = 1;
0165   // Implicitly sum over helicities (can be reset later).
0166   helMode = 1;
0167   // Restore print-out.
0168   cout.rdbuf (old);
0169
0170   return true;
0171
0172 }
0173
0174 //--------------------------------------------------------------------------
0175
0176 // Check if a matrix element is available.
0177
0178 bool ExternalMEsMadgraph::isAvailable(vector<int> idIn, vector<int> idOut) {
0179   set<int> req_s_channels;
0180   PY8MEs_namespace::PY8ME * query
0181     = libPtr->getProcess(idIn, idOut, req_s_channels);
0182   return (query != nullptr);
0183 }
0184
0185 bool ExternalMEsMadgraph::isAvailable(const Event& event) {
0186
0187   vector <int> in, out;
0188   fillIds(event, in, out);
0189   set<int> req_s_channels;
0190
0191   PY8MEs_namespace::PY8ME* query
0192     = libPtr->getProcess(in, out, req_s_channels);
0193   return (query != nullptr);
0194 }
0195
0196 bool ExternalMEsMadgraph::isAvailable(const Event& event, const int iBeg) {
0197
0198   vector <int> in, out;
0199   fillIds(event, in, out, iBeg);
0200   set<int> req_s_channels;
0201
0202   PY8MEs_namespace::PY8ME* query
0203     = libPtr->getProcess(in, out, req_s_channels);
0204   return (query != nullptr);
0205 }
0206
0207 bool ExternalMEsMadgraph::isAvailable(const vector<Particle>& state) {
0208
0209   vector <int> idIn, idOut;
0210   for (const Particle& ptcl : state) {
0211     if (ptcl.isFinal()) idOut.push_back(ptcl.id());
0212     else idIn.push_back(ptcl.id());
0213   }
0214   set<int> req_s_channels;
0215
0216   PY8MEs_namespace::PY8ME* query
0217     = libPtr->getProcess(idIn, idOut, req_s_channels);
0218   return (query != nullptr);
0219 }
0220
0221 //--------------------------------------------------------------------------
0222
0223 // Calcuate the squared matrix element.
0224
0225 double ExternalMEsMadgraph::calcME2(const vector<Particle>& state) {
0226
0227   // Prepare lists.
0228   vector<int> idIn, idOut, hels, cols;
0229   vector<vector<double>> pn;
0230   fillLists(state, idIn, idOut, hels, cols, pn);
0231   int nIn = idIn.size();
0232   if (nIn <= 0 || state.size() - nIn < 1) return -1.;
0233
0234   return calcME2(idIn, idOut, pn, hels, cols);
0235
0236 }
0237
0238 double ExternalMEsMadgraph::calcME2(const Pythia8::Event& event,
0239   const int iBeg) {
0240
0241   // Prepare lists.
0242   vector<int> in, out;
0243   fillIds(event, in, out, iBeg);
0244   vector<int> cols;
0245   fillCols(event, cols, iBeg);
0246   vector< vector<double> > pvec = fillMoms(event, iBeg);
0247   vector<int> helicities;
0248
0249   return calcME2(in, out, pvec, helicities, cols);
0250
0251 }
0252
0253 double ExternalMEsMadgraph::calcME2(vector<int>& idIn, vector<int>& idOut,
0254   vector< vector<double> >& pn, vector<int>& hels, vector<int>& cols,
0255   set<int> sChannels) {
0256
0257   // Access the process.
0258   PY8MEs_namespace::process_specifier proc_spec =
0259     libPtr->getProcessSpecifier(idIn, idOut, sChannels);
0260   PY8MEs_namespace::process_accessor proc_handle =
0261     libPtr->getProcess(proc_spec);
0262   // Return right away if unavailable.
0263   if (proc_handle.second.second < 0) return 0.0;
0264   PY8MEs_namespace::PY8ME* proc_ptr = proc_handle.first;
0265   proc_ptr->setPermutation(proc_handle.second.first);
0266   int procID = proc_handle.second.second;
0267   proc_ptr->setProcID(procID);
0268   int nIn = idIn.size();
0269
0270   // Save current helicity configuration.
0271   vector< vector<int> > helConf;
0272   helConf.push_back(hels);
0273   // Check if we are doing an explicit helicity sum.
0274   vector<int> i9;
0275   if (helMode == 0) {
0276     // Save indices of unpolarised partons.
0277     for (int i(0); i<(int)hels.size(); ++i)
0278       if (hels[i] == 9) i9.push_back(i);
0279   }
0280   // Manually calculate helicity average.
0281   double helAvgNorm = i9.size()>0 ? 1:proc_ptr->getHelicityAveragingFactor();
0282   while (i9.size() > 0) {
0283     int i  = i9.back();
0284     int id = (i < nIn) ? idIn[i] : idOut[i-nIn];
0285     // How many spin states.
0286     int nS = particleDataPtr->spinType(id);
0287     // Massless particles have max 2 (physical) spin states.
0288     if (particleDataPtr->m0(id) == 0.0) nS=min(nS,2);
0289     if (i < nIn) helAvgNorm *= (double)nS;
0290     // Create nS copies of helConf, one for each spin state.
0291     int helConfSizeNow = helConf.size();
0292     for (int iCopy = 1; iCopy <= nS; ++iCopy) {
0293       // Set hel for this particle in this copy.
0294       // Start from -1, then 1, then 0 (if 3 states).
0295       int h = -1;
0296       if (nS == 1) h = 0;
0297       else if (iCopy == 2) h = 1;
0298       else if (iCopy == 3) h = 0;
0299       else if (iCopy == 4) h = -2;
0300       else if (iCopy == 5) h = 2;
0301       for (int iHelConf=0; iHelConf<helConfSizeNow; ++iHelConf) {
0302         vector<int> helNow = helConf[iHelConf];
0303         helNow[i] = h;
0304         // First copy: use existing.
0305         if (iCopy == 1) helConf[iHelConf] = helNow;
0306         // Subsequent copies: create new.
0307         else helConf.push_back(helNow);
0308       }
0309     }
0310     // Remove the particle whose helicities have been summed over.
0311     i9.pop_back();
0312   }
0313
0314   // Set colour configurations according to requested colour mode.
0315   vector<vector<int>> colConf;
0316   if (colMode >= 3) {
0317     // For FC sum use empty colour vector to communicate it with MG5.
0318     cols.clear();
0319     colConf.push_back(cols);
0320   } else if (colMode == 2) {
0321     // For LC sum, fetch all LC colour configurations.
0322     vector<vector<int>> colorConfigs = proc_ptr->getColorConfigs(procID);
0323     for (const auto& cc : colorConfigs) {
0324       int colID = proc_ptr->getColorIDForConfig(cc);
0325       if (proc_ptr->getColorFlowRelativeNCPower(colID) == 0)
0326         colConf.push_back(cc);
0327     }
0328   } else colConf.push_back(cols);
0329
0330   // Compute sum over colours and helicities.
0331   // (Save helicity components if needed later).
0332   double me2 = 0.0;
0333   me2hel.clear();
0334   proc_ptr->setMomenta(pn);
0335   for (const auto& cc : colConf) {
0336     proc_ptr->setColors(cc);
0337     for (int iHC(0); iHC<(int)helConf.size(); ++iHC) {
0338       proc_ptr->setHelicities(helConf[iHC]);
0339       double me2now = proc_ptr->sigmaKin();
0340       // MG5 may produce inf/nan for unphysical hel combinations.
0341       if (!isnan(me2now) && !isinf(me2now)) {
0342         // Save helicity matrix element for possible later use.
0343         me2hel[helConf[iHC]] = me2now;
0344         me2 += me2now;
0345       }
0346     }
0347   }
0348
0349   // Potentially include symmetry and averaging factors.
0350   if (inclSymFac) me2 /= proc_ptr->getSymmetryFactor();
0351   if (inclHelAvgFac) me2 /= helAvgNorm;
0352   if (inclColAvgFac) me2 /= proc_ptr->getColorAveragingFactor();
0353   return me2;
0354
0355 }
0356
0357 //--------------------------------------------------------------------------
0358
0359 // Fill lists.
0360
0361 void ExternalMEsMadgraph::fillLists(const vector<Particle>& state,
0362   vector<int>& idsIn, vector<int>& idsOut, vector<int>& hels,
0363   vector<int>& cols, vector<vector<double>>& pn) const {
0364   for (const Particle& ptcl : state) {
0365     if (ptcl.isFinal()) idsOut.push_back(ptcl.id());
0366     else idsIn.push_back(ptcl.id());
0367     vector<double> pNow = {ptcl.e(), ptcl.px(), ptcl.py(), ptcl.pz()};
0368     pn.push_back(pNow);
0369     hels.push_back(ptcl.pol());
0370     cols.push_back(ptcl.col());
0371     cols.push_back(ptcl.acol());
0372   }
0373 }
0374
0375 //--------------------------------------------------------------------------
0376
0377 // Declare the plugin.
0378
0379 PYTHIA8_PLUGIN_CLASS(ExternalMEs, ExternalMEsMadgraph, false, false, false)
0380 PYTHIA8_PLUGIN_VERSIONS(PYTHIA_VERSION_INTEGER)
0381
0382 //==========================================================================
0383
0384 } // end namespace Pythia8
0385
0386 #endif // end Pythia8_ExternalMEsMadgraph_H