PYTHIA  8.312
VinciaCommon.h
1 // VinciaCommon.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2024 Peter Skands, 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 // This file contains global constants for the Vincia antenna shower as
7 // well as various auxiliary classes used by the Vincia shower model.
8 
9 #ifndef Pythia8_VinciaCommon_H
10 #define Pythia8_VinciaCommon_H
11 
12 // Maths headers.
13 #include <limits>
14 
15 // Include Pythia 8 headers.
16 #include "Pythia8/Event.h"
17 #include "Pythia8/Info.h"
18 #include "Pythia8/ParticleData.h"
19 #include "Pythia8/PartonSystems.h"
20 #include "Pythia8/PythiaStdlib.h"
21 #include "Pythia8/StandardModel.h"
22 
23 // Global Vincia constants, defined to live in a Vincia-specific
24 // namespace to avoid potential clashes with anything else defined in Pythia8.
25 namespace VinciaConstants {
26 
27 // Global numerical precision targets (should be kept within a
28 // reasonable margin of machine precision). Large values may result
29 // in strange results due to non-zero terms being artificially forced
30 // to zero.
31 const double UNITY = 1.0;
32 const double DECI = 1.0e-1;
33 const double CENTI = 1.0e-2;
34 const double MILLI = 1.0e-3;
35 const double MICRO = 1.0e-6;
36 const double NANO = 1.0e-9;
37 const double PICO = 1.0e-12;
38 
39 // Colour factors in Vincia normalization.
40 const double CA = 3.0;
41 const double CF = 8.0/3.0;
42 const double TR = 1.0;
43 const double NC = 3.0;
44 
45 // Mathematical constants (Euler–Mascheroni constant).
46 const double gammaE = 0.577215664901532860606512090082402431042;
47 
48 // Verbosity levels. Vincia has one more level (debug) beyond report.
49 const int DEBUG = 4;
50 
51 // Padding length for dashes in standardised Vincia verbose output.
52 const int DASHLEN = 80;
53 
54 }
55 
56 // Everything else lives in the Pythia 8 namespace.
57 
58 namespace Pythia8 {
59 
60 // Forward declaration.
61 class VinciaCommon;
62 
63 //==========================================================================
64 
65 // Enumerator for antenna function types, with "void" member NoFun.
66 enum AntFunType { NoFun,
67  QQEmitFF, QGEmitFF, GQEmitFF, GGEmitFF, GXSplitFF,
68  QQEmitRF, QGEmitRF, XGSplitRF,
69  QQEmitII, GQEmitII, GGEmitII, QXConvII, GXConvII,
70  QQEmitIF, QGEmitIF, GQEmitIF, GGEmitIF, QXConvIF,
71  GXConvIF, XGSplitIF };
72 
73 //==========================================================================
74 
75 // Convenient typedef for unsigned integers.
76 
77 typedef unsigned int uint;
78 
79 //==========================================================================
80 
81 // Print a method name using the appropritae pre-processor macro.
82 
83 // The following method was modified from
84 // boost/current_function.hpp - BOOST_CURRENT_FUNCTION
85 //
86 // Copyright (c) 2002 Peter Dimov and Multi Media Ltd.
87 //
88 // Distributed under the Boost Software License, Version 1.0.
89 // Boost Software License - Version 1.0 - August 17th, 2003
90 //
91 // Permission is hereby granted, free of charge, to any person or
92 // organization obtaining a copy of the software and accompanying
93 // documentation covered by this license (the "Software") to use,
94 // reproduce, display, distribute, execute, and transmit the
95 // Software, and to prepare derivative works of the Software, and to
96 // permit third-parties to whom the Software is furnished to do so,
97 // all subject to the following:
98 //
99 // The copyright notices in the Software and this entire statement,
100 // including the above license grant, this restriction and the
101 // following disclaimer, must be included in all copies of the
102 // Software, in whole or in part, and all derivative works of the
103 // Software, unless such copies or derivative works are solely in the
104 // form of machine-executable object code generated by a source
105 // language processor.
106 //
107 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
108 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
109 // MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND
110 // NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR
111 // ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR
112 // OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, ARISING
113 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
114 // OTHER DEALINGS IN THE SOFTWARE.
115 //
116 // http://www.boost.org/libs/utility/current_function.html
117 //
118 // Note that Boost Software License - Version 1.0 is fully compatible
119 // with GPLV2
120 // For more information see https://www.gnu.org/licenses/license-list.en.html
121 
122 #ifndef __METHOD_NAME__
123 
124 #ifndef VINCIA_FUNCTION
125 #if ( defined(__GNUC__) || (defined(__MWERKS__) && (__MWERKS__ >= 0x3000)) \
126 || (defined(__ICC) && (__ICC >= 600)) )
127 # define VINCIA_FUNCTION __PRETTY_FUNCTION__
128 #elif defined(__DMC__) && (__DMC__ >= 0x810)
129 # define VINCIA_FUNCTION __PRETTY_FUNCTION__
130 #elif defined(__FUNCSIG__)
131 # define VINCIA_FUNCTION __FUNCSIG__
132 #elif ( (defined(__INTEL_COMPILER) && (__INTEL_COMPILER >= 600)) \
133 || (defined(__IBMCPP__) && (__IBMCPP__ >= 500)) )
134 # define VINCIA_FUNCTION __FUNCTION__
135 #elif defined(__BORLANDC__) && (__BORLANDC__ >= 0x550)
136 # define VINCIA_FUNCTION __FUNC__
137 #elif defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901)
138 # define VINCIA_FUNCTION __func__
139 #else
140 # define VINCIA_FUNCTION "unknown"
141 #endif
142 #endif // end VINCIA_FUNCTION
143 
144 inline std::string methodName(const std::string& prettyFunction, bool
145  withPythia=false) {
146  size_t colons = prettyFunction.find("::");
147  // Include Pythia8:: or not.
148  size_t begin = colons + 2;
149  if (withPythia) begin = prettyFunction.substr(0,colons).rfind(" ") + 1;
150  size_t end = prettyFunction.rfind("(") - begin;
151  return prettyFunction.substr(begin,end) + "()";
152 }
153 
154 #define __METHOD_NAME__ methodName(VINCIA_FUNCTION)
155 #endif // end __METHOD_NAME__
156 
157 //==========================================================================
158 
159 // Global functions accessible in the Vincia namespace.
160 
161 // Utilities for printing out VINCIA Messages.
162 void printOut(string,string,int nPad=0,char padChar='-');
163 
164 // String utilities.
165 string num2str(int,int width=4) ;
166 string num2str(double,int width=9) ;
167 string bool2str(bool, int width=3) ;
168 
169 // Search and replace a string.
170 inline string replaceString(string subject, const string& search,
171  const string& replace) {
172  string::size_type pos = 0;
173  while ((pos = subject.find(search, pos)) != string::npos) {
174  subject.replace(pos, search.length(), replace);
175  pos += replace.length();
176  }
177  return subject;
178 
179 }
180 
181 // Remove ":" and "/" from a file name.
182 inline string sanitizeFileName(string fileName) {
183  map<string, string> rep;
184  rep["/"] = "_div_";
185  rep[":"] = "_colon_";
186  string retVal = fileName;
187  for (map<string, string>::const_iterator it = rep.begin(); it != rep.end();
188  ++it) {
189  retVal = replaceString(retVal, it->first, it->second);
190  }
191  return retVal;
192 
193 }
194 
195 // Utility for checking if a file exists.
196 inline bool fileExists(const std::string& name) {
197  if (FILE *file = fopen(name.c_str(), "r")) {
198  fclose(file);
199  return true;
200  } else {
201  return false;
202  }
203 
204 }
205 
206 //==========================================================================
207 
208 // A class to store and process colour information, e.g. colour maps,
209 // reconnections, etc.
210 
212 
213 public:
214 
215  // Initialize pointers (must be done before init).
216  void initPtr(Info* infoPtrIn) {
217  infoPtr = infoPtrIn;
218  particleDataPtr = infoPtr->particleDataPtr;
219  settingsPtr = infoPtr->settingsPtr;
220  partonSystemsPtr = infoPtr->partonSystemsPtr;
221  rndmPtr = infoPtr->rndmPtr;
222  isInitPtr=true;
223  }
224 
225  // Initialize.
226  bool init();
227 
228  // Get translation map from ordinary Les Houches tags to antenna
229  // indices (tags always ending on [1-9] with gluons having
230  // non-identical colour and anticolour indices).
231  bool colourise(int iSys, Event& event);
232 
233  // Method to sort list of partons in Vincia colour order. Returns
234  // vector<int> with indices in following order:
235  // (1) colourless incoming particles
236  // (2) triplet-ngluon-antitriplet contractions
237  // (3) ngluon rings
238  // (4) colourless outgoing particles
239  // No specific order is imposed on the relative ordering inside each
240  // of the four classes.
241  vector<int> colourSort(vector<Particle*>);
242 
243  // Method to create LC colour map and list of LC antennae.
244  void makeColourMaps(const int iSysIn, const Event& event,
245  map<int,int>& indexOfAcol, map<int,int>& indexOfCol,
246  vector< pair<int,int> >& antLC, const bool findFF, const bool findIX);
247 
248  // Determine whether 01 or 12 inherit the colour tag from a parent.
249  bool inherit01(double s01, double s12);
250 
251  // Set verbose level.
252  void setVerbose(int verboseIn) {verbose = verboseIn;};
253 
254 private:
255 
256  // Internal parameters.
257  int inheritMode{};
258 
259  // Is initialized.
260  bool isInitPtr{false}, isInit{false};
261 
262  // Pointers to PYTHIA 8 objects.
263  Info* infoPtr;
264  ParticleData* particleDataPtr;
265  Settings* settingsPtr;
266  PartonSystems* partonSystemsPtr;
267  Rndm* rndmPtr;
268 
269  // Verbose level.
270  int verbose{};
271 
272 };
273 
274 //==========================================================================
275 
276 // Simple struct to store information about a 3 -> 2 clustering.
277 
279 
280  // Set information about daughters from current event and daughters indices.
281  void setDaughters(const Event& state, int dau1In, int dau2In, int dau3In);
282  void setDaughters(const vector<Particle>& state, int dau1In, int dau2In,
283  int dau3In);
284 
285  // Set mother particle ids.
286  void setMothers(int idMot1In, int idMot2In) {
287  idMot1 = idMot1In;
288  idMot2 = idMot2In;
289  }
290 
291  // Set antenna information.
292  void setAntenna(bool isFSRin, enum AntFunType antFunTypeIn) {
293  isFSR = isFSRin;
294  antFunType = antFunTypeIn;
295  }
296 
297  // Initialise vectors of invariants and masses.
298  bool initInvariantAndMassVecs();
299 
300  // Set invariants and masses.
301  void setInvariantsAndMasses(const Event& state);
302  void setInvariantsAndMasses(const vector<Particle>& state);
303 
304  // Swap 1 <-> 3, including information about parents.
305  void swap13() {
306  swap(dau1,dau3);
307  swap(idMot1,idMot2);
308  swap(saj,sjb);
309  if (mDau.size() == 3)
310  swap(mDau[0],mDau[2]);
311  if (mMot.size() == 2)
312  swap(mMot[0],mMot[1]);
313  if (invariants.size() == 3) {
314  swap(invariants[1],invariants[2]);
315  }
316  }
317 
318  // Methods to get antenna type.
319  bool isFF() const {
320  if (!isFSR) return false;
321  else if (antFunType >= QQEmitFF && antFunType < QQEmitRF) return true;
322  else return false;
323  }
324  bool isRF() const {
325  if (!isFSR) return false;
326  else if (antFunType >= QQEmitRF && antFunType < QQEmitII) return true;
327  else return false;
328  }
329  bool isII() const {
330  if (isFSR) return false;
331  else if (antFunType >= QQEmitII && antFunType < QQEmitIF) return true;
332  return false;
333  }
334  bool isIF() const {
335  if (isFSR) return false;
336  else if (antFunType >= QQEmitIF) return true;
337  else return false;
338  }
339  string getAntName() const;
340 
341  // Methods to get branching type (currently only 2 -> 3).
342  bool is2to3() const { return true; }
343 
344  // Branching daughter information (indices in event record).
345  int dau1{}, dau2{}, dau3{};
346 
347  // Antenna information.
348  bool isFSR{true};
349  AntFunType antFunType{NoFun};
350 
351  // Mother ids.
352  int idMot1{}, idMot2{};
353 
354  // Helicities.
355  vector<int> helDau = {9, 9, 9};
356  vector<int> helMot = {9, 9};
357 
358  // Masses.
359  vector<double> mDau;
360  vector<double> mMot;
361 
362  // Invariants.
363  double saj{}, sjb{}, sab{};
364  // Vector of invariants (stored as sAB, saj, sjb, sab).
365  vector<double> invariants;
366 
367  // Value of sector resolution variable that this clustering produces.
368  double q2res{};
369 
370  // Value of evolution variable that this clustering produces.
371  double q2evol{};
372 
373  // Kinematic map (only used for FF).
374  int kMapType{};
375 
376 };
377 
378 //==========================================================================
379 
380 // A simple class for containing evolution variable definitions.
381 
382 class Resolution {
383 
384 public:
385 
386  // Initialize pointers (must be done before init).
387  void initPtr(Settings* settingsPtrIn, Info* infoPtrIn,
388  VinciaCommon* vinComPtrIn) {
389  settingsPtr = settingsPtrIn;
390  infoPtr = infoPtrIn;
391  loggerPtr = infoPtrIn->loggerPtr;
392  vinComPtr = vinComPtrIn;
393  isInitPtr = true;
394  }
395 
396  // Initialize.
397  bool init();
398 
399  // Method to calculate (and set) evolution variable.
400  double q2evol(VinciaClustering& clus);
401 
402  // Method to calculate dimensionless evolution variable.
403  // Note: calls q2evol(), so will set dimensionful evolution
404  // variable in the VinciaClustering object.
405  double xTevol(VinciaClustering& clus);
406 
407  // Top-level function to calculate (and set) sector resolution.
408  double q2sector(VinciaClustering& clus);
409 
410  // Find sector with minimal resolution.
411  // Optionally resolve Born: avoid clusterings that would not lead
412  // to a specified (Born) configuration.
413  VinciaClustering findSector(vector<Particle>& state,
414  map<int,int> flavsBorn);
415  // Find sector with minimal resolution.
416  // Optionally avoid sectors that cluster beyond a minimal number
417  // of quark pairs or gluons.
418  VinciaClustering findSector(vector<Particle>& state,
419  int nqpMin = 0, int ngMin = 0);
420 
421  // Find sector with minimal q2sector in list of clusterings.
422  VinciaClustering getMinSector(vector<VinciaClustering>& clusterings);
423 
424  // Sector veto to check whether given value of resolution is minimal,
425  // given we want to preserve a certain Born configuration.
426  // Returns true = vetoed, and false = not vetoed.
427  bool sectorVeto(double q2In, vector<Particle>& state,
428  map<int,int> nFlavsBorn) {
429  VinciaClustering clusMin = findSector(state, nFlavsBorn);
430  if (q2In > clusMin.q2res) return true;
431  else return false;
432  }
433  bool sectorVeto(const VinciaClustering& clusMin,
434  const VinciaClustering& clus);
435 
436  // Set verbosity level.
437  void setVerbose(int verboseIn) {verbose = verboseIn;}
438 
439 private:
440 
441  // Member functions.
442 
443  // Sector resolution for 2 -> 3 branchings.
444  double q2sector2to3FF(VinciaClustering& clus);
445  double q2sector2to3RF(VinciaClustering& clus);
446  double q2sector2to3IF(VinciaClustering& clus);
447  double q2sector2to3II(VinciaClustering& clus);
448 
449  // Sector resolution function for 3->4 branchings (currently only
450  // used for gluon splitting, with m2qq as the measure).
451  //TODO: currently disabled.
452  // double q2sector3to4(const Particle*, const Particle*,
453  // const Particle* j1, const Particle* j2) {return -1.;}
454 
455  // Sector resolution function for 2->4 branchings (double emission).
456  // Assume j1 and j2 are colour connected, with a and b hard
457  // recoilers.
458  //TODO: currently disabled.
459  // double q2sector2to4(const Particle* a, const Particle* b,
460  // const Particle* j1, const Particle* j2) {return -1.;}
461 
462  // Sector resolution function for 3->5 branchings (emission +
463  // splitting).
464  //TODO: currently disabled.
465  // double q2sector3to5(Particle* a, Particle* b,
466  // Particle* j1, Particle* j2, Particle* j3) {return -1;}
467 
468  // Members.
469 
470  // Initialized.
471  bool isInitPtr{false}, isInit{false};
472 
473  // Pointer to PYTHIA 8 settings database.
474  Settings* settingsPtr{};
475  Info* infoPtr{};
476  Logger* loggerPtr{};
477 
478  // Pointer to VinciaCommon.
479  VinciaCommon* vinComPtr{};
480 
481  // Number of flavours to be treated as massless.
482  int nFlavZeroMassSav{};
483 
484  // Verbosity level.
485  int verbose{};
486 
487 };
488 
489 //==========================================================================
490 
491 // Class which contains functions and variables shared by the
492 // VinciaShower and VinciaMatching classes.
493 
495 
496 public:
497 
498  // Initialize pointers.
499  bool initPtr(Info* infoPtrIn) {
500  infoPtr = infoPtrIn;
501  particleDataPtr = infoPtr->particleDataPtr;
502  settingsPtr = infoPtr->settingsPtr;
503  loggerPtr = infoPtr->loggerPtr;
504  rndmPtr = infoPtr->rndmPtr;
505  partonSystemsPtr = infoPtr->partonSystemsPtr;
506  isInitPtr = true;
507  return true;
508  }
509 
510  // Initialize data members.
511  bool init();
512 
513  // Function to check for the hadronization cutoff for a colour
514  // connected parton pair.
515  double mHadMin(const int id1, const int id2);
516 
517  // Function to check the event after each branching. Added by NF to
518  // see if certain Pythia warnings/error are caused by the shower.
519  bool showerChecks(Event& event, bool ISR);
520 
521  // Function to reset counters (print once every event for now).
522  void resetCounters() {
523  nUnkownPDG = 0;
524  nIncorrectCol = 0;
525  nNAN = 0;
526  nVertex = 0;
527  nChargeCons = 0;
528  nMotDau = 0;
529  for (int i=0; i<2; i++) {
530  nUnmatchedMass[i] = 0;
531  nEPcons[i] = 0;
532  }
533  }
534 
535  // Get number of active flavors at a given Q scale.
536  int getNf(double q) {
537  if (q <= mc) return 3;
538  else if (q <= mb) return 4;
539  else if (q <= mt) return 5;
540  else return 6;
541  }
542 
543  // Get the shower starting scale.
544  double getShowerStartingScale(int iSys, const Event& event,
545  double sbbSav);
546 
547  // Method to find all possible clusterings for a given system,
548  // given we want to resolve a certain Born configuration, i.e.,
549  // not cluster more gluons or quark flavours as we had in the Born.
550  vector<VinciaClustering> findClusterings(const vector<Particle>& state,
551  map<int, int> nFlavsBorn);
552  // Method to find all possible clusterings while retaining a certain
553  // minimal number of quark pairs and gluons.
554  vector<VinciaClustering> findClusterings(const vector<Particle>& state,
555  int nqpMin = 0, int ngMin = 0);
556 
557  // Check if clustering is sensible, i.e., corresponds to an existing antenna.
558  bool isValidClustering(const VinciaClustering& clus,
559  const Event& event, int verboseIn);
560 
561  // Perform a clustering.
562  bool clus3to2(const VinciaClustering& clus, const Event& event,
563  vector<Particle>& pClustered);
564  bool clus3to2(const VinciaClustering& clus, const vector<Particle>& state,
565  vector<Particle>& pClustered);
566  // Helper functions to perform clustering.
567  bool getCols3to2(const Particle* a, const Particle* j, const Particle* b,
568  const VinciaClustering& clus, pair<int,int>& colsA, pair<int,int>& colsB);
569  bool getMomenta3to2(vector<Vec4>& momNow, vector<Vec4>& momClus,
570  const VinciaClustering& clus, int iOffset = 0);
571 
572  // 3->2 clustering maps.
573  bool map3to2FF(vector<Vec4>& pClu, const vector<Vec4> pIn,
574  int kMapType, int a=0, int r=1, int b=2, double mI=0.0, double mK=0.0) {
575  if (mI == 0. && mK == 0.)
576  return map3to2FFmassless(pClu, pIn, kMapType, a, r, b);
577  else
578  return map3to2FFmassive(pClu, pIn, kMapType, mI, mK, a, r, b);
579  }
580  bool map3to2RF(vector<Vec4>& pClu, const vector<Vec4>& pIn, int a=0,
581  int r=1, int b=2, double mK = 0.);
582  bool map3to2IF(vector<Vec4>& pClu, const vector<Vec4>& pIn,
583  int a = 0, int r = 1, int b = 2,
584  double mj = 0., double mk = 0., double mK = 0.);
585  bool map3to2II(vector<Vec4>& pClu, const vector<Vec4>& pIn, bool doBoost,
586  int a = 0, int r = 2, int b = 1, double mj = 0.);
587 
588  // 2->3 kinematics maps for FF branchings. Original implementations;
589  // massless by Skands, massive by Ritzmann.
590  bool map2to3FF(vector<Vec4>& pNew, const vector<Vec4>& pOld, int kMapType,
591  const vector<double>& invariants, double phi, vector<double> masses) {
592  if ( masses.size() <= 2 || ( masses[0] == 0.0 && masses[1] == 0.0
593  && masses[2] == 0.0 )) {
594  return map2to3FFmassless(pNew, pOld, kMapType, invariants, phi);
595  } else {
596  return map2to3FFmassive(pNew, pOld, kMapType, invariants, phi, masses);
597  }
598  }
599 
600  // 2->3 kinematics maps for II branchings. Original implementations:
601  // massless by Fischer, massive by Verheyen.
602  bool map2to3II(vector<Vec4>& pNew, vector<Vec4>& pRec,
603  vector<Vec4>& pOld, double sAB, double saj, double sjb, double sab,
604  double phi, double m2j = 0.0) {
605  if (m2j == 0.0)
606  return map2to3IImassless(pNew, pRec, pOld, sAB, saj, sjb, sab, phi);
607  else
608  return map2to3IImassive(pNew, pRec, pOld, sAB, saj, sjb, sab, phi, m2j);
609  }
610 
611  // 2->3 kinematics maps for IF branchings. General massive case
612  // implemented by Verheyen.
613  bool map2to3IFlocal(vector<Vec4>& pNew, const vector<Vec4>& pOld,
614  double sOldAK, double saj, double sjk, double sak, double phi,
615  double m2oldK, double m2j, double m2k);
616  bool map2to3IFglobal(vector<Vec4>& pNew, vector<Vec4>& pRec,
617  const vector<Vec4>& pOld, const Vec4 &pB,
618  double sAK, double saj, double sjk, double sak, double phi,
619  double mK2, double mj2, double mk2);
620 
621  // Resonance decay kinematic maps.
622  bool map2toNRF(vector<Vec4>& pAfter, const vector<Vec4> pBefore,
623  unsigned int posR, unsigned int posF,
624  const vector<double> invariants, double phi,
625  const vector<double> masses);
626 
627  // 1->2 decay map for (already offshell) resonance decay
628  bool map1to2RF(vector<Vec4>& pNew, const Vec4 pRes, double m1,
629  double m2, double theta, double phi);
630 
631  // Check if 2-particle system is on-shell and rescale if not.
632  bool onShellCM(Vec4& p1, Vec4& p2, double m1, double m2, double tol = 1e-6);
633 
634  // Force initial-state and light-flavour partons to be massless.
635  bool mapToMassless(int iSys, Event& event, bool makeNewCopies);
636 
637  // Map a massless antenna to equivalent massive one. Boolean
638  // returns true if a modification was made.
639  bool mapToMassive(Vec4& p1, Vec4& p2, double m1, double m2) {
640  return (!onShellCM(p1,p2,m1,m2,1e-9));
641  }
642 
643  // Make list of particles as vector<Particle>.
644  // First 1 or 2 entries : incoming particle(s).
645  // Subseqent entries : outgoing particles.
646  // The two last arguments are optional and allow to specify a list
647  // of indices to be ignored, and a set of particles to be added, e.g.
648  // in the context of setting up a trial state after a branching.
649  // The newly added particles are then at the end of the respective
650  // lists, i.e. a newly added incoming particle is the last incoming
651  // one and newly added outgoing ones are the last among the outgoing
652  // ones.
653  vector<Particle> makeParticleList(const int iSys, const Event& event,
654  const vector<Particle> &pNew = vector<Particle>(),
655  const vector<int> &iOld = vector<int>());
656 
657  // Method to find all antennae that can produce a branching.
658  // IN: indices of clustering in event, where i2 is the emission.
659  // OUT: vector of VinciaClusterings.
660  // Also swap daughters to match antenna function convention if needed
661  // (e.g. for GXSplitFF, when dau2 and dau3 form the quark pair).
662  vector<VinciaClustering> findAntennae(Event& state, int i1, int i2, int i3);
663 
664  // Check whether two particles are colour connected.
665  bool colourConnected(const Particle& ptcl1, const Particle& ptcl2);
666 
667  // Print a list of Particles.
668  void list(const vector<Particle>& state, string title = "",
669  bool footer = true);
670 
671  // Print a list of VinciaClusterings.
672  void list(const vector<VinciaClustering>& clusterings, string title = "",
673  bool footer = true);
674 
675  // Get/set verbose parameter.
676  int getVerbose() {return verbose; };
677  void setVerbose(int verboseIn) { verbose = verboseIn;};
678 
679  // Public data members: strong coupling in MSbar and CMW schemes,
680  // user and default choices,
681  AlphaStrong alphaStrong{}, alphaStrongCMW{}, alphaStrongDef{},
682  alphaStrongDefCMW{};
683 
684  // Couplings for use in merging.
685  AlphaStrong alphaS{};
686  AlphaEM alphaEM{};
687  double mu2freeze{}, mu2min{}, alphaSmax{};
688 
689  // Quark masses.
690  double ms{}, mc{}, mb{}, mt{};
691  int nFlavZeroMass{};
692 
693  // Checks.
694  double epTolErr{}, epTolWarn{}, mTolErr{}, mTolWarn{};
695 
696 private:
697 
698  // Functions.
699 
700  // Special cases of 3 -> 2 maps.
701  bool map3to2FFmassive(vector<Vec4>& pClu, const vector<Vec4> pIn,
702  int kMapType, double mI, double mK, int a=0, int r=1, int b=2);
703  bool map3to2FFmassless(vector<Vec4>& pClu, const vector<Vec4> pIn,
704  int kMapType, int a=0, int r=1, int b=2);
705 
706  // Special cases of 2 -> 3 maps.
707  bool map2to3FFmassive(vector<Vec4>& pNew, const vector<Vec4>& pOld,
708  int kMapType, const vector<double>& invariants, double phi,
709  vector<double> masses);
710  bool map2to3FFmassless(vector<Vec4>& pNew, const vector<Vec4>& pOld,
711  int kMapType, const vector<double>& invariants, double phi);
712  bool map2to3IImassive(vector<Vec4>& pNew, vector<Vec4>& pRec,
713  vector<Vec4>& pOld, double sAB, double saj, double sjb, double sab,
714  double phi, double m2j = 0.0);
715  bool map2to3IImassless(vector<Vec4>& pNew, vector<Vec4>& pRec,
716  vector<Vec4>& pOld, double sAB, double saj, double sjb, double sab,
717  double phi);
718  bool map2to3RF(vector<Vec4>& pThree, const vector<Vec4> pTwo,
719  const vector<double> invariants, double phi,
720  const vector<double> masses);
721 
722  // Members.
723 
724  // Pointers.
725  Info* infoPtr{};
726  Settings* settingsPtr{};
727  ParticleData* particleDataPtr{};
728  Rndm* rndmPtr{};
729  Logger* loggerPtr{};
730  PartonSystems* partonSystemsPtr{};
731 
732  // Counter for output control.
733  int nUnkownPDG{}, nIncorrectCol{}, nNAN{}, nVertex{}, nChargeCons{},
734  nMotDau{};
735  vector<int> nUnmatchedMass, nEPcons;
736 
737  // Internal flags and settings.
738  bool isInitPtr{false}, isInit{false};
739  int verbose{};
740 
741 };
742 
743 //==========================================================================
744 
745 } // end namespace Pythia8
746 
747 #endif // Pythia8_VinciaCommon_H
void swap13()
Swap 1 <-> 3, including information about parents.
Definition: VinciaCommon.h:305
Definition: StandardModel.h:23
void setVerbose(int verboseIn)
Set verbosity level.
Definition: VinciaCommon.h:437
const double UNITY
Definition: VinciaCommon.h:31
Definition: Info.h:45
The Event class holds all info on the generated event.
Definition: Event.h:453
bool is2to3() const
Methods to get branching type (currently only 2 -> 3).
Definition: VinciaCommon.h:342
Simple struct to store information about a 3 -> 2 clustering.
Definition: VinciaCommon.h:278
Logger * loggerPtr
Pointer to the logger.
Definition: Info.h:86
Definition: VinciaCommon.h:211
vector< double > invariants
Vector of invariants (stored as sAB, saj, sjb, sab).
Definition: VinciaCommon.h:365
void printOut(string, string, int nPad=0, char padChar='-')
end METHOD_NAME
Definition: VinciaCommon.cc:4964
int getNf(double q)
Get number of active flavors at a given Q scale.
Definition: VinciaCommon.h:536
string replaceString(string subject, const string &search, const string &replace)
Search and replace a string.
Definition: VinciaCommon.h:170
string num2str(int, int width=4)
String utilities.
Definition: VinciaCommon.cc:4924
void initPtr(Settings *settingsPtrIn, Info *infoPtrIn, VinciaCommon *vinComPtrIn)
Initialize pointers (must be done before init).
Definition: VinciaCommon.h:387
double q2res
Value of sector resolution variable that this clustering produces.
Definition: VinciaCommon.h:368
Definition: StandardModel.h:106
Definition: Logger.h:23
bool map2to3II(vector< Vec4 > &pNew, vector< Vec4 > &pRec, vector< Vec4 > &pOld, double sAB, double saj, double sjb, double sab, double phi, double m2j=0.0)
Definition: VinciaCommon.h:602
AntFunType
Enumerator for antenna function types, with "void" member NoFun.
Definition: VinciaCommon.h:66
string sanitizeFileName(string fileName)
Remove ":" and "/" from a file name.
Definition: VinciaCommon.h:182
Definition: Basics.h:386
Definition: VinciaCommon.h:494
void resetCounters()
Function to reset counters (print once every event for now).
Definition: VinciaCommon.h:522
void setAntenna(bool isFSRin, enum AntFunType antFunTypeIn)
Set antenna information.
Definition: VinciaCommon.h:292
std::string methodName(const std::string &prettyFunction, bool withNamespace=false)
end PYTHIA_FUNCTION
Definition: PythiaStdlib.h:283
const double gammaE
Mathematical constants (Euler–Mascheroni constant).
Definition: VinciaCommon.h:46
void setVerbose(int verboseIn)
Set verbose level.
Definition: VinciaCommon.h:252
Definition: Event.h:32
int getVerbose()
Get/set verbose parameter.
Definition: VinciaCommon.h:676
void initPtr(Info *infoPtrIn)
Initialize pointers (must be done before init).
Definition: VinciaCommon.h:216
A simple class for containing evolution variable definitions.
Definition: VinciaCommon.h:382
ParticleData * particleDataPtr
Pointer to the particle data table.
Definition: Info.h:83
bool map2to3FF(vector< Vec4 > &pNew, const vector< Vec4 > &pOld, int kMapType, const vector< double > &invariants, double phi, vector< double > masses)
Definition: VinciaCommon.h:590
Maths headers.
Definition: VinciaCommon.h:25
bool fileExists(const std::string &name)
Utility for checking if a file exists.
Definition: VinciaCommon.h:196
double m2(const Vec4 &v1)
The squared invariant mass of one or more four-vectors.
Definition: Basics.cc:605
bool initPtr(Info *infoPtrIn)
Initialize pointers.
Definition: VinciaCommon.h:499
vector< double > mDau
Masses.
Definition: VinciaCommon.h:359
const int DEBUG
Verbosity levels. Vincia has one more level (debug) beyond report.
Definition: VinciaCommon.h:49
double phi(const Vec4 &v1, const Vec4 &v2)
phi is azimuthal angle between v1 and v2 around z axis.
Definition: Basics.cc:693
The PartonSystems class describes the whole set of subcollisions.
Definition: PartonSystems.h:42
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
const int DASHLEN
Padding length for dashes in standardised Vincia verbose output.
Definition: VinciaCommon.h:52
bool sectorVeto(double q2In, vector< Particle > &state, map< int, int > nFlavsBorn)
Definition: VinciaCommon.h:427
unsigned int uint
Convenient typedef for unsigned integers.
Definition: VinciaCommon.h:77
double theta(const Vec4 &v1, const Vec4 &v2)
theta is polar angle between v1 and v2.
Definition: Basics.cc:662
bool mapToMassive(Vec4 &p1, Vec4 &p2, double m1, double m2)
Definition: VinciaCommon.h:639
const double CA
Colour factors in Vincia normalization.
Definition: VinciaCommon.h:40
bool map3to2FF(vector< Vec4 > &pClu, const vector< Vec4 > pIn, int kMapType, int a=0, int r=1, int b=2, double mI=0.0, double mK=0.0)
3->2 clustering maps.
Definition: VinciaCommon.h:573
bool isFF() const
Methods to get antenna type.
Definition: VinciaCommon.h:319
This class holds a map of all ParticleDataEntries.
Definition: ParticleData.h:422
Definition: Basics.h:32
Definition: Settings.h:195
void setMothers(int idMot1In, int idMot2In)
Set mother particle ids.
Definition: VinciaCommon.h:286