StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
VinciaAntennaFunctions.cc
1 // VinciaAntennaFunctions.cc 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 // Function definitions (not found in the header) for the AntennaSet,
7 // AntennaFunction, and MECs classes for both FSR and ISR, and some
8 // related global functions.
9 
10 #include "Pythia8/VinciaAntennaFunctions.h"
11 
12 namespace Pythia8 {
13 
14 //==========================================================================
15 
16 // DGLAP splitting functions.
17 
18 //--------------------------------------------------------------------------
19 
20 // DGLAP splitting kernels with helicity-dependence (from Larkoski,
21 // Peskin, http://arxiv.org/abs/arXiv:0908.2450). Note, Pg2gg is
22 // written with the standard factor 2 normalization convention and
23 // gives a factor 2 difference with respect PYTHIA FSR. The
24 // difference is related to whether P is interpreted as the
25 // probability for gluon A to branch to BC (without factor 2), or as
26 // the probability for gluon B to have originated from A (with factor
27 // 2). Both quarks and gluons have definite helicities:
28 // 9: unpolarized
29 // +1: right-handed
30 // -1: left-handed
31 // where z is momentum fraction taken by B in A -> BC.
32 
33 //--------------------------------------------------------------------------
34 
35 // Pg2gg, written with the standard factor 2 normalization convention.
36 
37 double DGLAP::Pg2gg(double z, int hA, int hB, int hC) {
38 
39  // If A unpolarized, treat as unpolarized.
40  if (hA == 9) return 2*pow2(1. - z*(1. - z))/z/(1. - z);
41 
42  // Expressions coded for hA = 1. Flip helicities if hA = -1.
43  if (hA == -1) {hB *= -1; hC *= -1;}
44  if (hB == 1 && hC == 1) return 1./z/(1. - z);
45  else if (hB == -1 && hC == 1) return pow3(1. - z)/z;
46  else if (hB == 1 && hC == -1) return pow3(z)/(1. - z);
47  return 0.;
48 
49 }
50 
51 //--------------------------------------------------------------------------
52 
53 // Pg2qq, B is quark, C is antiquark. Note, mass corrections not
54 // implemented.
55 
56 double DGLAP::Pg2qq(double z, int hA, int hB, int hC, double mu2) {
57 
58  // If A unpolarized, treat as unpolarized.
59  if (hA == 9) return pow2(z) + pow2(1. - z) + 2.0*mu2;
60 
61  // Preserve quark helicity.
62  if (hB != -hC || abs(hB) != 1) return 0.;
63 
64  // Expressions coded for hA = 1. Flip helicities if hA = -1.
65  if (hA == -1) {hB *= -1; hC *= -1;}
66  if (hB == -1 && hC == 1) return pow2(1. - z);
67  else if (hB == 1 && hC == -1) return pow2(z);
68  return 0.;
69 
70 }
71 
72 //--------------------------------------------------------------------------
73 
74 // Pq2qg, B is quark, C is gluon. Note, mass corrections not
75 // implemented for polarised.
76 
77 double DGLAP::Pq2qg(double z, int hA, int hB, int hC, double mu2) {
78 
79  // If A unpolarized, treat as unpolarized.
80  if (hA == 9) return (1. + pow2(z))/(1. - z) - 2.0*mu2;
81 
82  // Preserve quark helicity
83  if (hA != hB || abs(hB) != 1) return 0.;
84 
85  // Expressions coded for hA = +1. Flip helicities if hA = -1.
86  if (hA == -1) {hB *= -1; hC *= -1;}
87  if (hB == 1 && hC == -1) return pow2(z)/(1. - z);
88  else if (hB == 1 && hC == 1) return 1./(1. - z);
89  return 0.;
90 
91 }
92 
93 //--------------------------------------------------------------------------
94 
95 // Pq2qg, B is gluon, C is quark. Note, mass corrections not
96 // implemented for polarised.
97 
98 double DGLAP::Pq2gq(double z, int hA, int hB, int hC, double mu) {
99  return Pq2qg(1-z,hA,hC,hB,mu);}
100 
101 //--------------------------------------------------------------------------
102 
103 // DGLAP splitting kernels with polarization-dependence (from Ellis,
104 // Stirling, Webber) Gluons are plane polarized, quarks have
105 // helicities:
106 // 9: unpolarized
107 // +1: in plane / right-handed
108 // -1: out-of-plane / left-handed
109 // where z is momentum fraction taken by B in A -> BC.
110 
111 //--------------------------------------------------------------------------
112 
113 // Pg2ggLin, written without factor 2 normalization convention.
114 
115 double DGLAP::Pg2ggLin(double z, int polA, int polB, int polC) {
116 
117  // If A unpolarized, treat as unpolarized.
118  if (polA == 9) return (1. - z + pow2(z))/z/(1 - z);
119 
120  // A in-plane polarized.
121  else if (polA == 1) {
122  if (polB == 1 && polC == 1) return (1. - z)/z + z/(1. - z) + z*(1. - z);
123  else if (polB == -1 && polC == -1) return z*(1. - z);
124  }
125 
126  // A out-of-plane polarized.
127  else if (polA == -1) {
128  if (polB == 1 && polC == -1) return (1. - z)/z;
129  else if (polB == -1 && polC == 1) return z/(1. - z);
130  }
131  return 0.;
132 
133 }
134 
135 //--------------------------------------------------------------------------
136 
137 // Pg2qqLin, B is quark, C is antiquark. Note, mass corrections not
138 // implemented.
139 
140 double DGLAP::Pg2qqLin(double z, int polA, int hB, int hC, double mu) {
141 
142  // If A unpolarized, treat as unpolarized.
143  if (polA == 9) return Pg2qq(z,9,9,9,mu);
144 
145  // Preserve quark helicity.
146  if (hB != -hC || abs(hB) != 1) return 0.;
147 
148  // A in-plane polarized.
149  else if (polA == 1) return pow2(1. - 2*z);
150 
151  // A out-of-plane polarized.
152  else if (polA == -1) return 1.;
153  return 0.;
154 
155 }
156 
157 //--------------------------------------------------------------------------
158 
159 // Pq2qgLin, B is quark, C is gluon. Note, mass corrections not
160 // implemented.
161 
162 double DGLAP::Pq2qgLin(double z, int hA, int hB, int polC, double mu) {
163 
164  // If A unpolarized, treat as unpolarized.
165  if (hA == 9) return Pq2qg(z, 9, 9, 9, mu);
166 
167  // Preserve quark helicity.
168  if (hA != hB || abs(hB) != 1) return 0.;
169 
170  // Gluon polarized in plane.
171  else if (polC == 1) return pow2(1. + z)/(1. - z);
172 
173  // Gluon polarized out of plane.
174  else if (polC == -1) return 1. - z;
175  return 0.;
176 
177 }
178 
179 //--------------------------------------------------------------------------
180 
181 // Pq2qgLin, B is gluon, C is quark. Note, mass corrections not
182 // implemented.
183 
184 double DGLAP::Pq2gqLin(double z, int hA, int polB, int hC, double mu) {
185  return Pq2qgLin(1 - z, hA, hC, polB, mu);}
186 
187 //==========================================================================
188 
189 // The AntennaFunction base class.
190 
191 //--------------------------------------------------------------------------
192 
193 // Default initialization.
194 
195 bool AntennaFunction::init() {
196 
197  // Check whether pointers are initialized.
198  if (!isInitPtr) return false;
199 
200  // Verbosity level.
201  verbose = settingsPtr->mode("Vincia:verbose");
202  // Charge factor
203  // (Use same charge factor for QGEmit and GQEmit)
204  if (vinciaName() == "Vincia:GQEmitFF")
205  chargeFacSav = settingsPtr->parm("Vincia:QGEmitFF:chargeFactor");
206  else
207  chargeFacSav = settingsPtr->parm(vinciaName() + ":chargeFactor");
208  if (chargeFacSav < 0.) chargeFacSav = 0.0;
209 
210  // Subleading-colour treatment.
211  // modeSLC = 0: all gluon-emission antennae normalised to CA.
212  // modeSLC = 1: use colour factors as specified by user.
213  // modeSLC = 2: QQ gets CF, GG gets CA, QG gets interpolation.
214  modeSLC = settingsPtr->mode("Vincia:modeSLC");
215  if (modeSLC == 0 && id1() == 21) chargeFacSav = CA;
216  if (modeSLC == 2 && id1() == 21) {
217  if (idA() == 21 && idB() == 21) chargeFacSav = CA;
218  else if (idA() == 21 || idB() == 21) chargeFacSav = (CA + CF)/2.;
219  else chargeFacSav = CF;
220  }
221 
222  // Kinematics map type (first look for specific, else default to global).
223  if (settingsPtr->isMode(vinciaName() + ":kineMap"))
224  kineMapSav = settingsPtr->mode(vinciaName()+":kineMap");
225  else {
226  // Gluon emission antennae.
227  if (id1() == 21) {
228  kineMapSav = settingsPtr->mode("Vincia:kineMapFFemit");
229 
230  // Gluon splitting antennae.
231  } else {
232  kineMapSav = settingsPtr->mode("Vincia:kineMapFFsplit");
233  // For gluon splittings, parton I is always emitter, K is recoiler.
234  if (kineMapSav == 2) kineMapSav = -1;
235  }
236  }
237 
238  // Collinear partitioning (for global antennae).
239  alphaSav = settingsPtr->parm("Vincia:octetPartitioning");
240 
241  // Sector shower on/off and sectorDamp parameter.
242  sectorShower = settingsPtr->flag("Vincia:sectorShower");
243  sectorDampSav = settingsPtr->parm("Vincia:sectorDamp");
244 
245  // Return.
246  isInit = true;
247  return isInit;
248 
249 }
250 
251 //--------------------------------------------------------------------------
252 
253 // Method to initialise internal helicity variables. Return value =
254 // number of helicity configurations to average over.
255 
256 int AntennaFunction::initHel(vector<int>* helBef, vector<int>* helNew) {
257 
258  // Initialise as unpolarised.
259  hA = 9; hB = 9; hi = 9; hj = 9; hk = 9;
260 
261  // Check if one or more partons have explicit helicities.
262  if (helNew->size() >= 3) {
263  hi = helNew->at(0); hj = helNew->at(1); hk = helNew->at(2);
264  }
265  if (helBef->size() >= 2) {
266  hA = helBef->at(0); hB = helBef->at(1);
267  }
268 
269  // Check if helicity assignments make sense.
270  bool physHel = true;
271  if (hA != 1 && hA != -1 && hA != 9) physHel = false;
272  if (hB != 1 && hB != -1 && hB != 9) physHel = false;
273  if (hi != 1 && hi != -1 && hi != 9) physHel = false;
274  if (hj != 1 && hj != -1 && hj != 9) physHel = false;
275  if (hk != 1 && hk != -1 && hk != 9) physHel = false;
276  if (!physHel) {
277  if (verbose >= normal){
278  stringstream ss;
279  ss << hA << " " << hB << " -> " << hi << " " << hj << " " << hk;
280  infoPtr->errorMsg("Warning in "+__METHOD_NAME__+
281  ": unphysical helicity configuration.",ss.str());
282  }
283  return 0;
284  }
285 
286  // How many configurations are we averaging over.
287  int nAvg = 1;
288  if (hA == 9) nAvg *= 2;
289  if (hB == 9) nAvg *= 2;
290  return nAvg;
291 
292 }
293 
294 //--------------------------------------------------------------------------
295 
296 // Generic check of antenna function for use during initialisation.
297 
298 bool AntennaFunction::check() {
299 
300  // Check soft singularity against massless dipole factor.
301  bool isOK = true;
302  if (id1() == 21 || id1() == 22) {
303  vector<double> invariants;
304 
305  // Test one symmetric and one asymmetric soft phase-space point.
306  for (int iTest = 0; iTest <= 1; ++iTest) {
307  // Test invariants.
308  double yij = 1.0e-3 * pow(10.0,iTest);
309  double yjk = 1.0e-3 / pow(10.0,iTest);
310  // sIK is arbitrary, going to divide out in antFun anyway,
311  // inclusion is just for consistent notation.
312  double sIK = 1.0;
313  invariants.resize(3);
314  invariants[0] = sIK;
315  invariants[1] = yij*sIK;
316  invariants[2] = yjk*sIK;
317  // Compute eikonal.
318  double eik = 2.*(1. - yij - yjk)/yij/yjk*(1./sIK);
319  // Compute antenna (without subleading-colour corrections).
320  int modeSLCsave = modeSLC;
321  modeSLC = 1;
322  double ant = antFun(invariants);
323  modeSLC = modeSLCsave;
324  // Compare.
325  double ratio = ant/eik;
326  if (abs(ratio - 1.) >= 0.001) {
327  isOK = false;
328  if (verbose >= quiet) {
329  infoPtr->errorMsg("Warning in "+__METHOD_NAME__+": Failed eikonal.",
330  "("+num2str(iTest,1) + ")");
331  }
332  } else if (verbose >= verylouddebug) printOut(__METHOD_NAME__,
333  vinciaName() + " OK (eikonal " + num2str(iTest, 1) + ")");
334  }
335  }
336 
337  // Test all helicity configurations for collinearity and posivity.
338  string helString = " 9 9 -> 9 9 9";
339  for (int iHel = -1; iHel < 32; ++iHel) {
340  if (iHel >= 0) {
341  hj = 1 - 2*(iHel%2);
342  hi = 1 - 2*((iHel/2)%2);
343  hk = 1 - 2*((iHel/4)%2);
344  hA = 1 - 2*((iHel/8)%2);
345  hB = 1 - 2*((iHel/16)%2);
346  helString = num2str(hA,2) + num2str(hB,2) + " ->"
347  + num2str(hi,2) + num2str(hj,2) + num2str(hk,2);
348  }
349  else {hA = 9; hB = 9; hi = 9; hj = 9; hk = 9;}
350  vector<int> helBef, helNew;
351  helBef.push_back(hA);
352  helBef.push_back(hB);
353  helNew.push_back(hi);
354  helNew.push_back(hj);
355  helNew.push_back(hk);
356 
357  // Check collinear singularities against massless AP kernels.
358  for (int iTest = 0; iTest <= 4; ++iTest) {
359  // Test invariants, for a few points along collinear axis.
360  double y1 = 1.0e-5;
361  double y2 = 0.1 + iTest*0.2;
362  vector<double> invariants1;
363 
364  // Test collinear limit on side 1.
365  invariants1.resize(3);
366  invariants1[0] = 1.;
367  invariants1[1] = y1;
368  invariants1[2] = y2;
369 
370  // Test collinear limit on side 2 (swap invariants).
371  vector<double> invariants2 = invariants1;
372  invariants2[1] = y2;
373  invariants2[2] = y1;
374 
375  // Compute DGLAP kernels.
376  double AP1 = AltarelliParisi(invariants1, mDum, helBef, helNew);
377  double AP2 = AltarelliParisi(invariants2, mDum, helBef, helNew);
378 
379  // For g->qq antennae, there is only either an ij or jk collinear term.
380  if (abs(id1()) <= 9) {
381  // Keep whichever is larger (more singular), kill the other
382  if (AP1 > AP2) AP2 = -1;
383  else AP1 = -1;
384  }
385  // There is only an IJ singularity if hB is conserved.
386  if (hB != hk) AP1 = -1.;
387  if (hA != hi) AP2 = -1.;
388  // Gluon emission: no external helicity flips for global antennae.
389  if (!sectorShower) {
390  if (id1() == 21) {
391  if (hA != hi || hB != hk) {
392  AP1 = -1.;
393  AP2 = -1.;
394  }
395  }
396  }
397  // Are there any limits to check?
398  if (AP1 < 0. && AP2 < 0.) continue;
399 
400  // Compute antennae (without subleading colour corrections).
401  int modeSLCsave = modeSLC;
402  modeSLC = 1;
403  double ant1 = antFun(invariants1,mDum,helBef,helNew);
404  double ant2 = antFun(invariants2,mDum,helBef,helNew);
405  // For global antennae, add "swapped" terms to represent neighbor.
406  if (!sectorShower && idA() == 21) {
407  invariants1[1] = y1;
408  invariants1[2] = 1. - y1 -y2;
409  vector<int> helComp;
410  helComp.push_back(hj);
411  helComp.push_back(hi);
412  helComp.push_back(hk);
413  ant1 += antFun(invariants1, mDum, helBef, helComp);
414  }
415  if (!sectorShower && idB() == 21) {
416  invariants2[1] = 1 - y1 - y2;
417  invariants2[2] = y1;
418  vector<int> helComp;
419  helComp.push_back(hi);
420  helComp.push_back(hk);
421  helComp.push_back(hj);
422  ant2 += antFun(invariants2, mDum, helBef, helComp);
423  }
424  // Restore subleading-colour level.
425  modeSLC = modeSLCsave;
426 
427  // Compare.
428  if (AP1 > 0.) {
429  double ratio = ant1/AP1;
430  // Require better than 5% agreement unless dominated by nonsingular.
431  if (abs(ratio-1.) >= 0.05 && abs(ant1 - AP1) > 10.) {
432  isOK = false;
433  if (verbose >= normal){
434  printOut(__METHOD_NAME__, "WARNING:" + vinciaName() +
435  " FAILED (collinear ij " + num2str(iTest,1) + " " +
436  helString+" )");
437  }
438  if (verbose >= quiteloud) {
439  cout << setprecision(6);
440  printOut(__METHOD_NAME__, " ant = " + num2str(ant1, 9) +
441  " y1 = " + num2str(y1, 9) +" y2 = " + num2str(y2, 9));
442  printOut(__METHOD_NAME__, " P(z) = "+num2str(AP1, 9) +
443  " zi = " + num2str(zA(invariants1), 9));
444  }
445  } else if (verbose >= verylouddebug) {
446  printOut(__METHOD_NAME__, vinciaName() + " OK (collinear ij " +
447  num2str(iTest, 1) + " " + helString + " )");
448  }
449  }
450 
451  // Require better than 5% agreement unless dominated by nonsingular.
452  if (AP2 > 0.) {
453  double ratio = ant2/AP2;
454  if (abs(ratio - 1.) >= 0.05 && abs(ant2 - AP2) > 10.) {
455  isOK = false;
456  if (verbose >= quiet) {
457  printOut(__METHOD_NAME__, "WARNING:" + vinciaName() +
458  " FAILED (collinear jk " + num2str(iTest, 1) + " " +
459  helString + " )");
460  }
461  if (verbose >= quiteloud) {
462  cout << setprecision(6);
463  printOut(__METHOD_NAME__, " ant = " + num2str(ant2, 9) +
464  " y1 = "+num2str(y2, 9) + " y2 = " + num2str(y1, 9));
465  printOut(__METHOD_NAME__, " P(z) = " + num2str(AP2, 9) +
466  " zk = " + num2str(zB(invariants1), 9));
467  }
468  } else if (verbose >= verylouddebug) {
469  printOut(__METHOD_NAME__, vinciaName() + " OK (collinear jk "
470  + num2str(iTest, 1) + " " + helString+" )");
471  }
472  }
473  }
474 
475  // Check positivity.
476  int nPointsCheck = settingsPtr->mode("Vincia:nPointsCheck");
477  bool isPositive = true;
478  bool isZero = true;
479  bool hasDeadZone = false;
480  double deadZoneAvoidance = settingsPtr->parm("Vincia:deadZoneAvoidance");
481  for (int iTest = 0; iTest < nPointsCheck; ++iTest) {
482 
483  // Generate a random point in phase-space triangle.
484  double y1 = rndmPtr->flat();
485  double y2 = rndmPtr->flat();
486  if (y1 + y2 > 1.0) {
487  y1 = 1 - y1;
488  y2 = 1 - y2;
489  }
490  vector<double> invariants;
491  invariants.resize(3);
492  invariants[0] = 1.;
493  invariants[1] = y1;
494  invariants[2] = y2;
495  double ant = antFun(invariants,mDum,helBef,helNew);
496 
497  // Check positivity (strict).
498  if (ant < 0.) {
499  isPositive = false;
500  if (verbose >= quiteloud){
501  printOut(__METHOD_NAME__, "ERROR:" + vinciaName() + " ant("
502  + num2str(y1, 9) + "," + num2str(y2, 9) + " ; " + helString
503  + ") = " + num2str(ant) + " < 0");
504  }
505  } else if (ant > 0.) isZero = false;
506  // Check for dead zones away from phase-space boundary.
507  if (1-y1-y2 > 0.05 && y1 > 0.05 && y2 > 0.05
508  && ant < deadZoneAvoidance) hasDeadZone = true;
509 
510  } // End nPointsCheck sum for positivity checks.
511  isOK = isOK && isPositive;
512 
513  // Verbose output.
514  if (!isPositive && verbose >= quiet)
515  printOut(__METHOD_NAME__, "ERROR" + vinciaName() +
516  " (ant < 0 encountered " + helString + " )");
517  else if (isPositive && !isZero && verbose >= verylouddebug)
518  printOut(__METHOD_NAME__, vinciaName() + " OK (is positive "
519  + helString + " )");
520  if (!isZero) {
521  if (hasDeadZone && verbose >= quiet)
522  printOut(__METHOD_NAME__, "WARNING" + vinciaName()
523  + " (dead zone encountered " + helString+" )");
524  else if (!hasDeadZone && verbose >= verylouddebug)
525  printOut(__METHOD_NAME__, vinciaName() + " OK (no dead zones "
526  + helString + " )");
527  }
528  } // End helicity sum.
529  return isOK;
530 
531 }
532 
533 //--------------------------------------------------------------------------
534 
535 // Initialize pointers. ParticleData: required for masses, colour
536 // charges, etc. Rndm: required to generate random test points to
537 // check antenna. Settings: required for all coefficients, switches,
538 // etc.
539 
540 void AntennaFunction::initPtr(Info* infoPtrIn, DGLAP* dglapPtrIn) {
541 
542  infoPtr = infoPtrIn;
543  particleDataPtr = infoPtr->particleDataPtr;
544  settingsPtr = infoPtr->settingsPtr;
545  rndmPtr = infoPtr->rndmPtr;
546  dglapPtr = dglapPtrIn;
547  isInitPtr = true;
548  isInit = false;
549 
550 }
551 
552 //--------------------------------------------------------------------------
553 
554 // Auxiliary function to translate an ID code to a string, for
555 // constructing antenna names.
556 
557 string AntennaFunction::id2str(int id) const {
558 
559  // Vector bosons.
560  if (id == 21) return "g";
561  else if (id == 22) return "gamma";
562  else if (id == 23) return "Z";
563  else if (abs(id) == 24) return "W";
564 
565  // Fermions.
566  else if (id >= 1 && id <= 4) return "q";
567  else if (-id >= 1 && -id <= 4) return "qbar";
568  else if (id == 5) return "b";
569  else if (id == -5) return "bbar";
570  else if (id == 6) return "t";
571  else if (id == -6) return "tbar";
572  else if (id >= 11 && id <= 20 && id%2 == 1) return "l-";
573  else if (id >= 11 && id <= 20 && id%2 == 0) return "nu";
574  else if (-id >= 11 && -id <= 20 && id%2 == 1) return "l+";
575  else if (-id >= 11 && -id <= 20 && id%2 == 0) return "nubar";
576 
577  // Octet fermion: use gluino.
578  else if (id == 1000021) return "~g";
579  // Charged scalar: use H+.
580  else if (id == 37) return "H+";
581  else if (id == -37) return "H-";
582  // Coloured scalars: use squarks.
583  else if (id >= 1000000 && id <= 1000010) return "~q";
584  else if (-id >= 1000000 && -id <= 1000010) return "~q*";
585  // Unknown particle.
586  else return "X";
587 
588 }
589 
590 //==========================================================================
591 
592 // Class QQEmitFF, final-final antenna function.
593 
594 //--------------------------------------------------------------------------
595 
596 // The antenna function [GeV^-2].
597 
598 double QQEmitFF::antFun(vector<double> invariants, vector<double> masses,
599  vector<int> helBef, vector<int> helNew) {
600 
601  // Make sure we have enough invariants.
602  if (invariants.size() <= 2) return 0.;
603  double sIK = invariants[0];
604  double sij = invariants[1];
605  double sjk = invariants[2];
606 
607  // Initialise masses and helicities. Return 0 for unphysical helicities.
608  initMasses(&masses);
609  int nAvg = initHel(&helBef, &helNew);
610  if (nAvg <= 0) return 0.;
611 
612  // Check helicity conservation for massless partons.
613  if (mi <= 0. && hi == -hA) return 0.0;
614  else if (mk <= 0. && hk == -hB) return 0.0;
615 
616  // Dimensionless parameters: yij = 2pi.pj / sIK.
617  double yij = sij/sIK;
618  double yjk = sjk/sIK;
619 
620  // Shorthands.
621  double eik = 1./yij/yjk;
622  double mTermI = (mi > 0.) ? pow2(mi)/sij/yij : 0.0;
623  double mTermK = (mk > 0.) ? pow2(mk)/sjk/yjk : 0.0;
624 
625  // Do helicity sum.
626  double hSum(0.);
627 
628  // (++ and --) parents.
629  if (hA * hB > 0 || hA == 9 || hB == 9) {
630 
631  // ++ > +++ && -- > --- antennae (MHV).
632  term = eik - mTermI/(1. - yjk) - mTermK/(1. - yij);
633  if (RH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
634  if (LH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
635 
636  // ++ > +-+ && -- > -+- antennae (NMHV).
637  term = eik * pow2(1. - yij - yjk) - mTermI*(1. - yjk) - mTermK*(1. - yij);
638  if (RH[hA] && RH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += term;
639  if (LH[hA] && LH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += term;
640 
641  // ++ > -++ && -- > +-- antennae (mI != 0).
642  if (mi != 0) {
643  term = mTermI * pow2(yjk)/(1. - yjk);
644  if (RH[hA] && RH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
645  if (LH[hA] && LH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
646  }
647 
648  // ++ > ++- && -- > --+ antennae (mK != 0).
649  if (mk != 0) {
650  term = mTermK * pow2(yij)/(1. - yij);
651  if (RH[hA] && RH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
652  if (LH[hA] && LH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
653  }
654  }
655 
656  // (+- and -+) parents.
657  if (hA * hB < 0 || hA == 9 || hB == 9) {
658 
659  // +- > ++- && -+ > --+ antennae.
660  term = eik * pow2(1. - yij) - mTermI/(1. - yjk) - mTermK*(1. - yij);
661  if (RH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
662  if (LH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
663 
664  // +- > +-- && -+ > -++ antennae.
665  term = eik * pow2(1. - yjk) - mTermI*(1. - yjk) - mTermK/(1. - yij);
666  if (RH[hA] && LH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
667  if (LH[hA] && RH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
668 
669  // +- > -+- && -+ > +-+ antennae (mI != 0).
670  if (mi != 0) {
671  term = mTermI * pow2(yjk)/(1. - yjk);
672  if (RH[hA] && LH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += term;
673  if (LH[hA] && RH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += term;
674  }
675 
676  // +- > +-+ && -+ > -+- antennae (mK != 0).
677  if (mk != 0) {
678  term = mTermK * pow2(yij)/(1. - yij);
679  if (RH[hA] && LH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += term;
680  if (LH[hA] && RH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += term;
681  }
682  }
683 
684  // Return helicity sum, averaged over initial helicities.
685  return hSum/nAvg/sIK;
686 
687 }
688 
689 //--------------------------------------------------------------------------
690 
691 // Function to give Altarelli-Parisi limits of this antenna.
692 // Defined as PI/sij + PK/sjk, i.e. equivalent to antennae.
693 
694 double QQEmitFF::AltarelliParisi(vector<double> invariants,
695  vector<double>, vector<int> helBef, vector<int> helNew) {
696 
697  int h0Now = helNew[0];
698  int h1Now = helNew[1];
699  int h2Now = helNew[2];
700  int hANow = helBef[0];
701  int hBNow = helBef[1];
702 
703  // Compute (sum of) DGLAP kernel(s)/Q^2
704  if (hANow != h0Now || hBNow != h2Now) return -1.;
705  else return dglapPtr->Pq2qg(zA(invariants),hANow,h0Now,h1Now)/invariants[1]
706  + dglapPtr->Pq2qg(zB(invariants),hBNow,h2Now,h1Now)/invariants[2];
707 
708 }
709 
710 //==========================================================================
711 
712 // Class QGEmitFF, final-final antenna function.
713 
714 //--------------------------------------------------------------------------
715 
716 // The antenna function [GeV^-2].
717 
718 double QGEmitFF::antFun(vector<double> invariants, vector<double> masses,
719  vector<int> helBef, vector<int> helNew) {
720 
721  // Make sure we have enough invariants.
722  if (invariants.size() <= 2) return 0.;
723  double sIK = invariants[0];
724  double sij = invariants[1];
725  double sjk = invariants[2];
726 
727  // Initialise masses and helicities. Return 0 for unphysical helicities.
728  initMasses(&masses);
729  int nAvg = initHel(&helBef, &helNew);
730  if (nAvg <= 0) return 0.;
731 
732  // Check helicity conservation for massless partons.
733  if (mi <= 0. && hi == -hA) return 0.0;
734  else if (hk == -hB) return 0.0;
735 
736  // Dimensionless parameters: yij = 2pi.pj / sIK.
737  double yij = sij/sIK;
738  double yjk = sjk/sIK;
739 
740  // Shorthands.
741  double yik = max(0., 1. - yij - yjk);
742  double eik = 1./yij/yjk;
743  double mTerm = pow2(mi)/sij/yij;
744  double a = 1. - alphaSav;
745  if (alphaSav == 0.0) a = 1.;
746  else if (alphaSav == 1.0) a = 0.;
747 
748  // Do helicity sum.
749  double hSum(0.);
750 
751  // (++ and --) parents.
752  if (hA * hB > 0 || hA == 9 || hB == 9) {
753 
754  // ++ > +++ && -- > --- antennae (MHV).
755  term = eik - mTerm/(1. - yjk);
756  if (a != 0.) term += a * (1. - yjk) * ( 1. - 2*yij - yjk )/yjk;
757  if (RH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
758  if (LH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
759 
760  // ++ > +-+ && -- > -+- antennae (NMHV).
761  term = eik * pow2(yik) * (1. - yij) - mTerm*(1. - yjk);
762  if (RH[hA] && RH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += term;
763  if (LH[hA] && LH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += term;
764 
765  // ++ > -++ && -- > +-- antennae (mI != 0).
766  if (mi != 0) {
767  term = mTerm*pow2(yjk)/(1. - yjk);
768  if (RH[hA] && RH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
769  if (LH[hA] && LH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
770  }
771  }
772 
773  // (+- and -+) parents
774  if (hA * hB < 0 || hA == 9 || hB == 9) {
775 
776  // +- > ++- && -+ > --+ antennae.
777  term = eik * pow3(1. - yij) - mTerm/(1. - yjk);
778  if (RH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
779  if (LH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
780 
781  // +- > +-- && -+ > -++ antennae.
782  term = eik * pow2(1. - yjk) - mTerm*(1. - yjk);
783  if (a != 0.) term += a * (1. - yjk) * ( 1. - 2*yij - yjk )/yjk;
784  if (RH[hA] && LH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
785  if (LH[hA] && RH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
786 
787  // +- > -+- && -+ > +-+ antennae (mI != 0).
788  if (mi != 0.) {
789  term = mTerm * pow2(yjk)/(1. - yjk);
790  if (RH[hA] && LH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += term;
791  if (LH[hA] && RH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += term;
792  }
793  }
794 
795  // Subleading colour correction.
796  if (modeSLC >= 2) hSum *= CF/chargeFacSav * (1 - yij)/(2 - yij - yjk)
797  + CA/chargeFacSav * (1 - yjk)/(2 - yij - yjk);
798 
799  // Return helicity sum, averaged over initial helicities.
800  return hSum/nAvg/sIK;
801 
802 }
803 
804 //--------------------------------------------------------------------------
805 
806 // Function to give Altarelli-Parisi limits of this antenna.
807 
808 double QGEmitFF::AltarelliParisi(vector<double> invariants,
809  vector<double>, vector<int> helBef, vector<int> helNew) {
810 
811  int h0Now = helNew[0];
812  int h1Now = helNew[1];
813  int h2Now = helNew[2];
814  int hANow = helBef[0];
815  int hBNow = helBef[1];
816  if (h0Now != hANow) return -1;
817 
818  // Compute (sum of) DGLAP kernel(s)/Q^2
819  double sum(0.);
820  if (h2Now == hBNow)
821  sum += dglapPtr->Pq2qg(zA(invariants), hANow, h0Now, h1Now)/invariants[1];
822  return sum + dglapPtr->Pg2gg(zB(invariants), hBNow, h2Now, h1Now)
823  /invariants[2];
824 
825 }
826 
827 //==========================================================================
828 
829 // Class GQEmitFF, final-final antenna function.
830 
831 //--------------------------------------------------------------------------
832 
833 // The antenna function [GeV^-2] (derived from QGEmit by swapping).
834 
835 double GQEmitFF::antFun(vector<double> invariants,
836  vector<double> mNew, vector<int> helBef, vector<int> helNew) {
837 
838  swap(invariants[1], invariants[2]);
839  swap(mNew[0], mNew[2]);
840  swap(helBef[0], helBef[1]);
841  swap(helNew[0], helNew[2]);
842  return QGEmitFF::antFun(invariants, mNew, helBef, helNew);
843 
844 }
845 
846 //--------------------------------------------------------------------------
847 
848 // Function to give Altarelli-Parisi limits of this antenna.
849 
850 double GQEmitFF::AltarelliParisi(vector<double> invariants,
851  vector<double>, vector<int> helBef, vector<int> helNew) {
852 
853  int h0Now = helNew[0];
854  int h1Now = helNew[1];
855  int h2Now = helNew[2];
856  int hANow = helBef[0];
857  int hBNow = helBef[1];
858  if (h2Now != hBNow) return -1;
859 
860  // Compute (sum of) DGLAP kernel(s)/Q^2
861  double sum(0.);
862  if (h0Now == hANow)
863  sum += dglapPtr->Pq2qg(zB(invariants), hBNow, h2Now, h1Now)/invariants[2];
864  return sum + dglapPtr->Pg2gg(zA(invariants), hANow, h0Now, h1Now)
865  /invariants[1];
866 
867 }
868 
869 //==========================================================================
870 
871 // Class GGEmitFF, final-final antenna function.
872 
873 //--------------------------------------------------------------------------
874 
875 // The antenna function [GeV^-2].
876 
877 double GGEmitFF::antFun(vector<double> invariants, vector<double>,
878  vector<int> helBef, vector<int> helNew) {
879 
880  // Make sure we have enough invariants.
881  if (invariants.size() <= 2) return 0.;
882  double sIK = invariants[0];
883  double sij = invariants[1];
884  double sjk = invariants[2];
885 
886  // Dimensionless parameters: yij = 2pi.pj / sIK.
887  double yij = sij/sIK;
888  double yjk = sjk/sIK;
889 
890  // Initialise helicities. Return 0 for unphysical helicities.
891  int nAvg = initHel(&helBef, &helNew);
892  if (nAvg <= 0) return 0.;
893 
894  // Check physical helicity combination.
895  if (hi == -hA) return 0.0;
896  else if (hk == -hB) return 0.0;
897 
898  // Shorthands.
899  double eik = 1./yij/yjk;
900  double yik = max(0.,1. - yij - yjk);
901  double a;
902  if (alphaSav == 0.0) a = 1.;
903  else if (alphaSav == 1.0) a = 0.;
904  else a = 1. - alphaSav;
905 
906  // Do helicity sum.
907  double hSum(0.);
908 
909  // (++ and --) parents.
910  if (hA * hB > 0 || hA == 9 || hB == 9) {
911 
912  // ++ > +++ && -- > --- antennae (MHV).
913  term = eik;
914  if (a != 0.0) term += a * ((1. - yjk)*(1. - 2*yij - yjk)/yjk
915  + (1. - yij)*(1. - 2*yjk - yij)/yij);
916  if (RH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
917  if (LH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
918 
919  // ++ > +-+ && -- > -+- antennae (NMHV).
920  term = eik * pow3(yik);
921  if (RH[hA] && RH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += term;
922  if (LH[hA] && LH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += term;
923 
924  }
925 
926  // (+- and -+) parents.
927  if (hA * hB < 0 || hA == 9 || hB == 9) {
928 
929  // +- > ++- && -+ > --+ antennae.
930  term = eik * pow3(1 - yij);
931  if (a != 0.) term += a * (1. - yij)*(1. - 2*yjk)/yij;
932  if (RH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
933  if (LH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
934 
935  // +- > +-- && -+ > -++ antennae.
936  term = eik * pow3(1 - yjk);
937  if (a != 0.) term += a * (1. - yjk)*(1. - 2*yij)/yjk;
938  if (RH[hA] && LH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
939  if (LH[hA] && RH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
940  }
941 
942  // Return helicity sum, averaged over initial helicities.
943  return hSum/nAvg/sIK;
944 
945 }
946 
947 //--------------------------------------------------------------------------
948 
949 // Function to give Altarelli-Parisi limits of this antenna.
950 
951 double GGEmitFF::AltarelliParisi(vector<double> invariants,
952  vector<double>, vector<int> helBef, vector<int> helNew) {
953 
954  int h0Now = helNew[0];
955  int h1Now = helNew[1];
956  int h2Now = helNew[2];
957  int hANow = helBef[0];
958  int hBNow = helBef[1];
959  double sum(0.);
960 
961  // Compute (sum of) DGLAP kernel(s)/Q^2
962  if (hBNow == h2Now)
963  sum += dglapPtr->Pg2gg(zA(invariants), hANow, h0Now, h1Now)/invariants[1];
964  if (hANow == h0Now)
965  sum += dglapPtr->Pg2gg(zB(invariants), hBNow, h2Now, h1Now)/invariants[2];
966  return sum;
967 
968 }
969 
970 //==========================================================================
971 
972 // Class GXSplitFF, final-final antenna function.
973 
974 //--------------------------------------------------------------------------
975 
976 // The antenna function [GeV^-2].
977 
978 double GXSplitFF::antFun(vector<double> invariants, vector<double> masses,
979  vector<int> helBef, vector<int> helNew) {
980 
981  // Make sure we have enough invariants.
982  if (invariants.size() <= 2) return 0.;
983  double sIK = invariants[0];
984  double sij = invariants[1];
985  double sjk = invariants[2];
986 
987  // Initialise masses and helicities. Return 0 for unphysical helicities.
988  initMasses(&masses);
989  int nAvg = initHel(&helBef,&helNew);
990  if (nAvg <= 0) return 0.;
991 
992  // Dimensionless parameters: yij = 2pi.pj / sIK.
993  double yij = sij/sIK;
994  double yjk = sjk/sIK;
995  double yik = 1. - yij - yjk - pow2(mi)/sIK - pow2(mj)/sIK;
996 
997  // Check phase space ok.
998  if (yij <= 0.0 || yjk <= 0.0 || yik <= 0.0) return 0.0;
999 
1000  // Shorthands.
1001  double mu2q = mj*mi/sIK;
1002  double mu2ij = yij + 2 * mu2q;
1003  double termQ = 0.5 * (pow2(yik) - mu2q/mu2ij * yik/(1. - yik)) / mu2ij;
1004  double termA = 0.5 * (pow2(yjk) - mu2q/mu2ij * yjk/(1. - yjk)) / mu2ij;
1005  double termF = (mu2q <= 0.) ? 0
1006  : 0.5 * mu2q/pow2(mu2ij)*(yik/(1. - yik) + yjk/(1. - yjk) + 2.);
1007 
1008  // Do helicity sum.
1009  double hSum(0.);
1010 
1011  // (++ and --) parents.
1012  if (hA * hB > 0 || hA == 9 || hB == 9) {
1013 
1014  // Quark (i) takes gluon helicity: ++ > +-+ && -- > -+- antennae.
1015  if (RH[hA] && RH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += termQ;
1016  if (LH[hA] && LH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += termQ;
1017 
1018  // Antiquark (j) takes gluon helicity: ++ > -++ && -- > +-- antennae.
1019  if (RH[hA] && RH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += termA;
1020  if (LH[hA] && LH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += termA;
1021 
1022  // Splitting with helicity flip (mQ != 0): ++ > +++ && -- > ---.
1023  if (mu2q > 0.) {
1024  if (RH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += termF;
1025  if (LH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += termF;
1026  }
1027  }
1028 
1029  // (+- and -+) parents.
1030  if ( hA * hB < 0 || hA == 9 || hB == 9) {
1031 
1032  // Quark (i) takes gluon helicity: +- > +-- && -+ > -++ antennae.
1033  if (RH[hA] && LH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += termQ;
1034  if (LH[hA] && RH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += termQ;
1035 
1036  // Antiquark (j) takes gluon helicity: +- > -+- && -+ > +-+ antennae.
1037  if (RH[hA] && LH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += termA;
1038  if (LH[hA] && RH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += termA;
1039 
1040  // Splitting with helicity flip (mQ != 0): +- > ++- && -+ > --+.
1041  if (mu2q > 0.) {
1042  if (RH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += termF;
1043  if (LH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += termF;
1044  }
1045  }
1046 
1047  // Return helicity sum, averaged over initial helicities.
1048  return hSum/nAvg/sIK;
1049 
1050 }
1051 
1052 //--------------------------------------------------------------------------
1053 
1054 // Function to give Altarelli-Parisi limits of this antenna.
1055 
1056 double GXSplitFF::AltarelliParisi(vector<double> invariants,
1057  vector<double>, vector<int> helBef, vector<int> helNew) {
1058 
1059  int h0Now = helNew[0];
1060  int h1Now = helNew[1];
1061  int h2Now = helNew[2];
1062  int hANow = helBef[0];
1063  int hBNow = helBef[1];
1064 
1065  // Compute DGLAP kernel/Q^2
1066  return hBNow != h2Now ? 0.0 : dglapPtr->Pg2qq(zA(invariants), hANow, h0Now,
1067  h1Now)/invariants[1];
1068 
1069 }
1070 
1071 //==========================================================================
1072 
1073 // Class QGEmitFFsec, sector final-final antenna function, explicit
1074 // symmetrisation of QGEmitFF.
1075 
1076 //--------------------------------------------------------------------------
1077 
1078 // The antenna function [GeV^-2].
1079 
1080 double QGEmitFFsec::antFun(vector<double> invariants,
1081  vector<double> mNew, vector<int> helBef, vector<int> helNew) {
1082 
1083  // Check if helicity vectors empty.
1084  double ant = QGEmitFF::antFun(invariants, mNew, helBef, helNew);
1085  if (helBef.size() < 2) {helBef.push_back(9); helBef.push_back(9);}
1086  if (helNew.size() < 3) {
1087  helNew.push_back(9); helNew.push_back(9); helNew.push_back(9);}
1088 
1089  // Check if j has same helicity as parent gluon.
1090  int hG = helBef[1];
1091  int hjNow = helNew[1];
1092  if ( hG == hjNow ) {
1093  // Define j<->k symmetrisation term.
1094  vector<double> invariantsSym = invariants;
1095  double s02 = invariants[0] - invariants[1] - invariants[2];
1096  vector<int> helSym = helNew;
1097  invariantsSym[1] = s02 + sectorDampSav * invariants[2];
1098  helSym[1] = helNew[2];
1099  helSym[2] = helNew[1];
1100  ant += QGEmitFF::antFun(invariantsSym, mNew, helBef, helSym);
1101  }
1102  return ant;
1103 }
1104 
1105 //==========================================================================
1106 
1107 // Class GQEmitFFsec, sector final-final antenna function, explicit
1108 // symmetrisation of GQEmitFF.
1109 
1110 //--------------------------------------------------------------------------
1111 
1112 // The antenna function [GeV^-2] (derived from QGEmitFFsec by swapping).
1113 
1114 double GQEmitFFsec::antFun(vector<double> invariants,
1115  vector<double> mNew, vector<int> helBef, vector<int> helNew) {
1116 
1117  swap(invariants[1], invariants[2]);
1118  swap(mNew[0], mNew[2]);
1119  swap(helBef[0], helBef[1]);
1120  swap(helNew[0], helNew[2]);
1121  return QGEmitFFsec::antFun(invariants, mNew, helBef, helNew);
1122 
1123 }
1124 
1125 //--------------------------------------------------------------------------
1126 
1127 // Function to give Altarelli-Parisi limits of this antenna.
1128 
1129 double GQEmitFFsec::AltarelliParisi(vector<double> invariants,
1130  vector<double>, vector<int> helBef, vector<int> helNew) {
1131 
1132  int h0Now = helNew[0];
1133  int h1Now = helNew[1];
1134  int h2Now = helNew[2];
1135  int hANow = helBef[0];
1136  int hBNow = helBef[1];
1137  if (h2Now != hBNow) return -1;
1138 
1139  // Compute (sum of) DGLAP kernel(s)/Q^2.
1140  double sum(0.);
1141  if (h0Now == hANow)
1142  sum += dglapPtr->Pq2qg(zB(invariants), hBNow, h2Now, h1Now)/invariants[2];
1143  return sum + dglapPtr->Pg2gg(zA(invariants), hANow, h0Now, h1Now)
1144  /invariants[1];
1145 
1146 }
1147 
1148 //==========================================================================
1149 
1150 // Class GGEmitFFsec, sector final-final antenna function, explicit
1151 // symmetrisation of GGEmitFF.
1152 
1153 //--------------------------------------------------------------------------
1154 
1155 // The antenna function [GeV^-2].
1156 
1157 double GGEmitFFsec::antFun(vector<double> invariants, vector<double> mNew,
1158  vector<int> helBef, vector<int> helNew) {
1159 
1160  // Check if helicity vectors empty
1161  double ant = GGEmitFF::antFun(invariants, mNew, helBef, helNew);
1162  if (helBef.size() < 2) {helBef.push_back(9); helBef.push_back(9);}
1163  if (helNew.size() < 3) {
1164  helNew.push_back(9); helNew.push_back(9); helNew.push_back(9);}
1165 
1166  // Check if j has same helicity as parent gluon 0.
1167  int hjNow = helNew[1];
1168  if (helBef[0] == hjNow) {
1169  // Define i<->j symmetrisation term.
1170  vector<double> invariantsSym = invariants;
1171  double s02 = invariants[0] - invariants[1] - invariants[2];
1172  vector<int> helSym = helNew;
1173  helSym[0] = helNew[1];
1174  helSym[1] = helNew[0];
1175  invariantsSym[2] = s02 + sectorDampSav * invariants[1];
1176  ant += GGEmitFF::antFun(invariantsSym, mNew, helBef, helSym);
1177  }
1178 
1179  // Check if j has same helicity as parent gluon 1.
1180  if (helBef[1] == hjNow) {
1181  // Define j<->k symmetrisation term.
1182  vector<double> invariantsSym = invariants;
1183  double s02 = invariants[0] - invariants[1] - invariants[2];
1184  vector<int> helSym = helNew;
1185  helSym[1] = helNew[2];
1186  helSym[2] = helNew[1];
1187  invariantsSym[1] = s02 + sectorDampSav * invariants[2];
1188  ant += GGEmitFF::antFun(invariantsSym, mNew, helBef, helSym);
1189  }
1190  return ant;
1191 
1192 }
1193 
1194 //==========================================================================
1195 
1196 // Class GXSplitFFsec, sector final-final antenna function, explicit
1197 // symmetrisation of GXSplitFF.
1198 
1199 //--------------------------------------------------------------------------
1200 
1201 // The antenna function [GeV^-2] (just 2*global).
1202 
1203 double GXSplitFFsec::antFun(vector<double> invariants, vector<double> mNew,
1204  vector<int> helBef, vector<int> helNew) {
1205  return 2*GXSplitFF::antFun(invariants,mNew,helBef,helNew);}
1206 
1207 //==========================================================================
1208 
1209 // AntennaFunctionIX class.
1210 
1211 //--------------------------------------------------------------------------
1212 
1213 // Method to initialise (can be different than that of the base class).
1214 
1215 bool AntennaFunctionIX::init() {
1216 
1217  // Check whether pointers are initialized.
1218  if (!isInitPtr) return false;
1219  verbose = settingsPtr->mode("Vincia:verbose");
1220  chargeFacSav = settingsPtr->parm(vinciaName()+":chargeFactor");
1221  if (chargeFacSav < 0.) chargeFacSav = 0.0;
1222 
1223  // Subleading-colour treatment.
1224  // modeSLC = 0: all gluon-emission antennae normalised to CA.
1225  // modeSLC = 1: use colour factors as specified by user.
1226  // modeSLC = 2: QQ gets CF, GG gets CA, QG gets interpolation.
1227  modeSLC = settingsPtr->mode("Vincia:modeSLC");
1228  if (modeSLC == 0 && id1() == 21) chargeFacSav = CA;
1229  if (modeSLC == 2 && id1() == 21) {
1230  if (idA() == 21 && idB() == 21 ) chargeFacSav = CA;
1231  else if (idA() == 21 || idB() == 21) chargeFacSav = (CA + CF)/2.;
1232  else chargeFacSav = CF;
1233  }
1234 
1235  // Collinear partitioning (for global antennae).
1236  alphaSav = settingsPtr->parm("Vincia:octetPartitioning");
1237 
1238  // Sector shower on/off and sectorDamp parameter.
1239  sectorShower = settingsPtr->flag("Vincia:sectorShower");
1240  sectorDampSav = settingsPtr->parm("Vincia:sectorDamp");
1241 
1242  // Return OK.
1243  isInit = true;
1244  return isInit;
1245 
1246 }
1247 
1248 //--------------------------------------------------------------------------
1249 
1250 // Function to check singularities, positivity, etc. For
1251 // initial-initial: sab = sAB + saj + sjb.
1252 
1253 bool AntennaFunctionIX::check() {
1254 
1255  // For check for positivity.
1256  bool isOK = true;
1257  int nPointsCheck = settingsPtr->mode("Vincia:nPointsCheck");
1258  double shhMax = pow2(14000.0);
1259  double deadZoneAvoidance = settingsPtr->parm("Vincia:deadZoneAvoidance");
1260 
1261  // Z and Higgs as testing case.
1262  double sAB[2] = {pow2(91.0), pow2(125.0)};
1263 
1264  // Check soft singularity against dipole factor.
1265  if (id1() == 21 || id1() == 22) {
1266  // Test one symmetric and one asymmetric soft phase-space point.
1267  for (int iTest = 0; iTest <= 1; ++iTest) {
1268  // Test two sAB values.
1269  for (int isAB = 0; isAB <= 1; ++isAB) {
1270  // Test invariants.
1271  double yaj = 1.0e-3 * pow(10.0,iTest); // saj/sab
1272  double yjb = 1.0e-3 / pow(10.0,iTest); // sjb/sab
1273  double sab = sAB[isAB]/(1.0-yaj-yjb); // 1 = sAB/sab + yaj + yjb.
1274  double saj = yaj*sab;
1275  double sjb = yjb*sab;
1276 
1277  // Compute eikonal.
1278  double eik = 2.0*sab/saj/sjb;
1279 
1280  // Compute antenna (without subleading-colour corrections).
1281  vector<double> invariants;
1282  invariants.push_back(sAB[isAB]);
1283  invariants.push_back(saj);
1284  invariants.push_back(sjb);
1285  int modeSLCsave = modeSLC;
1286  modeSLC = 1;
1287  double ant = antFun(invariants);
1288  modeSLC = modeSLCsave;
1289 
1290  // Compare.
1291  double ratio = ant/eik;
1292  if (abs(ratio - 1.0) >= 0.001) {
1293  isOK = false;
1294  if (verbose >= quiet)
1295  printOut(vinciaName() + ":check","WARNING: FAILED "
1296  "(eikonal " + num2str(iTest, 1) + " and sAB = " +
1297  num2str((int)sqrt(sAB[isAB])) + "^2)");
1298  } else if (verbose >= verylouddebug) {
1299  printOut(vinciaName() + ":check", "OK (eikonal " + num2str(iTest, 1)
1300  + " and sAB = " + num2str((int)sqrt(sAB[isAB])) + "^2)");
1301  }
1302  }
1303  }
1304  }
1305 
1306  // Check for collinearity (no helicity APs so far).
1307  // Test invariants, for a few points along collinear axis.
1308  for (int iTest = 0; iTest <= 2; ++iTest) {
1309  // Test two sAB values.
1310  for (int isAB = 0; isAB <= 1; ++isAB) {
1311  // Test invariants.
1312  double y1 = 1.0e-5; // saj/sab
1313  double y2 = 0.2 + iTest*0.3; // sjb/sab
1314  double sab = sAB[isAB]/(1.0-y1-y2); // 1 = sAB/sab + yaj + yjb
1315  double s1 = y1*sab;
1316  double s2 = y2*sab;
1317  vector<double> invariants1 = {sAB[isAB], s1, s2};
1318  vector<double> invariants2 = {sAB[isAB], s2, s1};
1319 
1320  // Compute AP kernels.
1321  double AP1 = AltarelliParisi(invariants1);
1322  double AP2 = AltarelliParisi(invariants2);
1323 
1324  // Compute antennae (without subleading colour corrections).
1325  int modeSLCsave = modeSLC;
1326  modeSLC = 1;
1327  double ant1 = antFun(invariants1);
1328  double ant2 = antFun(invariants2);
1329  modeSLC = modeSLCsave;
1330 
1331  // Check if only a singularity on one side exists (-> z = -1).
1332  double zs1 = zA(invariants1);
1333  double zs2 = zB(invariants1);
1334  if (zs1 != -1.0) {
1335  double ratio = ant1/AP1;
1336  if (abs(ratio - 1.0) >= 0.01) {
1337  isOK = false;
1338  if (verbose >= quiet)
1339  printOut(vinciaName() + ":check","WARNING: FAILED "
1340  "(collinear aj " + num2str(iTest, 1) + " and sAB = " +
1341  num2str((int)sqrt(sAB[isAB])) + "^2)");
1342  if (verbose >= 3)
1343  cout << setprecision(6) << " ant = " << num2str(ant1, 9)
1344  << " y1 = " << num2str(y1,9) << " y2 = " << num2str(y2, 9)
1345  << " " << endl << " P(z) = " << num2str(AP1, 9) << " z = "
1346  << num2str(zs1, 9) << endl;
1347  } else if (verbose >= 6)
1348  printOut(vinciaName() + ":check", "OK (collinear aj " +
1349  num2str(iTest, 1) + " and sAB = " +
1350  num2str((int)sqrt(sAB[isAB])) + "^2)");
1351  }
1352  if (zs2 != -1.0) {
1353  double ratio = ant2/AP2;
1354  if (abs(ratio - 1.0) >= 0.01) {
1355  isOK = false;
1356  if (verbose >= 1)
1357  printOut(vinciaName() + ":check","WARNING: FAILED "
1358  "(collinear jb " + num2str(iTest, 1) + " and sAB = " +
1359  num2str((int)sqrt(sAB[isAB])) + "^2)");
1360  if (verbose >= 3)
1361  cout << setprecision(6) << " ant = " << num2str(ant1, 9)
1362  << " y1 = " << num2str(y1, 9) << " y2 = " << num2str(y2, 9)
1363  << " " << endl << " P(z) = " << num2str(AP2,9) << " z = "
1364  << num2str(zs2, 9) << endl;
1365  } else if (verbose >= 6)
1366  printOut(vinciaName() + ":check", "OK (collinear jb" +
1367  num2str(iTest, 1) + " and sAB = " +
1368  num2str((int)sqrt(sAB[isAB])) + "^2)");
1369  }
1370  }
1371  }
1372 
1373  // Test all helicity configurations for posivity and dead zones.
1374  string helString = " 9 9 -> 9 9 9";
1375  for (int iHel = -1; iHel < 32; ++iHel) {
1376  bool isPositive = true;
1377  bool isZero = true;
1378  bool hasDeadZone = false;
1379  if (iHel >= 0) {
1380  hj = 1 - 2*(iHel%2);
1381  hi = 1 - 2*((iHel/2)%2);
1382  hk = 1 - 2*((iHel/4)%2);
1383  hA = 1 - 2*((iHel/8)%2);
1384  hB = 1 - 2*((iHel/16)%2);
1385  helString = num2str(hA,2) + num2str(hB,2) + " ->" +
1386  num2str(hi,2) + num2str(hj,2) + num2str(hk,2);
1387  }
1388  vector<int> helBef = {hA, hB};
1389  vector<int> helNew = {hi, hj, hk};
1390 
1391  // Gluon emission.
1392  if (id1() == 21) {
1393  // Massless: quark helicities must be conserved.
1394  if (idA() != 21 && hi != hA) continue;
1395  if (idB() != 21 && hk != hB) continue;
1396  // For gluons, sector terms can give helicity flips
1397  if (idA() == 21 && (hi != hA && hj != -hA)) continue;
1398  if (idB() == 21 && (hk != hB && hj != -hB)) continue;
1399  // Do not allow helicity flips on both sides simultaneously
1400  if (hi != hA && hk != hB) continue;
1401  }
1402  // Gluon conversion.
1403  else if (idA() == 21) {
1404  // Helicity of recoiler must be conserved.
1405  if (hk != hB) continue;
1406  // Massless: quark helicity must be conserved.
1407  if (hi != hj) continue;
1408  }
1409  // Gluon splitting.
1410  else {
1411  // Helicity of recoiler must be conserved.
1412  if (hk != hB) continue;
1413  // Massless: quark helicity must be conserved.
1414  if (hA == hj) continue;
1415  }
1416 
1417  // Check for positivity and dead zones.
1418  for (int iTest = 0; iTest <= nPointsCheck; ++iTest) {
1419  // sAB^0.5 between 1 and 1000 GeV.
1420  double sABnow = pow2(1.0+rndmPtr->flat()*999.0);
1421  // Generate a random point in phase-space triangle.
1422  double yaj = rndmPtr->flat(); // saj/sab
1423  double yjb = rndmPtr->flat(); // sjb/sab
1424  if (yaj + yjb > 1.0) {
1425  yaj = 1.0 - yaj;
1426  yjb = 1.0 - yjb;
1427  }
1428 
1429  // Compute sab from 1 = sAB/sab + yaj + yjb.
1430  double sab = sABnow/(1.0-yaj-yjb);
1431 
1432  // Check that sab < shhMax.
1433  if (sab > shhMax) continue;
1434  double saj = yaj*sab;
1435  double sjb = yjb*sab;
1436  vector<double> invariants = {sABnow, saj, sjb};
1437 
1438  // Compute antenna.
1439  double ant = antFun(invariants, mDum, helBef, helNew);
1440 
1441  // Check positivity (strict).
1442  if (ant <= 0.0) {
1443  isPositive = false;
1444  if (verbose >= 3)
1445  printOut(vinciaName() + ":check", "ERROR ant(" +
1446  num2str((int)sqrt(saj)) + "^2," + num2str((int)sqrt(sjb)) + "^2," +
1447  num2str((int)sqrt(sABnow)) + "^2 ; " + helString + ") < 0");
1448  } else if (ant > 0.0) isZero = false;
1449 
1450  // Check for dead zones away from phase-space boundary
1451  if ((1.0 - yaj - yjb > 0.05) && (yaj > 0.05) && (yjb > 0.05)
1452  && (ant*sABnow < deadZoneAvoidance)) hasDeadZone = true;
1453  }
1454  isOK = isOK && isPositive;
1455 
1456  // Verbose output
1457  if (!isPositive && verbose >= 1) printOut(vinciaName() + ":check",
1458  "ERROR (ant < 0 encountered " + helString + " )");
1459  else if (isPositive && !isZero && verbose >= 6)
1460  printOut(vinciaName() + ":check", "OK (is positive " + helString + " )");
1461  if (!isZero) {
1462  if (hasDeadZone && verbose >= 1) printOut(vinciaName() + ":check",
1463  "WARNING (dead zone encountered " + helString + " )");
1464  else if (!hasDeadZone && verbose >= 6) printOut(vinciaName() + ":check",
1465  "OK (no dead zones " + helString+" )");
1466  }
1467  } // End loop over helicities.
1468  return isOK;
1469 
1470 }
1471 
1472 //==========================================================================
1473 
1474 // Class QQEmitII, initial-initial antenna function.
1475 
1476 //--------------------------------------------------------------------------
1477 
1478 // The antenna function [GeV^-2].
1479 
1480 double QQEmitII::antFun(vector<double> invariants, vector<double> masses,
1481  vector<int> helBef, vector<int> helNew) {
1482 
1483  // Invariants.
1484  double sAB = invariants[0];
1485  double saj = invariants[1];
1486  double sjb = invariants[2];
1487 
1488  // Sanity check. Require positive invariants.
1489  if ((saj <= 0.0) || (sjb <= 0.0) || (sAB <= 0.0)) return 0.0;
1490 
1491  // Initialise masses and helicities. Return 0 for unphysical helicities.
1492  initMasses(&masses);
1493  int nAvg = initHel(&helBef,&helNew);
1494  if (nAvg <= 0) return 0.;
1495 
1496  // Shorthands.
1497  double sab = saj+sjb+sAB;
1498  double yaj = saj/sab;
1499  double yjb = sjb/sab;
1500  double yAB = sAB/sab;
1501  double eikJ = 1./(sAB*yaj*yjb);
1502 
1503  // Mass corrections.
1504  double facMa = (mi != 0.) ? pow2(mi) / sab / pow2(yaj) / sAB : 0.;
1505  double facMb = (mk != 0.) ? pow2(mk) / sab / pow2(yjb) / sAB : 0.;
1506 
1507  // Do helicity sum.
1508  double hSum = 0.0;
1509 
1510  // (++ and --) parents.
1511  if (hA * hB > 0 || hA == 9 || hB == 9) {
1512  // ++ > +++ && -- > --- antennae (MHV).
1513  term = eikJ - facMa - facMb;
1514  if (RH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
1515  if (LH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
1516 
1517  // ++ > +-+ && -- > -+- antennae (NMHV).
1518  term = pow2(yAB)*eikJ - pow2(1-yjb)*facMa - pow2(1-yaj)*facMb;
1519  if (RH[hA] && RH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += term;
1520  if (LH[hA] && LH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += term;
1521 
1522  // ++ > --+ && -- > ++- antennae (massive helicity flip).
1523  if (mi != 0.0) {
1524  term = facMa * pow2(yjb);
1525  if (RH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
1526  if (LH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
1527  }
1528 
1529  // ++ > +-- && -- > -++ antennae (massive helicity flip).
1530  if (mk != 0.0) {
1531  term = facMb * pow2(yaj);
1532  if (RH[hA] && RH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
1533  if (LH[hA] && LH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
1534  }
1535  }
1536 
1537  // (+- and -+) parents.
1538  if (hA * hB < 0 || hA == 9 || hB == 9) {
1539  // +- > ++- && -+ > --+ antennae.
1540  term = pow2(1. - yaj)*eikJ;
1541  if (RH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
1542  if (LH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
1543 
1544  // +- > +-- && -+ > -++ antennae.
1545  term = pow2(1. - yjb)*eikJ;
1546  if (RH[hA] && LH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
1547  if (LH[hA] && RH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
1548 
1549  // +- > --- && -+ > +++ antennae (massive helicity flip)
1550  if (mi != 0.0) {
1551  term = facMa * pow2(yjb);
1552  if (RH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
1553  if (LH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
1554  }
1555 
1556  // +- > -++ && -+ > +-- antennae (massive helicity flip)
1557  if (mk != 0.0) {
1558  term = facMb * pow2(yaj);
1559  if (RH[hA] && LH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
1560  if (LH[hA] && RH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
1561  }
1562  }
1563 
1564  // Return helicity sum, averaged over initial helicities.
1565  return hSum/nAvg;
1566 
1567 }
1568 
1569 //--------------------------------------------------------------------------
1570 
1571 // AP splitting kernel for collinear limit checks, P(z)/Q2.
1572 
1573 double QQEmitII::AltarelliParisi(vector<double> invariants, vector<double>,
1574  vector<int>, vector<int>) {
1575 
1576  // Sanity check. Require positive invariants.
1577  double sAB = invariants[0];
1578  double saj = invariants[1];
1579  double sjb = invariants[2];
1580  if ((saj <= 0.0) || (sjb <= 0.0) || (sAB <= 0.0)) return 0.0;
1581 
1582  // Using smaller invariant for collinear limit.
1583  double z = ( saj < sjb ? zA(invariants) : zB(invariants) );
1584  double Q2 = min(saj,sjb);
1585  double AP = 1.0/z * (1.0 + pow2(z))/(1.0 - z);
1586  // AP uses alpha/2pi * CF, ant uses alpha/4pi * 2CF, with CF = 4/3.
1587  AP *= 1.0;
1588  // ant with quark reproduces full AP.
1589  AP *= 1.0;
1590  return AP/Q2;
1591 
1592 }
1593 
1594 //==========================================================================
1595 
1596 // Class GQEmitII, initial-initial antenna function.
1597 
1598 //--------------------------------------------------------------------------
1599 
1600 // The antenna function [GeV^-2].
1601 
1602 double GQEmitII::antFun(vector<double> invariants, vector<double> masses,
1603  vector<int> helBef, vector<int> helNew) {
1604 
1605  // Invariants and helicities.
1606  double sAB = invariants[0];
1607  double saj = invariants[1];
1608  double sjb = invariants[2];
1609 
1610  // Sanity check. Require positive invariants.
1611  if ((saj <= 0.0) || (sjb <= 0.0) || (sAB <= 0.0)) return 0.0;
1612 
1613  // Initialise masses and helicities. Return 0 for unphysical helicities.
1614  initMasses(&masses);
1615  int nAvg = initHel(&helBef, &helNew);
1616  if (nAvg <= 0) return 0.;
1617 
1618  // Shorthands
1619  double sab = saj+sjb+sAB;
1620  double yaj = saj/sab;
1621  double yjb = sjb/sab;
1622  double yAB = sAB/sab;
1623  double eikJ = 1./(sAB*yaj*yjb);
1624  double collA = 1./(sAB*yaj*(1. - yjb));
1625 
1626  // Mass corrections.
1627  double facMb = (mk != 0.) ? pow2(mk) / sab / pow2(yjb) / sAB : 0.;
1628 
1629  // Do helicity sum.
1630  double hSum = 0.0;
1631 
1632  // (++ and --) parents.
1633  if (hA * hB > 0 || hA == 9 || hB == 9) {
1634  // ++ > +++ && -- > --- antennae (MHV).
1635  term = eikJ + collA - facMb;
1636  if (RH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
1637  if (LH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
1638 
1639  // ++ > +-+ && -- > -+- antennae (NMHV).
1640  term = (1. - yjb)*pow2(yAB)*eikJ - pow2(1-yaj) * facMb;
1641  if (RH[hA] && RH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += term;
1642  if (LH[hA] && LH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += term;
1643 
1644  // ++ > --+ && -- > ++- antennae.
1645  term = pow3(yjb)*collA;
1646  if (RH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
1647  if (LH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
1648 
1649  // ++ > +-- && -- > -++ antennae (massive helicity flip).
1650  if (mk != 0.) {
1651  term = facMb * pow2(yaj);
1652  if (RH[hA] && RH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
1653  if (LH[hA] && LH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
1654  }
1655  }
1656 
1657  // (+- and -+) parents.
1658  if (hA * hB < 0 || hA == 9 || hB == 9) {
1659  // +- > ++- && -+ > --+ antennae.
1660  term = pow2(1. - yaj)*eikJ + collA - pow2(1-yaj)*facMb;
1661  if (RH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
1662  if (LH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
1663 
1664  // +- > +-- && -+ > -++ antennae.
1665  term = pow3(1. - yjb)*eikJ - facMb;
1666  if (RH[hA] && LH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
1667  if (LH[hA] && RH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
1668 
1669  // +- > --- && -+ > +++ antennae.
1670  term = pow3(yjb)*collA;
1671  if (RH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
1672  if (LH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
1673 
1674  // +- > +++ && -+ > --- antennae (massive helicity flip).
1675  if (mk != 0.) {
1676  term = facMb * pow2(yaj);
1677  if (RH[hA] && LH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
1678  if (LH[hA] && RH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
1679  }
1680  }
1681 
1682  // Subleading colour correction.
1683  if (modeSLC >= 2) hSum *= CF/chargeFacSav * (1 - yjb)/(2 - yaj - yjb)
1684  + CA/chargeFacSav * (1 - yaj)/(2 - yaj - yjb);
1685 
1686  // Return helicity sum, averaged over initial helicities.
1687  return hSum/nAvg;
1688 
1689 }
1690 
1691 //--------------------------------------------------------------------------
1692 
1693 // AP splitting kernel for collinear limit checks, P(z)/Q2.
1694 
1695 double GQEmitII::AltarelliParisi(vector<double> invariants,
1696  vector<double>, vector<int>, vector<int>) {
1697 
1698  // Sanity check. Require positive invariants.
1699  double sAB = invariants[0];
1700  double saj = invariants[1];
1701  double sjb = invariants[2];
1702 
1703  // Sanity check. Require positive invariants.
1704  if ((saj <= 0.0) || (sjb <= 0.0) || (sAB <= 0.0)) return 0.0;
1705 
1706  // Using smaller invariant for collinear limit
1707  double z = ( saj < sjb ? zA(invariants) : zB(invariants));
1708  double Q2 = min(saj,sjb);
1709  double AP = 0.0;
1710 
1711  // Collinear to initial gluon.
1712  if (saj < sjb) {
1713  AP = 1.0/z*(pow(z, 4.0) + 1.0 + pow(1 - z, 4.0))/z/(1.0 - z);
1714  // AP uses alpha/2pi * CA, ant uses alpha/4pi * CA, with CA = 3.
1715  AP *= 2.0;
1716  // gluon appears in 2 ants, hence ant reproduces half of AP.
1717  AP *= 0.5;
1718 
1719  // Collinear to initial quark.
1720  } else {
1721  AP = 1.0/z * (1.0+pow2(z))/(1.0-z);
1722  // AP uses alpha/2pi * CF, ant uses alpha/4pi * 2CF, with CF = 4/3.
1723  AP *= 1.0;
1724  // ant with quark reproduces full AP.
1725  AP *= 1.0;
1726  }
1727  return AP/Q2;
1728 
1729 }
1730 
1731 //==========================================================================
1732 
1733 // Class GGEmitII, initial-initial antenna function.
1734 
1735 //--------------------------------------------------------------------------
1736 
1737 // The antenna function [GeV^-2].
1738 
1739 double GGEmitII::antFun(vector<double> invariants, vector<double>,
1740  vector<int> helBef, vector<int> helNew) {
1741 
1742  // Invariants and helicities.
1743  double sAB = invariants[0];
1744  double saj = invariants[1];
1745  double sjb = invariants[2];
1746 
1747  // Sanity check. Require positive invariants.
1748  if ((saj <= 0.0) || (sjb <= 0.0) || (sAB <= 0.0)) return 0.0;
1749 
1750  // Initialise helicities. Return 0 for unphysical helicities.
1751  int nAvg = initHel(&helBef, &helNew);
1752  if (nAvg <= 0) return 0.;
1753 
1754  // Shorthands.
1755  double sab = saj+sjb+sAB;
1756  double yaj = saj/sab;
1757  double yjb = sjb/sab;
1758  double yAB = sAB/sab;
1759  double eikJ = 1./(sAB*yaj*yjb);
1760  double collA = 1./(sAB*yaj*(1.-yjb));
1761  double collB = 1./(sAB*yjb*(1.-yaj));
1762 
1763  // Do helicity sum.
1764  double hSum = 0.0;
1765 
1766  // (++ and --) parents.
1767  if ( hA * hB > 0 || hA == 9 || hB == 9) {
1768  // ++ > +++ && -- > --- antennae (MHV).
1769  term = eikJ + collA + collB;
1770  if (RH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
1771  if (LH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
1772 
1773  // ++ > +-+ && -- > -+- antennae (NMHV).
1774  term = pow3(yAB)*eikJ;
1775  if (RH[hA] && RH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += term;
1776  if (LH[hA] && LH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += term;
1777 
1778  // ++ > --+ && -- > ++- antennae.
1779  term = pow3(yjb)*collA;
1780  if (RH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
1781  if (LH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
1782 
1783  // ++ > +-- && -- > -++ antennae.
1784  term = pow3(yaj)*collB;
1785  if (RH[hA] && RH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
1786  if (LH[hA] && LH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
1787  }
1788 
1789  // (+- and -+) parents.
1790  if (hA * hB < 0 || hA == 9 || hB == 9) {
1791  // +- > ++- && -+ > --+ antennae.
1792  term = pow3(1.-yaj)*eikJ + collA;
1793  if (RH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
1794  if (LH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
1795 
1796  // +- > +-- && -+ > -++ antennae.
1797  term = pow3(1 - yjb)*eikJ + collB;
1798  if (RH[hA] && LH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
1799  if (LH[hA] && RH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
1800 
1801  // +- > +++ && -+ > --- antennae.
1802  term = pow3(yaj)*collB;
1803  if (RH[hA] && LH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
1804  if (LH[hA] && RH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
1805 
1806  // +- > --- && -+ > +++ antennae.
1807  term = pow3(yjb)*collA;
1808  if (RH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
1809  if (LH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
1810  }
1811 
1812  // Return helicity sum, averaged over initial helicities.
1813  return hSum/nAvg;
1814 
1815 }
1816 
1817 //--------------------------------------------------------------------------
1818 
1819 // AP splitting kernel, P(z)/Q2.
1820 
1821 double GGEmitII::AltarelliParisi(vector<double> invariants, vector<double>,
1822  vector<int>, vector<int>) {
1823 
1824  // Sanity check. Require positive invariants.
1825  double sAB = invariants[0];
1826  double saj = invariants[1];
1827  double sjb = invariants[2];
1828 
1829  // Sanity check. Require positive invariants.
1830  if ((saj <= 0.0) || (sjb <= 0.0) || (sAB <= 0.0)) return 0.0;
1831 
1832  // Using smaller invariant for collinear limit.
1833  double z = ( saj < sjb ? zA(invariants) : zB(invariants) );
1834  double Q2 = min(saj,sjb);
1835  double AP = 1.0/z*(pow(z, 4.0) + 1.0+pow(1 - z, 4.0))/z/(1.0 - z);
1836  // AP uses alpha/2pi * CA, ant uses alpha/4pi * CA, with CA = 3.
1837  AP *= 2.0;
1838  // gluon appears in 2 ants, hence ant reproduces half of AP.
1839  AP *= 0.5;
1840  return AP/Q2;
1841 
1842 }
1843 
1844 //==========================================================================
1845 
1846 // Class QXSplitII, initial-initial antenna function.
1847 
1848 //--------------------------------------------------------------------------
1849 
1850 // The antenna function [GeV^-2].
1851 
1852 double QXSplitII::antFun(vector<double> invariants, vector<double> masses,
1853  vector<int> helBef, vector<int> helNew) {
1854 
1855  // Invariants and helicities.
1856  double sAB = invariants[0];
1857  double saj = invariants[1];
1858  double sjb = invariants[2];
1859 
1860  // Sanity check. Require positive invariants.
1861  if ((saj <= 0.0) || (sjb <= 0.0) || (sAB <= 0.0)) return 0.0;
1862 
1863  // Initialise masses and helicities. Return 0 for unphysical helicities.
1864  initMasses(&masses);
1865  int nAvg = initHel(&helBef,&helNew);
1866  if (nAvg <= 0) return 0.;
1867 
1868  // Shorthands.
1869  double sab = sAB + saj + sjb;
1870  double yaj = saj/sab;
1871  double yAB = sAB/sab;
1872  double splitA = 1.0/sAB/yaj;
1873  double facMj = (mj != 0. ) ? pow2(mj) / sab / pow2(yaj) / sAB : 0.0;
1874 
1875  // Do helicity sum.
1876  double hSum = 0.0;
1877 
1878  // (++ and --) parents.
1879  if (hA * hB > 0 || hA == 9 || hB == 9) {
1880  // ++ > +-+ && -- > -+- antennae.
1881  term = pow2(yAB) * splitA - pow2(yAB)/(1. - yAB) * facMj;
1882  if (RH[hA] && RH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += term;
1883  if (LH[hA] && LH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += term;
1884 
1885  // ++ > --+ && -- > ++- antennae.
1886  term = pow2(1. - yAB) * splitA - (1. - yAB) * facMj;
1887  if (RH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
1888  if (LH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
1889 
1890  // ++ > +++ && -- > --- antennae (massive helicity flip).
1891  if (mj != 0.0) {
1892  term = facMj / (1. - yAB);
1893  if (RH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
1894  if (LH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
1895  }
1896  }
1897 
1898  // (+- and -+) parents.
1899  if (hA * hB < 0 || hA == 9 || hB == 9) {
1900  // +- > +-- && -+ > -++ antennae.
1901  term = pow2(yAB) * splitA - pow2(yAB)/(1. - yAB) * facMj;
1902  if (RH[hA] && LH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
1903  if (LH[hA] && RH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
1904 
1905  // +- > --- && -+ > +++ antennae.
1906  term = pow2(1. - yAB) * splitA - (1. - yAB) * facMj;
1907  if (RH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
1908  if (LH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
1909 
1910  // +- > ++- && -+ > --+ antennae.
1911  if (mj != 0.0) {
1912  term = facMj / (1. - yAB);
1913  if (RH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
1914  if (LH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
1915  }
1916  }
1917 
1918  // Return helicity sum, averaged over initial helicities.
1919  return hSum/nAvg;
1920 
1921 }
1922 
1923 //--------------------------------------------------------------------------
1924 
1925 // AP splitting kernel, P(z)/Q2.
1926 
1927 double QXSplitII::AltarelliParisi(vector<double> invariants, vector<double>,
1928  vector<int>, vector<int>) {
1929 
1930  // Sanity check. Require positive invariants.
1931  double sAB = invariants[0];
1932  double saj = invariants[1];
1933  double sjb = invariants[2];
1934 
1935  // Sanity check. Require positive invariants.
1936  if ((saj <= 0.0) || (sjb <= 0.0) || (sAB <= 0.0)) return 0.0;
1937 
1938  // There is only the saj collinear limit.
1939  double z = zA(invariants);
1940  double Q2 = saj;
1941  double AP = 1.0/z * (pow2(z)+pow2(1 - z));
1942  // AP uses alpha/2pi * TR/2, ant uses alpha/4pi * TR, with TR = 1.
1943  AP *= 1.0;
1944  // Ant with quark reproduces full AP.
1945  AP *= 1.0;
1946  return AP/Q2;
1947 
1948 }
1949 
1950 //==========================================================================
1951 
1952 // Class GXConvII, initial-initial antenna function.
1953 
1954 //--------------------------------------------------------------------------
1955 
1956 // The antenna function [GeV^-2].
1957 
1958 double GXConvII::antFun(vector<double> invariants, vector<double> masses,
1959  vector<int> helBef, vector<int> helNew) {
1960 
1961  // Invariants and helicities.
1962  double sAB = invariants[0];
1963  double saj = invariants[1];
1964  double sjb = invariants[2];
1965 
1966  // Sanity check. Require positive invariants.
1967  if ((saj <= 0.0) || (sjb <= 0.0) || (sAB <= 0.0)) return 0.0;
1968 
1969  // Initialise masses and helicities. Return 0 for unphysical helicities.
1970  initMasses(&masses);
1971  int nAvg = initHel(&helBef, &helNew);
1972  if (nAvg <= 0) return 0.;
1973 
1974  // Shorthands.
1975  double sab = sAB + saj + sjb - 2*pow2(mj);
1976  double yaj = saj/sab;
1977  double yAB = sAB/sab;
1978  double mu2j = (mj != 0.) ? pow2(mj)/sab : 0.;
1979  double convA = 1.0/(2.*sAB*(yaj - 2.*mu2j)*yAB);
1980  double facMj = (mj != 0.) ? mu2j/(2.*sAB*pow2(yaj - 2.*mu2j)) : 0.;
1981 
1982  // Do helicity sum.
1983  double hSum = 0.0;
1984 
1985  // (++ and --) parents.
1986  if (hA * hB > 0 || hA == 9 || hB == 9) {
1987  // ++ > +++ && -- > --- antennae.
1988  term = convA - facMj*yAB/(1. - yAB);
1989  if (RH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
1990  if (LH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
1991 
1992  // ++ > --+ && -- > ++- antennae.
1993  term = pow2(1 - yAB)*convA - facMj*yAB*(1. - yAB);
1994  if (RH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
1995  if (LH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
1996 
1997  // ++ > +-+ && -- > -+- antennae (massive helicity flip).
1998  if (mj != 0.) {
1999  term = facMj*pow3(yAB)/(1. - yAB);
2000  if (RH[hA] && RH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += term;
2001  if (LH[hA] && LH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += term;
2002  }
2003  }
2004 
2005  // (+- and -+) parents.
2006  if (hA * hB < 0 || hA == 9 || hB == 9) {
2007  // +- > ++- && -+ > --+ antennae.
2008  term = convA - facMj*yAB/(1. - yAB);
2009  if (RH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
2010  if (LH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
2011 
2012  // +- > --- && -+ > +++ antennae.
2013  term = pow2(1 - yAB)*convA - facMj*yAB*(1. - yAB);
2014  if (RH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
2015  if (LH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
2016 
2017  // +- > +-- && -+ > -++ antennae (massive helicity flip).
2018  term = facMj*pow3(yAB)/(1. - yAB);
2019  if (RH[hA] && LH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
2020  if (LH[hA] && RH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
2021  }
2022 
2023  // Return helicity sum, averaged over initial helicities.
2024  return hSum/nAvg;
2025 
2026 }
2027 
2028 //--------------------------------------------------------------------------
2029 
2030 // AP splitting kernel, P(z)/Q2.
2031 
2032 double GXConvII::AltarelliParisi(vector<double> invariants, vector<double>,
2033  vector<int>, vector<int>) {
2034 
2035  // Sanity check. Require positive invariants.
2036  double sAB = invariants[0];
2037  double saj = invariants[1];
2038  double sjb = invariants[2];
2039 
2040  // Sanity check. Require positive invariants.
2041  if ((saj <= 0.0) || (sjb <= 0.0) || (sAB <= 0.0)) return 0.0;
2042 
2043  // There is only the saj collinear limit.
2044  double z = zA(invariants);
2045  double Q2 = saj;
2046  double AP = 1.0/z * (1.0+pow2(1.0-z))/z;
2047  // AP uses alpha/2pi * CF, ant uses alpha/4pi * 2CF, with CF = 4/3.
2048  AP *= 1.0;
2049  // gluon appears in 2 ants, hence ant reproduces half of AP.
2050  AP *= 0.5;
2051  return AP/Q2;
2052 
2053 }
2054 
2055 //==========================================================================
2056 
2057 // Class AntennaFunctionIF, base class for IF antenna functions which
2058 // implements QQEmitIF.
2059 
2060 //--------------------------------------------------------------------------
2061 
2062 // Method to initialise (can be different than that of the base class).
2063 
2064 bool AntennaFunctionIF::init() {
2065 
2066  // Check whether pointers are initialized.
2067  if (!isInitPtr) return false;
2068  verbose = settingsPtr->mode("Vincia:verbose");
2069  chargeFacSav = settingsPtr->parm(vinciaName()+":chargeFactor");
2070  if (chargeFacSav < 0.) chargeFacSav = 0.0;
2071 
2072  // Subleading-colour treatment.
2073  // modeSLC = 0: all gluon-emission antennae normalised to CA.
2074  // modeSLC = 1: use colour factors as specified by user.
2075  // modeSLC = 2: QQ gets CF, GG gets CA, QG gets interpolation.
2076  modeSLC = settingsPtr->mode("Vincia:modeSLC");
2077  if (modeSLC == 0 && id1() == 21) chargeFacSav = CA;
2078  if (modeSLC == 2 && id1() == 21) {
2079  if (idA() == 21 && idB() == 21 ) chargeFacSav = CA;
2080  else if (idA() == 21 || idB() == 21) chargeFacSav = (CA + CF)/2.;
2081  else chargeFacSav = CF;
2082  }
2083 
2084  // Set kinematics map, depending on whether this is an IF or RF antenna.
2085  if (settingsPtr->isMode(vinciaName() + ":kineMap"))
2086  kineMapSav = settingsPtr->mode(vinciaName()+":kineMap");
2087  else if (isRFant()) {
2088  // Gluon emission antennae.
2089  if (id1() == 21) kineMapSav = settingsPtr->mode("Vincia:kineMapRFemit");
2090  // Gluon splitting antennae.
2091  else kineMapSav = settingsPtr->mode("Vincia:kineMapRFsplit");
2092  }
2093  // IF antennae.
2094  else kineMapSav = settingsPtr->mode("Vincia:kineMapIF");
2095 
2096  // Collinear partitioning (for global antennae).
2097  alphaSav = settingsPtr->parm("Vincia:octetPartitioning");
2098 
2099  // Sector shower on/off and sectorDamp parameter.
2100  sectorShower = settingsPtr->flag("Vincia:sectorShower");
2101  sectorDampSav = settingsPtr->parm("Vincia:sectorDamp");
2102 
2103  // Return OK.
2104  isInit = true;
2105  return isInit;
2106 
2107 }
2108 
2109 //--------------------------------------------------------------------------
2110 
2111 // Function to check singularities, positivity, etc.
2112 
2113 bool AntennaFunctionIF::check() {
2114 
2115  // For check for positivity
2116  bool isOK = true;
2117  int nPointsCheck = settingsPtr->mode("Vincia:nPointsCheck");
2118  double shhMax = pow2(14000.0);
2119  double deadZoneAvoidance = settingsPtr->parm("Vincia:deadZoneAvoidance");
2120 
2121  // Resonance-final.
2122  if (isRFant()) return checkRes();
2123 
2124  // Initial-Final: sAK + sjk = sak + saj.
2125  else {
2126 
2127  // Two values for sAK.
2128  double sAK[2] = {pow2(50.0), pow2(150.0)};
2129 
2130  // Check soft singularity against dipole factor.
2131  if (id1() == 21 || id1() == 22) {
2132  // Test one symmetric and one asymmetric soft phase-space point.
2133  for (int iTest = 0; iTest <= 1; ++iTest) {
2134  // Test two sAK values.
2135  for (int isAK = 0; isAK <= 1; ++isAK) {
2136  // Test invariants.
2137  double yaj = 1.0e-3 * pow(10.0, iTest); // saj/sAK
2138  double yjk = 1.0e-3 / pow(10.0, iTest); // sjk/sAK
2139  double sak = sAK[isAK]*(1.0 + yjk - yaj); // 1 + yjk = sak/sAK + yaj
2140  double saj = yaj*sAK[isAK];
2141  double sjk = yjk*sAK[isAK];
2142  vector<double> invariants ={sAK[isAK], saj, sjk};
2143  // Compute eikonal.
2144  double eik = 2.0*sak/saj/sjk;
2145  // Compute antenna (without subleading-colour corrections).
2146  int modeSLCsave = modeSLC;
2147  modeSLC = 1;
2148  double ant = antFun(invariants);
2149  modeSLC = modeSLCsave;
2150  // Compare.
2151  double ratio = ant/eik;
2152  if (abs(ratio - 1.0) >= 0.001) {
2153  isOK = false;
2154  if (verbose >= 1) printOut(vinciaName() + ":check",
2155  "WARNING: FAILED (eikonal " + num2str(iTest, 1) +
2156  " and sAK = " + num2str((int)sqrt(sAK[isAK])) + "^2)");
2157  } else if (verbose >= 6) printOut(vinciaName() + ":check",
2158  "OK (eikonal " + num2str(iTest, 1) + " and sAK = " +
2159  num2str((int)sqrt(sAK[isAK])) + "^2)");
2160  }
2161  }
2162  }
2163 
2164  // Check for collinearity (no helicity APs so far).
2165  // Test invariants, for a few points along collinear axis.
2166  for (int iTest = 0; iTest <= 2; ++iTest) {
2167  // Test two sAK values.
2168  for (int isAK = 0; isAK <= 1; ++isAK) {
2169  // Test invariants.
2170  double y1 = 1.0e-5; // saj/sAK
2171  double y2 = 0.2 + iTest*0.3; // sjk/sAK
2172  double s1 = y1*sAK[isAK];
2173  double s2 = y2*sAK[isAK];
2174  vector<double> invariants1 = {sAK[isAK], s1, s2};
2175  vector<double> invariants2 = invariants1;
2176  invariants2[1] = invariants1[2];
2177  invariants2[2] = invariants1[1];
2178  // Compute AP kernels.
2179  double AP1 = AltarelliParisi(invariants1);
2180  double AP2 = AltarelliParisi(invariants2);
2181  // Compute antennae (without subleading colour corrections).
2182  int modeSLCsave = modeSLC;
2183  modeSLC = 1;
2184  double ant1 = antFun(invariants1);
2185  double ant2 = antFun(invariants2);
2186  modeSLC = modeSLCsave;
2187  // Check if only a singularity on one side exists (-> z = -1).
2188  double zs1 = zA(invariants1);
2189  double zs2 = zB(invariants1);
2190  if (zs1 != -1.0) {
2191  double ratio = ant1/AP1;
2192  if (abs(ratio - 1.0) >= 0.01) {
2193  isOK = false;
2194  if (verbose >= 1)
2195  printOut(vinciaName() + ":check","WARNING: FAILED "
2196  "(collinear aj " + num2str(iTest, 1) + " and sAK = " +
2197  num2str((int)sqrt(sAK[isAK])) + "^2)");
2198  if (verbose >= 3) cout << setprecision(6) << " ant = "
2199  << num2str(ant1, 9) << " y1 = " << num2str(y1, 9)
2200  << " y2 = " << num2str(y2, 9) << " " <<endl
2201  << " P(z) = " << num2str(AP1, 9)
2202  << " z = " << num2str(zs1, 9) << endl;
2203  } else if (verbose >= 6) printOut(vinciaName() + ":check",
2204  "OK (collinear aj " + num2str(iTest, 1)
2205  + " and sAK = " + num2str((int)sqrt(sAK[isAK])) + "^2)");
2206  }
2207  if (zs2 != -1.0) {
2208  double ratio = ant2/AP2;
2209  if (abs(ratio - 1.0) >= 0.01) {
2210  isOK = false;
2211  if (verbose >= 1) printOut(vinciaName() + ":check",
2212  "WARNING: FAILED (collinear jk " + num2str(iTest, 1) +
2213  " and sAK = " + num2str((int)sqrt(sAK[isAK]))+"^2)");
2214  if (verbose >= 3) cout << setprecision(6) << " ant = "
2215  << num2str(ant1, 9) << " y1 = " << num2str(y1, 9) << " y2 = "
2216  << num2str(y2, 9) << " " << endl << " P(z) = "
2217  << num2str(AP2,9) << " z = " << num2str(zs2, 9) << endl;
2218  } else if (verbose >= 6) printOut(vinciaName() + ":check",
2219  "OK (collinear jk " + num2str(iTest, 1) +
2220  " and sAK = " + num2str((int)sqrt(sAK[isAK])) + "^2)");
2221  }
2222  }
2223  }
2224 
2225  // Test all helicity configurations for posivity and dead zones.
2226  string helString = " 9 9 -> 9 9 9";
2227  for (int iHel = -1; iHel < 32; ++iHel) {
2228  bool isPositive = true;
2229  bool isZero = true;
2230  bool hasDeadZone = false;
2231  if (iHel >= 0) {
2232  hj = 1 - 2*(iHel%2);
2233  hi = 1 - 2*((iHel/2)%2);
2234  hk = 1 - 2*((iHel/4)%2);
2235  hA = 1 - 2*((iHel/8)%2);
2236  hB = 1 - 2*((iHel/16)%2);
2237  helString = num2str(hA,2) + num2str(hB,2) + " ->" +
2238  num2str(hi,2) + num2str(hj,2) + num2str(hk,2);
2239  }
2240  vector<int> helBef = {hA, hB};
2241  vector<int> helNew = {hi, hj, hk};
2242 
2243  // Gluon emission: do not allow external helicity flips.
2244  if (id1() == 21 || id1() == 22)
2245  if (hA != hi || hB != hk) continue;
2246 
2247  // Check for positivity and dead zones.
2248  for (int iTest = 0; iTest <= nPointsCheck; ++iTest) {
2249  // sAK^0.5 between 1 and 1000 GeV.
2250  double sAKnow = pow2(1.0+rndmPtr->flat()*999.0);
2251  // Pick random xA.
2252  double xA = rndmPtr->flat();
2253  // sjk between 0 and sAK(1-xA)/xA, saj between 0 and sAK+sjk.
2254  double sjk = rndmPtr->flat()*sAKnow*(1.0-xA)/xA;
2255  double saj = rndmPtr->flat()*(sAKnow+sjk);
2256  double sak = sAKnow+sjk-saj;
2257  vector<double> invariants = {sAKnow, saj, sjk};
2258  // Check that sxy < shhMax.
2259  if (saj > shhMax || sjk > shhMax || sak > shhMax || sAKnow > shhMax)
2260  continue;
2261  // Restrict checks to ordered region; require pT2 < sAK
2262  double pT2 = saj * sjk / (sAKnow + sjk);
2263  if (pT2 > sAKnow) continue;
2264  // Compute antenna.
2265  double ant = antFun(invariants, mDum, helBef, helNew);
2266  // Check positivity (strict).
2267  if (ant < 0.0) {
2268  cout << "ant = "<<ant * sAKnow << endl;
2269  isPositive = false;
2270  if (verbose >= 3) printOut(vinciaName() + ":check",
2271  "ERROR ant(" + num2str((int)sqrt(saj))
2272  + "^2," + num2str((int)sqrt(sjk)) + "^2,"
2273  + num2str((int)sqrt(sAKnow)) + "^2 ; " + helString + ") < 0");
2274  } else if (ant > 0.0) isZero = false;
2275 
2276  // Check for dead zones away from phase-space boundary.
2277  double yjk = sjk / (sAKnow + sjk);
2278  double yAK = sAKnow / (sAKnow + sjk);
2279  double yaj = saj / (sAKnow + sjk);
2280  double yak = sak / (sAKnow + sjk);
2281  double yMin = 5*sqrt(deadZoneAvoidance);
2282  if ((sjk < 0.95*sAKnow*(1.0-xA)/xA)
2283  && yjk > yMin && yaj > yMin && yak > yMin && yAK > yMin
2284  && (ant*(sAKnow+sjk) < deadZoneAvoidance)) {
2285  hasDeadZone = true;
2286  }
2287  }
2288  isOK = isOK && isPositive;
2289 
2290  // Verbose output.
2291  if (!isPositive && verbose >= 1) printOut(vinciaName() + ":check",
2292  "ERROR (ant < 0 encountered " + helString + " )");
2293  else if (isPositive && !isZero && verbose >= 6) printOut(vinciaName() +
2294  ":check", "OK (is positive " + helString + " )");
2295  if (!isZero) {
2296  if (hasDeadZone && verbose >= 1) printOut(vinciaName() + ":check",
2297  "WARNING (dead zone encountered " + helString + " )");
2298  else if (!hasDeadZone && verbose >= 6) printOut(vinciaName()
2299  + ":check", "OK (no dead zones " + helString+" )");
2300  }
2301  } // End loop over helicities.
2302  } // End loop over initi-final state.
2303  return isOK;
2304 
2305 }
2306 
2307 //--------------------------------------------------------------------------
2308 
2309 // Check for resonances.
2310 
2311 bool AntennaFunctionIF::checkRes() {
2312 
2313  // Get test masses.
2314  vector<double> masses;
2315  getTestMasses(masses);
2316 
2317  // Test soft singularity for gluon emissions. Only test asymmetric
2318  // soft phase-space point (masses limit available phase space).
2319  if (id1() == 21 || id1() == 22) {
2320  // Get scaled invariants.
2321  double yaj = 0.01;
2322  double yjk = 0.0001;
2323 
2324  // Get dimensionful invariants.
2325  vector<double> invariants;
2326  if (!getTestInvariants(invariants, masses, yaj, yjk)) return false;
2327 
2328  // Get massive eikonal.
2329  double antSoft = massiveEikonal(invariants, masses);
2330 
2331  // Get full antenna.
2332  double antNow = antFun(invariants, masses);
2333 
2334  // Check ratio.
2335  double ratio= antNow/antSoft;
2336  if (abs(ratio - 1.) >= 0.001) {
2337  if (verbose >= quiet) {
2338  stringstream ss;
2339  ss << "WARNING:" + vinciaName() << " FAILED soft eikonal: ratio to "
2340  << "soft = " << ratio;
2341  printOut(__METHOD_NAME__, ss.str());
2342  }
2343  return false;
2344  }
2345  else if (verbose >= verylouddebug)
2346  printOut(__METHOD_NAME__, vinciaName() + " OK (soft eikonal)");
2347  }
2348 
2349  // Check collinear singularity (only check for final state parton,
2350  // quasi-collinear limit not relevant for res here). Test
2351  // invariants, for a few points along collinear axis.
2352  for (int iTest = 0; iTest < 4; ++iTest) {
2353  // Get scaled invariants.
2354  double yaj = 0.2 + double(iTest)*0.2;
2355  double yjk = 0.01;
2356 
2357  // Get dimensionful invariants.
2358  vector<double> invariants;
2359  if (!getTestInvariants(invariants, masses, yaj, yjk)) {
2360  if (verbose>=normal) infoPtr->errorMsg("Error in "+__METHOD_NAME__
2361  +": Failed to get test invariants!");
2362  return false;
2363  }
2364 
2365  // Are we still in the available phase space?
2366  double gramdet = gramDet(invariants, masses);
2367  if (gramdet<0.) {
2368  if (verbose >= verylouddebug) printOut(__METHOD_NAME__,
2369  vinciaName() + " not in phase space. Continue.");
2370  break;
2371  }
2372 
2373  // Get the antenna (with sum over flipped invariants where
2374  // appropriate).
2375  double antNow = antFunCollLimit(invariants, masses);
2376 
2377  // Get AP function.
2378  double AP = AltarelliParisi(invariants, masses);
2379  if (AP > 0.) {
2380  double ratio = antNow/AP;
2381  // Require better than 1% agreement unless dominated by nonsingular.
2382  if (abs(ratio - 1.) >= 0.01 && abs(antNow - AP) > 10.) {
2383  if (verbose >= 1) printOut(__METHOD_NAME__, "WARNING:" + vinciaName() +
2384  "Failed (collinear ij " + num2str(iTest, 1)+" )");
2385  if (verbose >= 3) cout << setprecision(6) << " ant = "
2386  << num2str(antNow, 9) << " yaj = " << num2str(yaj, 9)
2387  << " yjk = " << num2str(yjk, 9) << " " << endl << " P(z) = "
2388  << num2str(AP, 9) << endl;
2389  return false;
2390  }
2391  else if (verbose >= verylouddebug) printOut(__METHOD_NAME__, vinciaName()
2392  + " OK (collinear ij " + num2str(iTest, 1) + " )");
2393  }
2394  }
2395  return true;
2396 
2397 }
2398 
2399 //--------------------------------------------------------------------------
2400 
2401 // Create the test invariants for the checkRes method.
2402 
2403 bool AntennaFunctionIF::getTestInvariants(vector<double> &invariants,
2404  vector<double> masses, double yaj, double yjk) {
2405 
2406  if (masses.size() != 4) return false;
2407  double mA = masses[0];
2408  double mK = masses[2];
2409  double mAK = masses[3];
2410  double sAK = mA*mA + mK*mK - mAK*mAK;
2411  double sjk = sAK*yjk/(1. - yjk);
2412  if((sjk+sAK)==0.0) return false;
2413  double saj = yaj*(sAK+sjk);
2414  double sak = sAK + sjk -saj;
2415  double det = saj*sjk*sak - saj*saj*mK*mK - sjk*sjk*mA*mA;
2416  if (det <0.) return false;
2417  invariants = {sAK, saj, sjk, sak};
2418  return true;
2419 
2420 }
2421 
2422 //--------------------------------------------------------------------------
2423 
2424 // Wrapper for comparing to AP functions, sums over flipped
2425 // invariants where appropriate.
2426 
2427 double AntennaFunctionIF::antFunCollLimit(vector<double> invariants,
2428  vector<double> masses){
2429 
2430  double ant = antFun(invariants,masses);
2431  if (idB() == 21) {
2432  vector<double> flipped_invariants = {
2433  invariants[0], invariants[3], invariants[2], invariants[1]};
2434  ant += antFun(flipped_invariants,masses);
2435  }
2436  return ant;
2437 
2438 }
2439 
2440 //==========================================================================
2441 
2442 // Class QQEmitIF, initial-final antenna function.
2443 
2444 //--------------------------------------------------------------------------
2445 
2446 // The antenna function [GeV^-2].
2447 
2448 double QQEmitIF::antFun(vector<double> invariants, vector<double> masses,
2449  vector<int> helBef, vector<int> helNew) {
2450 
2451  // Invariants and helicities.
2452  double sAK = invariants[0];
2453  double saj = invariants[1];
2454  double sjk = invariants[2];
2455 
2456  // Sanity check. Require positive invariants.
2457  if ((saj <= 0.0) || (sjk <= 0.0) || (sAK <= 0.0)) return 0.0;
2458 
2459  // Initialise masses and helicities. Return 0 for unphysical helicities.
2460  initMasses(&masses);
2461  int nAvg = initHel(&helBef, &helNew);
2462  if (nAvg <= 0) return 0.;
2463 
2464  // Shorthands.
2465  double s = sAK + sjk;
2466  double yaj = saj/s;
2467  double yjk = sjk/s;
2468  double eikJ = 1./(sAK*yaj*yjk);
2469  double facMa = (mi != 0.) ? pow2(mi)/s/sAK/pow2(yaj) : 0.;
2470  double facMk = (mk != 0.) ? pow2(mk)/s/sAK/pow2(yjk) : 0.;
2471 
2472  // Do helicity sum.
2473  double hSum = 0.0;
2474 
2475  // (++ and --) parents.
2476  if (hA * hB > 0 || hA == 9 || hB == 9) {
2477  // ++ > +++ && -- > --- antennae (MHV).
2478  term = eikJ - facMa - facMk/(1. - yaj);
2479  if (RH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
2480  if (LH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
2481 
2482  // ++ > +-+ && -- > -+- antennae (NMHV).
2483  term = (pow2(1. - yaj) + (pow2(1. - yjk) - 1.)*pow2(1. - yaj)) * eikJ
2484  - facMa*pow2(1. - yjk-yaj) - facMk*(1. - yaj)*pow2(1. - yjk);
2485  if (RH[hA] && RH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += term;
2486  if (LH[hA] && LH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += term;
2487 
2488  // ++ > --+ && -- > ++- antennae (massive helicity flip).
2489  if (mi != 0.) {
2490  term = facMa * pow2(yjk);
2491  if (RH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
2492  if (LH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
2493  }
2494 
2495  // ++ > ++- && -- > --+ antennae (massive helicity flip).
2496  if (mk != 0.) {
2497  term = facMk * pow2(yaj)/(1. - yaj);
2498  if (RH[hA] && RH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
2499  if (LH[hA] && LH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
2500  }
2501  }
2502 
2503  // (+- and -+) parents.
2504  if (hA * hB < 0 || hA == 9 || hB == 9) {
2505  term = pow2(1.-yaj) * eikJ - facMa*(1. - yaj) - facMk*(1. - yaj);
2506 
2507  // +- > ++- && -+ > --+ antennae.
2508  if (RH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
2509  if (LH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
2510 
2511  // +- > +-- && -+ > -++ antennae.
2512  term = pow2(1.-yjk) * eikJ - facMa*pow2(1. - yjk)
2513  - facMk*pow2(1. - yjk)/(1. - yaj);
2514  if (RH[hA] && LH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
2515  if (LH[hA] && RH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
2516 
2517  // +- > --- && -+ > +++ antennae (massive helicity flip).
2518  if (mi != 0.) {
2519  term = facMa * pow2(yjk);
2520  if (RH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
2521  if (LH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
2522  }
2523 
2524  // +- > +-+ && -+ > -+- antennae (massive helicity flip).
2525  if (mk != 0.) {
2526  term = facMk * pow2(yaj) / (1. - yaj);
2527  if (RH[hA] && LH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += term;
2528  if (LH[hA] && RH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += term;
2529  }
2530  }
2531 
2532  // Return helicity sum, averaged over initial helicities.
2533  return hSum/nAvg;
2534 
2535 }
2536 
2537 //--------------------------------------------------------------------------
2538 
2539 // The AP kernel, P(z)/Q2.
2540 
2541 double QQEmitIF::AltarelliParisi(vector<double> invariants, vector<double>,
2542  vector<int>, vector<int>) {
2543 
2544  // Sanity check. Require positive invariants.
2545  double sAK = invariants[0];
2546  double saj = invariants[1];
2547  double sjk = invariants[2];
2548  if ((saj <= 0.0) || (sjk <= 0.0) || (sAK <= 0.0)) return 0.0;
2549 
2550  // Using smaller invariant for collinear limit.
2551  double z = ( saj < sjk ? zA(invariants) : zB(invariants) );
2552  double Q2 = min(saj,sjk);
2553  double AP = 0.0;
2554 
2555  // Collinear to initial quark.
2556  if (saj < sjk) {
2557  AP = 1.0/z * (1.0 + pow2(z))/(1.0 - z);
2558  // AP uses alpha/2pi * CF, ant uses alpha/4pi * 2CF, with CF = 4/3.
2559  AP *= 1.0;
2560  // Ant with quark reproduces full AP.
2561  AP *= 1.0;
2562 
2563  // Collinear to final quark.
2564  } else {
2565  AP = (1.0 + pow2(z))/(1.0 - z);
2566  // AP uses alpha/2pi * CF, ant uses alpha/4pi * 2CF, with CF = 4/3.
2567  AP *= 1.0;
2568  // ant with quark reproduces full AP.
2569  AP *= 1.0;
2570  }
2571  return AP/Q2;
2572 
2573 }
2574 
2575 //==========================================================================
2576 
2577 // Class QGEmitIF, initial-final antenna function.
2578 
2579 //--------------------------------------------------------------------------
2580 
2581 // The antenna function [GeV^-2].
2582 
2583 double QGEmitIF::antFun(vector<double> invariants, vector<double> masses,
2584  vector<int> helBef, vector<int> helNew) {
2585 
2586  // Invariants and helicities.
2587  double sAK = invariants[0];
2588  double saj = invariants[1];
2589  double sjk = invariants[2];
2590 
2591  // Sanity check. Require positive invariants.
2592  if ((saj <= 0.0) || (sjk <= 0.0) || (sAK <= 0.0)) return 0.0;
2593 
2594  // Initialise masses and helicities. Return 0 for unphysical helicities.
2595  initMasses(&masses);
2596  int nAvg = initHel(&helBef, &helNew);
2597  if (nAvg <= 0) return 0.;
2598 
2599  // Shorthands.
2600  double s = sAK + sjk;
2601  double yaj = saj/s;
2602  double yjk = sjk/s;
2603  double eikJ = 1./(sAK*yaj*yjk);
2604  double colK = (1. - alphaSav) * (1. - 2.*yaj) / (sAK * yjk);
2605  double facMa = (mi != 0.) ? pow2(mi)/s/sAK/pow2(yaj) : 0.;
2606 
2607  // Do helicity sum.
2608  double hSum = 0.0;
2609 
2610  // (++ and --) parents.
2611  if (hA * hB > 0 || hA == 9 || hB == 9) {
2612  // ++ > +++ && -- > --- antennae (MHV).
2613  term = eikJ + colK - facMa;
2614  if (RH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
2615  if (LH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
2616 
2617  // ++ > +-+ && -- > -+- antennae (NMHV).
2618  term = (pow3(1. - yaj)+pow2(1. - yjk)-1) * eikJ
2619  - facMa*pow2(1. - yjk - yaj)*(1. - yaj) + (3. - pow2(yaj))/sAK;
2620  if (RH[hA] && RH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += term;
2621  if (LH[hA] && LH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += term;
2622 
2623  // ++ > --+ && -- > ++- antennae (massive helicity flip).
2624  if (mi != 0.) {
2625  term = facMa*pow2(yjk);
2626  if (RH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
2627  if (LH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
2628  }
2629  }
2630 
2631  // (+- and -+) parents.
2632  if (hA * hB < 0 || hA == 9 || hB == 9) {
2633  // +- > ++- && -+ > --+ antennae.
2634  term = pow3(1. - yaj) * eikJ - facMa*pow2(1. - yaj);
2635  if (RH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
2636  if (LH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
2637 
2638  // +- > +-- && -+ > -++ antennae.
2639  term = pow2(1. - yjk) * eikJ + colK - facMa*pow2(1. - yjk)
2640  + (2*yaj - yjk)/sAK;
2641  if (RH[hA] && LH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
2642  if (LH[hA] && RH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
2643 
2644  // +- > --- && -+ > +++ antennae (massive helicity flip).
2645  if (mi != 0.) {
2646  term = facMa * pow2(yjk);
2647  if (RH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
2648  if (LH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
2649  }
2650  }
2651 
2652  // Subleading colour correction.
2653  if (modeSLC >= 2) hSum *= CF/chargeFacSav * (1 - yaj)/(2 - yaj - yjk)
2654  + CA/chargeFacSav * (1 - yjk)/(2 - yaj - yjk);
2655 
2656  // Return helicity sum, averaged over initial helicities.
2657  return hSum/nAvg;
2658 
2659 }
2660 
2661 //--------------------------------------------------------------------------
2662 
2663 // The AP kernel, P(z)/Q2.
2664 
2665 double QGEmitIF::AltarelliParisi(vector<double> invariants, vector<double>,
2666  vector<int>, vector<int>) {
2667 
2668  // Sanity check. Require positive invariants.
2669  double sAK = invariants[0];
2670  double saj = invariants[1];
2671  double sjk = invariants[2];
2672  if ((saj <= 0.0) || (sjk <= 0.0) || (sAK <= 0.0)) return 0.0;
2673 
2674  // Using smaller invariant for collinear limit.
2675  double z = ( saj < sjk ? zA(invariants) : zB(invariants) );
2676  double Q2 = min(saj, sjk);
2677  double AP = 0.0;
2678 
2679  // Collinear to initial quark.
2680  if (saj < sjk) {
2681  AP = 1.0/z * (1.0 + pow2(z))/(1.0 - z);
2682  // AP uses alpha/2pi * CF, ant uses alpha/4pi * 2CF, with CF = 4/3.
2683  AP *= 1.0;
2684  // ant with quark reproduces full AP.
2685  AP *= 1.0;
2686 
2687  // Collinear to final gluon.
2688  } else {
2689  AP = (2.0*z/(1.0 - z)+z*(1.0 - z));
2690  // AP uses alpha/2pi * CA, ant uses alpha/4pi * CA, with CA = 3.
2691  AP *= 2.0;
2692  // gluon appears in 2 ants, hence ant reproduces half of AP.
2693  AP *= 0.5;
2694  }
2695  return AP/Q2;
2696 
2697 }
2698 
2699 //==========================================================================
2700 
2701 // Class GQEmitIF, initial-final antenna function.
2702 
2703 //--------------------------------------------------------------------------
2704 
2705 // The antenna function [GeV^-2].
2706 
2707 double GQEmitIF::antFun(vector<double> invariants, vector<double> masses,
2708  vector<int> helBef, vector<int> helNew) {
2709 
2710  // Invariants and helicities.
2711  double sAK = invariants[0];
2712  double saj = invariants[1];
2713  double sjk = invariants[2];
2714 
2715  // Sanity check. Require positive invariants.
2716  if ((saj <= 0.0) || (sjk <= 0.0) || (sAK <= 0.0)) return 0.0;
2717 
2718  // Initialise masses and helicities. Return 0 for unphysical helicities.
2719  initMasses(&masses);
2720  int nAvg = initHel(&helBef, &helNew);
2721  if (nAvg <= 0) return 0.;
2722 
2723  // Shorthands
2724  double s = sAK + sjk;
2725  double yaj = saj/s;
2726  double yjk = sjk/s;
2727  double eikJ = 1./(sAK*yaj*yjk);
2728  double colA = 1./(sAK*yaj*(1. - yjk));
2729  double facMk = (mk != 0.) ? pow2(mk)/s/sAK/pow2(yjk) : 0.;
2730 
2731  // Do helicity sum.
2732  double hSum = 0.0;
2733 
2734  // (++ and --) parents.
2735  if (hA * hB > 0 || hA == 9 || hB == 9) {
2736  // ++ > +++ && -- > --- antennae (MHV).
2737  term = eikJ + colA - facMk/(1. - yaj);
2738  if (RH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
2739  if (LH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
2740 
2741  // ++ > +-+ && -- > -+- antennae (NMHV).
2742  term = (pow2(1. - yaj)+(pow3(1. - yjk) - 1)*(pow2(1. - yaj)))* eikJ
2743  - facMk*(1. - yaj)*pow3(1 - yjk);
2744  if (RH[hA] && RH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += term;
2745  if (LH[hA] && LH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += term;
2746 
2747  // ++ > --+ && -- > ++- antennae.
2748  term = pow3(yjk) * colA;
2749  if (RH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
2750  if (LH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
2751 
2752  // ++ > ++- && -- > --+ antennae (massive helicity flip).
2753  if (mk != 0.) {
2754  term = facMk * pow2(yaj)/(1. - yaj);
2755  if (RH[hA] && RH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
2756  if (LH[hA] && LH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
2757  }
2758  }
2759 
2760  // (+- and -+) parents.
2761  if (hA * hB < 0 || hA == 9 || hB == 9) {
2762  // +- > ++- && -+ > --+ antennae.
2763  term = pow2(1 - yaj) * eikJ + colA - facMk*(1. - yaj);
2764  if (RH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
2765  if (LH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
2766 
2767  // +- > +-- && -+ > -++ antennae.
2768  term = pow3(1. - yjk) * eikJ - facMk*pow3(1. - yjk)/(1. - yaj);
2769  if (RH[hA] && LH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
2770  if (LH[hA] && RH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
2771 
2772  // +- > --- && -+ > +++ antennae.
2773  term = pow3(yjk) * colA;
2774  if (RH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
2775  if (LH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
2776 
2777  // +- > +-+ && -+ > -+- antennae (massive helicity flip).
2778  if (mk != 0.) {
2779  term = facMk * pow2(yaj)/(1. - yaj);
2780  if (RH[hA] && LH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += term;
2781  if (LH[hA] && RH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += term;
2782  }
2783  }
2784 
2785  // Subleading colour correction.
2786  if (modeSLC >= 2) hSum *= CF/chargeFacSav * (1 - yjk)/(2 - yaj - yjk)
2787  + CA/chargeFacSav * (1 - yaj)/(2 - yaj - yjk);
2788 
2789  // Return helicity sum, averaged over initial helicities.
2790  return hSum/nAvg;
2791 
2792 }
2793 
2794 //--------------------------------------------------------------------------
2795 
2796 // The AP kernel, P(z)/Q2.
2797 
2798 double GQEmitIF::AltarelliParisi(vector<double> invariants, vector<double>,
2799  vector<int>, vector<int>) {
2800 
2801  // Sanity check. Require positive invariants.
2802  double sAK = invariants[0];
2803  double saj = invariants[1];
2804  double sjk = invariants[2];
2805  if ((saj <= 0.0) || (sjk <= 0.0) || (sAK <= 0.0)) return 0.0;
2806 
2807  // Using smaller invariant for collinear limit.
2808  double z = ( saj < sjk ? zA(invariants) : zB(invariants) );
2809  double Q2 = min(saj,sjk);
2810  double AP = 0.0;
2811 
2812  // Collinear to initial gluon.
2813  if (saj < sjk) {
2814  AP = 1.0/z*(pow(z, 4.0) + 1.0 + pow(1 - z, 4.0))/z/(1.0 - z);
2815  // AP uses alpha/2pi * CA, ant uses alpha/4pi * CA, with CA = 3.
2816  AP *= 2.0;
2817  // gluon appears in 2 ants, hence ant reproduces half of AP.
2818  AP *= 0.5;
2819 
2820  // Collinear to final quark.
2821  } else {
2822  AP = (1.0 + pow2(z))/(1.0 - z);
2823  // AP uses alpha/2pi * CF, ant uses alpha/4pi * 2CF, with CF = 4/3.
2824  AP *= 1.0;
2825  // ant with quark reproduces full AP.
2826  AP *= 1.0;
2827  }
2828  return AP/Q2;
2829 
2830 }
2831 
2832 //==========================================================================
2833 
2834 // Class GGEmitIF, initial-final antenna function.
2835 
2836 //--------------------------------------------------------------------------
2837 
2838 // The antenna function [GeV^-2].
2839 
2840 double GGEmitIF::antFun(vector<double> invariants, vector<double>,
2841  vector<int> helBef, vector<int> helNew) {
2842 
2843  // Invariants and helicities.
2844  double sAK = invariants[0];
2845  double saj = invariants[1];
2846  double sjk = invariants[2];
2847 
2848  // Sanity check. Require positive invariants.
2849  if ((saj <= 0.0) || (sjk <= 0.0) || (sAK <= 0.0)) return 0.0;
2850 
2851  // Initialise helicities.
2852  int nAvg = initHel(&helBef, &helNew);
2853 
2854  // Shorthands.
2855  double s = sAK + sjk;
2856  double yaj = saj/s;
2857  double yjk = sjk/s;
2858  double yAK = sAK/s;
2859  double eikJ = 1./(sAK*yaj*yjk);
2860  double colA = 1./(sAK*yaj*yAK);
2861  double colK = (1.-alphaSav) * (1.-2.*yaj) / (sAK * yjk);
2862 
2863  // Do helicity sum.
2864  double hSum = 0.0;
2865 
2866  // (++ and --) parents.
2867  if (hA * hB > 0 || hA == 9 || hB == 9) {
2868  // ++ > +++ && -- > --- antennae (MHV).
2869  term = eikJ + colA + colK;
2870  if (RH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
2871  if (LH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
2872 
2873  // ++ > +-+ && -- > -+- antennae (NMHV).
2874  term = (pow3(1. - yaj) + pow3(1. - yjk) - 1.) * eikJ
2875  + ( 6. - 3.*(yaj+yjk) + yaj*yjk )/sAK;
2876  if (RH[hA] && RH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += term;
2877  if (LH[hA] && LH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += term;
2878 
2879  // ++ > --+ && -- > ++- antennae.
2880  term = pow3(yjk) * colA;
2881  if (RH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
2882  if (LH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
2883  }
2884 
2885  // (+- and -+) parents.
2886  if (hA * hB < 0 || hA == 9 || hB == 9) {
2887  // +- > ++- && -+ > --+ antennae.
2888  term = pow3(1. - yaj) * eikJ + colA;
2889  if (RH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
2890  if (LH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
2891 
2892  // +- > +-- && -+ > -++ antennae.
2893  term = pow3(1. - yjk) * eikJ + colK + (3.*yaj - yjk - yaj*yjk)/sAK;
2894  if (RH[hA] && LH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
2895  if (LH[hA] && RH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
2896 
2897  // +- > --- && -+ > +++ antennae.
2898  term = pow3(yjk) * colA;
2899  if (RH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
2900  if (LH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
2901  }
2902 
2903  // Return helicity sum, averaged over initial helicities.
2904  return hSum/nAvg;
2905 }
2906 
2907 //--------------------------------------------------------------------------
2908 
2909 // The AP kernel, P(z)/Q2.
2910 
2911 double GGEmitIF::AltarelliParisi(vector<double> invariants, vector<double>,
2912  vector<int>, vector<int>) {
2913 
2914  // Sanity check. Require positive invariants.
2915  double sAK = invariants[0];
2916  double saj = invariants[1];
2917  double sjk = invariants[2];
2918  if ((saj <= 0.0) || (sjk <= 0.0) || (sAK <= 0.0)) return 0.0;
2919 
2920  // Using smaller invariant for collinear limit.
2921  double z = ( saj < sjk ? zA(invariants) : zB(invariants) );
2922  double Q2 = min(saj,sjk);
2923  double AP = 0.0;
2924 
2925  // Collinear to initial gluon.
2926  if (saj < sjk) {
2927  AP = 1.0/z*(pow(z, 4.0) + 1.0 + pow(1 - z,4.0))/z/(1.0 - z);
2928  // AP uses alpha/2pi * CA, ant uses alpha/4pi * CA, with CA = 3.
2929  AP *= 2.0;
2930  // gluon appears in 2 ants, hence ant reproduces half of AP.
2931  AP *= 0.5;
2932 
2933  // Collinear to final gluon.
2934  } else {
2935  AP = (2.0*z/(1.0 - z) + z*(1.0 - z));
2936  // AP uses alpha/2pi * CA, ant uses alpha/4pi * CA, with CA = 3.
2937  AP *= 2.0;
2938  // gluon appears in 2 ants, hence ant reproduces half of AP.
2939  AP *= 0.5;
2940  }
2941  return AP/Q2;
2942 
2943 }
2944 
2945 //==========================================================================
2946 
2947 // Class QXSplitIF, initial-final antenna function.
2948 
2949 //--------------------------------------------------------------------------
2950 
2951 // The antenna function [GeV^-2].
2952 
2953 double QXSplitIF::antFun(vector<double> invariants, vector<double> masses,
2954  vector<int> helBef, vector<int> helNew) {
2955 
2956  // Invariants and helicities.
2957  double sAK = invariants[0];
2958  double saj = invariants[1];
2959  double sjk = invariants[2];
2960 
2961  // Sanity check. Require positive invariants.
2962  if ((saj <= 0.0) || (sjk <= 0.0) || (sAK <= 0.0)) return 0.0;
2963 
2964  // Initialise masses and helicities. Return 0 for unphysical helicities.
2965  initMasses(&masses);
2966  int nAvg = initHel(&helBef, &helNew);
2967  if (nAvg <= 0) return 0.;
2968 
2969  // Shorthands
2970  double s = sAK + sjk;
2971  double yAK = sAK/s;
2972  double yaj = saj/s;
2973  double colA = 1./(sAK * yaj);
2974  double facMj = (mj != 0.) ? pow2(mj)/s/sAK/pow2(yaj) : 0.;
2975 
2976  // Do helicity sum.
2977  double hSum = 0.0;
2978 
2979  // (++ and --) parents.
2980  if (hA * hB > 0 || hA == 9 || hB == 9) {
2981  // ++ > +-+ && -- > -+- antennae.
2982  term = pow2(yAK) * colA - facMj*pow2(yAK)/(1. - yAK);
2983  if (RH[hA] && RH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += term;
2984  if (LH[hA] && LH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += term;
2985 
2986  // ++ > --+ && -- > ++- antennae.
2987  term = pow2(1. - yAK) * colA - facMj*(1. - yAK);
2988  if (RH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
2989  if (LH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
2990 
2991  // ++ > +++ && -- > --- antennae (massive helicity flip).
2992  if (mj != 0.) {
2993  term = facMj / (1.-yAK);
2994  if (RH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
2995  if (LH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
2996  }
2997  }
2998 
2999  // (+- and -+) parents.
3000  if (hA * hB < 0 || hA == 9 || hB == 9) {
3001  // +- > +-- && -+ > -++ antennae.
3002  term = pow2(yAK) * colA - facMj*pow2(yAK)/(1. - yAK);
3003  if (RH[hA] && LH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
3004  if (LH[hA] && RH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
3005 
3006  // +- > --- && -+ > +++ antennae.
3007  term = pow2(1. - yAK) * colA - facMj*(1. - yAK);
3008  if (RH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
3009  if (LH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
3010 
3011  // +- > ++- && -+ > --+ antennae (massive helicity flip).
3012  if (mj != 0.) {
3013  term = facMj/(1. - yAK);
3014  if (RH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
3015  if (LH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
3016  }
3017  }
3018 
3019  // Return helicity sum, averaged over initial helicities.
3020  return hSum/nAvg;
3021 
3022 }
3023 
3024 //--------------------------------------------------------------------------
3025 
3026 // The AP kernel, P(z)/Q2.
3027 
3028 double QXSplitIF::AltarelliParisi(vector<double> invariants, vector<double>,
3029  vector<int>, vector<int>) {
3030 
3031  // Sanity check. Require positive invariants.
3032  double sAK = invariants[0];
3033  double saj = invariants[1];
3034  double sjk = invariants[2];
3035  if ((saj <= 0.0) || (sjk <= 0.0) || (sAK <= 0.0)) return 0.0;
3036 
3037  // There is only the saj collinear limit.
3038  double z = zA(invariants);
3039  double Q2 = saj;
3040  double AP = 1.0/z * (pow2(z) + pow2(1 - z));
3041  // AP uses alpha/2pi * TR/2, ant uses alpha/4pi * TR, with TR = 1.
3042  AP *= 1.0;
3043  // Ant with quark reproduces full AP.
3044  AP *= 1.0;
3045  return AP/Q2;
3046 
3047 }
3048 
3049 //==========================================================================
3050 
3051 // Class GXConvIF, initial-final antenna function.
3052 
3053 //--------------------------------------------------------------------------
3054 
3055 // The antenna function [GeV^-2].
3056 
3057 double GXConvIF::antFun(vector<double> invariants, vector<double> masses,
3058  vector<int> helBef, vector<int> helNew) {
3059 
3060  // Invariants and helicities.
3061  double sAK = invariants[0];
3062  double saj = invariants[1];
3063  double sjk = invariants[2];
3064 
3065  // Sanity check. Require positive invariants.
3066  if ((saj <= 0.0) || (sjk <= 0.0) || (sAK <= 0.0)) return 0.0;
3067 
3068  // Initialise masses and helicities. Return 0 for unphysical helicities.
3069  initMasses(&masses);
3070  int nAvg = initHel(&helBef, &helNew);
3071  if (nAvg <= 0) return 0.;
3072 
3073  // Shorthands
3074  double s = sAK + sjk + 2.*pow2(mj);
3075  double yaj = saj / s;
3076  double yAK = sAK / s;
3077  double mu2j = (mj != 0.) ? pow2(mj)/s : 0.;
3078  double colA = 1./(2.*sAK*yAK*(yaj - 2.*mu2j));
3079  double facMj = (mj != 0.) ? mu2j/(2.*sAK)/pow2(yaj - 2*mu2j) : 0.;
3080 
3081  // Do helicity sum.
3082  double hSum = 0.0;
3083 
3084  // (++ and --) parents.
3085  if (hA * hB > 0 || hA == 9 || hB == 9) {
3086  // ++ > +++ && -- > --- antennae.
3087  term = colA - facMj*yAK/(1. - yAK);
3088  if (RH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
3089  if (LH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
3090 
3091  // ++ > --+ && -- > ++- antennae.
3092  term = pow2(1. - yAK) * colA - facMj*yAK*(1. - yAK);
3093  if (RH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
3094  if (LH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
3095 
3096  // ++ > +-+ && -- > -+- antennae (massive helicity flip).
3097  if (mj != 0.) {
3098  term = facMj * pow3(yAK)/(1. - yAK);
3099  if (RH[hA] && RH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += term;
3100  if (LH[hA] && LH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += term;
3101  }
3102  }
3103 
3104  // (+- and -+) parents.
3105  if (hA * hB < 0 || hA == 9 || hB == 9) {
3106  // +- > ++- && -+ > --+ antennae.
3107  term = colA - facMj*yAK/(1. - yAK);
3108  if (RH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
3109  if (LH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
3110 
3111  // +- > --- && -+ > +++ antennae.
3112  term = pow2(1. - yAK) * colA - facMj*yAK*(1. - yAK) ;
3113  if (RH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
3114  if (LH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
3115 
3116  // +- > +-- && -+ > -++ antennae (massive helicity flip).
3117  if (mj != 0.) {
3118  term = facMj * pow3(yAK)/(1. - yAK);
3119  if (RH[hA] && LH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
3120  if (LH[hA] && RH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
3121  }
3122  }
3123 
3124  // Return helicity sum, averaged over initial helicities.
3125  return hSum/nAvg;
3126 
3127 }
3128 
3129 //--------------------------------------------------------------------------
3130 
3131 // The AP kernel, P(z)/Q2.
3132 
3133 double GXConvIF::AltarelliParisi(vector<double> invariants, vector<double>,
3134  vector<int>, vector<int>) {
3135 
3136  // Sanity check. Require positive invariants.
3137  double sAK = invariants[0];
3138  double saj = invariants[1];
3139  double sjk = invariants[2];
3140  if ((saj <= 0.0) || (sjk <= 0.0) || (sAK <= 0.0)) return 0.0;
3141 
3142  // There is only the saj collinear limit.
3143  double z = zA(invariants);
3144  double Q2 = saj;
3145  double AP = 1.0/z * (1.0 + pow2(1.0 - z))/z;
3146  // AP uses alpha/2pi * CF, ant uses alpha/4pi * 2CF, with CF = 4/3.
3147  AP *= 1.0;
3148  // gluon appears in 2 ants, hence ant reproduces half of AP.
3149  AP *= 0.5;
3150  return AP/Q2;
3151 
3152 }
3153 
3154 //==========================================================================
3155 
3156 // Class XGSplitIF, initial-final antenna function. Gluon splitting in
3157 // the final state.
3158 
3159 //--------------------------------------------------------------------------
3160 
3161 // The antenna function [GeV^-2].
3162 
3163 double XGSplitIF::antFun(vector<double> invariants, vector<double> masses,
3164  vector<int> helBef, vector<int> helNew) {
3165 
3166  // Invariants and helicities
3167  double sAK = invariants[0];
3168  double saj = invariants[1];
3169  double sjk = invariants[2];
3170 
3171  // Sanity check. Require positive invariants.
3172  if ((saj <= 0.0) || (sjk <= 0.0) || (sAK <= 0.0)) return 0.0;
3173 
3174  // Initialise masses and helicities. Return 0 for unphysical helicities.
3175  initMasses(&masses);
3176  int nAvg = initHel(&helBef, &helNew);
3177  if (nAvg <= 0) return 0.;
3178 
3179  // Shorthands.
3180  double s = sAK + sjk + 2*pow2(mj);
3181  double yaj = saj / s;
3182  double yak = 1. - yaj;
3183  double m2jk = sjk + 2*pow2(mj);
3184  double colK = 1./(2.*m2jk);
3185  double facMj = pow2(mj)/(2.*pow2(m2jk));
3186 
3187  // Do helicity sum.
3188  double hSum = 0.0;
3189 
3190  // (++ and --) parents.
3191  if (hA * hB > 0 || hA == 9 || hB == 9) {
3192  // ++ > +-+ && -- > -+- antennae.
3193  term = pow2(yak) * colK - facMj*yak/(1. - yak);
3194  if (RH[hA] && RH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += term;
3195  if (LH[hA] && LH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += term;
3196 
3197  // ++ > ++- && -- > --+ antennae.
3198  term = pow2(yaj) * colK - facMj*yaj/(1. - yaj);
3199  if (RH[hA] && RH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
3200  if (LH[hA] && LH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
3201 
3202  // ++ > +++ && -- > --- antennae (massive helicity flip).
3203  if (mj != 0.) {
3204  term = facMj*(yaj/(1. - yaj) + yak/(1. - yak) + 2);
3205  if (RH[hA] && RH[hB] && RH[hi] && RH[hj] && RH[hk]) hSum += term;
3206  if (LH[hA] && LH[hB] && LH[hi] && LH[hj] && LH[hk]) hSum += term;
3207  }
3208  }
3209 
3210  // (+- and -+) parents.
3211  if (hA * hB < 0 || hA == 9 || hB == 9) {
3212  // +- > ++- && -+ > --+ antennae.
3213  term = pow2(yak) * colK - facMj*yak/(1. - yak);
3214  if (RH[hA] && LH[hB] && RH[hi] && RH[hj] && LH[hk]) hSum += term;
3215  if (LH[hA] && RH[hB] && LH[hi] && LH[hj] && RH[hk]) hSum += term;
3216 
3217  // +- > +-+ && -+ > -+- antennae.
3218  term = pow2(yaj) * colK - facMj*yaj/(1. - yaj);
3219  if (RH[hA] && LH[hB] && RH[hi] && LH[hj] && RH[hk]) hSum += term;
3220  if (LH[hA] && RH[hB] && LH[hi] && RH[hj] && LH[hk]) hSum += term;
3221 
3222  // +- > +-- && -+ > -++ antennae (massive helicity flip).
3223  if (mj != 0.) {
3224  term = facMj*(yaj/(1. - yaj) + yak/(1. - yak) + 2);
3225  if (RH[hA] && LH[hB] && RH[hi] && LH[hj] && LH[hk]) hSum += term;
3226  if (LH[hA] && RH[hB] && LH[hi] && RH[hj] && RH[hk]) hSum += term;
3227  }
3228  }
3229 
3230  // Return helicity sum, averaged over initial helicities.
3231  return hSum/nAvg;
3232 
3233 }
3234 
3235 //--------------------------------------------------------------------------
3236 
3237 // The AP kernel, P(z)/Q2.
3238 
3239 double XGSplitIF::AltarelliParisi(vector<double> invariants, vector<double>,
3240  vector<int>, vector<int>) {
3241 
3242  // Sanity check. Require positive invariants.
3243  double sAK = invariants[0];
3244  double saj = invariants[1];
3245  double sjk = invariants[2];
3246  if ((saj <= 0.0) || (sjk <= 0.0) || (sAK <= 0.0)) return 0.0;
3247 
3248  // There is only the sjk collinear limit
3249  double z = zB(invariants);
3250  double Q2 = sjk;
3251  double AP = (pow2(z)+pow2(1 - z));
3252  // AP uses alpha/2pi * TR/2, ant uses alpha/4pi * TR, with TR = 1.
3253  AP *= 1.0;
3254  // Gluon appears in 2 ants, hence ant reproduces half of AP.
3255  AP *= 0.5;
3256  return AP/Q2;
3257 
3258 }
3259 
3260 //==========================================================================
3261 
3262 // Class QGEmitIFsec, derived class for sector initial-final antenna
3263 // functions.
3264 
3265 //--------------------------------------------------------------------------
3266 
3267 // The antenna function [GeV^-2].
3268 
3269 double QGEmitIFsec::antFun(vector<double> invariants,
3270  vector<double> mNew, vector<int> helBef, vector<int> helNew) {
3271 
3272  // Check if helicity vectors empty.
3273  double ant = QGEmitIF::antFun(invariants, mNew, helBef, helNew);
3274  if (helBef.size() < 2) {helBef.push_back(9); helBef.push_back(9);}
3275  if (helNew.size() < 3) {
3276  helNew.push_back(9); helNew.push_back(9); helNew.push_back(9);}
3277 
3278  // Check if j has same helicity as parent final-state gluon.
3279  int hG = helBef[1];
3280  int hjNow = helNew[1];
3281  if (hG == hjNow) {
3282  // Invariants are sAK, saj, sjk.
3283  // Define j<->k symmetrisation term, sak = sAK - saj + sjk.
3284  vector<double> invariantsSym = invariants;
3285  double s02 = invariants[0] - invariants[1] + invariants[2];
3286  vector<int> helSym = helNew;
3287  invariantsSym[1] = s02 + sectorDampSav * invariants[2];
3288  helSym[1] = helNew[2];
3289  helSym[2] = helNew[1];
3290  ant += QGEmitIF::antFun(invariantsSym, mNew, helBef, helSym);
3291  }
3292  return ant;
3293 
3294 }
3295 
3296 //==========================================================================
3297 
3298 // Class GGEmitIFsec, sector initial-final antenna function.
3299 
3300 //--------------------------------------------------------------------------
3301 
3302 // The antenna function [GeV^-2].
3303 
3304 double GGEmitIFsec::antFun(vector<double> invariants,
3305  vector<double> mNew, vector<int> helBef, vector<int> helNew) {
3306 
3307  // Check if helicity vectors empty.
3308  double ant = GGEmitIF::antFun(invariants, mNew, helBef, helNew);
3309  if (helBef.size() < 2) {helBef.push_back(9); helBef.push_back(9);}
3310  if (helNew.size() < 3) {
3311  helNew.push_back(9); helNew.push_back(9); helNew.push_back(9);}
3312 
3313  // Check if j has same helicity as parent final-state gluon.
3314  int hG = helBef[1];
3315  int hjNow = helNew[1];
3316  if ( hG == hjNow ) {
3317  // Invariants are sAK, saj, sjk
3318  // Define j<->k symmetrisation term: sak = sAK - saj + sjk.
3319  vector<double> invariantsSym = invariants;
3320  double s02 = invariants[0] - invariants[1] + invariants[2];
3321  vector<int> helSym = helNew;
3322  invariantsSym[1] = s02 + sectorDampSav * invariants[2];
3323  helSym[1] = helNew[2];
3324  helSym[2] = helNew[1];
3325  ant += GGEmitIF::antFun(invariantsSym, mNew, helBef, helSym);
3326  }
3327  return ant;
3328 
3329 }
3330 
3331 //==========================================================================
3332 
3333 // Class XGSplitIFsec, sector initial-final antenna function. Gluon
3334 // splitting in the final state.
3335 
3336 //--------------------------------------------------------------------------
3337 
3338 // The antenna function, just 2*global [GeV^-2].
3339 
3340 double XGSplitIFsec::antFun(vector<double> invariants, vector<double> mNew,
3341  vector<int> helBef, vector<int> helNew) {
3342  return 2*XGSplitIF::antFun(invariants,mNew,helBef,helNew);}
3343 
3344 //==========================================================================
3345 
3346 // The AntennaSetFSR class. Simple container of FF and RF antenna functions.
3347 
3348 //--------------------------------------------------------------------------
3349 
3350 // Initialize pointers.
3351 
3352 void AntennaSetFSR::initPtr(Info* infoPtrIn, DGLAP* dglapPtrIn) {
3353 
3354  infoPtr = infoPtrIn;
3355  particleDataPtr = infoPtr->particleDataPtr;
3356  settingsPtr = infoPtr->settingsPtr;
3357  rndmPtr = infoPtr->rndmPtr;
3358  dglapPtr = dglapPtrIn;
3359  isInitPtr = true;
3360 
3361 }
3362 
3363 //--------------------------------------------------------------------------
3364 
3365 // Initialize antenna set (optionally with min or max variation).
3366 
3367 void AntennaSetFSR::init() {
3368 
3369  // Check pointers and initialization.
3370  if (!isInitPtr) {
3371  printOut(__METHOD_NAME__, "Cannot initialize, pointers not set.");
3372  return;
3373  }
3374  verbose = settingsPtr->mode("Vincia:verbose");
3375  if (isInit) {
3376  if (verbose >= debug)
3377  printOut(__METHOD_NAME__, "Already initialized antenna set.");
3378  return;
3379  }
3380 
3381  // Create antenna objects (in map, so indices don't have to be consecutive).
3382  antFunPtrs.clear();
3383  bool sectorShower = settingsPtr->flag("Vincia:sectorShower");
3384  antFunPtrs[iQQemitFF] = sectorShower ? new QQEmitFFsec() : new QQEmitFF();
3385  antFunPtrs[iQGemitFF] = sectorShower ? new QGEmitFFsec() : new QGEmitFF();
3386  antFunPtrs[iGQemitFF] = sectorShower ? (AntennaFunction*)new GQEmitFFsec() :
3387  new GQEmitFF();
3388  antFunPtrs[iGGemitFF] = sectorShower ? new GGEmitFFsec() : new GGEmitFF();
3389  antFunPtrs[iGXsplitFF]= sectorShower ? new GXSplitFFsec() : new GXSplitFF();
3390  // Add RF antenna functions (no sector versions defined yet)
3391  antFunPtrs[iQQemitRF] = new QQEmitRF();
3392  antFunPtrs[iQGemitRF] = new QGEmitRF();
3393  antFunPtrs[iXGsplitRF]= new XGSplitRF();
3394  if (verbose >= quiteloud)
3395  printOut(__METHOD_NAME__, "Defined new antFunPtrs");
3396 
3397  // Loop through antFunPtrs and initialize.
3398  for (map<int,AntennaFunction*>::iterator antFunPtrIt = antFunPtrs.begin();
3399  antFunPtrIt != antFunPtrs.end(); ++antFunPtrIt) {
3400 
3401  // Initialize antenna pointers, antenna.
3402  AntennaFunction* antFunPtr = antFunPtrIt->second;
3403  antFunPtr->initPtr(infoPtr, dglapPtr);
3404  bool pass = antFunPtr->init();
3405 
3406  // Optionally check antenna (singularities, positivity, etc).
3407  if (settingsPtr->flag("Vincia:checkAntennae"))
3408  pass = pass && antFunPtr->check();
3409 
3410  // Everything OK?
3411  if (pass) {
3412  if (verbose > normal) printOut(__METHOD_NAME__, "Added to antenna list: "
3413  + antFunPtr->humanName());
3414  } else if (verbose >= quiet) infoPtr->errorMsg("Warning in "+
3415  __METHOD_NAME__+": one or more consistency checks failed.");
3416  }
3417  isInit = true;
3418 
3419 }
3420 
3421 //--------------------------------------------------------------------------
3422 
3423 // Method to return all iAntPhys values that are defined in antFunPtr.
3424 
3425 vector<int> AntennaSetFSR::getIant() {
3426 
3427  vector<int> iAnt;
3428  map<int,AntennaFunction*>::iterator it;
3429  for (it = antFunPtrs.begin(); it != antFunPtrs.end(); ++it)
3430  iAnt.push_back(it->first);
3431  return iAnt;
3432 
3433 }
3434 
3435 //==========================================================================
3436 
3437 // The AntennaSetISR class. Simple container of II and IF antenna functions.
3438 
3439 //--------------------------------------------------------------------------
3440 
3441 // Initialize pointers.
3442 
3443 void AntennaSetISR::initPtr(Info* infoPtrIn, DGLAP* dglapPtrIn) {
3444 
3445  infoPtr = infoPtrIn;
3446  particleDataPtr = infoPtr->particleDataPtr;
3447  settingsPtr = infoPtr->settingsPtr;
3448  rndmPtr = infoPtr->rndmPtr;
3449  dglapPtr = dglapPtrIn;
3450  isInitPtr = true;
3451 
3452 }
3453 
3454 //--------------------------------------------------------------------------
3455 
3456 // Initialize antenna set.
3457 
3458 void AntennaSetISR::init() {
3459 
3460  // Check pointers and initialization.
3461  if (!isInitPtr) {
3462  printOut(__METHOD_NAME__, "Cannot initialize, pointers not set.");
3463  return;
3464  }
3465  verbose = settingsPtr->mode("Vincia:verbose");
3466  if (isInit) {
3467  if (verbose >= debug)
3468  printOut(__METHOD_NAME__, "Already initialized antenna set.");
3469  return;
3470  }
3471 
3472  // Create antenna objects.
3473  bool sectorShower = settingsPtr->flag("Vincia:sectorShower");
3474  antFunPtrs[iQQemitII] = new QQEmitII();
3475  antFunPtrs[iGQemitII] = new GQEmitII();
3476  antFunPtrs[iGGemitII] = new GGEmitII();
3477  antFunPtrs[iQXsplitII] = new QXSplitII();
3478  antFunPtrs[iGXconvII] = new GXConvII();
3479  antFunPtrs[iQQemitIF] = new QQEmitIF();
3480  antFunPtrs[iQGemitIF] = sectorShower ? new QGEmitIFsec() : new QGEmitIF();
3481  antFunPtrs[iGQemitIF] = new GQEmitIF();
3482  antFunPtrs[iGGemitIF] = sectorShower ? new GGEmitIFsec() : new GGEmitIF();
3483  antFunPtrs[iQXsplitIF] = new QXSplitIF();
3484  antFunPtrs[iGXconvIF] = new GXConvIF();
3485  antFunPtrs[iXGsplitIF] = sectorShower ? new XGSplitIFsec() : new XGSplitIF();
3486 
3487  // Loop through antFunPtrs and initialize.
3488  for (map<int,AntennaFunctionIX*>::iterator antFunPtrIt = antFunPtrs.begin();
3489  antFunPtrIt != antFunPtrs.end(); ++antFunPtrIt) {
3490 
3491  // Initialize antenna pointers, antenna.
3492  AntennaFunction* antFunPtr = antFunPtrIt->second;
3493  antFunPtr->initPtr(infoPtr, dglapPtr);
3494  bool pass = antFunPtr->init();
3495 
3496  // Optionally check antenna (singularities, positivity, etc.).
3497  if (settingsPtr->flag("Vincia:checkAntennae"))
3498  pass = pass && antFunPtr->check();
3499 
3500  // Debug info.
3501  if (pass) {
3502  if (verbose > normal) printOut(__METHOD_NAME__, "Added to antenna list: "
3503  + antFunPtr->vinciaName());
3504  } else if (verbose >= quiet)
3505  infoPtr->errorMsg("Warning in "+__METHOD_NAME__
3506  +": one or more consistency checks failed.");
3507  }
3508  isInit = true;
3509 
3510 }
3511 
3512 //--------------------------------------------------------------------------
3513 
3514 // Method to return all iAntPhys values that are defined in antFunPtr.
3515 
3516 vector<int> AntennaSetISR::getIant() {
3517 
3518  vector<int> iAnt;
3519  map<int,AntennaFunctionIX*>::iterator it;
3520  for (it = antFunPtrs.begin(); it != antFunPtrs.end(); ++it)
3521  iAnt.push_back(it->first);
3522  return iAnt;
3523 
3524 }
3525 
3526 //==========================================================================
3527 
3528 // Class MECs, for computing matrix-element corrections to antenna
3529 // functions.
3530 
3531 //--------------------------------------------------------------------------
3532 
3533 // Initialize pointers.
3534 
3535 void MECs::initPtr(Info* infoPtrIn, ShowerMEs* mg5mesPtrIn,
3536  VinciaCommon* vinComPtrIn) {
3537 
3538  infoPtr = infoPtrIn;
3539  settingsPtr = infoPtr->settingsPtr;
3540  rndmPtr = infoPtr->rndmPtr;
3541  partonSystemsPtr = infoPtr->partonSystemsPtr;
3542  mg5mesPtr = mg5mesPtrIn;
3543  vinComPtr = vinComPtrIn;
3544  isInitPtr = true;
3545 
3546 }
3547 
3548 //--------------------------------------------------------------------------
3549 
3550 // Initialize.
3551 
3552 void MECs::init() {
3553 
3554  // MEC settings.
3555  verbose = settingsPtr->mode("Vincia:verbose");
3556  maxMECs2to1 = settingsPtr->mode("Vincia:maxMECs2to1");
3557  maxMECs2to2 = settingsPtr->mode("Vincia:maxMECs2to2");
3558  maxMECs2toN = settingsPtr->mode("Vincia:maxMECs2toN");
3559  maxMECsResDec = settingsPtr->mode("Vincia:maxMECsResDec");
3560  maxMECsMPI = settingsPtr->mode("Vincia:maxMECsMPI");
3561  matchingFullColour = settingsPtr->flag("Vincia:matchingFullColour");
3562  nFlavZeroMass = settingsPtr->mode("Vincia:nFlavZeroMass");
3563  if (maxMECs2to2 == 0) maxMECsMPI = 0;
3564  sizeOutBornSav.clear();
3565 
3566  // Initialise MG5 interface
3567  if (mg5mesPtr->initVincia()) {
3568  mg5mesPtr->setColourDepthVincia(matchingFullColour);
3569  } else {
3570  if (verbose >= 3) printOut("Vincia::MECs",
3571  "Could not initialise ShowerMEs interface.");
3572  maxMECs2to1 = -1;
3573  maxMECs2to2 = -1;
3574  maxMECs2toN = -1;
3575  maxMECsResDec = -1;
3576  maxMECsMPI = -1;
3577  }
3578  isInit = true;
3579 
3580 }
3581 
3582 //--------------------------------------------------------------------------
3583 
3584 // Function to return ME class (Born type) for a parton configuration.
3585 
3586 bool MECs::prepare(const int iSys, Event& event) {
3587 
3588  // Initialise for no MECs, then check if MECs should be applied.
3589  int nAll = partonSystemsPtr->sizeAll(iSys);
3590  int nOut = partonSystemsPtr->sizeOut(iSys);
3591  bool isHard = (iSys == 0 && nAll - nOut == 2);
3592  bool isMPI = (iSys >= 1 && nAll - nOut == 2);
3593  sizeOutBornSav[iSys] = nOut;
3594 
3595  // Check if MECs are switched on for this type of Born system.
3596  if (isHard) {
3597  if (nOut == 1 && maxMECs2to1 < 0) return false;
3598  else if (nOut == 2 && maxMECs2to2 < 0) return false;
3599  else if (nOut >= 3 && maxMECs2toN < 0) return false;
3600  } else if (isMPI) {
3601  if (infoPtr->nMPI() > maxMECsMPI) return false;
3602  } else {
3603  if (maxMECsResDec < 0) return false;
3604  }
3605 
3606  // Make vectors of ID codes.
3607  vector<int> idIn, idOut;
3608  if (partonSystemsPtr->hasInAB(iSys)) {
3609  idIn.push_back(event[partonSystemsPtr->getInA(iSys)].id());
3610  idIn.push_back(event[partonSystemsPtr->getInB(iSys)].id());
3611  } else if (partonSystemsPtr->hasInRes(iSys))
3612  idIn.push_back(event[partonSystemsPtr->getInRes(iSys)].id());
3613  for (int i = 0; i < partonSystemsPtr->sizeOut(iSys); ++i)
3614  idOut.push_back(event[partonSystemsPtr->getOut(iSys, i)].id());
3615 
3616  // Check whether MG5MEs interface has the process.
3617  set<int> sChan;
3618  return mg5mesPtr->hasProcessVincia(idIn, idOut, sChan);
3619 
3620 }
3621 
3622 //--------------------------------------------------------------------------
3623 
3624 // Function to assign helicities to particles (using MEs).
3625 
3626 bool MECs::polarise(const int iSys, Event& event) {
3627 
3628  // First check if we should be doing anything at all.
3629  if (partonSystemsPtr->hasInAB(iSys)) {
3630  // Hard System.
3631  if (iSys == 0) {
3632  int nOut = partonSystemsPtr->sizeOut(iSys);
3633  if (nOut==1 && maxMECs2to1 <= -1) return false;
3634  if (nOut==2 && maxMECs2to2 <= -1) return false;
3635  if (nOut>=3 && maxMECs2toN <= -1) return false;
3636  // Hardest MPI system
3637  } else if (iSys == 1 && maxMECsMPI <= -1) return false;
3638  // Further MPI systems just use unpolarised showers.
3639  else return false;
3640  // Resonance-Decay systems
3641  } else {if (maxMECsResDec <= -1) return false;}
3642 
3643  // If state does not already have one or more assigned helicities,
3644  // check if we can polarise it.
3645  bool checkIncoming = true;
3646  if (!isPolarised(iSys, event, checkIncoming)) {
3647 
3648  // Extract particles to use for ME evaluations (incoming first)/
3649  vector<Particle> parts = makeParticleList(iSys, event);
3650  if (parts.size() <= 2) return false;
3651  int nIn = partonSystemsPtr->hasInRes(iSys) ? 1 : 2;
3652 
3653  // Verbose output.
3654  if (verbose >= 4) cout << " MECs::polarise(): system " << iSys << " nIn = "
3655  << nIn << endl;
3656  // Check if MG5MEs interface can do this.
3657  if (!mg5mesPtr->selectHelicitiesVincia(parts, nIn)) return false;
3658 
3659  // Update particles in event record: incoming.
3660  if (nIn == 1) event[partonSystemsPtr->getInRes(iSys)].pol(parts[0].pol());
3661  else {
3662  event[partonSystemsPtr->getInA(iSys)].pol(parts[0].pol());
3663  event[partonSystemsPtr->getInB(iSys)].pol(parts[1].pol());
3664  }
3665 
3666  // Update particles in event record: outgoing.
3667  for (int iOut = 0; iOut < partonSystemsPtr->sizeOut(iSys); ++iOut) {
3668  int iEvent = partonSystemsPtr->getOut(iSys,iOut);
3669  int iParts = iOut + nIn;
3670  event[iEvent].pol(parts[iParts].pol());
3671  }
3672  }
3673 
3674  // Verbose output (showing polarisations).
3675  if (verbose >= 9) event.list(true);
3676 
3677  // All is well: return true if any particles remain with assigned helicities.
3678  return isPolarised(iSys,event,true);
3679 
3680 }
3681 
3682 //--------------------------------------------------------------------------
3683 
3684 // Make list of particles as vector<Particle>.
3685 
3686 vector<Particle> MECs::makeParticleList(const int iSys, const Event& event,
3687  const vector<Particle> pNew, const vector<int> iOld) {
3688 
3689  // Put incoming ones (initial-state partons or decaying resonance) first.
3690  vector<Particle> state;
3691  if (partonSystemsPtr->hasInAB(iSys)) {
3692  int iA = partonSystemsPtr->getInA(iSys);
3693  int iB = partonSystemsPtr->getInB(iSys);
3694  for (int j = 0; j < (int)iOld.size(); ++j) {
3695  // Exclude any partons in old state that should be replaced.
3696  if (iOld[j] == iA) iA = -1;
3697  if (iOld[j] == iB) iB = -1;
3698  }
3699  if (iA >= 0) state.push_back(event[iA]);
3700  if (iB >= 0) state.push_back(event[iB]);
3701  } else if (partonSystemsPtr->hasInRes(iSys)) {
3702  int iRes = partonSystemsPtr->getInRes(iSys);
3703  for (int j = 0; j < (int)iOld.size(); ++j) {
3704  // Exclude any partons in old state that should be replaced.
3705  if (iOld[j] == iRes) iRes = -1;
3706  }
3707  if (iRes >= 0) state.push_back(event[iRes]);
3708  }
3709  // Add any post-branching incoming particles.
3710  for (int j = 0; j < (int)pNew.size(); ++j) {
3711  if (!pNew[j].isFinal()) state.push_back(pNew[j]);
3712  }
3713 
3714  // Sanity check; must have at least 1 incoming particle.
3715  if (state.size() == 0) {
3716  if (verbose >= 5) cout << "Vincia::MECs::makeParticleList(): problem: "
3717  "could not identify incoming or decaying particle." << endl;
3718  if (verbose >= 9) event.list();
3719  return state;
3720  }
3721 
3722  // Then put outgoing ones.
3723  for (int i = 0; i < partonSystemsPtr->sizeOut(iSys); ++i) {
3724  int i1 = partonSystemsPtr->getOut(iSys, i);
3725  // Do not add any that are marked as branched.
3726  for (int j = 0; j < (int)iOld.size(); ++j) {if (iOld[j] == i1) i1 = -1;}
3727  if (i1 >= 0) state.push_back(event[i1]);
3728  }
3729  // Add any post-branching outgoing partons.
3730  for (int j=0; j<(int)pNew.size(); ++j)
3731  if (pNew[j].isFinal()) state.push_back(pNew[j]);
3732 
3733  // Return the state.
3734  return state;
3735 
3736 }
3737 
3738 //--------------------------------------------------------------------------
3739 
3740 // Check if state already has helicities.
3741 
3742 bool MECs::isPolarised(int iSys, Event& event, bool checkIncoming) {
3743 
3744  for (int i = 0; i < partonSystemsPtr->sizeAll(iSys); ++i) {
3745  int ip = partonSystemsPtr->getAll(iSys,i);
3746  if (ip == 0) continue;
3747  if (event[ip].isFinal() || checkIncoming)
3748  if (event[ip].pol() != 9) return true;
3749  }
3750  return false;
3751 
3752 }
3753 
3754 //--------------------------------------------------------------------------
3755 
3756 // Function to determine if MECs are requested at this order for this system.
3757 
3758 bool MECs::doMEC(int iSys, int nBranch) {
3759 
3760  // MECs in resonance decays.
3761  bool isResDec = partonSystemsPtr->hasInRes(iSys);
3762  if (isResDec) {if (nBranch <= maxMECsResDec) return true;
3763 
3764  // MECs in Scattering Processes.
3765  } else {
3766  // Hard processes.
3767  if (iSys == 0) {
3768  if (sizeOutBorn(iSys) == 1 && nBranch <= maxMECs2to1) return true;
3769  if (sizeOutBorn(iSys) == 2 && nBranch <= maxMECs2to2) return true;
3770  if (sizeOutBorn(iSys) >= 3 && nBranch <= maxMECs2toN) return true;
3771 
3772  // MPI.
3773  } else if (iSys == 1) {if (nBranch <= maxMECsMPI) return true;}
3774  }
3775 
3776  // If nobody said yes by now, return the sad news.
3777  return false;
3778 
3779 }
3780 
3781 //--------------------------------------------------------------------------
3782 
3783 // Get squared matrix element.
3784 
3785 double MECs::getME2(const vector<Particle> state, int nIn) {
3786  return mg5mesPtr->me2Vincia(state, nIn);}
3787 
3788 double MECs::getME2(const int iSys, const Event& event) {
3789  vector<Particle> state = makeParticleList(iSys, event);
3790  bool isResDec = partonSystemsPtr->hasInRes(iSys);
3791  return (isResDec) ?
3792  mg5mesPtr->me2Vincia(state, 1 ) : mg5mesPtr->me2Vincia(state, 2);
3793 }
3794 
3795 //--------------------------------------------------------------------------
3796 
3797 // Print header information.
3798 
3799 void MECs::header() {
3800 
3801  // Front matter.
3802  bool doPolarise = (maxMECs2to1 >= 0) || (maxMECs2to2 >= 0)
3803  || (maxMECs2toN >= 0) || (maxMECsResDec >= 0);
3804  bool doMecs = (maxMECs2to1 >= 1) || (maxMECs2to2 >= 1)
3805  || (maxMECs2toN >= 1) || (maxMECsResDec >= 1);
3806  cout << " |\n | MECs (-1:off, 0:selectHelicities, >=1:nMECs): ";
3807  if (doPolarise) {
3808  cout <<"\n | maxMECs2to1 = "
3809  << num2str(maxMECs2to1, 9) << "\n"
3810  << " | maxMECs2to2 = "
3811  << num2str(maxMECs2to2, 9) << "\n"
3812  <<" | maxMECs2toN = "
3813  << num2str(maxMECs2toN, 9) << "\n"
3814  <<" | maxMECsResDec = "
3815  << num2str(maxMECsResDec, 9) <<"\n";
3816 
3817  // Setings.
3818  if (doMecs) {
3819  cout << " | matchingFullColour = "
3820  << bool2str(matchingFullColour, 9) << "\n";
3821  int matchingRegOrder = settingsPtr->mode("Vincia:matchingRegOrder");
3822  int matchingRegShape = settingsPtr->mode("Vincia:matchingRegShape");
3823  double matchingScale = settingsPtr->parm("Vincia:matchingRegScale");
3824  bool matchingScaleIsAbs =
3825  settingsPtr->flag("Vincia:matchingRegScaleIsAbsolute");
3826  double matchingScaleRatio =
3827  settingsPtr->parm("Vincia:matchingRegScaleRatio");
3828  double matchingIRcutoff = settingsPtr->parm("Vincia:matchingIRcutoff");
3829  cout << " | regOrder = "
3830  << num2str(matchingRegOrder, 9) << endl;
3831  if (matchingScaleIsAbs)
3832  cout << " | regScale = "
3833  << num2str(matchingScale, 9) << endl;
3834  else
3835  cout << " | regScaleRatio = "
3836  << num2str(matchingScaleRatio, 9) << endl;
3837  if (verbose >= 2)
3838  cout << " | regShape = "
3839  << num2str(matchingRegShape, 9) << endl;
3840  cout << " | IR cutoff = "
3841  << num2str(matchingIRcutoff, 9) << endl;
3842  }
3843 
3844  // Print MG5MEs reference.
3845  cout << " | The MADGRAPH Matrix Element interface relies on:" << endl
3846  << " | MADGRAPH 5 : Alwall et al., JHEP06(2011)128, "
3847  << "arXiv:1106.0522 " << endl;
3848  }
3849  else cout << bool2str(false, 9) << "\n";
3850 
3851 }
3852 
3853 //==========================================================================
3854 
3855 } // end namespace Pythia8
Event
Definition: AgUStep.h:26
Generated by   doxygen 1.8.5