PYTHIA  8.312
VinciaEW.h
1 // VinciaEW.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2024 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{}, *beamBPtr{};
948  Info* infoPtr{};
949  PartonSystems* partonSystemsPtr{};
950  Rndm* rndmPtr{};
951  Settings* settingsPtr{};
952  Logger* loggerPtr{};
953  VinciaCommon* vinComPtr{};
954  AlphaEM* al{};
955 
956  // Antennae.
957  vector<EWAntennaFF> antVecFinal;
958  vector<EWAntennaII> antVecInitial;
959  vector<EWAntennaFFres> antVecRes;
960 
961  // Trial information.
962  EWAntenna* antTrial{};
963  double q2Trial;
964  bool lastWasSplitSav, lastWasDecSav, lastWasInitialSav, lastWasBelowCut;
965  int ISav, KSav;
966 
967  // Map of branchings.
968  unordered_map<pair<int, int>, vector<EWBranching> >
969  *brMapFinal{}, *brMapInitial{}, *brMapResonance{};
970 
971  // Cluster maps for spectator selection.
972  unordered_map<pair<int, int>, vector<pair<int, int> > >
973  *cluMapFinal{}, *cluMapInitial{};
974 
975  // Amplitude calculator.
976  AmpCalculator* ampCalcPtr{};
977 
978  // Flags.
979  bool isInit, doVetoHardEmissions;
980  int verbose;
981  double vetoHardEmissionsDeltaR2;
982 
983 };
984 
985 //==========================================================================
986 
987 // Top-level class for the electroweak shower module.
988 
989 class VinciaEW : public VinciaModule {
990 
991 public:
992 
993  // Constructor.
994  VinciaEW(): isLoaded{false} {;}
995 
996  // Initialize pointers (called at construction time).
997  void initPtr(Info* infoPtrIn, VinciaCommon* vinComPtrIn) override {
998  infoPtr = infoPtrIn;
999  particleDataPtr = infoPtr->particleDataPtr;
1000  loggerPtr = infoPtr->loggerPtr;
1001  partonSystemsPtr = infoPtr->partonSystemsPtr;
1002  rndmPtr = infoPtr->rndmPtr;
1003  settingsPtr = infoPtr->settingsPtr;
1004  vinComPtr = vinComPtrIn;
1005  ampCalc.initPtr(infoPtr, &al, &vinComPtr->alphaStrong);
1006  isInitPtr = true;}
1007 
1008  // Initialise settings for current run (called as part of Pythia::init()).
1009  void init(BeamParticle* beamAPtrIn = 0, BeamParticle* beamBPtrIn = 0)
1010  override;
1011 
1012  // Select helicities for a resonance-decay system.
1013  bool polarise(vector<Particle> &state) override {
1014  if (isLoaded) return ampCalc.polarise(state);
1015  else return false;}
1016 
1017  // Prepare to shower a system.
1018  // (If isBelowHadIn = true, assume only resonance decays may be left to do.)
1019  bool prepare(int iSysIn, Event &event, bool isBelowHadIn=false) override;
1020 
1021  // Update EW shower system each time something has changed.
1022  void update(Event &event, int iSysIn) override {
1023  if (verbose >= VinciaConstants::DEBUG)
1024  printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
1025  if (iSysIn != ewSystem.system()) return;
1026  else ewSystem.buildSystem(event);
1027  if (verbose >= VinciaConstants::DEBUG)
1028  printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);}
1029 
1030  // Set verbose level.
1031  void setVerbose(int verboseIn) override {
1032  verbose = verboseIn; ampCalc.setVerbose(verboseIn);}
1033 
1034  // Save information on masses and widths, and load branchings.
1035  void load() override;
1036 
1037  // Generate a trial scale.
1038  double q2Next(Event&, double q2Start,double q2End) override;
1039 
1040  // Query last branching.
1041  bool lastIsSplitting() override { return ewSystem.lastWasSplitToFermions(); }
1042  bool lastIsInitial() override { return ewSystem.lastWasInitial(); }
1043  bool lastIsResonanceDecay() override {
1044  return ewSystem.lastWasResonanceDecay(); }
1045 
1046  // Check veto.
1047  bool acceptTrial(Event& event) override {
1048  if (verbose >= VinciaConstants::DEBUG)
1049  printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
1050  bool success = false;
1051  if (ewSystem.hasTrial()) success = ewSystem.acceptTrial(event);
1052  else loggerPtr->ERROR_MSG("trial doesn't exist!");
1053  if (verbose >= VinciaConstants::DEBUG)
1054  printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);
1055  return success;}
1056 
1057  // Update event after branching accepted.
1058  void updateEvent(Event& event) override {
1059  if (verbose >= VinciaConstants::DEBUG)
1060  printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
1061  if (ewSystem.hasTrial()) ewSystem.updateEvent(event);
1062  else loggerPtr->ERROR_MSG("trial doesn't exist!");
1063  if (verbose >=VinciaConstants::DEBUG) {
1064  printOut(__METHOD_NAME__,"Event after update:"); event.list();
1065  printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);}}
1066 
1067  // Update partonSystems after branching accepted.
1068  void updatePartonSystems(Event& event) override {
1069  if (verbose >= VinciaConstants::DEBUG)
1070  printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
1071  if (ewSystem.hasTrial()) ewSystem.updatePartonSystems(event);
1072  else loggerPtr->ERROR_MSG("trial doesn't exist!");
1073  if (verbose >= VinciaConstants::DEBUG)
1074  printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);}
1075 
1076  // Clear EW system.
1077  void clear(int) override {
1078  ewSystem.clearAntennae();
1079  ampCalc.eventWeight(1);
1080  }
1081 
1082  // Get number of EW branchers.
1083  unsigned int nBranchers() override {return ewSystem.nBranchers();}
1084  unsigned int nResDec() override {return ewSystem.nResDec();}
1085 
1086  // Print loaded data on masses and widths.
1087  void printData();
1088 
1089  // Print branchings.
1090  void printBranchings();
1091 
1092  // Override additional base class methods.
1093  double q2min() override {return q2minSav;}
1094  double q2minColoured() override {return q2minSav;}
1095  double eventWeight() {return ampCalc.eventWeight();}
1096 
1097  // Setter for q2minSav
1098  void q2min(double q2minSavIn) {q2minSav = q2minSavIn;}
1099 
1100  // We only have one system.
1101  int sysWin() override {return ewSystem.system();}
1102 
1103  // Check if a particle is a resonance according to EW shower database.
1104  bool isResonance(int pid) {return ewData.isRes(pid);}
1105 
1106  // Public members.
1107 
1108  // Cluster maps for spectator selection.
1109  unordered_map<pair<int, int>, vector<pair<int, int> > >
1110  cluMapFinal, cluMapInitial;
1111 
1112  // Map from (PID, polarisation) of I to all possible branchings where I is:
1113  // final, initial or resonance.
1114  unordered_map<pair<int, int>, vector<EWBranching> >
1115  brMapFinal, brMapInitial, brMapResonance;
1116 
1117  // Locally store data about EW particles.
1119 
1120  // Amplitude calculator.
1122 
1123 private:
1124 
1125  // Read in data for overestimate function coefficients for EW branchings
1126  // from XML and store.
1127  bool readFile(string file);
1128  bool readLine(string line);
1129  bool attributeValue(string line, string attribute, string& val);
1130  template <class T> bool attributeValue(
1131  string line, string attribute, T& val) {
1132  string valString;
1133  if (!attributeValue(line, attribute, valString)) return false;
1134  istringstream valStream(valString);
1135  if ( !(valStream >> val) ) {
1136  loggerPtr->ERROR_MSG("failed to store attribute "
1137  + attribute + " " + valString);
1138  return false;
1139  } else return true;
1140  }
1141  bool addBranching(string line, unordered_map< pair<int, int>,
1142  vector<EWBranching> > & branchings, unordered_map< pair<int, int>,
1143  vector<pair<int, int> > > & clusterings, double headroom, bool decay);
1144  bool addParticle(int idIn, int polIn, bool isRes);
1145 
1146  // Members.
1147 
1148  // Cutoff scale.
1149  double q2minSav;
1150 
1151  // AlphaEM.
1152  AlphaEM al;
1153 
1154  // System.
1155  EWSystem ewSystem;
1156 
1157  // Trial information.
1158  double q2Trial;
1159 
1160  // Switches.
1161  bool isLoaded, doEW, doFFbranchings, doIIbranchings, doRFbranchings,
1162  doBosonInterference;
1163 
1164  // Massless quark flavors.
1165  int nFlavZeroMass;
1166 
1167  // Overestimate headroom.
1168  double headroomFinal, headroomInitial;
1169 
1170 };
1171 
1172 //==========================================================================
1173 
1174 // Class to do the veto for overlapping QCD/EW shower phase space.
1175 
1176 class VinciaEWVetoHook : public UserHooks {
1177 
1178  public:
1179 
1180  // Constructor.
1181  VinciaEWVetoHook() = default;
1182 
1183  // Define to activate veto.
1184  bool canVetoISREmission() override {return mayVeto;}
1185  bool canVetoFSREmission() override {return mayVeto;}
1186 
1187  // Methods to veto EW/QCD emissions.
1188  bool doVetoISREmission(int sizeOld, const Event &event, int iSys) override;
1189  bool doVetoFSREmission(int sizeOld, const Event &event, int iSys,
1190  bool inResonance = false) override;
1191 
1192  // Initialize.
1193  void init(shared_ptr<VinciaEW> ewShowerPtrIn);
1194 
1195  private:
1196 
1197  // Universal method called by doVeto(FSR/ISR)Emission.
1198  bool doVetoEmission(int sizeOld, const Event &event, int iSys);
1199  // Find all QCD clusterings and evaluate their kt measure.
1200  double findQCDScale(int sizeOld, const Event& event, int iSys);
1201  // Find all EW clusterings and evaluate their kt measure.
1202  double findEWScale(int sizeOld, const Event& event, int iSys);
1203  // Evaluate Durham kt measure for two particles i and j.
1204  double ktMeasure(const Event& event, int indexi, int indexj, double mI2);
1205  // Evaluate a QCD clustering - returns -1 if not a valid clustering.
1206  double findktQCD(const Event& event, int indexi, int indexj);
1207  // Evaluate an EW clustering - returns -1 if not a valid clustering.
1208  double findktEW(const Event& event, int indexi, int indexj);
1209  // Look up last FSR emission info from event record.
1210  bool setLastFSREmission(int sizeOld, const Event& event);
1211  // Look up last ISR emission info from event record.
1212  bool setLastISREmission(int sizeOld, const Event& event);
1213 
1214  // Members.
1215 
1216  // Verbosity.
1217  int verbose;
1218 
1219  // Flag to control if vetoing occurs.
1220  bool mayVeto{true};
1221 
1222  // delta R parameter used in kt measure.
1223  double deltaR;
1224 
1225  // Electroweak scale
1226  double q2EW;
1227 
1228  // Last emission info.
1229  bool lastIsQCD;
1230  double lastkT2;
1231 
1232  // Pointer to the EW shower
1233  shared_ptr<VinciaEW> ewShowerPtr{};
1234 
1235 };
1236 
1237 
1238 //==========================================================================
1239 
1240 } // end namespace Pythia8
1241 
1242 #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:1184
constexpr double pow2(const double &x)
Powers of small integers - for balance speed/code clarity.
Definition: PythiaStdlib.h:182
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:1185
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:1110
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:453
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:1118
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:1013
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:1176
Base class for Vincia&#39;s QED and EW shower modules.
Definition: VinciaQED.h:444
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:4964
Top-level class for the electroweak shower module.
Definition: VinciaEW.h:989
void setVerbose(int verboseIn) override
Set verbose level.
Definition: VinciaEW.h:1031
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:1041
int sysWin() override
We only have one system.
Definition: VinciaEW.h:1101
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:1083
EWAntenna()
Constructor and destructor.
Definition: VinciaEW.h:563
Definition: Basics.h:386
double eventWeight()
EW event weight.
Definition: VinciaEW.h:459
Definition: VinciaCommon.h:494
double q2min() override
Override additional base class methods.
Definition: VinciaEW.h:1093
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:1104
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:1115
bool acceptTrial(Event &event) override
Check veto.
Definition: VinciaEW.h:1047
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:1094
AmpCalculator ampCalc
Amplitude calculator.
Definition: VinciaEW.h:1121
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:1068
void initPtr(Info *infoPtrIn, VinciaCommon *vinComPtrIn, AlphaEM *alIn)
Initialize pointers.
Definition: VinciaEW.h:798
void clear(int) override
Clear EW system.
Definition: VinciaEW.h:1077
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:997
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:595
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:1098
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:1022
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:1058
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:45
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:994
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:195
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