PYTHIA  8.303
VinciaISR.h
1 // VinciaISR.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 Peter Skands, Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // This file contains header information for the VinciaISR class for
7 // QCD initial-state antenna showers (II and IF), and auxiliary classes.
8 
9 #ifndef Pythia_VinciaISR_H
10 #define Pythia_VinciaISR_H
11 
12 #include "Pythia8/SpaceShower.h"
13 #include "Pythia8/VinciaAntennaFunctions.h"
14 #include "Pythia8/VinciaCommon.h"
15 #include "Pythia8/VinciaQED.h"
16 #include "Pythia8/VinciaWeights.h"
17 
18 namespace Pythia8 {
19 
20 // Forward declarations.
21 class VinciaFSR;
22 
23 //==========================================================================
24 
25 // Base class for initial-state trial generators. Note, base class is
26 // coded for a soft-eikonal trial function.
27 
29 
30 public:
31 
32  // Constructor.
33  TrialGeneratorISR() : isInit(false) {;}
34  virtual ~TrialGeneratorISR() {;}
35 
36  // Initialize pointers.
37  void initPtr(Info* infoPtrIn);
38 
39  // Name of trial generator.
40  virtual string name() {return "TrialGeneratorISR";}
41 
42  // Initialize.
43  virtual void init(double mcIn, double mbIn);
44 
45  // Trial antenna function. Convention for what is coded here:
46  // when using x*PDF ratio <= const : aTrial
47  // when using PDF ratio <= const : aTrial * sab / sAB
48  // Base class implements soft eikonal with PDF ratio as overestimate.
49  virtual double aTrial(double saj, double sjb, double sAB);
50 
51  // Evolution scale.
52  virtual double getQ2(double saj, double sjb, double sAB) {
53  return (saj*sjb/(saj + sjb + sAB));}
54  virtual double getQ2max(double sAB, double, double) {
55  return (0.25*pow2(shhSav - sAB)/shhSav);}
56 
57  // Generate new Q value, with first-order running alphaS.
58  virtual double genQ2run(double q2old, double sAB, double zMin, double zMax,
59  double colFac, double PDFratio, double b0, double kR, double Lambda,
60  double eA, double eB, double headroomFac = 1.0, double enhanceFac = 1.0);
61 
62  // Generate new Q value, with constant trial alphaS.
63  virtual double genQ2(double q2old, double sAB, double zMin, double zMax,
64  double colFac, double alphaSvalue, double PDFratio, double eA, double eB,
65  double headroomFac = 1.0, double enhanceFac = 1.0);
66 
67  // Generate new Q value, with running of the PDFs towards the mass
68  // threshold.
69  virtual double genQ2thres(double q2old, double sAB, double zMin,
70  double zMax, double colFac, double alphaSvalue, double PDFratio, int idA,
71  int idB, double eA, double eB, bool useMpdf, double headroomFac = 1.0,
72  double enhanceFac = 1.0);
73 
74  // Generate a new zeta value in [zMin,zMax].
75  virtual double genZ(double zMin, double zMax);
76 
77  // The zeta integral.
78  virtual double getIz(double zMin, double zMax);
79 
80  // The zeta boundaries, for a given value of the evolution scale.
81  virtual double getZmin(double Qt2, double sAB, double eA, double eBeamUsed);
82  virtual double getZmax(double Qt2, double sAB, double eA, double eBeamUsed);
83 
84  // Inverse transforms to obtain saj and sjb from Qt2 and zeta.
85  virtual double getS1j(double Qt2, double zeta, double sAB);
86  virtual double getSj2(double Qt2, double zeta, double sAB);
87 
88  // Compute trial PDF ratio.
89  virtual double trialPDFratio(BeamParticle* beamAPtr, BeamParticle* beamBPtr,
90  int iSys, int idA, int idB, double eA, double eB,
91  double Qt2A, double Qt2B);
92 
93  // Return last trial PDF ratio.
94  virtual double getTrialPDFratio() {return trialPDFratioSav;}
95 
96  // Check initialization.
97  bool checkInit();
98 
99  // Return the last trial flavor.
100  int trialFlav() {return trialFlavSav;}
101 
102  protected:
103 
104  // Pointers.
106  Rndm* rndmPtr;
107  Settings* settingsPtr;
108 
109  // Use m or pT evolution for collinear singularities.
111 
112  // s for hadron hadron.
113  double shhSav;
114 
115  // For conversion trial generators.
117  int nGtoQISRSav;
118 
119  // Masses.
120  double mbSav;
121  double mcSav;
122 
123  // Saved trial PDF ratio and trial tolerance.
125  double TINYPDFtrial;
126 
127  private:
128 
129  // Status.
130  bool isInit;
131 
132  // Verbosity level.
133  int verbose;
134 
135 };
136 
137 //==========================================================================
138 
139 // Soft-eikonal trial function for II (clone of base class but with
140 // name change).
141 
143 
144 public:
145 
146  // Name of trial generator.
147  virtual string name() override {return "TrialIISoft";}
148 
149 };
150 
151 //==========================================================================
152 
153 // A collinear trial function for initial-initial.
154 
156 
157 public:
158 
159  // Name of trial generator.
160  virtual string name() override {return "TrialIIGCollA";}
161 
162  // Trial antenna function for g->gg (collinear to beam A).
163  // Used with x*PDF <= const, so no extra x-factor.
164  virtual double aTrial(double saj, double sjb, double sAB) override;
165 
166  // Evolution scale
167  virtual double getQ2(double saj, double sjb, double sAB) override {
168  return (saj*sjb/(saj + sjb + sAB));}
169  virtual double getQ2max(double sAB, double, double) override {
170  return (0.25*pow2(shhSav - sAB)/shhSav);}
171 
172  // Generate a new Q value, with first-order running alphaS.
173  virtual double genQ2run(double q2old, double sAB, double zMin, double zMax,
174  double colFac, double PDFratio, double b0, double kR, double Lambda,
175  double eA, double eB, double headroomFac = 1.0, double enhanceFac = 1.0)
176  override;
177 
178  // Generate a new Q value, with constant trial alphaS.
179  virtual double genQ2(double q2old, double sAB, double zMin, double zMax,
180  double colFac, double alphaSvalue, double PDFratio,
181  double eA, double eB, double headroomFac = 1.0, double enhanceFac = 1.0)
182  override;
183 
184  // Generate a new zeta value in [zMin,zMax].
185  virtual double genZ(double zMin, double zMax) override;
186 
187  // The zeta integral.
188  virtual double getIz(double zMin, double zMax) override;
189 
190  // The zeta boundaries, for a given value of the evolution scale
191  virtual double getZmin(double Qt2, double sAB, double eA, double eBeamUsed)
192  override;
193  virtual double getZmax(double Qt2, double sAB, double eA, double eBeamUsed)
194  override;
195 
196  // Inverse transforms to obtain saj and sjb from Qt2 and zeta.
197  virtual double getS1j(double Qt2, double zeta, double sAB) override;
198  virtual double getSj2(double Qt2, double zeta, double sAB) override;
199 
200 };
201 
202 //==========================================================================
203 
204 // B collinear trial function for initial-initial.
205 
206 class TrialIIGCollB : public TrialIIGCollA {
207 
208 public:
209 
210  // Name of trial generator.
211  virtual string name() override {return "TrialIIGCollB";}
212 
213  // Trial antenna function.
214  virtual double aTrial(double saj, double sjb, double sAB) override {
215  return TrialIIGCollA::aTrial(sjb, saj, sAB);}
216 
217  // Inverse transforms to obtain saj and sjb from Qt2 and zeta.
218  virtual double getS1j(double Qt2, double zeta, double sAB) override {
219  return TrialIIGCollA::getSj2(Qt2, zeta, sAB);}
220  virtual double getSj2(double Qt2, double zeta, double sAB) override {
221  return TrialIIGCollA::getS1j(Qt2, zeta, sAB);}
222 
223  // Trial PDF ratio.
224  virtual double trialPDFratio(BeamParticle* beamAPtr, BeamParticle* beamBPtr,
225  int iSys, int idA, int idB, double eA, double eB,
226  double Qt2A, double Qt2B) override {return TrialIIGCollA::trialPDFratio(
227  beamBPtr, beamAPtr, iSys, idB, idA, eB, eA, Qt2B, Qt2A);}
228 
229 };
230 
231 //==========================================================================
232 
233 // A splitting trial function for initial-initial, q -> gqbar.
234 
236 
237 public:
238 
239  // Name of trial generator.
240  virtual string name() override {return "TrialIISplitA";}
241 
242  // Trial antenna function. This trial used with PDF ratio <= const,
243  // so has sab/sAB prefactor.
244  virtual double aTrial(double saj, double sjb, double sAB) override;
245 
246  // Evolution scale.
247  virtual double getQ2(double saj, double sjb, double sAB) override {
248  return ((useMevolSav) ? saj : saj*sjb/(saj + sjb + sAB));}
249  virtual double getQ2max(double sAB, double, double) override {
250  return ((useMevolSav) ? shhSav-sAB : 0.25*pow2(shhSav-sAB)/shhSav);}
251 
252  // Generate new Q value, with first-order running alphaS. Same
253  // expression for QT and QA; Iz is different.
254  virtual double genQ2run(double q2old, double sAB, double zMin, double zMax,
255  double colFac, double PDFratio, double b0, double kR, double Lambda,
256  double eA, double eB, double headroomFac = 1.0, double enhanceFac = 1.0)
257  override;
258 
259  // Generate new Q value, with constant trial alphaS.
260  virtual double genQ2(double q2old, double sAB, double zMin, double zMax,
261  double colFac, double alphaSvalue, double PDFratio,
262  double eA, double eB, double headroomFac = 1.0, double enhanceFac = 1.0)
263  override;
264 
265  // Generate new Q value, with running of the PDFs towards the mass
266  // threshold.
267  virtual double genQ2thres(double q2old, double sAB,
268  double zMin, double zMax, double colFac, double alphaSvalue,
269  double PDFratio, int idA, int idB, double eA, double eB, bool useMpdf,
270  double headroomFac = 1.0, double enhanceFac = 1.0) override;
271 
272  // Generate a new zeta value in [zMin,zMax]. For QT 1/(1+z3)
273  // distribution; for QA 1/z4 distribution.
274  virtual double genZ(double zMin, double zMax) override;
275 
276  // The zeta integral (with alpha = 0). For QT integral dz4 1/(1+z3),
277  // for QA integral dz4 1/z4.
278  virtual double getIz(double zMin, double zMax) override;
279 
280  // The zeta boundaries, for a given value of the evolution scale.
281  virtual double getZmin(double Qt2, double sAB, double eA, double eBeamUsed)
282  override;
283  virtual double getZmax(double Qt2, double sAB, double eA, double eBeamUsed)
284  override;
285 
286  // Inverse transforms to obtain saj and sjb from Qt2 and zeta.
287  virtual double getS1j(double Qt2, double zeta, double sAB) override;
288  virtual double getSj2(double Qt2, double zeta, double sAB) override;
289 
290  // Trial PDF ratio.
291  virtual double trialPDFratio(BeamParticle* beamAPtr, BeamParticle* beamBPtr,
292  int iSys, int idA, int idB, double eA, double eB,
293  double Qt2A, double Qt2B) override;
294 
295 };
296 
297 //==========================================================================
298 
299 // B splitting trial function for initial-initial, q -> gqbar.
300 
301 class TrialIISplitB : public TrialIISplitA {
302 
303 public:
304 
305  // Name of trial generator.
306  virtual string name() override {return "TrialIISplitB";}
307 
308  // Trial antenna function.
309  virtual double aTrial(double saj, double sjb, double sAB) override {
310  return TrialIISplitA::aTrial(sjb, saj, sAB);}
311 
312  // Evolution scale.
313  virtual double getQ2(double saj, double sjb, double sAB) override {
314  return TrialIISplitA::getQ2(sjb, saj, sAB);}
315 
316  // Generate new Q value, with first-order running alphaS.
317  virtual double genQ2run(double q2old, double sAB, double zMin, double zMax,
318  double colFac, double PDFratio, double b0, double kR, double Lambda,
319  double eA, double eB, double headroomFac = 1.0, double enhanceFac = 1.0)
320  override {return TrialIISplitA::genQ2run(q2old, sAB, zMin, zMax, colFac,
321  PDFratio, b0, kR, Lambda, eB, eA, headroomFac, enhanceFac);}
322 
323  // Generate new Q value, with constant trial alphaS.
324  virtual double genQ2(double q2old, double sAB, double zMin, double zMax,
325  double colFac, double alphaSvalue, double PDFratio,
326  double eA, double eB, double headroomFac = 1.0, double enhanceFac = 1.0)
327  override {return TrialIISplitA::genQ2(q2old, sAB, zMin, zMax, colFac,
328  alphaSvalue, PDFratio, eB, eA, headroomFac, enhanceFac);}
329 
330  // Generate new Q value, with running of the PDFs towards the mass
331  // threshold.
332  virtual double genQ2thres(double q2old, double sAB,
333  double zMin, double zMax, double colFac, double alphaSvalue,
334  double PDFratio, int idA, int idB, double eA, double eB, bool useMpdf,
335  double headroomFac = 1.0, double enhanceFac = 1.0) override {
336  return TrialIISplitA::genQ2thres(q2old, sAB, zMin, zMax, colFac,
337  alphaSvalue, PDFratio, idB, idA, eB, eA, useMpdf, headroomFac,
338  enhanceFac);}
339 
340  // Inverse transforms to obtain saj and sjb from Qt2 and zeta.
341  virtual double getS1j(double Qt2, double zeta, double sAB) override {
342  return TrialIISplitA::getSj2(Qt2, zeta, sAB);}
343  virtual double getSj2(double Qt2, double zeta, double sAB) override {
344  return TrialIISplitA::getS1j(Qt2, zeta, sAB);}
345 
346  // Trial PDF ratio.
347  virtual double trialPDFratio(BeamParticle* beamAPtr, BeamParticle* beamBPtr,
348  int iSys, int idA, int idB, double eA, double eB, double Qt2A, double Qt2B)
349  override {return TrialIISplitA::trialPDFratio(beamBPtr, beamAPtr, iSys,
350  idB, idA, eB, eA, Qt2B, Qt2A);}
351 
352 };
353 
354 //==========================================================================
355 
356 // A conversion trial function for initial-initial, g -> qqbar.
357 
359 
360 public:
361 
362  // Name of trial generator
363  virtual string name() override {return "TrialIIConvA";}
364 
365  // Trial antenna function. Used with x*PDF ratio <= const, so no
366  // extra prefactor.
367  virtual double aTrial(double saj, double sjb, double sAB) override;
368 
369  // Evolution scale.
370  virtual double getQ2(double saj, double sjb, double sAB) override {
371  return ((useMevolSav) ? saj : (saj*sjb/(saj + sjb + sAB)));}
372  virtual double getQ2max(double sAB, double, double) override {
373  return ((useMevolSav) ? (shhSav - sAB) : 0.25*pow2(shhSav - sAB)/shhSav);}
374 
375  // Generate a new Q value, with first-order running alphaS.
376  virtual double genQ2run(double q2old, double sAB, double zMin, double zMax,
377  double colFac, double PDFratio, double b0, double kR, double Lambda,
378  double eA, double eB, double headroomFac = 1.0, double enhanceFac = 1.0)
379  override;
380 
381  // Generate a new Q value, with constant trial alphaS. Same
382  // expression for QT and QA; Iz is different.
383  virtual double genQ2(double q2old, double sAB, double zMin, double zMax,
384  double colFac, double alphaSvalue, double PDFratio,
385  double eA, double eB, double headroomFac = 1.0, double enhanceFac = 1.0)
386  override;
387 
388  // Generate a new zeta value in [zMin,zMax].
389  virtual double genZ(double zMin, double zMax) override;
390 
391  // The zeta integral.
392  virtual double getIz(double zMin, double zMax) override;
393 
394  // The zeta boundaries, for a given value of the evolution scale.
395  virtual double getZmin(double Qt2, double sAB, double eA, double eBeamUsed)
396  override;
397  virtual double getZmax(double Qt2, double sAB, double eA, double eBeamUsed)
398  override;
399 
400  // Inverse transforms to obtain saj and sjb from Qt2 and zeta.
401  virtual double getS1j(double Qt2, double zeta, double sAB) override;
402  virtual double getSj2(double Qt2, double zeta, double sAB) override;
403 
404  // Trial PDF ratio.
405  virtual double trialPDFratio(BeamParticle* beamAPtr, BeamParticle* beamBPtr,
406  int iSys, int idA, int idB, double eA, double eB,
407  double Qt2A, double Qt2B) override;
408 
409 };
410 
411 //==========================================================================
412 
413 // B conversion trial function for initial-initial, g -> qqbar.
414 
415 class TrialIIConvB : public TrialIIConvA {
416 
417 public:
418 
419  // Name of trial generator.
420  virtual string name() override {return "TrialIIConvB";}
421 
422  // Trial antenna function.
423  virtual double aTrial(double saj, double sjb, double sAB) override {
424  return TrialIIConvA::aTrial(sjb, saj, sAB);}
425 
426  // Evolution scale.
427  virtual double getQ2(double saj, double sjb, double sAB) override {
428  return TrialIIConvA::getQ2(sjb, saj, sAB);}
429 
430  // Generate a new Q value, with first-order running alphaS
431  virtual double genQ2run(double q2old, double sAB, double zMin, double zMax,
432  double colFac, double PDFratio, double b0, double kR, double Lambda,
433  double eA, double eB, double headroomFac = 1.0, double enhanceFac = 1.0)
434  override {return TrialIIConvA::genQ2run(q2old, sAB, zMin, zMax, colFac,
435  PDFratio, b0, kR, Lambda, eB, eA, headroomFac, enhanceFac);}
436 
437  // Generate a new Q value, with constant trial alphaS.
438  virtual double genQ2(double q2old, double sAB, double zMin, double zMax,
439  double colFac, double alphaSvalue, double PDFratio,
440  double eA, double eB, double headroomFac = 1.0, double enhanceFac = 1.0)
441  override {return TrialIIConvA::genQ2(q2old, sAB, zMin, zMax, colFac,
442  alphaSvalue, PDFratio, eB, eA, headroomFac, enhanceFac);}
443 
444  // Inverse transforms: to obtain saj and sjb from Qt2 and zeta.
445  virtual double getS1j(double Qt2, double zeta, double sAB) override {
446  return TrialIIConvA::getSj2(Qt2, zeta, sAB);}
447  virtual double getSj2(double Qt2, double zeta, double sAB) override {
448  return TrialIIConvA::getS1j(Qt2, zeta, sAB);}
449 
450  // Trial PDF ratio
451  virtual double trialPDFratio(BeamParticle* beamAPtr, BeamParticle* beamBPtr,
452  int iSys, int idA, int idB, double eA, double eB,
453  double Qt2A, double Qt2B) override {return TrialIIConvA::trialPDFratio(
454  beamBPtr, beamAPtr, iSys, idB, idA, eB, eA, Qt2B, Qt2A);}
455 
456 };
457 
458 //==========================================================================
459 
460 // Soft-eikonal trial function for initial-final.
461 
463 
464 public:
465 
466  // Name of trial generator.
467  virtual string name() override {return "TrialIFSoft";}
468 
469  // Trial antenna function. This trial generator uses x*PDF <= const
470  // as overestimate => no x-factor.
471  virtual double aTrial(double saj, double sjk, double sAK) override;
472 
473  // Evolution scale.
474  virtual double getQ2(double saj, double sjk, double sAK) override {
475  return (saj*sjk/(sAK+sjk));}
476  virtual double getQ2max(double sAK, double eA, double eBeamUsed) override {
477  double eAmax = ((sqrt(shhSav)/2.0) - (eBeamUsed-eA));
478  return (sAK*(eAmax-eA)/eA);}
479 
480  // Generate a new Q value, with first-order running alphaS.
481  virtual double genQ2run(double q2old, double sAK, double zMin, double zMax,
482  double colFac, double PDFratio, double b0, double kR, double Lambda,
483  double eA, double eK, double headroomFac = 1.0, double enhanceFac = 1.0)
484  override;
485 
486  // Generate a new Q value, with constant trial alphaS.
487  virtual double genQ2(double q2old, double sAK, double zMin, double zMax,
488  double colFac, double alphaSvalue, double PDFratio,
489  double eA, double eK, double headroomFac = 1.0, double enhanceFac = 1.0)
490  override;
491 
492  // Generate a new zeta value in [zMin,zMax].
493  virtual double genZ(double zMin, double zMax) override;
494 
495  // The zeta integral: dzeta/zeta/(zeta-1).
496  virtual double getIz(double zMin, double zMax) override;
497 
498  // The zeta boundaries, for a given value of the evolution scale.
499  virtual double getZmin(double Qt2, double sAK, double eA, double eBeamUsed)
500  override;
501  virtual double getZmax(double Qt2, double sAK, double eA, double eBeamUsed)
502  override;
503 
504  // Inverse transforms to obtain saj and sjk from Qt2 and zeta.
505  virtual double getS1j(double Qt2, double zeta, double sAK) override;
506  virtual double getSj2(double Qt2, double zeta, double sAK) override;
507 
508  // Trial PDF ratio.
509  virtual double trialPDFratio(BeamParticle* beamAPtr, BeamParticle* beamBPtr,
510  int iSys, int idA, int idK, double eA, double eK,
511  double Qt2A, double Qt2B) override;
512 
513 };
514 
515 //==========================================================================
516 
517 // Specialised soft-eikonal trial function for initial-final when
518 // initial-state parton is a valence quark.
519 
520 class TrialVFSoft : public TrialIFSoft {
521 
522 public:
523 
524  // Name of trial generator.
525  virtual string name() override {return "TrialVFSoft";}
526 
527  // Trial antenna function. This trial generator uses PDF <= const as
528  // overestimate => x-factor.
529  virtual double aTrial(double saj, double sjk, double sAK) override;
530 
531  // Generate a new zeta value in [zMin,zMax].
532  virtual double genZ(double zMin, double zMax) override;
533 
534  // The zeta integral: dzeta/(zeta-1).
535  virtual double getIz(double zMin, double zMax) override;
536 
537 };
538 
539 //==========================================================================
540 
541 // A gluon collinear trial function for initial-final.
542 
544 
545 public:
546 
547  // Name of trial generator.
548  virtual string name() override {return "TrialIFGCollA";}
549 
550  // Trial antenna function. This trial generator uses x*PDF ratio <=
551  // const as overestimate so the trial function has no extra
552  // x-factor.
553  virtual double aTrial(double saj, double sjk, double sAK) override;
554 
555  // Evolution scale.
556  virtual double getQ2(double saj, double sjk, double sAK) override {
557  return (saj*sjk/(sAK+sjk));}
558  virtual double getQ2max(double sAK, double eA, double eBeamUsed) override {
559  double eAmax = ( (sqrt(shhSav)/2.0) - (eBeamUsed-eA) );
560  return (sAK*(eAmax-eA)/eA);}
561 
562  // Generate a new Q value, with first-order running alphaS.
563  virtual double genQ2run(double q2old, double sAK, double zMin, double zMax,
564  double colFac, double PDFratio, double b0, double kR, double Lambda,
565  double eA, double eK, double headroomFac = 1.0, double enhanceFac = 1.0)
566  override;
567 
568  // Generate a new Q value, with constant trial alphaS.
569  virtual double genQ2(double q2old, double sAK, double zMin, double zMax,
570  double colFac, double alphaSvalue, double PDFratio,
571  double eA, double eK, double headroomFac = 1.0, double enhanceFac = 1.0)
572  override;
573 
574  // Generate a new zeta value in [zMin,zMax].
575  virtual double genZ(double zMin, double zMax) override;
576 
577  // The zeta integral.
578  virtual double getIz(double zMin, double zMax) override;
579 
580  // The zeta boundaries, for a given value of the evolution scale.
581  virtual double getZmin(double Qt2, double sAK, double eA, double eBeamUsed)
582  override;
583  virtual double getZmax(double Qt2, double sAK, double eA, double eBeamUsed)
584  override;
585 
586  // Inverse transforms to obtain saj and sjk from Qt2 and zeta.
587  virtual double getS1j(double Qt2, double zeta, double sAK) override;
588  virtual double getSj2(double Qt2, double zeta, double sAK) override;
589 
590  // Trial PDF ratio (= just a simple headroom factor).
591  virtual double trialPDFratio(BeamParticle* beamAPtr, BeamParticle* beamBPtr,
592  int iSys, int idA, int idK, double eA, double eK,
593  double Qt2A, double Qt2B) override;
594 
595 };
596 
597 
598 //==========================================================================
599 
600 // A splitting trial function for initial-final, q -> gqbar.
601 
603 
604 public:
605 
606  // Name of trial generator.
607  virtual string name() override {return "TrialIFSplitA";}
608 
609  // Trial antenna function. This trial generator uses the xf
610  // overestimate so no extra x-factor.
611  virtual double aTrial(double saj, double sjk, double sAK) override;
612 
613  // Evolution scale.
614  virtual double getQ2(double saj, double sjk, double sAK) override {
615  return ((useMevolSav) ? saj : saj*sjk/(sAK + sjk));}
616  virtual double getQ2max(double sAK, double eA, double eBeamUsed) override {
617  double xA = eA/(sqrt(shhSav)/2.0);
618  double eAmax = ((sqrt(shhSav)/2.0) - (eBeamUsed - eA));
619  return ((useMevolSav) ? sAK/xA : sAK*(eAmax - eA)/eA);}
620 
621  // Generate new Q value, with first-order running alphaS.
622  virtual double genQ2run(double q2old, double sAK, double zMin, double zMax,
623  double colFac, double PDFratio, double b0, double kR, double Lambda,
624  double eA, double eK, double headroomFac = 1.0, double enhanceFac = 1.0)
625  override;
626 
627  // Generate new Q value, with constant trial alphaS.
628  virtual double genQ2(double q2old, double sAK, double zMin, double zMax,
629  double colFac, double alphaSvalue, double PDFratio,
630  double eA, double eK, double headroomFac = 1.0, double enhanceFac = 1.0)
631  override;
632 
633  // Generate new Q value, with running of the PDFs towards the mass
634  // threshold.
635  virtual double genQ2thres(double q2old, double sAK,
636  double zMin, double zMax, double colFac, double alphaSvalue,
637  double PDFratio, int idA, int idK, double eA, double eK, bool useMpdf,
638  double headroomFac = 1.0, double enhanceFac = 1.0) override;
639 
640  // Generate a new zeta value in [zMin,zMax].
641  virtual double genZ(double zMin, double zMax) override;
642 
643  // The zeta integral.
644  virtual double getIz(double zMin, double zMax) override;
645 
646  // The zeta boundaries, for a given value of the evolution scale.
647  virtual double getZmin(double Qt2, double sAK, double eA, double eBeamUsed)
648  override;
649  virtual double getZmax(double Qt2, double sAK, double eA, double eBeamUsed)
650  override;
651 
652  // Inverse transforms to obtain saj and sjk from Qt2 and zeta.
653  virtual double getS1j(double Qt2, double zeta, double sAK) override;
654  virtual double getSj2(double Qt2, double zeta, double sAK) override;
655 
656  // Trial PDF ratio.
657  virtual double trialPDFratio(BeamParticle* beamAPtr, BeamParticle* beamBPtr,
658  int iSys, int idA, int idK, double eA, double eK,
659  double Qt2A, double Qt2B) override;
660 
661 };
662 
663 //==========================================================================
664 
665 // K splitting trial function for initial-final, g -> qqbar.
666 
668 
669  public:
670 
671  // Name of trial generator.
672  virtual string name() override {return "TrialIFSplitK";}
673 
674  // Trial antenna function. This trial uses the xf overestimate so no
675  // extra x factor.
676  virtual double aTrial(double saj, double sjk, double sAK) override;
677 
678  // Evolution scale.
679  virtual double getQ2(double saj, double sjk, double sAK) override {
680  return ((useMevolSav) ? sjk : saj*sjk/(sAK + sjk));}
681  virtual double getQ2max(double sAK, double eA, double eBeamUsed) override {
682  double xA = eA/(sqrt(shhSav)/2.0);
683  double eAmax = ((sqrt(shhSav)/2.0) - (eBeamUsed - eA));
684  return ((useMevolSav) ? sAK*(1.0 - xA)/xA : sAK*(eAmax - eA)/eA);}
685 
686  // Generate a new Q value, with first-order running alphaS.
687  virtual double genQ2run(double q2old, double sAK, double zMin, double zMax,
688  double colFac, double PDFratio, double b0, double kR, double Lambda,
689  double eA, double eK, double headroomFac = 1.0, double enhanceFac = 1.0)
690  override;
691 
692  // Generate a new Q value, with constant trial alphaS.
693  virtual double genQ2(double q2old, double sAK, double zMin, double zMax,
694  double colFac, double alphaSvalue, double PDFratio,
695  double eA, double eK, double headroomFac = 1.0, double enhanceFac = 1.0)
696  override;
697 
698  // Generate a new zeta value in [zMin,zMax].
699  virtual double genZ(double zMin, double zMax) override;
700 
701  // The zeta integral.
702  virtual double getIz(double zMin, double zMax) override;
703 
704  // The zeta boundaries, for a given value of the evolution scale.
705  virtual double getZmin(double Qt2, double sAK, double eA, double eBeamUsed)
706  override;
707  virtual double getZmax(double Qt2, double sAK, double eA, double eBeamUsed)
708  override;
709 
710  // Inverse transforms to obtain saj and sjk from Qt2 and zeta.
711  virtual double getS1j(double Qt2, double zeta, double sAK) override;
712  virtual double getSj2(double Qt2, double zeta, double sAK) override;
713 
714  // Trial PDF ratio.
715  virtual double trialPDFratio(BeamParticle* beamAPtr, BeamParticle* beamBPtr,
716  int iSys, int idA, int idK, double eA, double eK,
717  double Qt2A, double Qt2B) override;
718 
719 };
720 
721 //==========================================================================
722 
723 // A conversion trial function for initial-final, g -> qqbar.
724 
726 
727 public:
728 
729  // Name of trial generator.
730  virtual string name() override {return "TrialIFConvA";}
731 
732  // Trial antenna function. This trial currently uses the xf
733  // overestimate so no extra x-factor (but could be changed to use
734  // the f overestimate if many violations, and/or for valence
735  // flavours).
736  virtual double aTrial(double saj, double sjk, double sAK) override;
737 
738  // Evolution scale.
739  virtual double getQ2(double saj, double sjk, double sAK) override {
740  return ((useMevolSav) ? saj : saj*sjk/(sAK + sjk));}
741  virtual double getQ2max(double sAK, double eA, double eBeamUsed) override {
742  double xA = eA/(sqrt(shhSav)/2.0);
743  double eAmax = ((sqrt(shhSav)/2.0) - (eBeamUsed - eA));
744  return ((useMevolSav) ? sAK/xA : sAK*(eAmax-eA)/eA);}
745 
746  // Generate a new Q value, with first-order running alphaS.
747  virtual double genQ2run(double q2old, double sAK, double zMin, double zMax,
748  double colFac, double PDFratio, double b0, double kR, double Lambda,
749  double eA, double eK, double headroomFac = 1.0, double enhanceFac = 1.0)
750  override;
751 
752  // Generate a new Q value, with constant trial alphaS.
753  virtual double genQ2(double q2old, double sAK, double zMin, double zMax,
754  double colFac, double alphaSvalue, double PDFratio,
755  double eA, double eK, double headroomFac = 1.0, double enhanceFac = 1.0)
756  override;
757 
758  // Generate a new zeta value in [zMin,zMax].
759  virtual double genZ(double zMin, double zMax) override;
760 
761  // The zeta integral.
762  virtual double getIz(double zMin, double zMax) override;
763 
764  // The zeta boundaries, for a given value of the evolution scale.
765  virtual double getZmin(double Qt2, double sAK, double eA, double eBeamUsed)
766  override;
767  virtual double getZmax(double Qt2, double sAK, double eA, double eBeamUsed)
768  override;
769 
770  // Inverse transforms to obtain saj and sjk from Qt2 and zeta.
771  virtual double getS1j(double Qt2, double zeta, double sAK) override;
772  virtual double getSj2(double Qt2, double zeta, double sAK) override;
773 
774  // Trial PDF ratio.
775  virtual double trialPDFratio(BeamParticle* beamAPtr, BeamParticle* beamBPtr,
776  int iSys, int idA, int idK, double eA, double eK,
777  double Qt2A, double Qt2B) override;
778 
779 };
780 
781 //==========================================================================
782 
783 // The BranchElementalISR class, container for 2 -> 3 trial branchings.
784 
785 // Input: i1In carries colour (or incoming anticolour)
786 // i2In carries anticolour (or incoming colour)
787 // colIn colour tag
788 // isVal1In (true if i1In is incoming and is a valence parton on its
789 // side/system)
790 // isVal2In (true if i2In is incoming and is a valence parton on its
791 // side/system)
792 // Internal storage for II:
793 // i1sav: whichever of i1In, i2In has positive pz
794 // i2sav: whichever of i1In, i2In has negative pz
795 // is1A : always true (i1sav always has positive pz)
796 // isII : true
797 // Internal storage for IF:
798 // i1sav: whichever of i1In, i2In is the initial-state leg
799 // i2sav: whichever of i1In, i2In is the final-state leg
800 // is1A : true if i1sav has positive pz, false if it has negative pz
801 // isII : false
802 
804 
805 public:
806 
807  // Constructors.
808  BranchElementalISR() = default;
809  BranchElementalISR(int iSysIn, Event& event, int iOld1In,
810  int iOld2In, int colIn, bool isVal1In, bool isVal2In) {
811  reset(iSysIn, event, iOld1In, iOld2In, colIn, isVal1In, isVal2In);}
812 
813  // Main method to initialize/reset a BranchElemental.
814  void reset(int iSysIn, Event& event, int i1In, int i2In, int colIn,
815  bool isVal1In, bool isVal2In);
816 
817  // Antenna mass (negative if spacelike virtual).
818  double mAnt() const {return mAntSav;}
819  double m2Ant() const {return m2AntSav;}
820  double sAnt() const {return sAntSav;}
821  // Dot products of daughters.
822  double s12() const {return 2*new1.p()*new2.p();}
823  double s23() const {return 2*new2.p()*new3.p();}
824  double s13() const {return 2*new1.p()*new3.p();}
825 
826  // This is an II or an IF type.
827  bool isII() const {return isIIsav;}
828  // Is 1 a side A (p+) guy (need to know for pdfs in IF).
829  bool is1A() const {return is1Asav;}
830  // Valence.
831  bool isVal1() const {return isVal1sav;}
832  bool isVal2() const {return isVal2sav;}
833  int colType1() const {return colType1sav;}
834  int colType2() const {return colType2sav;}
835  int col() const {return colSav;}
836  int geti1() {return i1sav;}
837  int geti2() {return i2sav;}
838  int getSystem() {return system;}
839 
840  // Function to reset all trial generators for this branch elemental.
841  void clearTrialGenerators();
842 
843  // Add a trial generator to this BranchElemental.
844  void addTrialGenerator(int iAntPhysIn, bool swapIn,
845  TrialGeneratorISR* trialGenPtrIn);
846 
847  // Add to and get rescue levels.
848  void addRescue(int iTrial) {nShouldRescue[iTrial]++;}
849  int getNshouldRescue(int iTrial) {return nShouldRescue[iTrial];}
850  void resetRescue() {
851  for (int i=0; i<(int)nShouldRescue.size(); i++) nShouldRescue[i] = 0;}
852 
853  // Function to return number of trial generators for this antenna.
854  int nTrialGenerators() const {return trialGenPtrsSav.size();}
855 
856  // Save a generated trial branching.
857  void saveTrial(int iTrial, double qOld, double qTrial, double zMin=0.,
858  double zMax=0., double colFac=0.,double alphaEff=0., double pdfRatio=0.,
859  int trialFlav=0, double extraMpdf=0., double headroom = 1.0,
860  double enhanceFac = 1.0);
861 
862  // Add the physical pdf ratio.
863  void addPDF(int iTrial,double pdfRatio) {physPDFratioSav[iTrial] = pdfRatio;}
864 
865  // Generate invariants for saved branching.
866  bool genTrialInvariants(double& s1j, double& sj2,
867  double eBeamUsed, int iTrial = -1);
868 
869  // Get trial function index of winner.
870  int getTrialIndex() const;
871 
872  // Check if a saved trial exists for a particular trialGenerator.
873  bool hasTrial(int iTrial) const {
874  if (iTrial < int(hasSavedTrial.size())) return hasSavedTrial[iTrial];
875  else return false;}
876 
877  // Get whether physical antenna is swapped.
878  bool getIsSwapped(int iTrial = -1) const {
879  if (iTrial <= -1) iTrial = getTrialIndex();
880  return isSwappedSav[iTrial];}
881 
882  // Get physical antenna function index of winner.
883  int getPhysIndex(int iTrial = -1) const {
884  if (iTrial <= -1) iTrial = getTrialIndex();
885  return iAntPhysSav[iTrial];}
886 
887  // Get scale for a specific saved trial.
888  double getTrialScale(int iTrial) const {
889  if (iTrial < int(scaleSav.size())) return scaleSav[iTrial];
890  else return -1.0;}
891 
892  // Get scale of winner.
893  double getTrialScale() const;
894 
895  // Get colour factor.
896  double getColFac(int iTrial = -1) {
897  if (iTrial <= -1) iTrial = getTrialIndex();
898  return colFacSav[iTrial];}
899 
900  // Get headroom factor.
901  double getHeadroomFac(int iTrial = -1) {
902  if (iTrial <= -1) iTrial = getTrialIndex();
903  return headroomSav[iTrial];}
904 
905  // Get headroom factor.
906  double getEnhanceFac(int iTrial = -1) {
907  if (iTrial <= -1) iTrial = getTrialIndex();
908  return enhanceFacSav[iTrial];}
909 
910  // Get alpha S.
911  double getAlphaTrial(int iTrial = -1) {
912  if (iTrial <= -1) iTrial = getTrialIndex();
913  return alphaSav[iTrial];}
914 
915  // Get pdf ratio.
916  double getPDFratioTrial(int iTrial = -1) {
917  if (iTrial <= -1) iTrial = getTrialIndex();
918  return trialPDFratioSav[iTrial];}
919 
920  // For gluon conversions, get ID of quark flavour to convert to.
921  int getTrialFlav(int iTrial = -1) {
922  if (iTrial <= -1) iTrial = getTrialIndex();
923  return trialFlavSav[iTrial];}
924 
925  // Get pdf ratio.
926  double getPDFratioPhys(int iTrial = -1) {
927  if (iTrial <= -1) iTrial = getTrialIndex();
928  return physPDFratioSav[iTrial];}
929 
930  // Get the extra factor when getting rid of the heavy quarks.
931  double getExtraMassPDFfactor(int iTrial = -1) {
932  if (iTrial <= -1) iTrial = getTrialIndex();
933  return extraMassPDFfactorSav[iTrial];}
934 
935  // Flag to set if a saved trial should be ignored and a new one generated.
936  // Default value -1 : force all trials to renew.
937  void renewTrial(int iTrial = -1) {
938  if (iTrial >= 0) hasSavedTrial[iTrial] = false;
939  else for (iTrial = 0; iTrial < int(hasSavedTrial.size()); ++iTrial)
940  hasSavedTrial[iTrial] = false;}
941 
942  // List function.
943  void list(bool header = false, bool footer = false) const;
944 
945  // Data storage members.
946  int i1sav{}, i2sav{}, id1sav{}, id2sav{}, colType1sav{}, colType2sav{},
947  h1sav{}, h2sav{};
948  double e1sav{}, e2sav{};
949  bool isVal1sav{}, isVal2sav{}, isIIsav{}, is1Asav{};
950  Particle new1{}, new2{}, new3{};
951  // Colour, not obvious, since for e.g. gg -> H we have two II antennae.
952  int colSav{};
953  // System and counter for vetos.
954  int system{0}, nVeto{}, nHull{}, nHadr{};
955  // We have to force a splitting (heavy quarks).
956  bool forceSplittingSav{};
957 
958  // Trial Generators and properties of saved trials.
959  vector<TrialGeneratorISR*> trialGenPtrsSav{};
960  vector<double> zMinSav{}, zMaxSav{}, colFacSav{}, alphaSav{};
961  vector<double> physPDFratioSav{}, trialPDFratioSav{};
962  vector<double> extraMassPDFfactorSav{};
963  vector<double> scaleSav{}, scaleOldSav{}, headroomSav{}, enhanceFacSav{};
964  vector<bool> hasSavedTrial{}, isSwappedSav{};
965  vector<int> iAntPhysSav{}, nShouldRescue{}, trialFlavSav{};
966  // Note: isSwapped = true for II means physical antenna function is
967  // coded for side A but trial generator is for side B. For IF, is1A
968  // = true for 1 being on side A, false for 1 being on side B.
969 
970  private:
971 
972  // Saved antenna invariant mass value.
973  double m2AntSav{}, mAntSav{}, sAntSav{};
974 
975 };
976 
977 //==========================================================================
978 
979 // The VinciaISR class
980 // Main shower class for initial-state (II and IF) antenna showers
981 // Inherits from SpaceShower in Pythia 8 so can be used as alternative to
982 // SpaceShower.
983 // Methods that must replace ones in SpaceShower are marked with override.
984 
985 class VinciaISR : public SpaceShower {
986 
987  // Allow VinciaFSR to access private information.
988  friend class VinciaFSR;
989 
990 public:
991 
992  // Constructor.
993  VinciaISR() : isInit(false) {;}
994 
995  // Destructor.
996  virtual ~VinciaISR() {;}
997 
998  // Initialize shower. Possibility to force re-initialization by hand.
999  void init(BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn) override;
1000 
1001  // Possible limitation of first emission.
1002  bool limitPTmax(Event& event, double Q2Fac = 0., double Q2Ren = 0.) override;
1003 
1004  // Prepare system for evolution; identify ME.
1005  void prepare(int iSys, Event& event, bool limitPTmaxIn = false) override;
1006 
1007  // Update dipole list after each ISR emission.
1008  void update( int iSys, Event& event, bool hasWeakRad = false) override;
1009 
1010  // Select next pT in downwards evolution.
1011  double pTnext(Event& event, double pTbegAll, double pTendAll,
1012  int nRadIn = -1, bool doTrialIn = false) override;
1013 
1014  // Perform a branching (as defined by current "winner"). Returns
1015  // true of the branching was accepted and false if rejected.
1016  bool branch(Event& event) override;
1017 
1018  // Print a list of II and IF dipole-antennae.
1019  void list() const override;
1020 
1021  // Initialize data members for calculation of uncertainty bands. So
1022  // far purely dummy in Vincia.
1023  bool initUncertainties() override {return false;}
1024 
1025  // Flag for failure in branch(...) that will force a retry of parton
1026  // level. Returns false in base class, and rescatterFail in
1027  // simpleTimeShower. Not implemented in Vincia yet since we do not
1028  // do rescattering.
1029  bool doRestart() const override {return false;}
1030 
1031  // Flag for last ISR branching being gamma -> qqbar. Not implemented
1032  // in Vincia's QED evolution yet.
1033  bool wasGamma2qqbar() override {return false;}
1034 
1035  // Tell whether ISR has done a weak emission. Not implemented in Vincia yet.
1036  bool getHasWeaklyRadiated() override {return false;}
1037 
1038  // Tell which system was the last processed one.
1039  int system() const override {return iSysWin;}
1040 
1041  // Potential enhancement factor of pTmax scale for hardest emission.
1042  // Used if limitPTmax = true.
1043  double enhancePTmax() const override {return pTmaxFudge;}
1044 
1045  // Initialise pointers to Vincia objects.
1046  void initVinciaPtrs(Colour* colourPtrIn, shared_ptr<VinciaFSR> fsrPtrIn,
1047  QEDShower* qedPtrIn,MECs* mecsPtrIn,Resolution* resolutionPtrIn,
1048  VinciaCommon* vinComPtrIn,VinciaWeights* vinWeightsPtrIn);
1049 
1050  // Clear all containers.
1051  void clearContainers();
1052 
1053  // Initialize pointers to antenna sets.
1054  void initAntPtr(AntennaSetISR* antSetIn) {antSetPtr = antSetIn;}
1055 
1056  // Function to tell if VinciaISR::prepare() has treated this system.
1057  bool prepared(int iSys) {
1058  if (hasPrepared.find(iSys) != hasPrepared.end()) return hasPrepared[iSys];
1059  else return false;}
1060 
1061  // Wrapper function to return a specific antenna inside antenna set
1062  AntennaFunctionIX* getAnt(int iAnt) {return antSetPtr->getAnt(iAnt);}
1063 
1064  // Print final statistics information.
1065  void printInfo(bool pluginCall = false);
1066 
1067  // Evolution windows, phase space region boundaries.
1068  int getRegion(double q) {
1069  for (int i=1; i<(int)regMinScalesSav.size(); i++)
1070  if (q < regMinScalesSav[i]) return i-1;
1071  return (int)regMinScalesSav.size() - 1;}
1072 
1073  // Evolution window, phase space region boundaries.
1074  double getQmin(int iRegion) {
1075  iRegion = max(0,iRegion);
1076  iRegion = min(iRegion,(int)regMinScalesSav.size()-1);
1077  return regMinScalesSav[iRegion];}
1078 
1079  // Number of active flavors.
1080  int getNf(int iRegion) {
1081  if (iRegion <= 1) return 3;
1082  else if (iRegion <= 2) return 4;
1083  else if (iRegion <= 4) return 5;
1084  else return 6;}
1085 
1086  // Lambda value
1087  double getLambda(int nFin) {
1088  if (nFin <= 3) return alphaSptr->Lambda3();
1089  else if (nFin <= 4) return alphaSptr->Lambda4();
1090  else if (nFin <= 5) return alphaSptr->Lambda5();
1091  else return alphaSptr->Lambda6();}
1092 
1093  // Add trial functions to a BranchElemental.
1094  void resetTrialGenerators(BranchElementalISR* trial);
1095 
1096  // Method to check if a gluon splitting in the initial state (to get
1097  // rid of heavy quarks) is still possible after the current
1098  // branching.
1099  bool checkHeavyQuarkPhaseSpace(vector<Particle> parts, int iSyst);
1100 
1101  // Method to check if heavy quark left after passing the evolution window.
1102  bool heavyQuarkLeft(double qTrial);
1103 
1104  // Function to ask if a system is hard.
1105  bool isSysHard(int iSys) {
1106  if (!isInit) return false;
1107  if ((int)isHardSys.size() <= iSys) return false;
1108  return isHardSys[iSys];}
1109 
1110  // Return a vector of the masses.
1111  vector<double> getMasses() {return vector<double> {mt, mtb, mb, mc, ms};}
1112 
1113  // Get number of systems.
1114  int getNsys() {return nBranchISR.size();}
1115  // Get number of branchings in a system (return -1 if no such system).
1116  int getNbranch(int iSys = -1) {
1117  int n = 0;
1118  if (iSys < 0) for (int i = 0; i < (int)nBranchISR.size(); ++i)
1119  n += nBranchISR[iSys];
1120  else if (iSys < (int)nBranchISR.size()) n = nBranchISR[iSys];
1121  else n = -1;
1122  return n;}
1123 
1124  // Get scale of branchings; use (0,1) for first branching in 1st system.
1125  // Could be extended so (i,0) would return starting scale for system i.
1126  double getQbranch(int iSys, int iBranch) {
1127  if (iSys < 4 && iSys >= 0 && iBranch <= 10 && iBranch >= 1)
1128  return qBranch[iSys][iBranch];
1129  else return 0.;}
1130  double getPTphys(int iSys, int iBranch) {
1131  if (iSys < 4 && iSys >= 0 && iBranch <= 10 && iBranch >= 1)
1132  return pTphys[iSys][iBranch];
1133  else return 0.;}
1134 
1135  // Set verbosity level.
1136  void setVerbose(int verboseIn) {verbose = verboseIn;}
1137 
1138  // Check the antennae.
1139  bool checkAntennae(const Event& event);
1140 
1141  // Pointer to global AlphaStrong instance.
1143 
1144 private:
1145 
1146  // Set starting scale of shower (power vs wimpy) for system iSys.
1147  void setStartScale(int iSys, Event& event);
1148 
1149  // Function to return headroom factor.
1150  double getHeadroomFac(int iSys, int iAntPhys, double qMinNow);
1151 
1152  // Generate trial branching kinematics and check physical phase space
1153  bool generateKinematics(Event& event, BranchElementalISR* trialPtr,
1154  vector<Vec4>& pRec) {
1155  return ( trialPtr->isII()
1156  ? generateKinematicsII(event, trialPtr, pRec)
1157  : generateKinematicsIF(event, trialPtr, pRec) ); }
1158 
1159  // Generate kinematics (II) and set flavours and masses.
1160  bool generateKinematicsII(Event& event, BranchElementalISR* trialPtr,
1161  vector<Vec4>& pRec);
1162 
1163  // Generate kinematics (IF) and set flavours and masses.
1164  bool generateKinematicsIF(Event& event, BranchElementalISR* trialPtr,
1165  vector<Vec4>& pRec);
1166 
1167  // Main trial accept function.
1168  bool acceptTrial(const Event& event, BranchElementalISR* winnerPtr);
1169 
1170  // Method to assign colour flow.
1171  bool assignColourFlow(Event& event, BranchElementalISR* trialPtr);
1172 
1173  // Initialised.
1174  bool isInit;
1175  bool isPrepared;
1176 
1177  // Possibility to allow user veto of emission step.
1178  bool hasUserHooks, canVetoEmission;
1179 
1180  // Beams, saved as positive and negative pz respectively.
1181  int beamFrameType;
1182  double eBeamA, eBeamB, eCMBeamsSav, m2BeamsSav;
1183  double TINYPDF;
1184 
1185  // Main Vincia ISR on/off switches.
1186  bool doII, doIF;
1187 
1188  // Map of which systems ISR::prepare() has treated.
1189  map<int, bool> hasPrepared;
1190 
1191  // Shower parameters.
1192  bool helicityShower, sectorShower, convGluonToQuarkI, convQuarkToGluonI;
1193  bool kineMapIFretry;
1194  int nGluonToQuarkI, nGluonToQuarkF;
1195  double cutoffScaleII, cutoffScaleIF;
1196  int nFlavZeroMass;
1197 
1198  // Factorization scale and shower starting settings.
1199  int pTmaxMatch;
1200  double pTmaxFudge, pT2maxFudge, pT2maxFudgeMPI;
1201 
1202  // AlphaS parameters.
1203  bool useCMW;
1204  int alphaSorder;
1205  double alphaSvalue, alphaSmax, alphaSmuFreeze, alphaSmuMin;
1206  double aSkMu2EmitI, aSkMu2SplitI, aSkMu2SplitF, aSkMu2Conv;
1207  double mt, mtb, ms, mb, mc;
1208  // Calculated values.
1209  double mu2freeze, mu2min;
1210 
1211  // Trial generators.
1212  TrialIISoft trialIISoft;
1213  TrialIIGCollA trialIIGCollA;
1214  TrialIIGCollB trialIIGCollB;
1215  TrialIISplitA trialIISplitA;
1216  TrialIISplitB trialIISplitB;
1217  TrialIIConvA trialIIConvA;
1218  TrialIIConvB trialIIConvB;
1219  TrialIFSoft trialIFSoft;
1220  TrialVFSoft trialVFSoft;
1221  TrialIFGCollA trialIFGCollA;
1222  TrialIFSplitA trialIFSplitA;
1223  TrialIFSplitK trialIFSplitK;
1224  TrialIFConvA trialIFConvA;
1225 
1226  // Enhanceing switches and parameters.
1227  bool enhanceInHard, enhanceInResDec, enhanceInMPI;
1228  double enhanceAll, enhanceBottom, enhanceCharm, enhanceCutoff;
1229 
1230  // Pointer to VINCIA objects.
1231  AntennaSetISR* antSetPtr;
1232  MECs* mecsPtr;
1233  Colour* colourPtr;
1234  Resolution* resolutionPtr;
1235  QEDShower* qedShowerPtr;
1236  shared_ptr<VinciaFSR> fsrPtr;
1237  VinciaCommon* vinComPtr;
1238  VinciaWeights* weightsPtr;
1239 
1240  // Total and MEC accept probability.
1241  vector<double> Paccept;
1242 
1243  // Evolution windows.
1244  vector<double> regMinScalesMtSav;
1245  vector<double> regMinScalesSav;
1246  vector<double> regMinScalesNow;
1247 
1248  // Vector of dipoles (with trial branchings, 4 at most).
1249  vector<BranchElementalISR > branchElementals;
1250 
1251  // Current winner.
1252  BranchElementalISR* winnerPtr;
1253  int indxWin;
1254  int iSysWin;
1255 
1256  // Flags to tell a few basic properties of each parton system.
1257  map<int, bool> isHardSys, isResonanceSys, polarisedSys, doMECsSys;
1258 
1259  // Saved particle state and number in event record.
1260  map<int, vector< Particle > > partsSav;
1261  map<int, vector< int > > indexSav;
1262 
1263  // Save initial ISR starting scale system by system.
1264  map<int, double> Q2hat;
1265 
1266  // Count the number of branchings in the system.
1267  map<int, int> nBranch, nBranchISR;
1268 
1269  // Saved incoming guys.
1270  map<int, Particle> initialA;
1271  map<int, Particle> initialB;
1272  double eBeamAUsed, eBeamBUsed;
1273 
1274  // Count numbers of quarks and gluons.
1275  map<int, int> nG, nQQ;
1276 
1277  // Statistics.
1278  long nTrialsSum;
1279  vector<long> nTrials, nTrialsAccepted, nFailedVeto, nFailedHull, nFailedKine;
1280  vector<long> nFailedMass, nFailedCutoff, nClosePSforHQ, nSectorReject;
1281  vector<double> rFailedVetoPDF;
1282  int nForceSplit;
1283 
1284  // Rescue mechanism.
1285  bool doRescue;
1286  int nRescue;
1287  double rescueMin;
1288 
1289  // Verbose setting, P>1 warning settings, and version number.
1290  int verbose;
1291  double qBranch[4][11], pTphys[4][11];
1292 
1293  // Counter for numbers of events.
1294  long nAccepted, nSelected;
1295  int nVetoUserHooks, nFailHadLevel, nCallPythiaNext;
1296 
1297 };
1298 
1299 //==========================================================================
1300 
1301 } // end namespace Pythia8
1302 
1303 #endif // end Pythia8_VinciaISR_H
virtual double aTrial(double saj, double sjb, double sAB)
Trial antenna function.
Definition: VinciaISR.cc:58
virtual string name() override
Name of trial generator.
Definition: VinciaISR.h:240
virtual string name() override
Name of trial generator.
Definition: VinciaISR.h:147
virtual double trialPDFratio(BeamParticle *beamAPtr, BeamParticle *beamBPtr, int iSys, int idA, int idB, double eA, double eB, double Qt2A, double Qt2B) override
Trial PDF ratio.
Definition: VinciaISR.h:224
virtual double getQ2(double saj, double sjb, double sAB) override
Evolution scale.
Definition: VinciaISR.h:427
virtual string name() override
Name of trial generator.
Definition: VinciaISR.h:467
virtual double getS1j(double Qt2, double zeta, double sAB) override
Inverse transforms to obtain saj and sjb from Qt2 and zeta.
Definition: VinciaISR.h:341
virtual double genQ2(double q2old, double sAB, double zMin, double zMax, double colFac, double alphaSvalue, double PDFratio, double eA, double eB, double headroomFac=1.0, double enhanceFac=1.0)
Generate new Q value, with constant trial alphaS.
Definition: VinciaISR.cc:95
A conversion trial function for initial-initial, g -> qqbar.
Definition: VinciaISR.h:358
bool wasGamma2qqbar() override
Definition: VinciaISR.h:1033
constexpr double pow2(const double &x)
Powers of small integers - for balance speed/code clarity.
Definition: PythiaStdlib.h:153
Definition: VinciaISR.h:520
void initPtr(Info *infoPtrIn)
Initialize pointers.
Definition: VinciaISR.cc:22
double getLambda(int nFin)
Lambda value.
Definition: VinciaISR.h:1087
void addRescue(int iTrial)
Add to and get rescue levels.
Definition: VinciaISR.h:848
virtual double genQ2thres(double q2old, double sAB, double zMin, double zMax, double colFac, double alphaSvalue, double PDFratio, int idA, int idB, double eA, double eB, bool useMpdf, double headroomFac=1.0, double enhanceFac=1.0) override
Definition: VinciaISR.h:332
A splitting trial function for initial-final, q -> gqbar.
Definition: VinciaISR.h:602
double getExtraMassPDFfactor(int iTrial=-1)
Get the extra factor when getting rid of the heavy quarks.
Definition: VinciaISR.h:931
virtual double getQ2(double saj, double sjb, double sAB) override
Evolution scale.
Definition: VinciaISR.h:167
bool hasTrial(int iTrial) const
Check if a saved trial exists for a particular trialGenerator.
Definition: VinciaISR.h:873
virtual double getS1j(double Qt2, double zeta, double sAB) override
Inverse transforms to obtain saj and sjb from Qt2 and zeta.
Definition: VinciaISR.h:218
Definition: StandardModel.h:23
virtual double genQ2run(double q2old, double sAB, double zMin, double zMax, double colFac, double PDFratio, double b0, double kR, double Lambda, double eA, double eB, double headroomFac=1.0, double enhanceFac=1.0) override
Generate a new Q value, with first-order running alphaS.
Definition: VinciaISR.cc:544
virtual string name() override
Name of trial generator.
Definition: VinciaISR.h:730
virtual double genQ2run(double q2old, double sAB, double zMin, double zMax, double colFac, double PDFratio, double b0, double kR, double Lambda, double eA, double eB, double headroomFac=1.0, double enhanceFac=1.0) override
Generate a new Q value, with first-order running alphaS.
Definition: VinciaISR.h:431
VinciaISR()
Constructor.
Definition: VinciaISR.h:993
Definition: Info.h:42
int system() const override
Tell which system was the last processed one.
Definition: VinciaISR.h:1039
int getPhysIndex(int iTrial=-1) const
Get physical antenna function index of winner.
Definition: VinciaISR.h:883
int trialFlavSav
For conversion trial generators.
Definition: VinciaISR.h:116
The Event class holds all info on the generated event.
Definition: Event.h:379
double shhSav
s for hadron hadron.
Definition: VinciaISR.h:113
bool isVal1() const
Valence.
Definition: VinciaISR.h:831
bool isII() const
This is an II or an IF type.
Definition: VinciaISR.h:827
int getTrialFlav(int iTrial=-1)
For gluon conversions, get ID of quark flavour to convert to.
Definition: VinciaISR.h:921
int trialFlav()
Return the last trial flavor.
Definition: VinciaISR.h:100
Definition: BeamParticle.h:133
virtual double getQ2(double saj, double sjk, double sAK) override
Evolution scale.
Definition: VinciaISR.h:739
virtual double getSj2(double Qt2, double zeta, double sAB) override
Definition: VinciaISR.cc:331
virtual double getS1j(double Qt2, double zeta, double sAB)
Inverse transforms to obtain saj and sjb from Qt2 and zeta.
Definition: VinciaISR.cc:167
A collinear trial function for initial-initial.
Definition: VinciaISR.h:155
AntennaFunctionIX * getAnt(int iAnt)
Wrapper function to return a specific antenna inside antenna set.
Definition: VinciaISR.h:1062
K splitting trial function for initial-final, g -> qqbar.
Definition: VinciaISR.h:667
double s12() const
Dot products of daughters.
Definition: VinciaISR.h:822
B splitting trial function for initial-initial, q -> gqbar.
Definition: VinciaISR.h:301
double getPDFratioPhys(int iTrial=-1)
Get pdf ratio.
Definition: VinciaISR.h:926
Definition: VinciaAntennaFunctions.h:415
virtual string name() override
Name of trial generator.
Definition: VinciaISR.h:363
virtual string name() override
Name of trial generator.
Definition: VinciaISR.h:607
virtual double getQ2(double saj, double sjb, double sAB) override
Evolution scale.
Definition: VinciaISR.h:370
The BranchElementalISR class, container for 2 -> 3 trial branchings.
Definition: VinciaISR.h:803
Definition: VinciaISR.h:28
virtual string name()
Name of trial generator.
Definition: VinciaISR.h:40
double enhancePTmax() const override
Definition: VinciaISR.h:1043
Soft-eikonal trial function for initial-final.
Definition: VinciaISR.h:462
virtual double genQ2thres(double q2old, double sAB, double zMin, double zMax, double colFac, double alphaSvalue, double PDFratio, int idA, int idB, double eA, double eB, bool useMpdf, double headroomFac=1.0, double enhanceFac=1.0)
Definition: VinciaISR.cc:119
bool isSysHard(int iSys)
Function to ask if a system is hard.
Definition: VinciaISR.h:1105
A conversion trial function for initial-final, g -> qqbar.
Definition: VinciaISR.h:725
virtual double getS1j(double Qt2, double zeta, double sAB) override
Inverse transforms to obtain saj and sjb from Qt2 and zeta.
Definition: VinciaISR.cc:483
virtual double getSj2(double Qt2, double zeta, double sAB)
Definition: VinciaISR.cc:179
virtual double getS1j(double Qt2, double zeta, double sAB) override
Inverse transforms: to obtain saj and sjb from Qt2 and zeta.
Definition: VinciaISR.h:445
Definition: VinciaISR.h:142
virtual double getQ2(double saj, double sjk, double sAK) override
Evolution scale.
Definition: VinciaISR.h:679
bool useMevolSav
Use m or pT evolution for collinear singularities.
Definition: VinciaISR.h:110
virtual double aTrial(double saj, double sjb, double sAB) override
Trial antenna function.
Definition: VinciaISR.h:423
virtual double getTrialPDFratio()
Return last trial PDF ratio.
Definition: VinciaISR.h:94
virtual string name() override
Name of trial generator.
Definition: VinciaISR.h:548
virtual double genQ2run(double q2old, double sAB, double zMin, double zMax, double colFac, double PDFratio, double b0, double kR, double Lambda, double eA, double eB, double headroomFac=1.0, double enhanceFac=1.0) override
Generate a new Q value, with first-order running alphaS.
Definition: VinciaISR.cc:364
virtual double genQ2run(double q2old, double sAB, double zMin, double zMax, double colFac, double PDFratio, double b0, double kR, double Lambda, double eA, double eB, double headroomFac=1.0, double enhanceFac=1.0) override
Generate new Q value, with first-order running alphaS.
Definition: VinciaISR.h:317
void setVerbose(int verboseIn)
Set verbosity level.
Definition: VinciaISR.h:1136
double trialPDFratioSav
Saved trial PDF ratio and trial tolerance.
Definition: VinciaISR.h:124
double getHeadroomFac(int iTrial=-1)
Get headroom factor.
Definition: VinciaISR.h:901
double mAnt() const
Antenna mass (negative if spacelike virtual).
Definition: VinciaISR.h:818
B collinear trial function for initial-initial.
Definition: VinciaISR.h:206
virtual string name() override
Name of trial generator.
Definition: VinciaISR.h:160
virtual double trialPDFratio(BeamParticle *beamAPtr, BeamParticle *beamBPtr, int iSys, int idA, int idB, double eA, double eB, double Qt2A, double Qt2B) override
Trial PDF ratio.
Definition: VinciaISR.h:347
A splitting trial function for initial-initial, q -> gqbar.
Definition: VinciaISR.h:235
Definition: VinciaCommon.h:347
Info * infoPtr
Pointers.
Definition: VinciaISR.h:105
void renewTrial(int iTrial=-1)
Definition: VinciaISR.h:937
int nTrialGenerators() const
Function to return number of trial generators for this antenna.
Definition: VinciaISR.h:854
Definition: VinciaFSR.h:600
Definition: Basics.h:349
bool doRestart() const override
Definition: VinciaISR.h:1029
Definition: VinciaCommon.h:490
virtual double getQ2(double saj, double sjb, double sAB) override
Evolution scale.
Definition: VinciaISR.h:313
Definition: VinciaISR.h:985
virtual string name() override
Name of trial generator.
Definition: VinciaISR.h:525
int getNsys()
Get number of systems.
Definition: VinciaISR.h:1114
int getRegion(double q)
Evolution windows, phase space region boundaries.
Definition: VinciaISR.h:1068
TrialGeneratorISR()
Constructor.
Definition: VinciaISR.h:33
virtual double aTrial(double saj, double sjb, double sAB) override
Trial antenna function.
Definition: VinciaISR.h:214
bool prepared(int iSys)
Function to tell if VinciaISR::prepare() has treated this system.
Definition: VinciaISR.h:1057
virtual double aTrial(double saj, double sjb, double sAB) override
A conversion trial function for initial-initial, g -> qqbar.
Definition: VinciaISR.cc:533
double getTrialScale(int iTrial) const
Get scale for a specific saved trial.
Definition: VinciaISR.h:888
virtual double genQ2thres(double q2old, double sAB, double zMin, double zMax, double colFac, double alphaSvalue, double PDFratio, int idA, int idB, double eA, double eB, bool useMpdf, double headroomFac=1.0, double enhanceFac=1.0) override
Definition: VinciaISR.cc:413
void initAntPtr(AntennaSetISR *antSetIn)
Initialize pointers to antenna sets.
Definition: VinciaISR.h:1054
virtual double trialPDFratio(BeamParticle *beamAPtr, BeamParticle *beamBPtr, int iSys, int idA, int idB, double eA, double eB, double Qt2A, double Qt2B) override
Trial PDF ratio.
Definition: VinciaISR.h:451
virtual double genQ2(double q2old, double sAB, double zMin, double zMax, double colFac, double alphaSvalue, double PDFratio, double eA, double eB, double headroomFac=1.0, double enhanceFac=1.0) override
Generate new Q value, with constant trial alphaS.
Definition: VinciaISR.h:324
virtual string name() override
Name of trial generator.
Definition: VinciaISR.h:306
void addPDF(int iTrial, double pdfRatio)
Add the physical pdf ratio.
Definition: VinciaISR.h:863
virtual double trialPDFratio(BeamParticle *beamAPtr, BeamParticle *beamBPtr, int iSys, int idA, int idB, double eA, double eB, double Qt2A, double Qt2B) override
Trial PDF ratio.
Definition: VinciaISR.cc:669
virtual double getQ2(double saj, double sjk, double sAK) override
Evolution scale.
Definition: VinciaISR.h:614
int getNbranch(int iSys=-1)
Get number of branchings in a system (return -1 if no such system).
Definition: VinciaISR.h:1116
virtual double getQ2(double saj, double sjb, double sAB) override
Evolution scale.
Definition: VinciaISR.h:247
Definition: VinciaWeights.h:17
virtual double getIz(double zMin, double zMax)
The zeta integral.
Definition: VinciaISR.cc:138
bool checkInit()
Check initialization.
Definition: VinciaISR.cc:205
virtual double trialPDFratio(BeamParticle *beamAPtr, BeamParticle *beamBPtr, int iSys, int idA, int idB, double eA, double eB, double Qt2A, double Qt2B) override
Trial PDF ratio.
Definition: VinciaISR.cc:516
virtual double getQ2(double saj, double sjb, double sAB)
Evolution scale.
Definition: VinciaISR.h:52
Definition: Event.h:32
virtual double genQ2(double q2old, double sAB, double zMin, double zMax, double colFac, double alphaSvalue, double PDFratio, double eA, double eB, double headroomFac=1.0, double enhanceFac=1.0) override
Generate new Q value, with constant trial alphaS.
Definition: VinciaISR.cc:389
double getAlphaTrial(int iTrial=-1)
Get alpha S.
Definition: VinciaISR.h:911
virtual double getSj2(double Qt2, double zeta, double sAB) override
Definition: VinciaISR.cc:498
virtual string name() override
Name of trial generator.
Definition: VinciaISR.h:672
virtual double getS1j(double Qt2, double zeta, double sAB) override
Inverse transforms to obtain saj and sjb from Qt2 and zeta.
Definition: VinciaISR.cc:318
virtual double genZ(double zMin, double zMax)
Generate a new zeta value in [zMin,zMax].
Definition: VinciaISR.cc:127
vector< double > getMasses()
Return a vector of the masses.
Definition: VinciaISR.h:1111
A simple class for containing evolution variable definitions.
Definition: VinciaCommon.h:421
virtual double getZmin(double Qt2, double sAB, double eA, double eBeamUsed)
The zeta boundaries, for a given value of the evolution scale.
Definition: VinciaISR.cc:147
virtual double aTrial(double saj, double sjb, double sAB) override
A collinear trial function for initial-initial.
Definition: VinciaISR.cc:219
double getQbranch(int iSys, int iBranch)
Definition: VinciaISR.h:1126
virtual ~VinciaISR()
Destructor.
Definition: VinciaISR.h:996
int getNf(int iRegion)
Number of active flavors.
Definition: VinciaISR.h:1080
virtual double getQ2(double saj, double sjk, double sAK) override
Evolution scale.
Definition: VinciaISR.h:556
The AntennaSetISR class. Simple container of II and IF antenna functions.
Definition: VinciaAntennaFunctions.h:1107
virtual double getQ2(double saj, double sjk, double sAK) override
Evolution scale.
Definition: VinciaISR.h:474
virtual string name() override
Name of trial generator.
Definition: VinciaISR.h:211
virtual double trialPDFratio(BeamParticle *beamAPtr, BeamParticle *beamBPtr, int iSys, int idA, int idB, double eA, double eB, double Qt2A, double Qt2B)
Compute trial PDF ratio.
Definition: VinciaISR.cc:195
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
virtual double getSj2(double Qt2, double zeta, double sAB) override
Definition: VinciaISR.cc:650
B conversion trial function for initial-initial, g -> qqbar.
Definition: VinciaISR.h:415
A gluon collinear trial function for initial-final.
Definition: VinciaISR.h:543
virtual string name() override
Name of trial generator.
Definition: VinciaISR.h:420
bool getHasWeaklyRadiated() override
Tell whether ISR has done a weak emission. Not implemented in Vincia yet.
Definition: VinciaISR.h:1036
bool initUncertainties() override
Definition: VinciaISR.h:1023
Definition: VinciaAntennaFunctions.h:1170
double getEnhanceFac(int iTrial=-1)
Get headroom factor.
Definition: VinciaISR.h:906
bool is1A() const
Is 1 a side A (p+) guy (need to know for pdfs in IF).
Definition: VinciaISR.h:829
The SpaceShower class does spacelike showers.
Definition: SpaceShower.h:33
virtual double getS1j(double Qt2, double zeta, double sAB) override
Inverse transforms to obtain saj and sjb from Qt2 and zeta.
Definition: VinciaISR.cc:635
virtual double genQ2run(double q2old, double sAB, double zMin, double zMax, double colFac, double PDFratio, double b0, double kR, double Lambda, double eA, double eB, double headroomFac=1.0, double enhanceFac=1.0)
Generate new Q value, with first-order running alphaS.
Definition: VinciaISR.cc:70
virtual double aTrial(double saj, double sjb, double sAB) override
A splitting trial function for initial-initial, q -> gqbar.
Definition: VinciaISR.cc:352
double mbSav
Masses.
Definition: VinciaISR.h:120
Class for performing QED showers.
Definition: VinciaQED.h:475
virtual double genQ2(double q2old, double sAB, double zMin, double zMax, double colFac, double alphaSvalue, double PDFratio, double eA, double eB, double headroomFac=1.0, double enhanceFac=1.0) override
Generate a new Q value, with constant trial alphaS.
Definition: VinciaISR.h:438
double getPDFratioTrial(int iTrial=-1)
Get pdf ratio.
Definition: VinciaISR.h:916
virtual double aTrial(double saj, double sjb, double sAB) override
Trial antenna function.
Definition: VinciaISR.h:309
Definition: Settings.h:195
double getColFac(int iTrial=-1)
Get colour factor.
Definition: VinciaISR.h:896
virtual void init(double mcIn, double mbIn)
Initialize.
Definition: VinciaISR.cc:32
AlphaStrong * alphaSptr
Pointer to global AlphaStrong instance.
Definition: VinciaISR.h:1142
virtual double genQ2(double q2old, double sAB, double zMin, double zMax, double colFac, double alphaSvalue, double PDFratio, double eA, double eB, double headroomFac=1.0, double enhanceFac=1.0) override
Generate a new Q value, with constant trial alphaS.
Definition: VinciaISR.cc:569
double getQmin(int iRegion)
Evolution window, phase space region boundaries.
Definition: VinciaISR.h:1074
bool getIsSwapped(int iTrial=-1) const
Get whether physical antenna is swapped.
Definition: VinciaISR.h:878