PYTHIA  8.315
VinciaEW.h
1 // VinciaEW.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2025 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 the EW antenna-shower class and auxiliary
7 // classes. Main author is Rob Verheyen.
8 
9 #ifndef Pythia8_VinciaEW_H
10 #define Pythia8_VinciaEW_H
11 
12 // PYTHIA 8 headers.
13 #include "Pythia8/Event.h"
14 #include "Pythia8/StandardModel.h"
15 #include "Pythia8/PartonSystems.h"
16 #include "Pythia8/PythiaComplex.h"
17 #include "Pythia8/BeamParticle.h"
18 #include "Pythia8/UserHooks.h"
19 
20 // VINCIA functionality.
21 #include "Pythia8/VinciaCommon.h"
22 #include "Pythia8/VinciaQED.h"
23 
24 namespace Pythia8 {
25 
26 //==========================================================================
27 
28 // Simple class to save information about particles.
29 
30 class EWParticle {
31 
32  public:
33 
34  // Constructor.
35  EWParticle() = default;
36  EWParticle(double massIn, double widthIn, bool isResIn) :
37  mass(massIn), width(widthIn), isRes(isResIn) {};
38 
39  // Members.
40  double mass{0.}, width{0.};
41  bool isRes{false};
42 
43 };
44 
45 //==========================================================================
46 
47 // Locally saved particle data in glorified map.
48 
50 
51  public:
52 
53  // Find particle data.
54  bool find(int id, int pol){
55  return data.find(make_pair(id, pol)) != data.end();}
56 
57  // Add particle data.
58  void add(int id, int pol, double massIn, double widthIn, bool isResIn) {
59  data[make_pair(id, pol)] = EWParticle(massIn, widthIn, isResIn);}
60 
61  // Return particle mass.
62  double mass(int id, int pol) {
63  return find(id, pol) ? data[make_pair(id, pol)].mass : 0.;}
64 
65  // Return particle mass.
66  double mass(int id) {
67  // Every particle has either pol = 1 or pol = 0.
68  if (find(id, 1)) return data[make_pair(id, 1)].mass;
69  return find(id, 0) ? data[make_pair(id, 0)].mass : 0.;
70  }
71 
72  // Return particle width.
73  double width(int id, int pol) {
74  return find(id, pol) ? data[make_pair(id, pol)].width : 0.;}
75 
76  // Return if particle is resonance.
77  bool isRes(int id, int pol) {
78  return find(id, pol) && data[make_pair(id, pol)].isRes;}
79 
80  // Return if particle is resonance.
81  bool isRes(int id) {
82  // Every particle has either pol = 1 or pol = 0.
83  if (find(id, 1)) return data[make_pair(id, 1)].isRes;
84  else if (find(id, 0)) return data[make_pair(id, 0)].isRes;
85  else return false;
86  }
87 
88  // Data access.
89  EWParticle* get(int id, int pol) { return &data.at(make_pair(id, pol)); }
90  unordered_map<pair<int, int>, EWParticle>::iterator begin() {
91  return data.begin();}
92  unordered_map<pair<int, int>, EWParticle>::iterator end() {
93  return data.end();}
94 
95  // Member.
96  unordered_map<pair<int,int>, EWParticle> data;
97 
98 };
99 
100 //==========================================================================
101 
102 // Class to contain an antenna function and two outgoing polarizations.
103 
104 class AntWrapper {
105 
106 public:
107 
108  // Constructor.
109  AntWrapper(double valIn, int poliIn, int poljIn):
110  val(valIn), poli(poliIn), polj(poljIn) {}
111 
112  // Print.
113  void print() {cout << "(" << poli << ", " << polj << ") " << val;}
114 
115  // Members.
116  double val;
117  int poli, polj;
118 
119 };
120 
121 //==========================================================================
122 
123 // Class to contain an amplitude and two outgoing polarizations.
124 
125 class AmpWrapper {
126 
127  public:
128 
129  // Constructor.
130  AmpWrapper(complex valIn, int poliIn, int poljIn):
131  val(valIn), poli(poliIn), polj(poljIn) {}
132 
133  // Normalization.
134  AntWrapper norm() {return AntWrapper(std::norm(val), poli, polj);}
135 
136  // Operators.
137  AmpWrapper& operator+=(complex c) {this->val += c; return *this;}
138  AmpWrapper& operator*=(complex c) {this->val *= c; return *this;}
139  void print() {cout << "(" << poli << ", " << polj << ") " << val;}
140 
141  // Members.
143  int poli, polj;
144 
145 };
146 
147 //==========================================================================
148 
149 // Calculator class for amplitudes, antennae, and Breit-Wigners.
150 
152 
153 public:
154 
155  // Initialize the pointers.
156  void initPtr(Info* infoPtrIn, AlphaEM* alphaPtrIn,
157  AlphaStrong* alphaSPtrIn) {
158  infoPtr = infoPtrIn;
159  partonSystemsPtr = infoPtr->partonSystemsPtr;
160  rndmPtr = infoPtr->rndmPtr;
161  settingsPtr = infoPtr->settingsPtr;
162  loggerPtr = infoPtr->loggerPtr;
163  alphaPtr = alphaPtrIn;
164  alphaSPtr = alphaSPtrIn;
165  isInitPtr = true;
166  }
167 
168  // Initialize with maps.
169  void init(EWParticleData* dataIn, unordered_map< pair<int, int>,
170  vector<pair<int, int> > >* cluMapFinalIn,
171  unordered_map< pair<int, int>, vector<pair<int, int> > >*
172  cluMapInitialIn);
173 
174  // Set verbosity level.
175  void setVerbose(int verboseIn) {verbose = verboseIn;}
176 
177  // Spinor products used for the calculation of the amplitudes.
178  complex spinProd(int pol, const Vec4& ka, const Vec4& kb);
179  complex spinProd(int pol, const Vec4& ka, const Vec4& pa, const Vec4& kb);
180  complex spinProd(int pol, const Vec4& ka, const Vec4& pa, const Vec4& pb,
181  const Vec4& kb);
182  complex spinProd(int pol, const Vec4& ka, const Vec4& pa, const Vec4& pb,
183  const Vec4& pc, const Vec4& kb);
184  complex spinProd(int pol, const Vec4& ka, const Vec4& pa, const Vec4& pb,
185  const Vec4& pc, const Vec4& pd, const Vec4& kb);
186  Vec4 spinProdFlat(string method, const Vec4& ka, const Vec4& pa);
187 
188  // Initialize couplings.
189  void initCoup(bool va, int id1, int id2, int pol, bool m);
190 
191  // Initialize an FSR branching amplitude.
192  void initFSRAmp(bool va, int id1, int id2, int pol,
193  const Vec4& pi, const Vec4 &pj, const double& mMot, const double& widthQ2);
194 
195  // Check for zero denominator in an FSR amplitude.
196  bool zdenFSRAmp(const string& method, const Vec4& pi, const Vec4& pj,
197  bool check);
198 
199  // Initialize an ISR branching amplitude.
200  void initISRAmp(bool va, int id1, int id2, int pol,
201  const Vec4& pa, const Vec4 &pj, double& mA);
202 
203  // Check for zero denominator in an ISR amplitude.
204  bool zdenISRAmp(const string& method, const Vec4& pa, const Vec4& pj,
205  bool check);
206 
207  // Branching amplitude formalism.
208  // Naming scheme: f = fermion, v = vector boson, h = higgs.
209 
210  // Final-state branching amplitudes.
211  complex ftofvFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
212  int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
213  complex ftofhFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
214  int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
215  complex fbartofbarvFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
216  int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
217  complex fbartofbarhFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
218  int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
219  complex vTtoffbarFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
220  int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
221  complex vTtovhFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
222  int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
223  complex vTtovvFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
224  int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
225  complex vLtoffbarFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
226  int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
227  complex vLtovhFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
228  int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
229  complex vLtovvFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
230  int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
231  complex htoffbarFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
232  int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
233  complex htovvFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
234  int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
235  complex htohhFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
236  int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
237 
238  // Initial-state branching amplitudes.
239  complex ftofvISRAmp(const Vec4& pa, const Vec4& pj, int idA, int ida,
240  int idj, double mA, int polA, int pola, int polj);
241  complex ftofhISRAmp(const Vec4& pa, const Vec4& pj, int idA, int ida,
242  int idj, double mA, int polA, int pola, int polj);
243  complex fbartofbarvISRAmp(const Vec4& pa, const Vec4& pj, int idA, int ida,
244  int idj, double mA, int polA, int pola, int polj);
245  complex fbartofbarhISRAmp(const Vec4& pa, const Vec4& pj, int idA, int ida,
246  int idj, double mA, int polA, int pola, int polj);
247 
248  // Branching amplitude selector.
249  complex branchAmpFSR(const Vec4& pi, const Vec4& pj, int idMot, int idi,
250  int idj, double mMot, double widthQ2, int polMot, int poli=9, int polj=9);
251  complex branchAmpISR(const Vec4& pa, const Vec4& pj, int idA, int ida,
252  int idj, double mA, int polA, int pola=9, int polj=9);
253 
254  // Compute FF antenna function from amplitudes.
255  double branchKernelFF(const Vec4& pi, const Vec4& pj, int idMot, int idi,
256  int idj, double mMot, double widthQ2, int polMot, int poli, int polj) {
257  return norm(branchAmpFSR(pi, pj, idMot, idi, idj, mMot, widthQ2,
258  polMot, poli, polj));}
259 
260  // Compute FF antenna functions from amplitudes for all polarizations.
261  vector<AntWrapper> branchKernelFF(Vec4 pi, Vec4 pj, int idMot, int idi,
262  int idj, double mMot, double widthQ2, int polMot);
263 
264  // Compute II antenna function from amplitudes.
265  double branchKernelII(Vec4 pa, Vec4 pj, int idA, int ida,
266  int idj, double mA, int polA, int pola, int polj) {
267  return norm(branchAmpISR(pa, pj, idA, ida, idj, mA, polA, pola, polj));}
268 
269  // Compute II antenna functions from amplitudes for all polarizations.
270  vector<AntWrapper> branchKernelII(Vec4 pa, Vec4 pj, int idA, int ida,
271  int idj, double mA, int polA);
272 
273  // Initialize an FF antenna function.
274  void initFFAnt(bool va, int id1, int id2, int pol, const double& Q2,
275  const double& widthQ2, const double& xi, const double& xj,
276  const double& mMot, const double& miIn, const double& mjIn);
277 
278  // Report helicity combination error for an FF antenna function.
279  void hmsgFFAnt(int polMot, int poli, int polj) {
280  stringstream ss;
281  ss << "helicity combination was not found:\n "
282  << "polMot = " << polMot << " poli = " << poli << " polj = " << polj;
283  loggerPtr->ERROR_MSG(ss.str());}
284 
285  // Initialize an II antenna function.
286  void initIIAnt(int id1, int id2, int pol, const double& Q2,
287  const double& xA, const double& xj,
288  const double& mA, const double& maIn, const double& mjIn);
289 
290  // Report helicity combination error for an II antenna function.
291  void hmsgIIAnt(int polA, int pola, int polj) {
292  stringstream ss;
293  ss << "helicity combination was not found:\n "
294  << "polA = " << polA << " pola = " << pola << " polj = " << polj;
295  loggerPtr->ERROR_MSG(ss.str());}
296 
297  // FF Antenna functions for branching process I (K) -> i j (k).
298  // Q2 is the offshellness of I.
299  // xi and xj are energy fractions of pi and pj in the collinear limit.
300  double ftofvFFAnt(double Q2, double widthQ2, double xi, double xj,
301  int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
302  int polMot, int poli, int polj);
303  double ftofhFFAnt(double Q2, double widthQ2, double xi, double xj,
304  int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
305  int polMot, int poli, int polj);
306  double fbartofbarvFFAnt(double Q2, double widthQ2, double xi, double xj,
307  int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
308  int polMot, int poli, int polj);
309  double fbartofbarhFFAnt(double Q2, double widthQ2, double xi, double xj,
310  int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
311  int polMot, int poli, int polj);
312  double vtoffbarFFAnt(double Q2, double widthQ2, double xi, double xj,
313  int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
314  int polMot, int poli, int polj);
315  double vtovhFFAnt(double Q2, double widthQ2, double xi, double xj,
316  int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
317  int polMot, int poli, int polj);
318  double vtovvFFAnt(double Q2, double widthQ2, double xi, double xj,
319  int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
320  int polMot, int poli, int polj);
321  double htoffbarFFAnt(double Q2, double widthQ2, double xi, double xj,
322  int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
323  int polMot, int poli, int polj);
324  double htovvFFAnt(double Q2, double widthQ2, double xi, double xj,
325  int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
326  int polMot, int poli, int polj);
327  double htohhFFAnt(double Q2, double widthQ2, double xi, double xj,
328  int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
329  int polMot, int poli, int polj);
330 
331  // II Antenna functions for branching process A (B) -> a j (b).
332  // Q2 is the offshellness of A.
333  // xA and xj are energy fractions of pA and pj in the collinear limit.
334  double ftofvIIAnt(double Q2, double xA, double xj,
335  int idA, int ida, int idj, double mA, double maIn, double mjIn,
336  int polA, int pola, int polj);
337  double fbartofbarvIIAnt(double Q2, double xA, double xj,
338  int idA, int ida, int idj, double mA, double maIn, double mjIn,
339  int polA, int pola, int polj);
340 
341  // FF antenna function calculator.
342  double antFuncFF(double Q2, double widthQ2, double xi, double xj,
343  int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
344  int polMot, int poli, int polj);
345 
346  // FF antenna function calculator for all outgoing polarizations.
347  vector<AntWrapper> antFuncFF(double Q2, double widthQ2, double xi,
348  double xj, int idMot, int idi, int idj, double mMot, double miIn,
349  double mjIn, int polMot);
350 
351  // II antenna function calculator.
352  double antFuncII(double Q2, double xA, double xj,
353  int idA, int ida, int idj, double mA, double maIn, double mjIn,
354  int polA, int pola, int polj);
355 
356  // II antenna function calculator for all outgoing polarizations.
357  vector<AntWrapper> antFuncII(double Q2, double xA, double xj,
358  int idA, int ida, int idj, double mA, double maIn, double mjIn, int polA);
359 
360  // Initialize an FSR splitting kernel.
361  void initFSRSplit(bool va, int id1, int id2, int pol,
362  const double& mMot, const double& miIn, const double& mjIn) {
363  mi = miIn; mj = mjIn; mMot2 = pow2(mMot); mi2 = pow2(mi); mj2 = pow2(mj);
364  initCoup(va, id1, id2, pol, true);}
365 
366  // Check for zero denominator in an FSR splitting kernel.
367  bool zdenFSRSplit(const string& method, const double& Q2, const double& z,
368  bool check);
369 
370  // Report helicty combination error for an FSR splitting kernel.
371  void hmsgFSRSplit(int polMot, int poli, int polj) {
372  stringstream ss;
373  ss << "helicity combination was not found:\n "
374  << "polMot = " << polMot << " poli = " << poli << " polj = " << polj;
375  loggerPtr->ERROR_MSG(ss.str());}
376 
377  // Initialize an ISR splitting kernel.
378  void initISRSplit(bool va, int id1, int id2, int pol,
379  const double& mA, const double& maIn, const double& mjIn) {
380  ma = maIn; mj = mjIn; mA2 = pow2(mA); ma2 = pow2(ma); mj2 = pow2(mj);
381  initCoup(va, id1, id2, pol, mA > VinciaConstants::NANO);}
382 
383  // Check for zero denominator in an ISR splitting kernel.
384  bool zdenISRSplit(const string& method, const double& Q2, const double& z,
385  bool flip, bool check);
386 
387  // Report helicty combination error for an ISR splitting kernel.
388  void hmsgISRSplit(int polA, int pola, int polj) {
389  stringstream ss;
390  ss << "helicity combination was not found:\n "
391  << "polA = " << polA << " pola = " << pola << " polj = " << polj;
392  loggerPtr->ERROR_MSG(ss.str());}
393 
394  // Final-state splitting kernels.
395  double ftofvFSRSplit(double Q2, double z, int idMot, int idi, int idj,
396  double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
397  double ftofhFSRSplit(double Q2, double z, int idMot, int idi, int idj,
398  double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
399  double fbartofbarvFSRSplit(double Q2, double z, int idMot, int idi, int idj,
400  double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
401  double fbartofbarhFSRSplit(double Q2, double z, int idMot, int idi, int idj,
402  double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
403  double vTtoffbarFSRSplit(double Q2, double z, int idMot, int idi, int idj,
404  double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
405  double vTtovhFSRSplit(double Q2, double z, int idMot, int idi, int idj,
406  double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
407  double vTtovvFSRSplit(double Q2, double z, int idMot, int idi, int idj,
408  double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
409  double vLtoffbarFSRSplit(double Q2, double z, int idMot, int idi, int idj,
410  double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
411  double vLtovhFSRSplit(double Q2, double z, int idMot, int idi, int idj,
412  double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
413  double vLtovvFSRSplit(double Q2, double z, int idMot, int idi, int idj,
414  double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
415  double htoffbarFSRSplit(double Q2, double z, int idMot, int idi, int idj,
416  double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
417  double htovvFSRSplit(double Q2, double z, int idMot, int idi, int idj,
418  double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
419  double htohhFSRSplit(double Q2, double z, int idMot, int idi, int idj,
420  double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
421 
422  // Initial-state splitting kernels.
423  double ftofvISRSplit(double Q2, double z, int idA, int ida, int idj,
424  double mA, double maIn, double mjIn, int polA, int pola, int polj);
425  double ftofhISRSplit(double Q2, double z, int idA, int ida, int idj,
426  double mA, double maIn, double mjIn, int polA, int pola, int polj);
427  double fbartofbarvISRSplit(double Q2, double z, int idA, int ida, int idj,
428  double mA, double maIn, double mjIn, int polA, int pola, int polj);
429  double fbartofbarhISRSplit(double Q2, double z, int idA, int ida, int idj,
430  double mA, double maIn, double mjIn, int polA, int pola, int polj);
431 
432  // Splitting kernel caller.
433  double splitFuncFSR(double Q2, double z, int idMot, int idi, int idj,
434  double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
435  double splitFuncISR(double Q2, double z, int idA, int ida, int idj,
436  double mA, double maIn, double mjIn, int polA, int pola, int polj);
437 
438  // Compute partial decay width.
439  double getPartialWidth(int idMot, int idi, int idj, double mMot, int polMot);
440 
441  // Compute total decay width.
442  double getTotalWidth(int idMot, double mMot, int polMot);
443 
444  // Breit-Wigner calculators.
445  double getBreitWigner(int id, double m, int pol);
446  double getBreitWignerOverestimate(int id, double m, int pol);
447 
448  // Generate Breit-Wigner mass.
449  double sampleMass(int id, int pol);
450 
451  // Bosonic interference factor.
452  void applyBosonInterferenceFactor(Event &event, int XYEv, Vec4 pi, Vec4 pj,
453  int idi, int idj, int poli, int polj);
454 
455  // Polarise a resonance decay.
456  bool polarise(vector<Particle> &state);
457 
458  // EW event weight.
459  double eventWeight() {return eventWeightSave;}
460  void eventWeight(double eventWeightIn) {eventWeightSave = eventWeightIn;}
461 
462  // Public data members.
463 
464  // EW data.
465  EWParticleData* dataPtr{};
466 
467  // Maps of coupling constants.
468  // TODO: read this from data file (rather than hard code).
469  unordered_map<pair<int, int>, double> vMap, aMap, gMap, vCKM;
470 
471 private:
472 
473  // Data members.
474 
475  // Constants used for Breit-Wigner sampling.
476  unordered_map<int, vector<double> > cBW;
477  // BW normalizations.
478  unordered_map<pair<int,int>, double> nBW, np;
479 
480  // Event weight.
481  double eventWeightSave{1};
482 
483  // Electroweak constants.
484  double mw, mw2, sw, sw2;
485 
486  // Mode of matching to the Breit-Wigner.
487  int bwMatchMode;
488 
489  // Maps of EW clusterings.
490  unordered_map<pair<int, int>, vector<pair<int, int> > >* cluMapFinal{};
491  unordered_map<pair<int, int>, vector<pair<int, int> > >* cluMapInitial{};
492 
493  // Vectors with spin assignments.
494  vector<int> fermionPols, vectorPols, scalarPols;
495 
496  // Couplings.
497  double v, a, vPls, vMin, g;
498 
499  // Masses and scale (FSR and ISR).
500  double mMot2, mi, mi2, mj, mj2, mA2, ma, ma2, isrQ2;
501  complex M, fsrQ2;
502 
503  // Reference vectors (FSR and ISR).
504  Vec4 kij, ki, kj, pij, kaj, ka, paj;
505  double wij, wi, wj, wij2, wi2, wj2, waj, wa, waj2, wa2;
506 
507  // Antenna function members.
508  double Q4, Q4gam, Q2til, ant;
509 
510  // Pointers.
511  Info* infoPtr{};
512  PartonSystems* partonSystemsPtr{};
513  Rndm* rndmPtr{};
514  Settings* settingsPtr{};
515  Logger* loggerPtr{};
516  AlphaEM* alphaPtr{};
517  AlphaStrong* alphaSPtr{};
518 
519  // Initializations.
520  bool isInit{false};
521  bool isInitPtr{false};
522  int verbose{0};
523 };
524 
525 //==========================================================================
526 
527 // Class that contains an electroweak branching.
528 
529 class EWBranching {
530 
531  public:
532 
533  // Constructor.
534  EWBranching(int idMotIn, int idiIn, int idjIn, int polMotIn,
535  double c0In = 0, double c1In = 0, double c2In = 0, double c3In = 0):
536  idMot(idMotIn), idi(idiIn), idj(idjIn), polMot(polMotIn), c0(c0In),
537  c1(c1In), c2(c2In), c3(c3In),
538  isSplitToFermions(abs(idMot) > 20 && abs(idi) < 20 && abs(idj) < 20) {;}
539 
540  // Print.
541  void print() {cout <<" (" << idMot << ", " << polMot << ") -> " << idi <<
542  "," << idj << ": (" << c0 << ", " << c1 << ", " << c2 << ", " << c3 <<
543  ") \n";}
544 
545  // For a branching I->ij, particle IDs and polarisation of I.
546  int idMot, idi, idj, polMot;
547  // Overestimate constants used to sample antennae.
548  double c0, c1, c2, c3;
549  // Store if branching is a splitting.
551 
552 };
553 
554 //==========================================================================
555 
556 // Base class for an electroweak antenna.
557 
558 class EWAntenna {
559 
560  public:
561 
562  // Constructor and destructor.
563  EWAntenna(): iSys(-1), shat(0), doBosonInterference(false) {};
564  virtual ~EWAntenna() = default;
565 
566  // Print, must be implemented by base classes.
567  void print() {
568  stringstream ss;
569  ss << "Brancher = (" << iMot << ", " << polMot
570  << "), Recoiler = " << iRec;
571  printOut(__METHOD_NAME__, ss.str());
572  for (int i = 0; i < (int)brVec.size(); i++) brVec[i].print();}
573 
574  // Initialize pointers.
575  void initPtr(Info* infoPtrIn, VinciaCommon* vinComPtrIn, AlphaEM* alphaPtrIn,
576  AmpCalculator* ampCalcPtrIn) {
577  infoPtr = infoPtrIn;
578  rndmPtr = infoPtr->rndmPtr;
579  loggerPtr = infoPtr->loggerPtr;
580  partonSystemsPtr = infoPtr->partonSystemsPtr;
581  vinComPtr = vinComPtrIn;
582  alphaPtr = alphaPtrIn;
583  ampCalcPtr = ampCalcPtrIn;
584  };
585 
586  // Initialize, must be implemented by base classes.
587  virtual bool init(Event &event, int iMotIn, int iRecIn, int iSysIn,
588  vector<EWBranching> &branchings, Settings* settingsPtr) = 0;
589 
590  // Set verbosity level.
591  void setVerbose(int verboseIn) {verbose = verboseIn;}
592 
593  // Generate a trial scale to be compared against other QCD and EW trials.
594  virtual double generateTrial(double q2Start, double q2End,
595  double alphaIn) = 0;
596 
597  // Accept/reject step to accept a trial.
598  virtual bool acceptTrial(Event& event) = 0;
599 
600  // Update an event and parton system.
601  virtual void updateEvent(Event& event) = 0;
602  virtual void updatePartonSystems(Event& event);
603 
604  // Return index.
605  int getIndexMot() {return iMot;};
606  int getIndexRec() {return iRec;};
607 
608  // Check if splitting to fermions, inital, or resoance.
610  return brTrial != nullptr && brTrial->isSplitToFermions;}
611  virtual bool isInitial() {return false;}
612  virtual bool isResonanceDecay() {return false;}
613 
614  // Select a channel.
615  bool selectChannel(int idx, const double& cSum, const map<double, int>&
616  cSumSoFar, int& idi, int& idj, double& mi2, double& mj2);
617 
618 protected:
619 
620  // Indices, PID, and polarization of I, K in Pythia event record.
621  int iMot, iRec, idMot, idRec, polMot;
622  // Momenta of antenna constituents.
623  Vec4 pMot, pRec;
624  // Masses and invariants of antenna constituents.
625  double sAnt, mMot, mMot2, mRec, mRec2;
626  // Overestimate of QED coupling.
627  double alpha;
628  // Parton system this antenna is part of.
629  int iSys;
630 
631  // EW branchings.
632  vector<EWBranching> brVec;
633 
634  // Trial variables.
635  bool hasTrial;
636  double q2Trial, sijTrial, sjkTrial;
637  int poliTrial, poljTrial;
638 
639  // Outgoing momenta after branching.
640  vector<Vec4> pNew;
641 
642  // Info on coefficents.
643  double c0Sum, c1Sum, c2Sum, c3Sum;
644  map<double, int> c0SumSoFar, c1SumSoFar, c2SumSoFar, c3SumSoFar;
645 
646  // Matching scale.
647  double q2Match;
648 
649  // Information for partonSystems.
650  int jNew;
651  unordered_map<int,int> iReplace;
652  double shat;
653 
654  // Pointers.
655  EWBranching* brTrial{};
656  Info* infoPtr{};
657  Rndm* rndmPtr{};
658  Logger* loggerPtr{};
659  PartonSystems* partonSystemsPtr{};
660  VinciaCommon* vinComPtr{};
661  AlphaEM* alphaPtr{};
662  AmpCalculator* ampCalcPtr{};
663 
664  // Settings.
666 
667  // Verbosity.
668  int verbose;
669 
670 };
671 
672 //==========================================================================
673 
674 // Final-final electroweak antenna.
675 
676 class EWAntennaFF : public EWAntenna {
677  public:
678 
679  // Overridden virtual functions.
680  virtual bool init(Event &event, int iMotIn, int iRecIn, int iSysIn,
681  vector<EWBranching> &branchings, Settings* settingsPtr) override;
682  virtual double generateTrial(double q2Start, double q2End, double alphaIn)
683  override;
684  virtual bool acceptTrial(Event &event) override;
685  virtual void updateEvent(Event &event) override;
686 
687  private:
688 
689  // Data members.
690  double mAnt2, sqrtKallen;
691  // Kinematic map.
692  int kMapFinal;
693  // Controls if resonances with too large offshellness are vetoed.
694  bool vetoResonanceProduction;
695 
696 };
697 
698 //==========================================================================
699 
700 // Final-final electroweak resonance antenna.
701 
702 class EWAntennaFFres : public EWAntennaFF {
703 
704 public:
705 
706  // Overridden virtual functions.
707  bool init(Event &event, int iMotIn, int iRecIn, int iSysIn,
708  vector<EWBranching> &branchings, Settings* settingsPtr) override;
709  bool isResonanceDecay() override {return true;}
710  double generateTrial(double q2Start, double q2End, double alphaIn) override;
711  bool acceptTrial(Event &event) override;
712  void updateEvent(Event &event) override;
713 
714  // Generate the kinematics and channel for decays below matching scale.
715  bool genForceDecay(Event &event);
716 
717 private:
718 
719  // Check if the trial was a resonance decay.
720  bool trialIsResDecay;
721  // Matching mode.
722  int bwMatchMode;
723  // Offshellness of the resonance and EW scale.
724  double q2Dec, q2EW;
725  // Switch for the special case of a resonance without a recoiler.
726  bool doDecayOnly{false};
727 
728 };
729 
730 //==========================================================================
731 
732 // Initial-initial electroweak antenna.
733 
734 class EWAntennaII : public EWAntenna {
735 
736 public:
737 
738  // Constructor.
739  EWAntennaII(BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn):
740  beamAPtr(beamAPtrIn), beamBPtr(beamBPtrIn), shh(0), xMot(0), xRec(0),
741  vetoResonanceProduction(false), TINYPDFtrial(1e-10) {;}
742 
743  // Overridden virtual functions.
744  bool init(Event &event, int iMotIn, int iRecIn, int iSysIn,
745  vector<EWBranching> &branchings,Settings* settingsPtr) override;
746  bool isInitial() override {return true;}
747  double generateTrial(double q2Start, double q2End, double alphaIn) override;
748  bool acceptTrial(Event &event) override;
749  void updateEvent(Event &event) override;
750  void updatePartonSystems(Event &event) override;
751 
752 private:
753 
754  // Members.
755  // Beam pointers.
756  BeamParticle* beamAPtr{};
757  BeamParticle* beamBPtr{};
758  // Hadronic invariant mass.
759  double shh;
760  // Antenna hadronic momentum fractions.
761  double xMot, xRec;
762  // Controls if resonances with too large offshellness are vetoed.
763  bool vetoResonanceProduction;
764  // Global recoil momenta.
765  vector<Vec4> pRecVec;
766  vector<int> iRecVec;
767 
768  // Tolerance for PDFs.
769  double TINYPDFtrial;
770 
771 };
772 
773 //==========================================================================
774 
775 // Class that performs electroweak showers in a single parton system.
776 
777 class EWSystem {
778 
779 public:
780 
781  // Constructors.
782  EWSystem(): antTrial(nullptr) {clearLastTrial();}
783  EWSystem(unordered_map<pair<int, int>, vector<EWBranching> >* brMapFinalIn,
784  unordered_map<pair<int, int>, vector<EWBranching> >* brMapInitialIn,
785  unordered_map<pair<int, int>, vector<EWBranching> >* brMapResonanceIn,
786  unordered_map<pair<int, int>, vector<pair<int, int> > >* cluMapFinalIn,
787  unordered_map<pair<int, int>, vector<pair<int, int> > >* cluMapInitialIn,
788  AmpCalculator * ampCalcIn):
789  shh(0), iSysSav(0), resDecOnlySav(false), q2Cut(0), q2Trial(0),
790  lastWasSplitSav(false), lastWasDecSav(false), lastWasInitialSav(false),
791  lastWasBelowCut(false), ISav(0), KSav(0), brMapFinal(brMapFinalIn),
792  brMapInitial(brMapInitialIn), brMapResonance(brMapResonanceIn),
793  cluMapFinal(cluMapFinalIn), cluMapInitial(cluMapInitialIn),
794  ampCalcPtr(ampCalcIn), isInit(false), doVetoHardEmissions(false),
795  verbose(false), vetoHardEmissionsDeltaR2(0) {clearLastTrial();}
796 
797  // Initialize pointers.
798  void initPtr(Info* infoPtrIn, VinciaCommon* vinComPtrIn, AlphaEM* alIn) {
799  infoPtr = infoPtrIn;
800  partonSystemsPtr = infoPtr->partonSystemsPtr;
801  rndmPtr = infoPtr->rndmPtr;
802  settingsPtr = infoPtr->settingsPtr;
803  loggerPtr = infoPtr->loggerPtr;
804  vinComPtr = vinComPtrIn;
805  al = alIn;}
806 
807  // Initialize.
808  void init(BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn) {
809  beamAPtr = beamAPtrIn; beamBPtr = beamBPtrIn;
810  doVetoHardEmissions = settingsPtr->flag("Vincia:EWoverlapVeto");
811  vetoHardEmissionsDeltaR2 =
812  pow2(settingsPtr->parm("Vincia:EWoverlapVetoDeltaR"));
813  isInit = true;}
814 
815  // Set verbosity.
816  void setVerbose(int verboseIn) {verbose = verboseIn;}
817 
818  // Prepare an event.
819  bool prepare(Event &event, int iSysIn, double q2CutIn, bool resDecOnlyIn) {
820  iSysSav = iSysIn; resDecOnlySav = resDecOnlyIn; q2Cut = q2CutIn;
821  shh = infoPtr->s(); return buildSystem(event);}
822 
823  // Build a system.
824  bool buildSystem(Event &event);
825 
826  // Return which system we are handling.
827  int system() {return iSysSav;}
828 
829  // Generate the next Q2.
830  double q2Next(double q2Start, double q2End);
831 
832  // Add an antenna.
833  template <class T> void addAntenna(T ant, vector<T>& antVec,
834  Event &event, int iMot, int iRec,
835  unordered_map<pair<int, int>, vector<EWBranching> >* brMapPtr) {
836  if (iMot == 0) return; // Check mother.
837  int idA(event[iMot].id()), polA(event[iMot].pol());
838  if (idA == 21) return; // Skip gluons.
839  auto it = brMapPtr->find(make_pair(idA, polA));
840  if (it != brMapPtr->end()) {
841  // Found. Pass verbosity.
842  ant.setVerbose(verbose);
843  // Pass pointers.
844  ant.initPtr(infoPtr, vinComPtr, al, ampCalcPtr);
845  // Initialise and if success, store.
846  if (ant.init(event, iMot, iRec, iSysSav, it->second, settingsPtr)) {
847  antVec.push_back(std::move(ant));
848  if (verbose >= VinciaConstants::DEBUG) {
849  stringstream ss;
850  ss << "Added EW antenna with iEv = "
851  << iMot << " and iRec = "<< iRec<< " in system "<< iSysSav;
852  printOut(__METHOD_NAME__, ss.str());}}}}
853 
854  // Generate a trial.
855  template <class T> void generateTrial(vector<T> & antVec, double q2Start,
856  double q2End, double alphaIn) {
857  if (q2End > q2Start) return;
858  for (int i = 0; i < (int)antVec.size(); i++) {
859  // Generate a trial scale for the current antenna.
860  double q2New = antVec[i].generateTrial(q2Start,q2End,alphaIn);
861  // Current winner.
862  if (q2New > q2Trial && q2New > q2End) {
863  // Save trial information.
864  q2Trial = q2New;
865  antTrial = &(antVec[i]);
866  lastWasDecSav = antTrial->isResonanceDecay();
867  lastWasInitialSav = antTrial->isInitial();
868  // This is done to avoid issues with resonance antennae
869  // not deciding their channel until acceptTrial.
870  lastWasSplitSav = lastWasDecSav ? true : antTrial->isSplitToFermions();
871  lastWasBelowCut = (q2Trial < q2Cut)? true : false;
872  ISav = antTrial->getIndexMot(); KSav = antTrial->getIndexRec();}}}
873 
874  // Overloaded version passing event, just for resonance decays.
875  void generateTrial(Event &event, vector<EWAntenna> & antVec,double q2Start,
876  double q2End, double alphaIn);
877 
878  // Accept a trial.
879  bool acceptTrial(Event &event) {
880  bool passed = antTrial->acceptTrial(event);
881  if (verbose >= VinciaConstants::DEBUG)
882  printOut(__METHOD_NAME__, passed ? "Passed veto" : "Vetoed branching");
883  return passed;}
884 
885  // Update an event.
886  void updateEvent(Event &event) {
887  if (verbose >= VinciaConstants::DEBUG)
888  printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
889  if (antTrial != nullptr) antTrial->updateEvent(event);
890  else loggerPtr->ERROR_MSG("trial doesn't exist!");
891  if (verbose >= VinciaConstants::DEBUG)
892  printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);}
893 
894  // Update parton systems.
895  void updatePartonSystems(Event &event) {
896  if (verbose >= VinciaConstants::DEBUG)
897  printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
898  if (antTrial!=nullptr) antTrial->updatePartonSystems(event);
899  else loggerPtr->ERROR_MSG("trial doesn't exist!");
900  if (verbose >= VinciaConstants::DEBUG)
901  printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);}
902 
903  // Print the antennas.
904  void printAntennae() {
905  for (int i = 0; i < (int)antVecFinal.size(); i++) antVecFinal[i].print();
906  for (int i = 0; i < (int)antVecRes.size(); i++) antVecRes[i].print();
907  for (int i = 0; i < (int)antVecInitial.size(); i++)
908  antVecInitial[i].print();}
909 
910  // Get number of branchers (total and res-dec only)
911  unsigned int nBranchers() {return antVecFinal.size() +
912  antVecInitial.size() + antVecRes.size();}
913  unsigned int nResDec() {return antVecRes.size();}
914 
915  // Check if system has a trial pointer.
916  bool hasTrial() {return antTrial != nullptr;}
917 
918  // Check if split to fermions, initial, or resonance.
919  bool lastWasSplitToFermions() {return lastWasSplitSav;}
920  bool lastWasInitial() {return lastWasInitialSav;}
921  bool lastWasResonanceDecay() {return lastWasDecSav;}
922 
923  // Clear the antennae.
924  void clearAntennae() {antVecFinal.clear(); antVecInitial.clear();
925  antVecRes.clear(); clearLastTrial();}
926 
927  // Clear the last saved trial (but keep the antennae).
928  void clearLastTrial() {
929  q2Trial = 0; antTrial = nullptr; lastWasSplitSav = false;
930  lastWasDecSav = false; lastWasInitialSav = false; lastWasBelowCut = false;
931  ISav = 0; KSav = 0;}
932 
933 private:
934 
935  // Members.
936  // Hadronic invariant mass.
937  double shh;
938 
939  // System and whether to do only resonance decays or full EW shower.
940  int iSysSav;
941  bool resDecOnlySav;
942 
943  // Cutoff scale.
944  double q2Cut;
945 
946  // Pointers.
947  BeamParticle* beamAPtr{};
948  BeamParticle* beamBPtr{};
949  Info* infoPtr{};
950  PartonSystems* partonSystemsPtr{};
951  Rndm* rndmPtr{};
952  Settings* settingsPtr{};
953  Logger* loggerPtr{};
954  VinciaCommon* vinComPtr{};
955  AlphaEM* al{};
956 
957  // Antennae.
958  vector<EWAntennaFF> antVecFinal;
959  vector<EWAntennaII> antVecInitial;
960  vector<EWAntennaFFres> antVecRes;
961 
962  // Trial information.
963  EWAntenna* antTrial{};
964  double q2Trial;
965  bool lastWasSplitSav, lastWasDecSav, lastWasInitialSav, lastWasBelowCut;
966  int ISav, KSav;
967 
968  // Map of branchings.
969  unordered_map<pair<int, int>, vector<EWBranching> >
970  *brMapFinal{}, *brMapInitial{}, *brMapResonance{};
971 
972  // Cluster maps for spectator selection.
973  unordered_map<pair<int, int>, vector<pair<int, int> > >
974  *cluMapFinal{}, *cluMapInitial{};
975 
976  // Amplitude calculator.
977  AmpCalculator* ampCalcPtr{};
978 
979  // Flags.
980  bool isInit, doVetoHardEmissions;
981  int verbose;
982  double vetoHardEmissionsDeltaR2;
983 
984 };
985 
986 //==========================================================================
987 
988 // Top-level class for the electroweak shower module.
989 
990 class VinciaEW : public VinciaModule {
991 
992 public:
993 
994  // Constructor.
995  VinciaEW(): isLoaded{false} {;}
996 
997  // Initialize pointers (called at construction time).
998  void initPtr(Info* infoPtrIn, VinciaCommon* vinComPtrIn) override {
999  infoPtr = infoPtrIn;
1000  particleDataPtr = infoPtr->particleDataPtr;
1001  loggerPtr = infoPtr->loggerPtr;
1002  partonSystemsPtr = infoPtr->partonSystemsPtr;
1003  rndmPtr = infoPtr->rndmPtr;
1004  settingsPtr = infoPtr->settingsPtr;
1005  vinComPtr = vinComPtrIn;
1006  ampCalc.initPtr(infoPtr, &al, &vinComPtr->alphaStrong);
1007  isInitPtr = true;}
1008 
1009  // Initialise settings for current run (called as part of Pythia::init()).
1010  void init(BeamParticle* beamAPtrIn = nullptr,
1011  BeamParticle* beamBPtrIn = nullptr) override;
1012 
1013  // Select helicities for a resonance-decay system.
1014  bool polarise(vector<Particle> &state) override {
1015  if (isLoaded) return ampCalc.polarise(state);
1016  else return false;}
1017 
1018  // Prepare to shower a system.
1019  // (If isBelowHadIn = true, assume only resonance decays may be left to do.)
1020  bool prepare(int iSysIn, Event &event, int scaleRegionIn = 0) override;
1021 
1022  // Update EW shower system each time something has changed.
1023  void update(Event &event, int iSysIn) override {
1024  if (verbose >= VinciaConstants::DEBUG)
1025  printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
1026  if (iSysIn != ewSystem.system()) return;
1027  else ewSystem.buildSystem(event);
1028  if (verbose >= VinciaConstants::DEBUG)
1029  printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);}
1030 
1031  // Set verbose level.
1032  void setVerbose(int verboseIn) override {
1033  verbose = verboseIn; ampCalc.setVerbose(verboseIn);}
1034 
1035  // Save information on masses and widths, and load branchings.
1036  void load() override;
1037 
1038  // Generate a trial scale.
1039  double q2Next(Event&, double q2Start,double q2End) override;
1040 
1041  // Query last branching.
1042  bool lastIsSplitting() override { return ewSystem.lastWasSplitToFermions(); }
1043  bool lastIsInitial() override { return ewSystem.lastWasInitial(); }
1044  bool lastIsResonanceDecay() override {
1045  return ewSystem.lastWasResonanceDecay(); }
1046 
1047  // Check veto.
1048  bool acceptTrial(Event& event) override {
1049  if (verbose >= VinciaConstants::DEBUG)
1050  printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
1051  bool success = false;
1052  if (ewSystem.hasTrial()) success = ewSystem.acceptTrial(event);
1053  else loggerPtr->ERROR_MSG("trial doesn't exist!");
1054  if (verbose >= VinciaConstants::DEBUG)
1055  printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);
1056  return success;}
1057 
1058  // Update event after branching accepted.
1059  void updateEvent(Event& event) override {
1060  if (verbose >= VinciaConstants::DEBUG)
1061  printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
1062  if (ewSystem.hasTrial()) ewSystem.updateEvent(event);
1063  else loggerPtr->ERROR_MSG("trial doesn't exist!");
1064  if (verbose >=VinciaConstants::DEBUG) {
1065  printOut(__METHOD_NAME__,"Event after update:"); event.list();
1066  printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);}}
1067 
1068  // Update partonSystems after branching accepted.
1069  void updatePartonSystems(Event& event) override {
1070  if (verbose >= VinciaConstants::DEBUG)
1071  printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
1072  if (ewSystem.hasTrial()) ewSystem.updatePartonSystems(event);
1073  else loggerPtr->ERROR_MSG("trial doesn't exist!");
1074  if (verbose >= VinciaConstants::DEBUG)
1075  printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);}
1076 
1077  // Clear EW system.
1078  void clear(int) override {
1079  ewSystem.clearAntennae();
1080  ampCalc.eventWeight(1);
1081  }
1082 
1083  // Get number of EW branchers.
1084  unsigned int nBranchers() override {return ewSystem.nBranchers();}
1085  unsigned int nResDec() override {return ewSystem.nResDec();}
1086 
1087  // Print loaded data on masses and widths.
1088  void printData();
1089 
1090  // Print branchings.
1091  void printBranchings();
1092 
1093  // Override additional base class methods.
1094  double q2min() override {return q2minSav;}
1095  double q2minColoured() override {return q2minSav;}
1096  double eventWeight() {return ampCalc.eventWeight();}
1097 
1098  // Setter for q2minSav
1099  void q2min(double q2minSavIn) {q2minSav = q2minSavIn;}
1100 
1101  // We only have one system.
1102  int sysWin() override {return ewSystem.system();}
1103 
1104  // Check if a particle is a resonance according to EW shower database.
1105  bool isResonance(int pid) {return ewData.isRes(pid);}
1106 
1107  // Public members.
1108 
1109  // Cluster maps for spectator selection.
1110  unordered_map<pair<int, int>, vector<pair<int, int> > >
1111  cluMapFinal, cluMapInitial;
1112 
1113  // Map from (PID, polarisation) of I to all possible branchings where I is:
1114  // final, initial or resonance.
1115  unordered_map<pair<int, int>, vector<EWBranching> >
1116  brMapFinal, brMapInitial, brMapResonance;
1117 
1118  // Locally store data about EW particles.
1120 
1121  // Amplitude calculator.
1123 
1124 private:
1125 
1126  // Read in data for overestimate function coefficients for EW branchings
1127  // from XML and store.
1128  bool readFile(string file);
1129  bool readLine(string line);
1130  bool attributeValue(string line, string attribute, string& val);
1131  template <class T> bool attributeValue(
1132  string line, string attribute, T& val) {
1133  string valString;
1134  if (!attributeValue(line, attribute, valString)) return false;
1135  istringstream valStream(valString);
1136  if ( !(valStream >> val) ) {
1137  loggerPtr->ERROR_MSG("failed to store attribute "
1138  + attribute + " " + valString);
1139  return false;
1140  } else return true;
1141  }
1142  bool addBranching(string line, unordered_map< pair<int, int>,
1143  vector<EWBranching> > & branchings, unordered_map< pair<int, int>,
1144  vector<pair<int, int> > > & clusterings, double headroom, bool decay);
1145  bool addParticle(int idIn, int polIn, bool isRes);
1146 
1147  // Members.
1148 
1149  // Cutoff scale.
1150  double q2minSav;
1151 
1152  // AlphaEM.
1153  AlphaEM al;
1154 
1155  // System.
1156  EWSystem ewSystem;
1157 
1158  // Trial information.
1159  double q2Trial;
1160 
1161  // Switches.
1162  bool isLoaded, doEW, doFFbranchings, doIIbranchings, doRFbranchings,
1163  doBosonInterference;
1164 
1165  // Massless quark flavors.
1166  int nFlavZeroMass;
1167 
1168  // Overestimate headroom.
1169  double headroomFinal, headroomInitial;
1170 
1171 };
1172 
1173 //==========================================================================
1174 
1175 // Class to do the veto for overlapping QCD/EW shower phase space.
1176 
1177 class VinciaEWVetoHook : public UserHooks {
1178 
1179  public:
1180 
1181  // Constructor.
1182  VinciaEWVetoHook() = default;
1183 
1184  // Define to activate veto.
1185  bool canVetoISREmission() override {return mayVeto;}
1186  bool canVetoFSREmission() override {return mayVeto;}
1187 
1188  // Methods to veto EW/QCD emissions.
1189  bool doVetoISREmission(int sizeOld, const Event &event, int iSys) override;
1190  bool doVetoFSREmission(int sizeOld, const Event &event, int iSys,
1191  bool inResonance = false) override;
1192 
1193  // Initialize.
1194  void init(shared_ptr<VinciaEW> ewShowerPtrIn);
1195 
1196  private:
1197 
1198  // Universal method called by doVeto(FSR/ISR)Emission.
1199  bool doVetoEmission(int sizeOld, const Event &event, int iSys);
1200  // Find all QCD clusterings and evaluate their kt measure.
1201  double findQCDScale(int sizeOld, const Event& event, int iSys);
1202  // Find all EW clusterings and evaluate their kt measure.
1203  double findEWScale(int sizeOld, const Event& event, int iSys);
1204  // Evaluate Durham kt measure for two particles i and j.
1205  double ktMeasure(const Event& event, int indexi, int indexj, double mI2);
1206  // Evaluate a QCD clustering - returns -1 if not a valid clustering.
1207  double findktQCD(const Event& event, int indexi, int indexj);
1208  // Evaluate an EW clustering - returns -1 if not a valid clustering.
1209  double findktEW(const Event& event, int indexi, int indexj);
1210  // Look up last FSR emission info from event record.
1211  bool setLastFSREmission(int sizeOld, const Event& event);
1212  // Look up last ISR emission info from event record.
1213  bool setLastISREmission(int sizeOld, const Event& event);
1214 
1215  // Members.
1216 
1217  // Verbosity.
1218  int verbose;
1219 
1220  // Flag to control if vetoing occurs.
1221  bool mayVeto{true};
1222 
1223  // delta R parameter used in kt measure.
1224  double deltaR;
1225 
1226  // Electroweak scale
1227  double q2EW;
1228 
1229  // Last emission info.
1230  bool lastIsQCD;
1231  double lastkT2;
1232 
1233  // Pointer to the EW shower
1234  shared_ptr<VinciaEW> ewShowerPtr{};
1235 
1236 };
1237 
1238 
1239 //==========================================================================
1240 
1241 } // end namespace Pythia8
1242 
1243 #endif // Pythia8_VinciaEW_H
void init(BeamParticle *beamAPtrIn, BeamParticle *beamBPtrIn)
Initialize.
Definition: VinciaEW.h:808
AntWrapper norm()
Normalization.
Definition: VinciaEW.h:134
bool canVetoISREmission() override
Define to activate veto.
Definition: VinciaEW.h:1185
constexpr double pow2(const double &x)
Powers of small integers - for balance speed/code clarity.
Definition: PythiaStdlib.h:173
double width(int id, int pol)
Return particle width.
Definition: VinciaEW.h:73
bool canVetoFSREmission() override
Possibility to veto an emission in the FSR machinery.
Definition: VinciaEW.h:1186
int iMot
Indices, PID, and polarization of I, K in Pythia event record.
Definition: VinciaEW.h:621
Class to contain an amplitude and two outgoing polarizations.
Definition: VinciaEW.h:125
unordered_map< pair< int, int >, vector< pair< int, int > > > cluMapFinal
Public members.
Definition: VinciaEW.h:1111
double c0Sum
Info on coefficents.
Definition: VinciaEW.h:643
double val
Members.
Definition: VinciaEW.h:116
std::complex< double > complex
Convenient typedef for double precision complex numbers.
Definition: PythiaComplex.h:17
Rndm * rndmPtr
Pointer to the random number generator.
Definition: Info.h:89
EWBranching(int idMotIn, int idiIn, int idjIn, int polMotIn, double c0In=0, double c1In=0, double c2In=0, double c3In=0)
Constructor.
Definition: VinciaEW.h:534
Simple class to save information about particles.
Definition: VinciaEW.h:30
int system()
Return which system we are handling.
Definition: VinciaEW.h:827
void hmsgIIAnt(int polA, int pola, int polj)
Report helicity combination error for an II antenna function.
Definition: VinciaEW.h:291
Definition: StandardModel.h:23
Class that contains an electroweak branching.
Definition: VinciaEW.h:529
AmpWrapper & operator+=(complex c)
Operators.
Definition: VinciaEW.h:137
Definition: Info.h:45
The Event class holds all info on the generated event.
Definition: Event.h:408
unordered_map< pair< int, int >, EWParticle > data
Member.
Definition: VinciaEW.h:96
Definition: BeamParticle.h:133
EWParticleData ewData
Locally store data about EW particles.
Definition: VinciaEW.h:1119
Initial-initial electroweak antenna.
Definition: VinciaEW.h:734
void print()
Print.
Definition: VinciaEW.h:541
void clearAntennae()
Clear the antennae.
Definition: VinciaEW.h:924
bool polarise(vector< Particle > &state) override
Select helicities for a resonance-decay system.
Definition: VinciaEW.h:1014
int idMot
For a branching I->ij, particle IDs and polarisation of I.
Definition: VinciaEW.h:546
EWParticle()=default
Constructor.
Calculator class for amplitudes, antennae, and Breit-Wigners.
Definition: VinciaEW.h:151
double q2Match
Matching scale.
Definition: VinciaEW.h:647
void initPtr(Info *infoPtrIn, AlphaEM *alphaPtrIn, AlphaStrong *alphaSPtrIn)
Initialize the pointers.
Definition: VinciaEW.h:156
Class to do the veto for overlapping QCD/EW shower phase space.
Definition: VinciaEW.h:1177
Base class for Vincia&#39;s QED and EW shower modules.
Definition: VinciaQED.h:456
Final-final electroweak antenna.
Definition: VinciaEW.h:676
void print()
Print.
Definition: VinciaEW.h:113
double alpha
Overestimate of QED coupling.
Definition: VinciaEW.h:627
vector< Vec4 > pNew
Outgoing momenta after branching.
Definition: VinciaEW.h:640
void printOut(string, string, int nPad=0, char padChar='-')
end METHOD_NAME
Definition: VinciaCommon.cc:4979
Top-level class for the electroweak shower module.
Definition: VinciaEW.h:990
void setVerbose(int verboseIn) override
Set verbose level.
Definition: VinciaEW.h:1032
bool find(int id, int pol)
Find particle data.
Definition: VinciaEW.h:54
AmpWrapper(complex valIn, int poliIn, int poljIn)
Constructor.
Definition: VinciaEW.h:130
void hmsgFSRSplit(int polMot, int poli, int polj)
Report helicty combination error for an FSR splitting kernel.
Definition: VinciaEW.h:371
void hmsgISRSplit(int polA, int pola, int polj)
Report helicty combination error for an ISR splitting kernel.
Definition: VinciaEW.h:388
void hmsgFFAnt(int polMot, int poli, int polj)
Report helicity combination error for an FF antenna function.
Definition: VinciaEW.h:279
void print()
Print, must be implemented by base classes.
Definition: VinciaEW.h:567
Class that performs electroweak showers in a single parton system.
Definition: VinciaEW.h:777
Definition: StandardModel.h:106
Definition: Logger.h:23
bool lastIsSplitting() override
Query last branching.
Definition: VinciaEW.h:1042
int sysWin() override
We only have one system.
Definition: VinciaEW.h:1102
void setVerbose(int verboseIn)
Set verbosity.
Definition: VinciaEW.h:816
double mass(int id)
Return particle mass.
Definition: VinciaEW.h:66
double mass
Members.
Definition: VinciaEW.h:40
unsigned int nBranchers() override
Get number of EW branchers.
Definition: VinciaEW.h:1084
EWAntenna()
Constructor and destructor.
Definition: VinciaEW.h:563
Definition: Basics.h:388
double eventWeight()
EW event weight.
Definition: VinciaEW.h:459
Definition: VinciaCommon.h:494
double q2min() override
Override additional base class methods.
Definition: VinciaEW.h:1094
double branchKernelFF(const Vec4 &pi, const Vec4 &pj, int idMot, int idi, int idj, double mMot, double widthQ2, int polMot, int poli, int polj)
Compute FF antenna function from amplitudes.
Definition: VinciaEW.h:255
void updateEvent(Event &event)
Update an event.
Definition: VinciaEW.h:886
bool hasTrial
Trial variables.
Definition: VinciaEW.h:635
int jNew
Information for partonSystems.
Definition: VinciaEW.h:650
complex val
Members.
Definition: VinciaEW.h:142
int getIndexMot()
Return index.
Definition: VinciaEW.h:605
bool isResonance(int pid)
Check if a particle is a resonance according to EW shower database.
Definition: VinciaEW.h:1105
void initFSRSplit(bool va, int id1, int id2, int pol, const double &mMot, const double &miIn, const double &mjIn)
Initialize an FSR splitting kernel.
Definition: VinciaEW.h:361
unordered_map< pair< int, int >, vector< EWBranching > > brMapFinal
Definition: VinciaEW.h:1116
bool acceptTrial(Event &event) override
Check veto.
Definition: VinciaEW.h:1048
bool isRes(int id)
Return if particle is resonance.
Definition: VinciaEW.h:81
Final-final electroweak resonance antenna.
Definition: VinciaEW.h:702
double q2minColoured() override
End scales.
Definition: VinciaEW.h:1095
AmpCalculator ampCalc
Amplitude calculator.
Definition: VinciaEW.h:1122
void initISRSplit(bool va, int id1, int id2, int pol, const double &mA, const double &maIn, const double &mjIn)
Initialize an ISR splitting kernel.
Definition: VinciaEW.h:378
Locally saved particle data in glorified map.
Definition: VinciaEW.h:49
Class to contain an antenna function and two outgoing polarizations.
Definition: VinciaEW.h:104
void updatePartonSystems(Event &event) override
Update partonSystems after branching accepted.
Definition: VinciaEW.h:1069
void initPtr(Info *infoPtrIn, VinciaCommon *vinComPtrIn, AlphaEM *alIn)
Initialize pointers.
Definition: VinciaEW.h:798
void clear(int) override
Clear EW system.
Definition: VinciaEW.h:1078
int iSys
Parton system this antenna is part of.
Definition: VinciaEW.h:629
Base class for an electroweak antenna.
Definition: VinciaEW.h:558
ParticleData * particleDataPtr
Pointer to the particle data table.
Definition: Info.h:83
PartonSystems * partonSystemsPtr
Pointer to information on subcollision parton locations.
Definition: Info.h:99
Vec4 pMot
Momenta of antenna constituents.
Definition: VinciaEW.h:623
bool lastWasSplitToFermions()
Check if split to fermions, initial, or resonance.
Definition: VinciaEW.h:919
vector< EWBranching > brVec
EW branchings.
Definition: VinciaEW.h:632
int verbose
Verbosity.
Definition: VinciaEW.h:668
void initPtr(Info *infoPtrIn, VinciaCommon *vinComPtrIn) override
Initialize pointers (called at construction time).
Definition: VinciaEW.h:998
EWAntennaII(BeamParticle *beamAPtrIn, BeamParticle *beamBPtrIn)
Constructor.
Definition: VinciaEW.h:739
bool initPtr(Info *infoPtrIn)
Initialize pointers.
Definition: VinciaCommon.h:499
bool prepare(Event &event, int iSysIn, double q2CutIn, bool resDecOnlyIn)
Prepare an event.
Definition: VinciaEW.h:819
const int DEBUG
Verbosity levels. Vincia has one more level (debug) beyond report.
Definition: VinciaCommon.h:49
void addAntenna(T ant, vector< T > &antVec, Event &event, int iMot, int iRec, unordered_map< pair< int, int >, vector< EWBranching > > *brMapPtr)
Add an antenna.
Definition: VinciaEW.h:833
void initPtr(Info *infoPtrIn, VinciaCommon *vinComPtrIn, AlphaEM *alphaPtrIn, AmpCalculator *ampCalcPtrIn)
Initialize pointers.
Definition: VinciaEW.h:575
The PartonSystems class describes the whole set of subcollisions.
Definition: PartonSystems.h:42
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
double m(const Vec4 &v1)
Invariant mass and its square.
Definition: Basics.cc:588
double branchKernelII(Vec4 pa, Vec4 pj, int idA, int ida, int idj, double mA, int polA, int pola, int polj)
Compute II antenna function from amplitudes.
Definition: VinciaEW.h:265
void q2min(double q2minSavIn)
Setter for q2minSav.
Definition: VinciaEW.h:1099
AntWrapper(double valIn, int poliIn, int poljIn)
Constructor.
Definition: VinciaEW.h:109
UserHooks is base class for user access to program execution.
Definition: UserHooks.h:32
bool acceptTrial(Event &event)
Accept a trial.
Definition: VinciaEW.h:879
const int DASHLEN
Padding length for dashes in standardised Vincia verbose output.
Definition: VinciaCommon.h:52
bool isSplitToFermions
Store if branching is a splitting.
Definition: VinciaEW.h:550
void setVerbose(int verboseIn)
Set verbosity level.
Definition: VinciaEW.h:175
void update(Event &event, int iSysIn) override
Update EW shower system each time something has changed.
Definition: VinciaEW.h:1023
bool isSplitToFermions()
Check if splitting to fermions, inital, or resoance.
Definition: VinciaEW.h:609
void updateEvent(Event &event) override
Update event after branching accepted.
Definition: VinciaEW.h:1059
void updatePartonSystems(Event &event)
Update parton systems.
Definition: VinciaEW.h:895
EWSystem()
Constructors.
Definition: VinciaEW.h:782
void add(int id, int pol, double massIn, double widthIn, bool isResIn)
Add particle data.
Definition: VinciaEW.h:58
double c0
Overestimate constants used to sample antennae.
Definition: VinciaEW.h:548
void clearLastTrial()
Clear the last saved trial (but keep the antennae).
Definition: VinciaEW.h:928
void printAntennae()
Print the antennas.
Definition: VinciaEW.h:904
void setVerbose(int verboseIn)
Set verbosity level.
Definition: VinciaEW.h:591
bool flag(string key) const
Shorthand to read settings values.
Definition: PhysicsBase.h:41
bool doBosonInterference
Settings.
Definition: VinciaEW.h:665
void generateTrial(vector< T > &antVec, double q2Start, double q2End, double alphaIn)
Generate a trial.
Definition: VinciaEW.h:855
Definition: Basics.h:32
VinciaEW()
Constructor.
Definition: VinciaEW.h:995
double sAnt
Masses and invariants of antenna constituents.
Definition: VinciaEW.h:625
double mass(int id, int pol)
Return particle mass.
Definition: VinciaEW.h:62
Definition: Settings.h:196
unordered_map< pair< int, int >, double > vMap
Definition: VinciaEW.h:469
bool isRes(int id, int pol)
Return if particle is resonance.
Definition: VinciaEW.h:77
bool hasTrial()
Check if system has a trial pointer.
Definition: VinciaEW.h:916
unsigned int nBranchers()
Get number of branchers (total and res-dec only)
Definition: VinciaEW.h:911