StRoot  1
VinciaAntennaFunctions.h
1 // VinciaAntennaFunctions.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 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 header information for the AntennaFunction base
7 // class, its derived classes for FF, IF, and II antenna functions,
8 // the AntennaSetFSR and AntennaSetISR classes, and the MEC class for
9 // matrix- element corrections.
10
11 #ifndef Pythia8_VinciaAntennaFunctions_H
12 #define Pythia8_VinciaAntennaFunctions_H
13
15 #include "Pythia8/Basics.h"
16 #include "Pythia8/Event.h"
17 #include "Pythia8/PythiaStdlib.h"
18 #include "Pythia8/ShowerMEs.h"
19
21 #include "Pythia8/VinciaCommon.h"
22
23 namespace Pythia8 {
24
25 //==========================================================================
26
27 // A class containing DGLAP splitting functions for limit checking.
28
29 class DGLAP {
30
31 public:
32
33  // Constructor.
34  DGLAP() {;}
35
36  // Helicity-dependent Altarelli-Parisi kernels (mu = m/Q). Note,
37  // Pg2gg is written with standard factor 2 normalization convention.
38  double Pg2gg(double z, int hA=9, int hB=9, int hC=9);
39  double Pg2qq(double z, int hA=9, int hB=9, int hC=9, double mu=0.);
40  double Pq2qg(double z, int hA=9, int hB=9, int hC=9, double mu=0.);
41  double Pq2gq(double z, int hA=9, int hB=9, int hC=9, double mu=0.);
42
43  // Wrappers to get unpolarized massive Altarelli-Pariis kernels (mu = m/Q).
44  double Pg2qq(double z, double mu) {return Pg2qq(z, 9, 9, 9, mu);}
45  double Pq2qg(double z, double mu) {return Pq2qg(z, 9, 9, 9, mu);}
46  double Pq2gq(double z, double mu) {return Pq2gq(z, 9, 9, 9, mu);}
47
48  // Altarelli-Parisi kernels with linear in/out polarizations for
49  // gluons: pol = +1 for in-plane, -1 for out-of-plane.
50  double Pg2ggLin(double z, int polA = 9, int polB = 9, int polC = 9);
51  double Pg2qqLin(double z, int polA = 9, int hB = 9, int hC = 9,
52  double mu = 0.);
53  double Pq2qgLin(double z, int hA = 9, int hB=9, int polC = 9,
54  double mu = 0.);
55  double Pq2gqLin(double z, int hA = 9, int polB=9, int hC = 9,
56  double mu = 0.);
57
58 };
59
60 //==========================================================================
61
62 // The AntennaFunction base class. Base class implementation for all
63 // AntennaFunction objects.
64
66
67 public:
68
69  // Constructor.
70  AntennaFunction() = default;
71
72  // Destructor
73  virtual ~AntennaFunction() {};
74
75  // Names of this antenna, for VINCIA, and for humans.
76  virtual string vinciaName() const = 0;
77
78  // Parton types (AB -> 0i 1j 2k): needed by soft- and collinear-limit checks.
79  virtual int idA() const = 0;
80  virtual int idB() const = 0;
81  virtual int id1() const = 0;
82
83  // The antenna function [GeV^-2].
84  virtual double antFun(vector<double> invariants, vector<double> mNew,
85  vector<int> helBef, vector<int> helNew) = 0;
86
87  // Optional implementation of the DGLAP kernels for collinear-limit checks
88  // Defined as PI/sij + PK/sjk, i.e. equivalent to antennae.
89  virtual double AltarelliParisi(vector<double> invariants,
90  vector<double> mNew, vector<int> helBef, vector<int> helNew) = 0;
91
92  // Default initialization.
93  virtual bool init();
94
95  // Construct baseName from idA, idB, and id1.
96  virtual string baseName() const {
97  return id2str(id1()) + "/" + id2str(idA()) + id2str(idB());}
98
99  // Wrapper that can modify baseName to more human readable form if required.
100  virtual string humanName() const {return baseName();}
101
102  // Function to check singularities, positivity, etc.
103  virtual bool check();
104
105  // Method to intialise mass values.
106  virtual void initMasses(vector<double>* masses) {
107  if (masses->size() >= 3) {
108  mi = masses->at(0); mj = masses->at(1); mk = masses->at(2);
109  } else {mi = 0.0; mj = 0.0; mk = 0.0;}}
110
111  // Method to initialise internal helicity variables.
112  virtual int initHel(vector<int>* helBef, vector<int>* helNew);
113
114  // Wrapper for helicity-summed/averaged antenna function.
115  double antFun(vector<double> invariants, vector<double> masses) {
116  return antFun(invariants, masses, hDum, hDum);}
117
118  // Wrapper for massless, helicity-summed/averaged antenna function.
119  double antFun(vector<double> invariants) {
120  return antFun(invariants, mDum, hDum, hDum);}
121
122  // Wrapper without helicity assignments.
123  double AltarelliParisi(vector<double> invariants, vector<double> masses) {
124  return AltarelliParisi(invariants, masses, hDum, hDum);}
125
126  // Wrapper for massless helicity-summed/averaged DGLAP kernels.
127  double AltarelliParisi(vector<double> invariants) {
128  return AltarelliParisi(invariants, mDum, hDum, hDum);}
129
130  // Initialize pointers.
131  void initPtr(Info* infoPtrIn, DGLAP* dglapPtrIn);
132
133  // Get parameters.
134  double chargeFac() {return chargeFacSav;}
135  int kineMap() {return kineMapSav;}
136  double alpha() {return alphaSav;}
137  double sectorDamp() {return sectorDampSav;}
138
139  // Functions to get Altarelli-Parisi energy fractions from invariants.
140  double zA(vector<double> invariants) {
141  double yij = invariants[1]/invariants[0];
142  double yjk = invariants[2]/invariants[0];
143  return (1.-yjk)/(1.+yij);}
144  double zB(vector<double> invariants) {
145  double yij = invariants[1]/invariants[0];
146  double yjk = invariants[2]/invariants[0];
147  return (1.-yij)/(1.+yjk);}
148
149  // Auxiliary function to translate an ID code to a string.
150  string id2str(int id) const;
151
152 protected:
153
154  // Is initialized.
155  bool isInitPtr{false}, isInit{false};
156
157  // Charge factor, kinematics map, and subleading-colour treatment.
158  double chargeFacSav{0.0};
159  int kineMapSav{0}, modeSLC{-1};
160  bool sectorShower{false};
161
162  // The alpha collinear-partitioning parameter.
163  double alphaSav{0.0};
164
165  // The sector-shower collinear dampening parameter.
166  double sectorDampSav{0.0};
167
168  // Shorthand for commonly used variable(s).
169  double term{}, preFacFiniteTermSav{0.0}, antMinSav{0.0};
170  bool isMinVar{};
171
172  // Variables for internal storage of masses and helicities.
173  double mi{0.0}, mj{0.0}, mk{0.0};
174  int hA{9}, hB{9}, hi{9}, hj{9}, hk{9};
175
176  // Map to tell whether a given helicity value maps to L- and/or
177  // R-handed. Defined by constructor and not to be changed
178  // dynamically.
179  map<int, bool> LH{{9, true}, {1, false}, {-1, true}};
180  map<int, bool> RH{{9, true}, {1, true}, {-1, false}};
181
182  // Verbosity level.
183  int verbose{1};
184
185  // Pointers to Pythia8 classes.
186  Info* infoPtr{};
187  ParticleData* particleDataPtr{};
188  Settings* settingsPtr{};
189  Rndm* rndmPtr{};
190
191  // Pointer to VINCIA DGLAP class.
192  DGLAP* dglapPtr{};
193
194  // Dummy vectors.
195  vector<double> mDum{0, 0, 0, 0};
196  vector<int> hDum{9, 9, 9, 9};
197
198 };
199
200 //==========================================================================
201
202 // Class QQEmitFF, final-final antenna function.
203
204 class QQEmitFF : public AntennaFunction {
205
206 public:
207
208  // Names (remember to redefine both for each inherited class).
209  virtual string vinciaName() const {return "Vincia:QQEmitFF";};
210
211  // Functions needed by soft- and collinear-limit checks (AB -> 0i 1j 2k).
212  virtual int idA() const {return 1;}
213  virtual int idB() const {return -1;}
214  virtual int id1() const {return 21;}
215
216  // The antenna function [GeV^-2].
217  virtual double antFun(vector<double> invariants, vector<double> mNew,
218  vector<int> helBef, vector<int> helNew);
219
220  // Function to give Altarelli-Parisi limits of this antenna.
221  // Defined as PI/sij + PK/sjk, i.e. equivalent to antennae.
222  virtual double AltarelliParisi(vector<double> invariants,
223  vector<double>, vector<int> helBef, vector<int> helNew);
224
225 };
226
227 //==========================================================================
228
229 // Class QGEmitFF, final-final antenna function.
230
231 class QGEmitFF : public AntennaFunction {
232
233 public:
234
235  // Names (remember to redefine both for each inherited class).
236  virtual string vinciaName() const {return "Vincia:QGEmitFF";};
237
238  // Parton types (AB -> 0i 1j 2k): needed by soft- and collinear-limit checks.
239  virtual int idA() const {return 1;}
240  virtual int idB() const {return 21;}
241  virtual int id1() const {return 21;}
242
243  // The antenna function [GeV^-2].
244  virtual double antFun(vector<double> invariants,
245  vector<double> mNew, vector<int> helBef, vector<int> helNew);
246
247  // Function to give Altarelli-Parisi limits of this antenna.
248  virtual double AltarelliParisi(vector<double> invariants,
249  vector<double> /* mNew */, vector<int> helBef, vector<int> helNew);
250
251 };
252
253 //==========================================================================
254
255 // Class GQEmitFF, final-final antenna function.
256
257 class GQEmitFF : public QGEmitFF {
258
259 public:
260
261  // Names (remember to redefine both for each inherited class).
262  virtual string vinciaName() const {return "Vincia:GQEmitFF";};
263
264  // Parton types (AB -> 0i 1j 2k): needed by soft- and collinear-limit checks.
265  virtual int idA() const {return 21;}
266  virtual int idB() const {return -1;}
267  virtual int id1() const {return 21;}
268
269  // The antenna function [GeV^-2] (derived from QGEmit by swapping).
270  virtual double antFun(vector<double> invariants,
271  vector<double> mNew, vector<int> helBef, vector<int> helNew);
272
273  // Function to give Altarelli-Parisi limits of this antenna.
274  virtual double AltarelliParisi(vector<double> invariants,
275  vector<double>, vector<int> helBef, vector<int> helNew);
276
277 };
278
279 //==========================================================================
280
281 // Class GQEmitFF, final-final antenna function.
282
283 class GGEmitFF : public AntennaFunction {
284
285 public:
286
287  // Names (remember to redefine both for each inherited class).
288  virtual string vinciaName() const {return "Vincia:GGEmitFF";};
289
290  // Parton types (AB -> 0i 1j 2k): needed by soft- and collinear-limit checks.
291  virtual int idA() const {return 21;}
292  virtual int idB() const {return 21;}
293  virtual int id1() const {return 21;}
294
295  // The antenna function [GeV^-2].
296  virtual double antFun(vector<double> invariants,
297  vector<double> mNew, vector<int> helBef, vector<int> helNew);
298
299  // Function to give Altarelli-Parisi limits of this antenna.
300  virtual double AltarelliParisi(vector<double> invariants,
301  vector<double>, vector<int> helBef, vector<int> helNew);
302
303 };
304
305 //==========================================================================
306
307 // Class GXSplitFF, final-final antenna function.
308
309 class GXSplitFF : public AntennaFunction {
310
311 public:
312
313  // Names (remember to redefine both for each inherited class).
314  virtual string vinciaName() const {return "Vincia:GXSplitFF";};
315
316  // Parton types (AB -> 0i 1j 2k): needed by soft- and collinear-limit checks.
317  virtual int idA() const {return 21;}
318  virtual int idB() const {return 0;}
319  virtual int id1() const {return -1;}
320
321  // The antenna function [GeV^-2].
322  virtual double antFun(vector<double> invariants,
323  vector<double> mNew, vector<int> helBef, vector<int> helNew);
324
325  // Function to give Altarelli-Parisi limits of this antenna.
326  virtual double AltarelliParisi(vector<double> invariants,
327  vector<double>, vector<int> helBef, vector<int> helNew);
328
329 };
330
331 //==========================================================================
332
333 // Class QQEmitFFsec, sector final-final antenna function, identical
334 // to global one.
335
336 class QQEmitFFsec : public QQEmitFF {};
337
338 //==========================================================================
339
340 // Class QGEmitFFsec, sector final-final antenna function, explicit
341 // symmetrisation of QGEmitFF.
342
343 class QGEmitFFsec : public QGEmitFF {
344
345 public:
346
347  // The antenna function [GeV^-2].
348  virtual double antFun(vector<double> invariants,
349  vector<double> mNew, vector<int> helBef, vector<int> helNew);
350
351 };
352
353 //==========================================================================
354
355 // Class GQEmitFFsec, sector final-final antenna function, explicit
356 // symmetrisation of GQEmitFF.
357
358 class GQEmitFFsec : public QGEmitFFsec {
359
360 public:
361
362  // Parton types (AB -> 0i 1j 2k): needed by soft- and collinear-limit checks.
363  virtual int idA() const {return 21;}
364  virtual int idB() const {return -1;}
365  virtual int id1() const {return 21;}
366
367  // The antenna function [GeV^-2] (derived from QGEmitFFsec by swapping).
368  virtual double antFun(vector<double> invariants,
369  vector<double> mNew, vector<int> helBef, vector<int> helNew);
370
371  // Function to give Altarelli-Parisi limits of this antenna.
372  virtual double AltarelliParisi(vector<double> invariants,
373  vector<double>, vector<int> helBef, vector<int> helNew);
374
375 };
376
377 //==========================================================================
378
379 // Class GGEmitFFsec, sector final-final antenna function, explicit
380 // symmetrisation of GGEmitFF.
381
382 class GGEmitFFsec : public GGEmitFF {
383
384 public:
385
386  // The dimensionless antenna function.
387  virtual double antFun(vector<double> invariants, vector<double> mNew,
388  vector<int> helBef, vector<int> helNew);
389
390 };
391
392 //==========================================================================
393
394 // Class GXSplitFFsec, sector final-final antenna function, explicit
395 // symmetrisation of GXSplitFF.
396
397 class GXSplitFFsec : public GXSplitFF {
398
399  public:
400
401  // The antenna function [GeV^-2].
402  virtual double antFun(vector<double> invariants, vector<double> mNew,
403  vector<int> helBef, vector<int> helNew);
404
405 };
406
407 //==========================================================================
408
409 // Class AntennaFunctionIX, base class for initial-initial and
410 // initial-final antenna functions which implements II. All derived
411 // classes for initial-initial antenna functions are the same for
412 // global and sector cases since even the global ones include sector
413 // terms representing "emission into the initial state".
414
416
417 public:
418
419  // Method to initialise (can be different than that of the base class).
420  virtual bool init();
421
422  // Names (remember to redefine both for each inherited class).
423  virtual string vinciaName() const {return "Vincia:AntennaFunctionIX";}
424
425  // Parton types AB -> 0a 1j 2b with A,B,a,b initial and j final.
426  virtual int idA() const {return 0;}
427  virtual int idB() const {return 0;}
428  virtual int id0() const {return 0;}
429  virtual int id1() const {return 0;}
430  virtual int id2() const {return 0;}
431
432  // Functions to get Altarelli-Parisi energy fractions.
433  virtual double zA(vector<double> invariants) {double sAB = invariants[0];
434  double sjb = invariants[2]; return sAB/(sAB+sjb);}
435  virtual double zB(vector<double> invariants) {double sAB = invariants[0];
436  double saj = invariants[1]; return sAB/(sAB+saj);}
437
438  // Function to tell if this is an II antenna.
439  virtual bool isIIant() {return true;}
440
441  // Function to check singularities, positivity, etc.
442  virtual bool check();
443
444 };
445
446 //==========================================================================
447
448 // Class QQEmitII, initial-initial antenna function.
449
450 class QQEmitII : public AntennaFunctionIX {
451
452 public:
453
454  // Names (remember to redefine both for each inherited class).
455  virtual string vinciaName() const {return "Vincia:QQEmitII";}
456
457  // Parton types AB -> 0a 1j 2b with A,B,a,b initial and j final.
458  virtual int idA() const {return 1;}
459  virtual int idB() const {return -1;}
460  virtual int id0() const {return 1;}
461  virtual int id1() const {return 21;}
462  virtual int id2() const {return -1;}
463
464  // The antenna function [GeV^-2].
465  virtual double antFun(vector<double> invariants, vector<double> masses,
466  vector<int> helBef, vector<int> helNew);
467
468  // AP splitting kernel for collinear limit checks.
469  virtual double AltarelliParisi(vector<double> invariants,
470  vector<double>, vector<int> helBef, vector<int> helNew);
471
472 };
473
474 //==========================================================================
475
476 // Class GQEmitII, initial-initial antenna function.
477
478 class GQEmitII : public AntennaFunctionIX {
479
480 public:
481
482  // Names (remember to redefine both for each inherited class).
483  virtual string vinciaName() const {return "Vincia:GQEmitII";}
484
485  // Parton types AB -> 0a 1j 2b with A,B,a,b initial and j final.
486  virtual int idA() const {return 21;}
487  virtual int idB() const {return 1;}
488  virtual int id0() const {return 21;}
489  virtual int id1() const {return 21;}
490  virtual int id2() const {return 1;}
491
492  // The antenna function.
493  virtual double antFun(vector<double> invariants, vector<double> masses,
494  vector<int> helBef, vector<int> helNew);
495
496  // AP splitting kernel for collinear limit checks.
497  virtual double AltarelliParisi(vector<double> invariants,
498  vector<double>, vector<int> helBef, vector<int> helNew);
499
500 };
501
502 //==========================================================================
503
504 // Class GGEmitII, initial-initial antenna function.
505
506 class GGEmitII : public AntennaFunctionIX {
507
508 public:
509
510  // Names (remember to redefine both for each inherited class).
511  virtual string vinciaName() const {return "Vincia:GGEmitII";}
512
513  // Parton types AB -> 0a 1j 2b with A,B,a,b initial and j final.
514  virtual int idA() const {return 21;}
515  virtual int idB() const {return 21;}
516  virtual int id0() const {return 21;}
517  virtual int id1() const {return 21;}
518  virtual int id2() const {return 21;}
519
520  // The antenna function [GeV^-2].
521  virtual double antFun(vector<double> invariants, vector<double> masses,
522  vector<int> helBef, vector<int> helNew);
523
524  // AP splitting kernel, P(z)/Q2.
525  virtual double AltarelliParisi(vector<double> invariants,
526  vector<double>, vector<int> helBef, vector<int> helNew);
527
528 };
529
530 //==========================================================================
531
532 // Class QXSplitII, initial-initial antenna function. Splitting is in
533 // the forwards sense, i.e. quark backwards evolving to a gluon and
534 // emitting an antiquark in the final state.
535
536 class QXSplitII : public AntennaFunctionIX {
537
538 public:
539
540  // Names (remember to redefine both for each inherited class).
541  virtual string vinciaName() const { return "Vincia:QXSplitII";}
542
543  // Parton types AB -> 0a 1j 2b with A,B, a,b initial and j final.
544  virtual int idA() const {return 1;}
545  virtual int idB() const {return 0;}
546  virtual int id0() const {return 21;}
547  virtual int id1() const {return -1;}
548  virtual int id2() const {return 0;}
549
550  // The antenna function [GeV^-2].
551  virtual double antFun(vector<double> invariants, vector<double> masses,
552  vector<int> helBef, vector<int> helNew);
553
554  // AP splitting kernel, P(z)/Q2.
555  virtual double AltarelliParisi(vector<double> invariants,
556  vector<double>, vector<int> helBef, vector<int> helNew);
557
558  // Mark that this function has no zB collinear limit.
559  virtual double zB(vector<double>) {return -1.0;}
560
561 };
562
563 //==========================================================================
564
565 // Class GXConvII, initial-initial antenna function. Gluon evolves
566 // backwards into a quark and emits a quark in the final state.
567
568 class GXConvII : public AntennaFunctionIX {
569
570 public:
571
572  // Names (remember to redefine both for each inherited class).
573  virtual string vinciaName() const {return "Vincia:GXConvII";}
574
575  // Parton types AB -> 0a 1j 2b with A,B,a,b initial and j final.
576  virtual int idA() const {return 21;}
577  virtual int idB() const {return 0;}
578  virtual int id0() const {return 2;}
579  virtual int id1() const {return 2;}
580  virtual int id2() const {return 0;}
581
582  // The antenna function [GeV^-2].
583  virtual double antFun(vector<double> invariants, vector<double> masses,
584  vector<int> helBef, vector<int> helNew);
585
586  // AP splitting kernel, P(z)/Q2.
587  virtual double AltarelliParisi(vector<double> invariants,
588  vector<double>, vector<int> helBef, vector<int> helNew);
589
590  // Mark that this function has no zB collinear limit.
591  virtual double zB(vector<double>) {return -1.0;}
592
593 };
594
595 //==========================================================================
596
597 // Class AntennaFunctionIF, base class for IF/RF antenna functions
598 // which implements QQEmitIF. Derived classes are for global
599 // initial-final and resonance-final antenna functions. The method
600 // isRFant() distinguishes between the two.
601
603
604 public:
605
606  // Method to initialise (can be different than that of the base class).
607  virtual bool init();
608
609  // Function to check singularities, positivity, etc.
610  virtual bool check();
611
612  // Names (remember to redefine both for each inherited class).
613  virtual string vinciaName() const {return "Vincia:AntennaFunctionIF";}
614
615  // Parton types AB -> 0a 1j 2b with A,a initial and B,b,j final.
616  virtual int idA() const {return 1;}
617  virtual int idB() const {return -1;}
618  virtual int id0() const {return 1;}
619  virtual int id1() const {return 21;}
620  virtual int id2() const {return -1;}
621
622  // Functions to get Altarelli-Parisi energy fractions.
623  virtual double zA(vector<double> invariants) {double sAK = invariants[0];
624  double sjk = invariants[2]; return sAK/(sAK+sjk);}
625  virtual double zB(vector<double> invariants) {double sAK = invariants[0];
626  double saj = invariants[1]; return (sAK-saj)/sAK;}
627
628  // Methods to tell II, IF, and RF apart.
629  virtual bool isIIant() {return false;}
630  virtual bool isRFant() {return false;}
631
632  // Check for resonances.
633  virtual bool checkRes();
634
635  // Create the test masses for the checkRes method.
636  virtual void getTestMasses(vector<double> &masses) {masses.resize(4, 0.0);}
637
638  // Create the test invariants for the checkRes method.
639  virtual bool getTestInvariants(vector<double> &invariants,
640  vector<double> masses, double yaj, double yjk);
641
642 protected:
643
644  // Massive eikonal factor, n.b. is positive definite - proportional
645  // to gram determinant.
646  double massiveEikonal(double saj, double sjk, double sak,
647  double m_a, double m_k) {return 2.0*sak/(saj*sjk) - 2.0*m_a*m_a/(saj*saj)
648  - 2.0*m_k*m_k/(sjk*sjk);}
649
650  // Massive eikonal factor, given invariants and masses.
651  double massiveEikonal(vector<double> invariants, vector<double> masses) {
652  return massiveEikonal(invariants[1], invariants[2], invariants[3],
653  masses[0], masses[2]);}
654
655  // Return the Gram determinant.
656  double gramDet(vector<double> invariants, vector<double> masses) {
657  double saj(invariants[1]), sjk(invariants[2]), sak(invariants[3]),
658  mares(masses[0]), mjres(masses[1]), mkres(masses[2]);
659  return 0.25*(saj*sjk*sak - saj*saj*mkres*mkres -sak*sak*mjres*mjres
660  - sjk*sjk*mares*mares + 4.0*mares*mares*mjres*mjres*mkres*mkres);}
661
662  // Wrapper for comparing to AP functions, sums over flipped
663  // invariants where appropriate.
664  double antFunCollLimit(vector<double> invariants,vector<double> masses);
665
666 };
667
668 //==========================================================================
669
670 // Class QQEmitIF, initial-final antenna function.
671
672 class QQEmitIF : public AntennaFunctionIF {
673
674 public:
675
676  // Names (remember to redefine both for each inherited class).
677  virtual string vinciaName() const { return "Vincia:QQEmitIF";}
678
679  // Parton types AB -> 0a 1j 2b with A,a initial and B,b,j final.
680  virtual int idA() const {return 1;}
681  virtual int idB() const {return -1;}
682  virtual int id0() const {return 1;}
683  virtual int id1() const {return 21;}
684  virtual int id2() const {return -1;}
685
686  // The antenna function [GeV^-2].
687  virtual double antFun(vector<double> invariants, vector<double> masses,
688  vector<int> helBef, vector<int> helNew);
689
690  // The AP kernel, P(z)/Q2.
691  virtual double AltarelliParisi(vector<double> invariants,
692  vector<double>, vector<int> helBef, vector<int> helNew);
693
694  // Functions to get Altarelli-Parisi energy fractions.
695  virtual double zA(vector<double> invariants) {double sAK = invariants[0];
696  double sjk = invariants[2]; return sAK/(sAK+sjk);}
697  virtual double zB(vector<double> invariants) {double sAK = invariants[0];
698  double saj = invariants[1]; return (sAK-saj)/sAK;}
699
700 };
701
702 //==========================================================================
703
704 // Class QGEmitIF, initial-final antenna function.
705
706 class QGEmitIF : public AntennaFunctionIF {
707
708 public:
709
710  // Names (remember to redefine both for each inherited class).
711  virtual string vinciaName() const {return "Vincia:QGEmitIF";}
712
713  // Parton types AB -> 0a 1j 2b with A,a initial and B,b,j final.
714  virtual int idA() const {return 1;}
715  virtual int idB() const {return 21;}
716  virtual int id0() const {return 1;}
717  virtual int id1() const {return 21;}
718  virtual int id2() const {return 21;}
719
720  // The antenna function [GeV^-2].
721  virtual double antFun(vector<double> invariants, vector<double> masses,
722  vector<int> helBef, vector<int> helNew);
723
724  // The AP kernel, P(z)/Q2.
725  virtual double AltarelliParisi(vector<double> invariants,
726  vector<double>, vector<int> helBef, vector<int> helNew);
727
728 };
729
730 //==========================================================================
731
732 // Class GQEmitIF, initial-final antenna function.
733
734 class GQEmitIF : public AntennaFunctionIF {
735
736 public:
737
738  // Names (remember to redefine both for each inherited class).
739  virtual string vinciaName() const {return "Vincia:GQEmitIF";}
740
741  // Parton types AB -> 0a 1j 2b with A,a initial and B,b,j final.
742  virtual int idA() const {return 21;}
743  virtual int idB() const {return 1;}
744  virtual int id0() const {return 21;}
745  virtual int id1() const {return 21;}
746  virtual int id2() const {return 1;}
747
748  // The antenna function [GeV^-2].
749  virtual double antFun(vector<double> invariants, vector<double> masses,
750  vector<int> helBef, vector<int> helNew);
751
752  // The AP kernel, P(z)/Q2.
753  virtual double AltarelliParisi(vector<double> invariants,
754  vector<double>, vector<int> helBef, vector<int> helNew);
755
756 };
757
758 //==========================================================================
759
760 // Class GGEmitIF, initial-final antenna function.
761
762 class GGEmitIF : public AntennaFunctionIF {
763
764 public:
765
766  // Names (remember to redefine both for each inherited class).
767  virtual string vinciaName() const {return "Vincia:GGEmitIF";}
768
769  // Parton types AB -> 0a 1j 2b with A,a initial and B,b,j final.
770  virtual int idA() const {return 21;}
771  virtual int idB() const {return 21;}
772  virtual int id0() const {return 21;}
773  virtual int id1() const {return 21;}
774  virtual int id2() const {return 21;}
775
776  // The antenna function [GeV^-2].
777  virtual double antFun(vector<double> invariants, vector<double> masses,
778  vector<int> helBef, vector<int> helNew);
779
780  // The AP kernel, P(z)/Q2.
781  virtual double AltarelliParisi(vector<double> invariants,
782  vector<double>, vector<int> helBef, vector<int> helNew);
783
784 };
785
786 //==========================================================================
787
788 // Class QXSplitIF, initial-final antenna function. Splitting is in
789 // the forwards sense, i.e. quark backwards evolving to a gluon and
790 // emitting an antiquark in the final state.
791
792 class QXSplitIF : public AntennaFunctionIF {
793
794 public:
795
796  // Names (remember to redefine both for each inherited class).
797  virtual string vinciaName() const {return "Vincia:QXSplitIF";}
798
799  // Parton types AB -> 0a 1j 2b with A,a initial and B,b,j final.
800  virtual int idA() const {return 1;}
801  virtual int idB() const {return 0;}
802  virtual int id0() const {return 21;}
803  virtual int id1() const {return -1;}
804  virtual int id2() const {return 0;}
805
806  // The antenna function [GeV^-2].
807  virtual double antFun(vector<double> invariants, vector<double> masses,
808  vector<int> helBef, vector<int> helNew);
809
810  virtual double AltarelliParisi(vector<double> invariants,
811  vector<double> /* mNew */, vector<int> helBef, vector<int> helNew);
812
813  // Mark that this function does not have a zB collinear limit.
814  virtual double zB(vector<double>) {return -1.0;}
815
816 };
817
818 //==========================================================================
819
820 // Class GXConvIF, initial-final antenna function. Gluon evolves
821 // backwards into a quark and emits a quark in the final state.
822
823 class GXConvIF : public AntennaFunctionIF {
824
825 public:
826
827  // Names (remember to redefine both for each inherited class).
828  virtual string vinciaName() const {return "Vincia:GXConvIF";}
829
830  // Parton types AB -> 0a 1j 2b with A,a initial and B,b,j final.
831  virtual int idA() const {return 21;}
832  virtual int idB() const {return 0;}
833  virtual int id0() const {return 2;}
834  virtual int id1() const {return 2;}
835  virtual int id2() const {return 0;}
836
837  // The antenna function [GeV^-2].
838  virtual double antFun(vector<double> invariants, vector<double> masses,
839  vector<int> helBef, vector<int> helNew);
840
841  // The AP kernel, P(z)/Q2.
842  virtual double AltarelliParisi(vector<double> invariants,
843  vector<double>, vector<int> helBef, vector<int> helNew);
844
845  // Mark that this function does not have a zB collinear limit.
846  virtual double zB(vector<double>) {return -1.0;}
847
848 };
849
850 //==========================================================================
851
852 // Class XGSplitIF, initial-final antenna function. Gluon splitting in
853 // the final state.
854
855 class XGSplitIF : public AntennaFunctionIF {
856
857 public:
858
859  // Names (remember to redefine both for each inherited class).
860  virtual string vinciaName() const {return "Vincia:XGsplitIF";}
861
862  // Parton types AB -> 0a 1j 2b with A,a initial and B,b,j final.
863  virtual int idA() const {return 0;}
864  virtual int idB() const {return 21;}
865  virtual int id0() const {return 0;}
866  virtual int id1() const {return -1;}
867  virtual int id2() const {return 1;}
868
869  // The antenna function [GeV^-2].
870  virtual double antFun(vector<double> invariants, vector<double> masses,
871  vector<int> helBef, vector<int> helNew);
872
873  // The AP kernel, P(z)/Q2.
874  virtual double AltarelliParisi(vector<double> invariants,
875  vector<double>, vector<int> helBef, vector<int> helNew);
876
877  // Mark that this function does not have a zA collinear limit.
878  virtual double zA(vector<double>) {return -1.0;}
879
880 };
881
882 //==========================================================================
883
884 // Class QGEmitIFsec, derived class for sector initial-final antenna
885 // function. Note only the final-state leg needs to be symmetrised,
886 // as the global IF functions already contain sector terms on their
887 // initial-state legs to account for the absence of "emission into the
888 // initial state".
889
890 class QGEmitIFsec : public QGEmitIF {
891
892 public:
893
894  // The antenna function [GeV^-2].
895  virtual double antFun(vector<double> invariants,
896  vector<double> mNew, vector<int> helBef, vector<int> helNew);
897
898 };
899
900 //==========================================================================
901
902 // Class GGEmitIFsec, sector initial-final antenna function.
903
904 class GGEmitIFsec : public GGEmitIF {
905
906 public:
907
908  // The antenna function [GeV^-2].
909  virtual double antFun(vector<double> invariants,
910  vector<double> mNew, vector<int> helBef, vector<int> helNew);
911
912 };
913
914 //==========================================================================
915
916 // Class XGSplitIFsec, sector initial-final antenna function. Gluon
917 // splitting in the final state.
918
919 class XGSplitIFsec : public XGSplitIF {
920
921 public:
922
923  // The antenna function, just 2*global [GeV^-2].
924  virtual double antFun(vector<double> invariants, vector<double> mNew,
925  vector<int> helBef, vector<int> helNew);
926
927 };
928
929 //==========================================================================
930
931 // Class QQEmitRF, resonance-final antenna function.
932
933 class QQEmitRF : public QQEmitIF {
934
935 public:
936
937  // Names (remember to redefine both for each inherited class).
938  string vinciaName() const {return "Vincia:QQEmitRF";}
939
940  // Parton types AB -> ijk with A,i initial and B,k,j final.
941  int idA() const {return 6;}
942  int idB() const {return 5;}
943  int id0() const {return 6;}
944  int id1() const {return 21;}
945  int id2() const {return 5;}
946
947  // Mark that this function does not have a zA collinear limit.
948  double zA(vector<double>) {return -1;}
949
950  // Return this is a resonance-final antenna.
951  bool isRFant() {return true;}
952
953  // Test masses (top, gluon, bottom, and W mass).
954  void getTestMasses(vector<double> &masses) {masses = {particleDataPtr->m0(6),
955  0.0, particleDataPtr->m0(5), particleDataPtr->m0(24)};}
956
957  // AP with dummy helicities.
958  virtual double AltarelliParisi(vector<double> invariants,
959  vector<double> masses, vector<int>, vector<int>) {
960  double sjk(invariants[2]), mkres(masses[2]), z(zB(invariants)),
961  mu2(mkres*mkres/sjk), Pz(dglapPtr->Pq2gq(z,9,9,9,mu2));
962  return Pz/sjk;};
963
964 };
965
966 //==========================================================================
967
968 // Class QGEmitRF, resonance-final antenna function.
969
970 class QGEmitRF : public QGEmitIF {
971
972 public:
973
974  // Names (remember to redefine both for each inherited class).
975  string vinciaName() const {return "Vincia:QGEmitRF";}
976
977  // Parton types AB -> ijk with A,i initial and B,k,j final.
978  int idA() const {return 6;}
979  int idB() const {return 21;}
980  int ida() const {return 6;}
981  int idb() const {return 21;}
982  int id1() const {return 21;}
983
984  // Mark that this function does not have a zA collinear limit.
985  double zA(vector<double>) {return -1;}
986
987  // Return this is a resonance-final antenna.
988  bool isRFant() {return true;}
989
990  // Test masses (top, gluon, gluon, X).
991  void getTestMasses(vector<double> &masses) {masses = {particleDataPtr->m0(6),
992  0.0, 0.0, 0.6*particleDataPtr->m0(6)};}
993
994  // AP with dummy helicities and masses.
995  virtual double AltarelliParisi(vector<double> invariants, vector<double>,
996  vector<int>, vector<int>) {
997  double sjk(invariants[2]), z(zB(invariants)),
998  Pz(dglapPtr->Pg2gg(z, 9, 9, 9));
999  return Pz/sjk;}
1000
1001 };
1002
1003 //==========================================================================
1004
1005 // Class XGSplitRF, resonance-final antenna function.
1006
1007 class XGSplitRF : public XGSplitIF {
1008
1009 public:
1010
1011  // Names (remember to redefine both for each inherited class)
1012  string vinciaName() const {return "Vincia:XGSplitRF";}
1013
1014  // Mark that this function does not have a zA collinear limit.
1015  double zA(vector<double>){ return -1;}
1016
1017  // Return this is a resonance-final antenna.
1018  bool isRFant() {return true;}
1019
1020  // Parton types AB -> ijk with A,i initial and B,k,j final.
1021  int idA() const {return 6;}
1022  int idB() const {return 21;}
1023  int ida() const {return 6;}
1024  int idb() const {return 2;}
1025  int id1() const {return -2;}
1026
1027  // Test masses (top, gluon, gluon, X).
1028  void getTestMasses(vector<double> &masses) {masses = {particleDataPtr->m0(6),
1029  0.0, 0.0, 0.6*particleDataPtr->m0(6)};}
1030
1031  // AP with dummy helicities.
1032  double AltarelliParisi(vector<double> invariants, vector<double> masses,
1033  vector<int>, vector<int>) {
1034  double sAK(invariants[0]), saj(invariants[1]), sjk(invariants[2]),
1035  mkres(masses[2]), m2q(mkres*mkres), Q2(sjk + 2*m2q), mu2(m2q/Q2),
1036  z((sAK+saj-Q2)/sAK), Pz(dglapPtr->Pg2qq(z, 9, 9, 9, mu2));
1037  return Pz/Q2;}
1038
1039 };
1040
1041 //==========================================================================
1042
1043 // The AntennaSetFSR class. Simple container of FF and RF antenna functions.
1044
1046
1047 public:
1048
1049  // Default constructor.
1050  AntennaSetFSR() = default;
1051
1052  // Destructor, delete the antennae.
1053  virtual ~AntennaSetFSR() {
1054  for (map<int, AntennaFunction*>::iterator it = antFunPtrs.begin();
1055  it != antFunPtrs.end(); ++it) delete it->second;
1056  antFunPtrs.clear();}
1057
1058  // Initialize pointers.
1059  void initPtr(Info* infoPtrIn, DGLAP* dglapPtrIn);
1060
1061  // Initialize antenna set (optionally with min or max variation).
1062  void init();
1063
1064  // Function to chek if an antenna with the given index exists.
1065  bool exists(int iAntIn) {return antFunPtrs.count(iAntIn);}
1066
1067  // Gets an antenna iAntIn from the AntennaSet.
1068  AntennaFunction* getAnt(int iAntIn) {
1069  return (exists(iAntIn)) ? antFunPtrs[iAntIn] : nullptr;}
1070
1071  // Method to return all iAntPhys values that are defined in antFunPtr.
1072  vector<int> getIant();
1073
1074  // Get Vincia name, e.g. "Vincia:QQEmitFF".
1075  string vinciaName(int iAntIn) {
1076  return exists(iAntIn) ? antFunPtrs[iAntIn]->vinciaName() : "noVinciaName";}
1077
1078  // Get human name, e.g. "g/qq".
1079  string humanName(int iAntIn) {
1080  return exists(iAntIn) ? antFunPtrs[iAntIn]->humanName() : "noHumanName";}
1081
1082 private:
1083
1084  // Use a map of AntennaFunction pointers, create them with new on
1085  // initialization.
1086  map<int,AntennaFunction*> antFunPtrs{};
1087
1088  // Pointers to Pythia8 classes, needed to initialise antennae.
1089  bool isInitPtr{false}, isInit{false};
1090  Info* infoPtr{};
1091  ParticleData* particleDataPtr{};
1092  Settings* settingsPtr{};
1093  Rndm* rndmPtr{};
1094
1095  // Pointer to VINCIA DGLAP class.
1096  DGLAP* dglapPtr{};
1097
1098  // Verbosity level
1099  int verbose{};
1100
1101 };
1102
1103 //==========================================================================
1104
1105 // The AntennaSetISR class. Simple container of II and IF antenna functions.
1106
1108
1109  public:
1110
1111  // Default constructor.
1112  AntennaSetISR() = default;
1113
1114  // Destructor, delete the antennae.
1115  ~AntennaSetISR() {
1116  for (map<int, AntennaFunctionIX*>::iterator it = antFunPtrs.begin();
1117  it != antFunPtrs.end(); ++it) delete it->second;
1118  antFunPtrs.clear();}
1119
1120  // Initialize pointers.
1121  void initPtr(Info* infoPtrIn, DGLAP* dglapPtrIn);
1122
1123  // Initialize antenna set.
1124  void init();
1125
1126  // Function to chek if an antenna with the given index exists.
1127  bool exists(int iAntIn) {return antFunPtrs.count(iAntIn);}
1128
1129  // Gets an antenna iAntIn from the AntennaSetISR.
1130  AntennaFunctionIX* getAnt(vector<AntennaFunctionIX*>::size_type iAntIn) {
1131  return (exists(iAntIn)) ? antFunPtrs[iAntIn] : nullptr;}
1132
1133  // Method to return all iAntPhys values that are defined in antFunPtr.
1134  vector<int> getIant();
1135
1136  // Get Vincia name, e.g. "Vincia:QQEmitII".
1137  string vinciaName(vector<AntennaFunctionIX*>::size_type iAntIn) {
1138  return exists(iAntIn) ? antFunPtrs[iAntIn]->vinciaName() : "noVinciaName";}
1139
1140  // Get human name, e.g. "g/qq".
1141  string humanName(int iAntIn) {
1142  return exists(iAntIn) ? antFunPtrs[iAntIn]->humanName() : "noHumanName";}
1143
1144 private:
1145
1146  // Use a map of AntennaFunction pointers, create them with new on
1147  // initialization.
1148  map<int,AntennaFunctionIX*> antFunPtrs{};
1149
1150  // Pointers to Pythia 8 classes, needed to initialise antennae.
1151  bool isInitPtr{false}, isInit{false};
1152  Info* infoPtr{};
1153  ParticleData* particleDataPtr{};
1154  Settings* settingsPtr{};
1155  Rndm* rndmPtr{};
1156
1157  // Pointer to VINCIA DGLAP class
1158  DGLAP* dglapPtr{};
1159
1160  // Verbosity level
1161  int verbose{};
1162
1163 };
1164
1165 //==========================================================================
1166
1167 // Class MECs, for computing matrix-element corrections to antenna
1168 // functions.
1169
1170 class MECs {
1171
1172 public:
1173
1174  // Constructor.
1175  MECs() {isInitPtr = false; isInit = false;}
1176
1177  // Destructor.
1178  virtual ~MECs() {};
1179
1180  // Initialize pointers.
1181  void initPtr(Info* infoPtrIn, ShowerMEs* mg5mesPtrIn,
1182  VinciaCommon* vinComPtrIn);
1183
1184  // Initialize pointers to antenna sets.
1185  void initAntPtr(AntennaSetFSR* antFSRusr, AntennaSetISR* antISRusr) {
1186  antSetFSR = antFSRusr; antSetISR = antISRusr;}
1187
1188  // Initialize.
1189  void init();
1190
1191  // Function to return ME class (Born type) for a parton
1192  // configuration. Can be called from either of the ISR::prepare() or
1193  // FSR::prepare() functions, or from the ISR::branch() or
1194  // FSR::branch() functions. Returns >= 0 if there an ME for this
1195  // configuration, associated with the (arbitrary) integer code label
1196  // denoted by the return value. If return < 0 we don't have an ME /
1197  // no ME should be used for this system.
1198  bool prepare(const int iSys, Event& event);
1199
1200  // Function to assign helicities to particles (using MEs).
1201  bool polarise(const int iSys, Event& event);
1202
1203  // Make list of particles as vector<Particle>.
1204  // First 1 or 2 entries : incoming particle(s).
1205  // Subseqent entries : outgoing particles.
1206  // The two last arguments are optional and allow to specify a list
1207  // of indices to be ignored, and a set of particles to be added, e.g.
1208  // in the context of setting up a trial state after a branching.
1209  // The newly added particles are then at the end of the respective
1210  // lists, i.e. a newly added incoming particle is the last incoming
1211  // one and newly added outgoing ones are the last among the outgoing
1212  // ones.
1213  vector<Particle> makeParticleList(const int iSys, const Event& event,
1214  const vector<Particle> pNew = vector<Particle>(),
1215  const vector<int> iOld = vector<int>());
1216
1217  // Check if state already has helicities. Return true if any
1218  // particle in state already has a specified helicity.
1219  bool isPolarised(int iSys, Event& event, bool checkIncoming);
1220
1221  // Wrapper function to return a specific antenna function.
1222  AntennaFunction* getAntFSR(const int iAnt) {
1223  return antSetFSR->getAnt(iAnt); }
1224  AntennaFunctionIX* getAntISR(const int iAnt) {
1225  return antSetISR->getAnt(iAnt); }
1226
1227  // Function to determine if MECs are requested at this order for this system.
1228  bool doMEC(const int iSys, const int nBranch);
1229
1230  // Get squared matrix element.
1231  double getME2(const vector<Particle> parts, const int nIn);
1232  double getME2(const int iSys, const Event& event);
1233
1234  // Get matrix element correction factor (TODO: dummy implementation).
1235  double getMEC(const int, const Event&) {return 1.;}
1236
1237  // Return number of partons added since Born (as defined by prepare).
1238  int sizeOutBorn(const int iSys) {return sizeOutBornSav[iSys];}
1239
1240  // Function to set level of verbose output.
1241  void setVerbose(const int verboseIn) {verbose = verboseIn;}
1242
1245
1246 private:
1247
1248  // Verbosity level.
1249  int verbose;
1250
1251  // Is initialized.
1252  bool isInitPtr, isInit;
1253
1254  // Pointers to PYTHIA objects.
1255  Info* infoPtr;
1256  Settings* settingsPtr;
1257  Rndm* rndmPtr;
1258  PartonSystems* partonSystemsPtr;
1259  ShowerMEs* mg5mesPtr;
1260
1261  // Pointers to VINCIA objects.
1262  VinciaCommon* vinComPtr;
1263
1264  // Antenna sets.
1265  AntennaSetFSR* antSetFSR;
1266  AntennaSetISR* antSetISR;
1267
1268  // Matching settings.
1269  bool matchingFullColour;
1270  int maxMECs2to1, maxMECs2to2, maxMECs2toN, maxMECsResDec, maxMECsMPI;
1271  int nFlavZeroMass;
1272
1273  // Map from iSys to Born multiplicity and ME class; set in prepare().
1274  map<int,int> sizeOutBornSav;
1275
1276 };
1277
1278 //==========================================================================
1279
1280 } // end namespace Pythia8
1281
1282 #endif // end Pythia8_VinciaAntennaFunctions_H