11 #ifndef Pythia8_ExternalMEsMadgraph_H 12 #define Pythia8_ExternalMEsMadgraph_H 15 #include "Pythia8/ExternalMEs.h" 16 #include "Pythia8/Plugins.h" 34 isInitPtr =
false; isInit =
false; libPtr =
nullptr; modelPtr =
nullptr;}
38 if (modelPtr !=
nullptr)
delete modelPtr;}
44 bool isAvailable(vector<int> idIn, vector<int> idOut)
override;
47 bool isAvailable(
const vector<Particle>& state)
override;
50 double calcME2(
const vector<Particle>& state)
override;
51 double calcME2(
const Event& event,
const int iBeg = 3)
override;
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;
61 double calcME2(vector<int>& idIn, vector<int>& idOut,
62 vector< vector<double> >& pn, vector<int>& hels, vector<int>& cols,
63 set<int> sChannels = {});
65 PY8MEs_namespace::PY8MEs* libPtr;
78 int verbose = settingsPtr->mode(
"Vincia:verbose");
80 cout <<
" (ExternalMEsMadgraph::init()) begin -------" << endl;
82 loggerPtr->ERROR_MSG(
"cannot initialize, pointers not set");
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();
95 modelPtr->printIndependentParameters();
96 modelPtr->printIndependentCouplings();
111 double alpS = 1.0 / ( 4 * M_PI );
112 modelPtr->setDependentParameters(particleDataPtr, coupSMPtr, slhaPtr,
114 modelPtr->setDependentCouplings();
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);
122 libPtr->seProcessesIncludeSymmetryFactors(
false);
123 libPtr->seProcessesIncludeHelicityAveragingFactors(
false);
124 libPtr->seProcessesIncludeColorAveragingFactors(
false);
127 inclHelAvgFac =
true;
128 inclColAvgFac =
true;
143 set<int> req_s_channels;
144 PY8MEs_namespace::PY8ME * query
145 = libPtr->getProcess(idIn, idOut, req_s_channels);
146 return (query !=
nullptr);
151 vector <int> in, out;
153 set<int> req_s_channels;
155 PY8MEs_namespace::PY8ME* query
156 = libPtr->getProcess(in, out, req_s_channels);
157 return (query !=
nullptr);
162 vector <int> in, out;
164 set<int> req_s_channels;
166 PY8MEs_namespace::PY8ME* query
167 = libPtr->getProcess(in, out, req_s_channels);
168 return (query !=
nullptr);
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());
178 set<int> req_s_channels;
180 PY8MEs_namespace::PY8ME* query
181 = libPtr->getProcess(idIn, idOut, req_s_channels);
182 return (query !=
nullptr);
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.;
198 return calcME2(idIn, idOut, pn, hels, cols);
210 vector< vector<double> > pvec =
fillMoms(event, iBeg);
211 vector<int> helicities;
213 return calcME2(in, out, pvec, helicities, cols);
218 vector< vector<double> >& pn, vector<int>& hels, vector<int>& cols,
219 set<int> sChannels) {
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);
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();
235 vector< vector<int> > helConf;
236 helConf.push_back(hels);
241 for (
int i(0); i<(int)hels.size(); ++i)
242 if (hels[i] == 9) i9.push_back(i);
245 double helAvgNorm = i9.size()>0 ? 1:proc_ptr->getHelicityAveragingFactor();
246 while (i9.size() > 0) {
248 int id = (i < nIn) ? idIn[i] : idOut[i-nIn];
250 int nS = particleDataPtr->spinType(
id);
252 if (particleDataPtr->m0(
id) == 0.0) nS=min(nS,2);
253 if (i < nIn) helAvgNorm *= (double)nS;
255 int helConfSizeNow = helConf.size();
256 for (
int iCopy = 1; iCopy <= nS; ++iCopy) {
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];
269 if (iCopy == 1) helConf[iHelConf] = helNow;
271 else helConf.push_back(helNow);
279 vector<vector<int>> colConf;
283 colConf.push_back(cols);
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);
292 }
else colConf.push_back(cols);
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();
305 if (!isnan(me2now) && !isinf(me2now)) {
307 me2hel[helConf[iHC]] = me2now;
314 if (
inclSymFac) me2 /= proc_ptr->getSymmetryFactor();
315 if (inclHelAvgFac) me2 /= helAvgNorm;
316 if (inclColAvgFac) me2 /= proc_ptr->getColorAveragingFactor();
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()};
333 hels.push_back(ptcl.pol());
334 cols.push_back(ptcl.col());
335 cols.push_back(ptcl.acol());
345 PYTHIA8_PLUGIN_VERSIONS(PYTHIA_VERSION_INTEGER)
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
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
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
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