StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 
14 // Pythia headers.
15 #include "Pythia8/Basics.h"
16 #include "Pythia8/Event.h"
17 #include "Pythia8/PythiaStdlib.h"
18 #include "Pythia8/ShowerMEs.h"
19 
20 // Vincia headers.
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 
65 class AntennaFunction {
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 
415 class AntennaFunctionIX : public AntennaFunction {
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 
602 class AntennaFunctionIF : public AntennaFunctionIX {
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 
1045 class AntennaSetFSR {
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 
1107 class AntennaSetISR {
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 
1243  // Header.
1244  void header();
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
Generated by   doxygen 1.8.5