Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-09-15 09:06:52

0001 // VinciaEW.h is a part of the PYTHIA event generator.
0002 // Copyright (C) 2025 Peter Skands, Torbjorn Sjostrand.
0003 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
0004 // Please respect the MCnet Guidelines, see GUIDELINES for details.
0005 
0006 // This file contains the EW antenna-shower class and auxiliary
0007 // classes. Main author is Rob Verheyen.
0008 
0009 #ifndef Pythia8_VinciaEW_H
0010 #define Pythia8_VinciaEW_H
0011 
0012 // PYTHIA 8 headers.
0013 #include "Pythia8/Event.h"
0014 #include "Pythia8/StandardModel.h"
0015 #include "Pythia8/PartonSystems.h"
0016 #include "Pythia8/PythiaComplex.h"
0017 #include "Pythia8/BeamParticle.h"
0018 #include "Pythia8/UserHooks.h"
0019 
0020 // VINCIA functionality.
0021 #include "Pythia8/VinciaCommon.h"
0022 #include "Pythia8/VinciaQED.h"
0023 
0024 namespace Pythia8 {
0025 
0026 //==========================================================================
0027 
0028 // Simple class to save information about particles.
0029 
0030 class EWParticle {
0031 
0032  public:
0033 
0034   // Constructor.
0035   EWParticle() = default;
0036   EWParticle(double massIn, double widthIn, bool isResIn) :
0037     mass(massIn), width(widthIn), isRes(isResIn) {};
0038 
0039   // Members.
0040   double mass{0.}, width{0.};
0041   bool isRes{false};
0042 
0043 };
0044 
0045 //==========================================================================
0046 
0047 // Locally saved particle data in glorified map.
0048 
0049 class EWParticleData {
0050 
0051  public:
0052 
0053   // Find particle data.
0054   bool find(int id, int pol){
0055     return data.find(make_pair(id, pol)) != data.end();}
0056 
0057   // Add particle data.
0058   void add(int id, int pol, double massIn, double widthIn, bool isResIn) {
0059     data[make_pair(id, pol)] = EWParticle(massIn, widthIn, isResIn);}
0060 
0061   // Return particle mass.
0062   double mass(int id, int pol) {
0063     return find(id, pol) ? data[make_pair(id, pol)].mass : 0.;}
0064 
0065   // Return particle mass.
0066   double mass(int id) {
0067     // Every particle has either pol = 1 or pol = 0.
0068     if (find(id, 1)) return data[make_pair(id, 1)].mass;
0069     return find(id, 0) ? data[make_pair(id, 0)].mass : 0.;
0070   }
0071 
0072   // Return particle width.
0073   double width(int id, int pol) {
0074     return find(id, pol) ? data[make_pair(id, pol)].width : 0.;}
0075 
0076   // Return if particle is resonance.
0077   bool isRes(int id, int pol) {
0078     return find(id, pol) && data[make_pair(id, pol)].isRes;}
0079 
0080   // Return if particle is resonance.
0081   bool isRes(int id) {
0082     // Every particle has either pol = 1 or pol = 0.
0083     if (find(id, 1)) return data[make_pair(id, 1)].isRes;
0084     else if (find(id, 0)) return data[make_pair(id, 0)].isRes;
0085     else return false;
0086   }
0087 
0088   // Data access.
0089   EWParticle* get(int id, int pol) { return &data.at(make_pair(id, pol)); }
0090   unordered_map<pair<int, int>, EWParticle>::iterator begin() {
0091     return data.begin();}
0092   unordered_map<pair<int, int>, EWParticle>::iterator end() {
0093     return data.end();}
0094 
0095   // Member.
0096   unordered_map<pair<int,int>, EWParticle> data;
0097 
0098 };
0099 
0100 //==========================================================================
0101 
0102 // Class to contain an antenna function and two outgoing polarizations.
0103 
0104 class AntWrapper {
0105 
0106 public:
0107 
0108   // Constructor.
0109   AntWrapper(double valIn, int poliIn, int poljIn):
0110     val(valIn), poli(poliIn), polj(poljIn) {}
0111 
0112   // Print.
0113   void print() {cout << "(" << poli << ", " << polj << ") " << val;}
0114 
0115   // Members.
0116   double val;
0117   int poli, polj;
0118 
0119 };
0120 
0121 //==========================================================================
0122 
0123 // Class to contain an amplitude and two outgoing polarizations.
0124 
0125 class AmpWrapper {
0126 
0127  public:
0128 
0129   // Constructor.
0130   AmpWrapper(complex valIn, int poliIn, int poljIn):
0131     val(valIn), poli(poliIn), polj(poljIn) {}
0132 
0133   // Normalization.
0134   AntWrapper norm() {return AntWrapper(std::norm(val), poli, polj);}
0135 
0136   // Operators.
0137   AmpWrapper& operator+=(complex c) {this->val += c; return *this;}
0138   AmpWrapper& operator*=(complex c) {this->val *= c; return *this;}
0139   void print() {cout << "(" << poli << ", " << polj << ") " << val;}
0140 
0141   // Members.
0142   complex val;
0143   int poli, polj;
0144 
0145 };
0146 
0147 //==========================================================================
0148 
0149 // Calculator class for amplitudes, antennae, and Breit-Wigners.
0150 
0151 class AmpCalculator {
0152 
0153 public:
0154 
0155   // Initialize the pointers.
0156   void initPtr(Info* infoPtrIn, AlphaEM* alphaPtrIn,
0157     AlphaStrong* alphaSPtrIn) {
0158     infoPtr          = infoPtrIn;
0159     partonSystemsPtr = infoPtr->partonSystemsPtr;
0160     rndmPtr          = infoPtr->rndmPtr;
0161     settingsPtr      = infoPtr->settingsPtr;
0162     loggerPtr        = infoPtr->loggerPtr;
0163     alphaPtr         = alphaPtrIn;
0164     alphaSPtr        = alphaSPtrIn;
0165     isInitPtr        = true;
0166   }
0167 
0168   // Initialize with maps.
0169   void init(EWParticleData* dataIn, unordered_map< pair<int, int>,
0170             vector<pair<int, int> > >* cluMapFinalIn,
0171             unordered_map< pair<int, int>, vector<pair<int, int> > >*
0172             cluMapInitialIn);
0173 
0174   // Set verbosity level.
0175   void setVerbose(int verboseIn) {verbose = verboseIn;}
0176 
0177   // Spinor products used for the calculation of the amplitudes.
0178   complex spinProd(int pol, const Vec4& ka, const Vec4& kb);
0179   complex spinProd(int pol, const Vec4& ka, const Vec4& pa, const Vec4& kb);
0180   complex spinProd(int pol, const Vec4& ka, const Vec4& pa, const Vec4& pb,
0181                    const Vec4& kb);
0182   complex spinProd(int pol, const Vec4& ka, const Vec4& pa, const Vec4& pb,
0183                    const Vec4& pc, const Vec4& kb);
0184   complex spinProd(int pol, const Vec4& ka, const Vec4& pa, const Vec4& pb,
0185                    const Vec4& pc, const Vec4& pd, const Vec4& kb);
0186   Vec4    spinProdFlat(string method, const Vec4& ka, const Vec4& pa);
0187 
0188   // Initialize couplings.
0189   void initCoup(bool va, int id1, int id2, int pol, bool m);
0190 
0191   // Initialize an FSR branching amplitude.
0192   void initFSRAmp(bool va, int id1, int id2, int pol,
0193     const Vec4& pi, const Vec4 &pj, const double& mMot, const double& widthQ2);
0194 
0195   // Check for zero denominator in an FSR amplitude.
0196   bool zdenFSRAmp(const string& method, const Vec4& pi, const Vec4& pj,
0197     bool check);
0198 
0199   // Initialize an ISR branching amplitude.
0200   void initISRAmp(bool va, int id1, int id2, int pol,
0201     const Vec4& pa, const Vec4 &pj, double& mA);
0202 
0203   // Check for zero denominator in an ISR amplitude.
0204   bool zdenISRAmp(const string& method, const Vec4& pa, const Vec4& pj,
0205     bool check);
0206 
0207   // Branching amplitude formalism.
0208   // Naming scheme: f = fermion, v = vector boson, h = higgs.
0209 
0210   // Final-state branching amplitudes.
0211   complex ftofvFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0212     int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0213   complex ftofhFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0214     int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0215   complex fbartofbarvFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0216     int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0217   complex fbartofbarhFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0218     int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0219   complex vTtoffbarFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0220     int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0221   complex vTtovhFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0222     int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0223   complex vTtovvFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0224     int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0225   complex vLtoffbarFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0226     int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0227   complex vLtovhFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0228     int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0229   complex vLtovvFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0230     int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0231   complex htoffbarFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0232     int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0233   complex htovvFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0234     int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0235   complex htohhFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0236     int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0237 
0238   // Initial-state branching amplitudes.
0239   complex ftofvISRAmp(const Vec4& pa, const Vec4& pj, int idA, int ida,
0240     int idj, double mA, int polA, int pola, int polj);
0241   complex ftofhISRAmp(const Vec4& pa, const Vec4& pj, int idA, int ida,
0242     int idj, double mA, int polA, int pola, int polj);
0243   complex fbartofbarvISRAmp(const Vec4& pa, const Vec4& pj, int idA, int ida,
0244     int idj, double mA, int polA, int pola, int polj);
0245   complex fbartofbarhISRAmp(const Vec4& pa, const Vec4& pj, int idA, int ida,
0246     int idj, double mA, int polA, int pola, int polj);
0247 
0248   // Branching amplitude selector.
0249   complex branchAmpFSR(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0250     int idj, double mMot, double widthQ2, int polMot, int poli=9, int polj=9);
0251   complex branchAmpISR(const Vec4& pa, const Vec4& pj, int idA, int ida,
0252     int idj, double mA, int polA, int pola=9, int polj=9);
0253 
0254   // Compute FF antenna function from amplitudes.
0255   double branchKernelFF(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0256     int idj, double mMot, double widthQ2, int polMot, int poli, int polj) {
0257     return norm(branchAmpFSR(pi, pj, idMot, idi, idj, mMot, widthQ2,
0258                              polMot, poli, polj));}
0259 
0260   // Compute FF antenna functions from amplitudes for all polarizations.
0261   vector<AntWrapper> branchKernelFF(Vec4 pi, Vec4 pj, int idMot, int idi,
0262     int idj, double mMot, double widthQ2, int polMot);
0263 
0264   // Compute II antenna function from amplitudes.
0265   double branchKernelII(Vec4 pa, Vec4 pj, int idA, int ida,
0266     int idj, double mA, int polA, int pola, int polj) {
0267     return norm(branchAmpISR(pa, pj, idA, ida, idj, mA, polA, pola, polj));}
0268 
0269   // Compute II antenna functions from amplitudes for all polarizations.
0270   vector<AntWrapper> branchKernelII(Vec4 pa, Vec4 pj, int idA, int ida,
0271     int idj, double mA, int polA);
0272 
0273   // Initialize an FF antenna function.
0274   void initFFAnt(bool va, int id1, int id2, int pol, const double& Q2,
0275     const double& widthQ2, const double& xi, const double& xj,
0276     const double& mMot, const double& miIn, const double& mjIn);
0277 
0278   // Report helicity combination error for an FF antenna function.
0279   void hmsgFFAnt(int polMot, int poli, int polj) {
0280     stringstream ss;
0281     ss << "helicity combination was not found:\n    "
0282        << "polMot = " << polMot << " poli = " << poli << " polj = " << polj;
0283     loggerPtr->ERROR_MSG(ss.str());}
0284 
0285   // Initialize an II antenna function.
0286   void initIIAnt(int id1, int id2, int pol, const double& Q2,
0287     const double& xA, const double& xj,
0288     const double& mA, const double& maIn, const double& mjIn);
0289 
0290   // Report helicity combination error for an II antenna function.
0291   void hmsgIIAnt(int polA, int pola, int polj) {
0292     stringstream ss;
0293     ss << "helicity combination was not found:\n    "
0294        << "polA = " << polA << " pola = " << pola << " polj = " << polj;
0295     loggerPtr->ERROR_MSG(ss.str());}
0296 
0297   // FF Antenna functions for branching process I (K) -> i j (k).
0298   // Q2 is the offshellness of I.
0299   // xi and xj are energy fractions of pi and pj in the collinear limit.
0300   double ftofvFFAnt(double Q2, double widthQ2, double xi, double xj,
0301     int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0302     int polMot, int poli, int polj);
0303   double ftofhFFAnt(double Q2, double widthQ2, double xi, double xj,
0304     int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0305     int polMot, int poli, int polj);
0306   double fbartofbarvFFAnt(double Q2, double widthQ2, double xi, double xj,
0307     int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0308     int polMot, int poli, int polj);
0309   double fbartofbarhFFAnt(double Q2, double widthQ2, double xi, double xj,
0310     int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0311     int polMot, int poli, int polj);
0312   double vtoffbarFFAnt(double Q2, double widthQ2, double xi, double xj,
0313     int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0314     int polMot, int poli, int polj);
0315   double vtovhFFAnt(double Q2, double widthQ2, double xi, double xj,
0316     int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0317     int polMot, int poli, int polj);
0318   double vtovvFFAnt(double Q2, double widthQ2, double xi, double xj,
0319     int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0320     int polMot, int poli, int polj);
0321   double htoffbarFFAnt(double Q2, double widthQ2, double xi, double xj,
0322     int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0323     int polMot, int poli, int polj);
0324   double htovvFFAnt(double Q2, double widthQ2, double xi, double xj,
0325     int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0326     int polMot, int poli, int polj);
0327   double htohhFFAnt(double Q2, double widthQ2, double xi, double xj,
0328     int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0329     int polMot, int poli, int polj);
0330 
0331   // II Antenna functions for branching process A (B) -> a j (b).
0332   // Q2 is the offshellness of A.
0333   // xA and xj are energy fractions of pA and pj in the collinear limit.
0334   double ftofvIIAnt(double Q2, double xA, double xj,
0335     int idA, int ida, int idj, double mA, double maIn, double mjIn,
0336     int polA, int pola, int polj);
0337   double fbartofbarvIIAnt(double Q2, double xA, double xj,
0338     int idA, int ida, int idj, double mA, double maIn, double mjIn,
0339     int polA, int pola, int polj);
0340 
0341   // FF antenna function calculator.
0342   double antFuncFF(double Q2, double widthQ2, double xi, double xj,
0343     int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0344     int polMot, int poli, int polj);
0345 
0346   // FF antenna function calculator for all outgoing polarizations.
0347   vector<AntWrapper> antFuncFF(double Q2, double widthQ2, double xi,
0348     double xj, int idMot, int idi, int idj, double mMot, double miIn,
0349     double mjIn, int polMot);
0350 
0351   // II antenna function calculator.
0352   double antFuncII(double Q2, double xA, double xj,
0353     int idA, int ida, int idj, double mA, double maIn, double mjIn,
0354     int polA, int pola, int polj);
0355 
0356   // II antenna function calculator for all outgoing polarizations.
0357   vector<AntWrapper> antFuncII(double Q2, double xA, double xj,
0358     int idA, int ida, int idj, double mA, double maIn, double mjIn, int polA);
0359 
0360   // Initialize an FSR splitting kernel.
0361   void initFSRSplit(bool va, int id1, int id2, int pol,
0362     const double& mMot, const double& miIn, const double& mjIn) {
0363     mi = miIn; mj = mjIn; mMot2 = pow2(mMot); mi2 = pow2(mi); mj2 = pow2(mj);
0364     initCoup(va, id1, id2, pol, true);}
0365 
0366   // Check for zero denominator in an FSR splitting kernel.
0367   bool zdenFSRSplit(const string& method, const double& Q2, const double& z,
0368     bool check);
0369 
0370   // Report helicty combination error for an FSR splitting kernel.
0371   void hmsgFSRSplit(int polMot, int poli, int polj) {
0372     stringstream ss;
0373     ss << "helicity combination was not found:\n    "
0374        << "polMot = " << polMot << " poli = " << poli << " polj = " << polj;
0375     loggerPtr->ERROR_MSG(ss.str());}
0376 
0377   // Initialize an ISR splitting kernel.
0378   void initISRSplit(bool va, int id1, int id2, int pol,
0379     const double& mA, const double& maIn, const double& mjIn) {
0380     ma = maIn; mj = mjIn; mA2 = pow2(mA); ma2 = pow2(ma); mj2 = pow2(mj);
0381     initCoup(va, id1, id2, pol, mA > VinciaConstants::NANO);}
0382 
0383   // Check for zero denominator in an ISR splitting kernel.
0384   bool zdenISRSplit(const string& method, const double& Q2, const double& z,
0385     bool flip, bool check);
0386 
0387   // Report helicty combination error for an ISR splitting kernel.
0388   void hmsgISRSplit(int polA, int pola, int polj) {
0389     stringstream ss;
0390     ss << "helicity combination was not found:\n    "
0391        << "polA = " << polA << " pola = " << pola << " polj = " << polj;
0392     loggerPtr->ERROR_MSG(ss.str());}
0393 
0394   // Final-state splitting kernels.
0395   double ftofvFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0396     double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0397   double ftofhFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0398     double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0399   double fbartofbarvFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0400     double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0401   double fbartofbarhFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0402     double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0403   double vTtoffbarFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0404     double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0405   double vTtovhFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0406     double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0407   double vTtovvFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0408     double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0409   double vLtoffbarFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0410     double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0411   double vLtovhFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0412     double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0413   double vLtovvFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0414     double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0415   double htoffbarFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0416     double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0417   double htovvFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0418     double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0419   double htohhFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0420     double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0421 
0422   // Initial-state splitting kernels.
0423   double ftofvISRSplit(double Q2, double z, int idA, int ida, int idj,
0424     double mA, double maIn, double mjIn, int polA, int pola, int polj);
0425   double ftofhISRSplit(double Q2, double z, int idA, int ida, int idj,
0426     double mA, double maIn, double mjIn, int polA, int pola, int polj);
0427   double fbartofbarvISRSplit(double Q2, double z, int idA, int ida, int idj,
0428     double mA, double maIn, double mjIn, int polA, int pola, int polj);
0429   double fbartofbarhISRSplit(double Q2, double z, int idA, int ida, int idj,
0430     double mA, double maIn, double mjIn, int polA, int pola, int polj);
0431 
0432   // Splitting kernel caller.
0433   double splitFuncFSR(double Q2, double z, int idMot, int idi, int idj,
0434     double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0435   double splitFuncISR(double Q2, double z, int idA, int ida, int idj,
0436     double mA, double maIn, double mjIn, int polA, int pola, int polj);
0437 
0438   // Compute partial decay width.
0439   double getPartialWidth(int idMot, int idi, int idj, double mMot, int polMot);
0440 
0441   // Compute total decay width.
0442   double getTotalWidth(int idMot, double mMot, int polMot);
0443 
0444   // Breit-Wigner calculators.
0445   double getBreitWigner(int id, double m, int pol);
0446   double getBreitWignerOverestimate(int id, double m, int pol);
0447 
0448   // Generate Breit-Wigner mass.
0449   double sampleMass(int id, int pol);
0450 
0451   // Bosonic interference factor.
0452   void applyBosonInterferenceFactor(Event &event, int XYEv, Vec4 pi, Vec4 pj,
0453     int idi, int idj, int poli, int polj);
0454 
0455   // Polarise a resonance decay.
0456   bool polarise(vector<Particle> &state);
0457 
0458   // EW event weight.
0459   double eventWeight() {return eventWeightSave;}
0460   void   eventWeight(double eventWeightIn) {eventWeightSave = eventWeightIn;}
0461 
0462   // Public data members.
0463 
0464   // EW data.
0465   EWParticleData* dataPtr{};
0466 
0467   // Maps of coupling constants.
0468   // TODO: read this from data file (rather than hard code).
0469   unordered_map<pair<int, int>, double> vMap, aMap, gMap, vCKM;
0470 
0471 private:
0472 
0473   // Data members.
0474 
0475   // Constants used for Breit-Wigner sampling.
0476   unordered_map<int, vector<double> > cBW;
0477   // BW normalizations.
0478   unordered_map<pair<int,int>, double> nBW, np;
0479 
0480   // Event weight.
0481   double eventWeightSave{1};
0482 
0483   // Electroweak constants.
0484   double mw, mw2, sw, sw2;
0485 
0486   // Mode of matching to the Breit-Wigner.
0487   int bwMatchMode;
0488 
0489   // Maps of EW clusterings.
0490   unordered_map<pair<int, int>, vector<pair<int, int> > >* cluMapFinal{};
0491   unordered_map<pair<int, int>, vector<pair<int, int> > >* cluMapInitial{};
0492 
0493   // Vectors with spin assignments.
0494   vector<int> fermionPols, vectorPols, scalarPols;
0495 
0496   // Couplings.
0497   double v, a, vPls, vMin, g;
0498 
0499   // Masses and scale (FSR and ISR).
0500   double mMot2, mi, mi2, mj, mj2, mA2, ma, ma2, isrQ2;
0501   complex M, fsrQ2;
0502 
0503   // Reference vectors (FSR and ISR).
0504   Vec4 kij, ki, kj, pij, kaj, ka, paj;
0505   double wij, wi, wj, wij2, wi2, wj2, waj, wa, waj2, wa2;
0506 
0507   // Antenna function members.
0508   double Q4, Q4gam, Q2til, ant;
0509 
0510   // Pointers.
0511   Info* infoPtr{};
0512   PartonSystems* partonSystemsPtr{};
0513   Rndm* rndmPtr{};
0514   Settings* settingsPtr{};
0515   Logger* loggerPtr{};
0516   AlphaEM* alphaPtr{};
0517   AlphaStrong* alphaSPtr{};
0518 
0519   // Initializations.
0520   bool isInit{false};
0521   bool isInitPtr{false};
0522   int verbose{0};
0523 };
0524 
0525 //==========================================================================
0526 
0527 // Class that contains an electroweak branching.
0528 
0529 class EWBranching {
0530 
0531  public:
0532 
0533   // Constructor.
0534   EWBranching(int idMotIn, int idiIn, int idjIn, int polMotIn,
0535     double c0In = 0, double c1In = 0, double c2In = 0, double c3In = 0):
0536     idMot(idMotIn), idi(idiIn), idj(idjIn), polMot(polMotIn), c0(c0In),
0537     c1(c1In), c2(c2In), c3(c3In),
0538     isSplitToFermions(abs(idMot) > 20 && abs(idi) < 20 && abs(idj) < 20) {;}
0539 
0540   // Print.
0541   void print() {cout <<"    (" << idMot << ", " << polMot << ") -> " << idi <<
0542     "," << idj << ": (" << c0 << ", " << c1 << ", " << c2 << ", " << c3 <<
0543     ") \n";}
0544 
0545   // For a branching I->ij, particle IDs and polarisation of I.
0546   int idMot, idi, idj, polMot;
0547   // Overestimate constants used to sample antennae.
0548   double c0, c1, c2, c3;
0549   // Store if branching is a splitting.
0550   bool isSplitToFermions;
0551 
0552 };
0553 
0554 //==========================================================================
0555 
0556 // Base class for an electroweak antenna.
0557 
0558 class EWAntenna {
0559 
0560  public:
0561 
0562   // Constructor and destructor.
0563   EWAntenna(): iSys(-1), shat(0), doBosonInterference(false) {};
0564   virtual ~EWAntenna() = default;
0565 
0566   // Print, must be implemented by base classes.
0567   void print() {
0568     stringstream ss;
0569     ss << "Brancher = (" << iMot << ", " << polMot
0570        << "), Recoiler = " << iRec;
0571     printOut(__METHOD_NAME__, ss.str());
0572     for (int i = 0; i < (int)brVec.size(); i++) brVec[i].print();}
0573 
0574   // Initialize pointers.
0575   void initPtr(Info* infoPtrIn, VinciaCommon* vinComPtrIn, AlphaEM* alphaPtrIn,
0576     AmpCalculator* ampCalcPtrIn) {
0577     infoPtr          = infoPtrIn;
0578     rndmPtr          = infoPtr->rndmPtr;
0579     loggerPtr        = infoPtr->loggerPtr;
0580     partonSystemsPtr = infoPtr->partonSystemsPtr;
0581     vinComPtr        = vinComPtrIn;
0582     alphaPtr         = alphaPtrIn;
0583     ampCalcPtr       = ampCalcPtrIn;
0584   };
0585 
0586   // Initialize, must be implemented by base classes.
0587   virtual bool init(Event &event, int iMotIn, int iRecIn, int iSysIn,
0588     vector<EWBranching> &branchings, Settings* settingsPtr) = 0;
0589 
0590   // Set verbosity level.
0591   void setVerbose(int verboseIn) {verbose = verboseIn;}
0592 
0593   // Generate a trial scale to be compared against other QCD and EW trials.
0594   virtual double generateTrial(double q2Start, double q2End,
0595     double alphaIn) = 0;
0596 
0597   // Accept/reject step to accept a trial.
0598   virtual bool acceptTrial(Event& event) = 0;
0599 
0600   // Update an event and parton system.
0601   virtual void updateEvent(Event& event) = 0;
0602   virtual void updatePartonSystems(Event& event);
0603 
0604   // Return index.
0605   int getIndexMot() {return iMot;};
0606   int getIndexRec() {return iRec;};
0607 
0608   // Check if splitting to fermions, inital, or resoance.
0609   bool isSplitToFermions() {
0610     return brTrial != nullptr && brTrial->isSplitToFermions;}
0611   virtual bool isInitial() {return false;}
0612   virtual bool isResonanceDecay() {return false;}
0613 
0614   // Select a channel.
0615   bool selectChannel(int idx, const double& cSum, const map<double, int>&
0616     cSumSoFar, int& idi, int& idj, double& mi2, double& mj2);
0617 
0618 protected:
0619 
0620   // Indices, PID, and polarization of I, K in Pythia event record.
0621   int iMot, iRec, idMot, idRec, polMot;
0622   // Momenta of antenna constituents.
0623   Vec4 pMot, pRec;
0624   // Masses and invariants of antenna constituents.
0625   double sAnt, mMot, mMot2, mRec, mRec2;
0626   // Overestimate of QED coupling.
0627   double alpha;
0628   // Parton system this antenna is part of.
0629   int iSys;
0630 
0631   // EW branchings.
0632   vector<EWBranching> brVec;
0633 
0634   // Trial variables.
0635   bool hasTrial;
0636   double q2Trial, sijTrial, sjkTrial;
0637   int poliTrial, poljTrial;
0638 
0639   // Outgoing momenta after branching.
0640   vector<Vec4> pNew;
0641 
0642   // Info on coefficents.
0643   double c0Sum, c1Sum, c2Sum, c3Sum;
0644   map<double, int> c0SumSoFar, c1SumSoFar, c2SumSoFar, c3SumSoFar;
0645 
0646   // Matching scale.
0647   double q2Match;
0648 
0649   // Information for partonSystems.
0650   int jNew;
0651   unordered_map<int,int> iReplace;
0652   double shat;
0653 
0654   // Pointers.
0655   EWBranching* brTrial{};
0656   Info* infoPtr{};
0657   Rndm* rndmPtr{};
0658   Logger* loggerPtr{};
0659   PartonSystems* partonSystemsPtr{};
0660   VinciaCommon* vinComPtr{};
0661   AlphaEM* alphaPtr{};
0662   AmpCalculator* ampCalcPtr{};
0663 
0664   // Settings.
0665   bool doBosonInterference;
0666 
0667   // Verbosity.
0668   int verbose;
0669 
0670 };
0671 
0672 //==========================================================================
0673 
0674 // Final-final electroweak antenna.
0675 
0676 class EWAntennaFF : public EWAntenna {
0677  public:
0678 
0679   // Overridden virtual functions.
0680   virtual bool init(Event &event, int iMotIn, int iRecIn, int iSysIn,
0681     vector<EWBranching> &branchings, Settings* settingsPtr) override;
0682   virtual double generateTrial(double q2Start, double q2End, double alphaIn)
0683     override;
0684   virtual bool acceptTrial(Event &event) override;
0685   virtual void updateEvent(Event &event) override;
0686 
0687  private:
0688 
0689   // Data members.
0690   double mAnt2, sqrtKallen;
0691   // Kinematic map.
0692   int kMapFinal;
0693   // Controls if resonances with too large offshellness are vetoed.
0694   bool vetoResonanceProduction;
0695 
0696 };
0697 
0698 //==========================================================================
0699 
0700 // Final-final electroweak resonance antenna.
0701 
0702 class EWAntennaFFres : public EWAntennaFF {
0703 
0704 public:
0705 
0706   // Overridden virtual functions.
0707   bool init(Event &event, int iMotIn, int iRecIn, int iSysIn,
0708     vector<EWBranching> &branchings, Settings* settingsPtr) override;
0709   bool isResonanceDecay() override {return true;}
0710   double generateTrial(double q2Start, double q2End, double alphaIn) override;
0711   bool acceptTrial(Event &event) override;
0712   void updateEvent(Event &event) override;
0713 
0714   // Generate the kinematics and channel for decays below matching scale.
0715   bool genForceDecay(Event &event);
0716 
0717 private:
0718 
0719   // Check if the trial was a resonance decay.
0720   bool trialIsResDecay;
0721   // Matching mode.
0722   int bwMatchMode;
0723   // Offshellness of the resonance and EW scale.
0724   double q2Dec, q2EW;
0725   // Switch for the special case of a resonance without a recoiler.
0726   bool doDecayOnly{false};
0727 
0728 };
0729 
0730 //==========================================================================
0731 
0732 // Initial-initial electroweak antenna.
0733 
0734 class EWAntennaII : public EWAntenna {
0735 
0736 public:
0737 
0738   // Constructor.
0739   EWAntennaII(BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn):
0740     beamAPtr(beamAPtrIn), beamBPtr(beamBPtrIn), shh(0), xMot(0), xRec(0),
0741       vetoResonanceProduction(false), TINYPDFtrial(1e-10) {;}
0742 
0743   // Overridden virtual functions.
0744   bool init(Event &event, int iMotIn, int iRecIn, int iSysIn,
0745     vector<EWBranching> &branchings,Settings* settingsPtr) override;
0746   bool isInitial() override {return true;}
0747   double generateTrial(double q2Start, double q2End, double alphaIn) override;
0748   bool acceptTrial(Event &event) override;
0749   void updateEvent(Event &event) override;
0750   void updatePartonSystems(Event &event) override;
0751 
0752 private:
0753 
0754   // Members.
0755   // Beam pointers.
0756   BeamParticle* beamAPtr{};
0757   BeamParticle* beamBPtr{};
0758   // Hadronic invariant mass.
0759   double shh;
0760   // Antenna hadronic momentum fractions.
0761   double xMot, xRec;
0762   // Controls if resonances with too large offshellness are vetoed.
0763   bool vetoResonanceProduction;
0764   // Global recoil momenta.
0765   vector<Vec4> pRecVec;
0766   vector<int>  iRecVec;
0767 
0768   // Tolerance for PDFs.
0769   double TINYPDFtrial;
0770 
0771 };
0772 
0773 //==========================================================================
0774 
0775 // Class that performs electroweak showers in a single parton system.
0776 
0777 class EWSystem {
0778 
0779 public:
0780 
0781   // Constructors.
0782   EWSystem(): antTrial(nullptr) {clearLastTrial();}
0783   EWSystem(unordered_map<pair<int, int>, vector<EWBranching> >* brMapFinalIn,
0784     unordered_map<pair<int, int>, vector<EWBranching> >* brMapInitialIn,
0785     unordered_map<pair<int, int>, vector<EWBranching> >* brMapResonanceIn,
0786     unordered_map<pair<int, int>, vector<pair<int, int> > >* cluMapFinalIn,
0787     unordered_map<pair<int, int>, vector<pair<int, int> > >* cluMapInitialIn,
0788     AmpCalculator * ampCalcIn):
0789     shh(0), iSysSav(0), resDecOnlySav(false), q2Cut(0), q2Trial(0),
0790     lastWasSplitSav(false), lastWasDecSav(false), lastWasInitialSav(false),
0791     lastWasBelowCut(false), ISav(0), KSav(0), brMapFinal(brMapFinalIn),
0792     brMapInitial(brMapInitialIn), brMapResonance(brMapResonanceIn),
0793     cluMapFinal(cluMapFinalIn), cluMapInitial(cluMapInitialIn),
0794     ampCalcPtr(ampCalcIn), isInit(false), doVetoHardEmissions(false),
0795     verbose(false), vetoHardEmissionsDeltaR2(0) {clearLastTrial();}
0796 
0797   // Initialize pointers.
0798   void initPtr(Info* infoPtrIn, VinciaCommon* vinComPtrIn, AlphaEM* alIn) {
0799     infoPtr = infoPtrIn;
0800     partonSystemsPtr = infoPtr->partonSystemsPtr;
0801     rndmPtr          = infoPtr->rndmPtr;
0802     settingsPtr      = infoPtr->settingsPtr;
0803     loggerPtr        = infoPtr->loggerPtr;
0804     vinComPtr = vinComPtrIn;
0805     al = alIn;}
0806 
0807   // Initialize.
0808   void init(BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn)  {
0809   beamAPtr = beamAPtrIn; beamBPtr = beamBPtrIn;
0810   doVetoHardEmissions = settingsPtr->flag("Vincia:EWoverlapVeto");
0811   vetoHardEmissionsDeltaR2 =
0812     pow2(settingsPtr->parm("Vincia:EWoverlapVetoDeltaR"));
0813   isInit = true;}
0814 
0815   // Set verbosity.
0816   void setVerbose(int verboseIn) {verbose = verboseIn;}
0817 
0818   // Prepare an event.
0819   bool prepare(Event &event, int iSysIn, double q2CutIn, bool resDecOnlyIn) {
0820     iSysSav = iSysIn; resDecOnlySav = resDecOnlyIn; q2Cut = q2CutIn;
0821     shh = infoPtr->s(); return buildSystem(event);}
0822 
0823   // Build a system.
0824   bool buildSystem(Event &event);
0825 
0826   // Return which system we are handling.
0827   int system() {return iSysSav;}
0828 
0829   // Generate the next Q2.
0830   double q2Next(double q2Start, double q2End);
0831 
0832   // Add an antenna.
0833   template <class T> void addAntenna(T ant, vector<T>& antVec,
0834     Event &event, int iMot, int iRec,
0835     unordered_map<pair<int, int>, vector<EWBranching> >* brMapPtr) {
0836     if (iMot == 0) return; // Check mother.
0837     int idA(event[iMot].id()), polA(event[iMot].pol());
0838     if (idA == 21) return; // Skip gluons.
0839     auto it = brMapPtr->find(make_pair(idA, polA));
0840     if (it != brMapPtr->end()) {
0841       // Found. Pass verbosity.
0842       ant.setVerbose(verbose);
0843       // Pass pointers.
0844       ant.initPtr(infoPtr, vinComPtr, al, ampCalcPtr);
0845       // Initialise and if success, store.
0846       if (ant.init(event, iMot, iRec, iSysSav, it->second, settingsPtr)) {
0847         antVec.push_back(std::move(ant));
0848         if (verbose >= VinciaConstants::DEBUG) {
0849           stringstream ss;
0850           ss << "Added EW antenna with iEv = "
0851              << iMot << " and iRec = "<< iRec<< " in system "<< iSysSav;
0852           printOut(__METHOD_NAME__, ss.str());}}}}
0853 
0854   // Generate a trial.
0855   template <class T> void generateTrial(vector<T> & antVec, double q2Start,
0856     double q2End, double alphaIn) {
0857     if (q2End > q2Start) return;
0858     for (int i = 0; i < (int)antVec.size(); i++) {
0859       // Generate a trial scale for the current antenna.
0860       double q2New = antVec[i].generateTrial(q2Start,q2End,alphaIn);
0861       // Current winner.
0862       if (q2New > q2Trial && q2New > q2End) {
0863         // Save trial information.
0864         q2Trial = q2New;
0865         antTrial = &(antVec[i]);
0866         lastWasDecSav = antTrial->isResonanceDecay();
0867         lastWasInitialSav = antTrial->isInitial();
0868         // This is done to avoid issues with resonance antennae
0869         // not deciding their channel until acceptTrial.
0870         lastWasSplitSav = lastWasDecSav ? true : antTrial->isSplitToFermions();
0871         lastWasBelowCut = (q2Trial < q2Cut)? true : false;
0872         ISav = antTrial->getIndexMot(); KSav = antTrial->getIndexRec();}}}
0873 
0874   // Overloaded version passing event, just for resonance decays.
0875   void generateTrial(Event &event, vector<EWAntenna> & antVec,double q2Start,
0876     double q2End, double alphaIn);
0877 
0878   // Accept a trial.
0879   bool acceptTrial(Event &event) {
0880     bool passed = antTrial->acceptTrial(event);
0881     if (verbose >= VinciaConstants::DEBUG)
0882       printOut(__METHOD_NAME__, passed ? "Passed veto" : "Vetoed branching");
0883     return passed;}
0884 
0885   // Update an event.
0886   void updateEvent(Event &event) {
0887     if (verbose >= VinciaConstants::DEBUG)
0888       printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
0889     if (antTrial != nullptr) antTrial->updateEvent(event);
0890     else loggerPtr->ERROR_MSG("trial doesn't exist!");
0891     if (verbose >= VinciaConstants::DEBUG)
0892       printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);}
0893 
0894   // Update parton systems.
0895   void updatePartonSystems(Event &event) {
0896     if (verbose >= VinciaConstants::DEBUG)
0897       printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
0898     if (antTrial!=nullptr) antTrial->updatePartonSystems(event);
0899     else loggerPtr->ERROR_MSG("trial doesn't exist!");
0900     if (verbose >= VinciaConstants::DEBUG)
0901       printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);}
0902 
0903   // Print the antennas.
0904   void printAntennae() {
0905     for (int i = 0; i < (int)antVecFinal.size(); i++) antVecFinal[i].print();
0906     for (int i = 0; i < (int)antVecRes.size(); i++) antVecRes[i].print();
0907     for (int i = 0; i < (int)antVecInitial.size(); i++)
0908       antVecInitial[i].print();}
0909 
0910   // Get number of branchers (total and res-dec only)
0911   unsigned int nBranchers() {return antVecFinal.size() +
0912       antVecInitial.size() + antVecRes.size();}
0913   unsigned int nResDec() {return antVecRes.size();}
0914 
0915   // Check if system has a trial pointer.
0916   bool hasTrial() {return antTrial != nullptr;}
0917 
0918   // Check if split to fermions, initial, or resonance.
0919   bool lastWasSplitToFermions() {return lastWasSplitSav;}
0920   bool lastWasInitial() {return lastWasInitialSav;}
0921   bool lastWasResonanceDecay() {return lastWasDecSav;}
0922 
0923   // Clear the antennae.
0924   void clearAntennae() {antVecFinal.clear(); antVecInitial.clear();
0925     antVecRes.clear(); clearLastTrial();}
0926 
0927   // Clear the last saved trial (but keep the antennae).
0928   void clearLastTrial() {
0929     q2Trial = 0; antTrial = nullptr; lastWasSplitSav = false;
0930     lastWasDecSav = false; lastWasInitialSav = false; lastWasBelowCut = false;
0931     ISav = 0; KSav = 0;}
0932 
0933 private:
0934 
0935   // Members.
0936   // Hadronic invariant mass.
0937   double shh;
0938 
0939   // System and whether to do only resonance decays or full EW shower.
0940   int iSysSav;
0941   bool resDecOnlySav;
0942 
0943   // Cutoff scale.
0944   double q2Cut;
0945 
0946   // Pointers.
0947   BeamParticle* beamAPtr{};
0948   BeamParticle* beamBPtr{};
0949   Info* infoPtr{};
0950   PartonSystems* partonSystemsPtr{};
0951   Rndm* rndmPtr{};
0952   Settings* settingsPtr{};
0953   Logger* loggerPtr{};
0954   VinciaCommon* vinComPtr{};
0955   AlphaEM* al{};
0956 
0957   // Antennae.
0958   vector<EWAntennaFF> antVecFinal;
0959   vector<EWAntennaII> antVecInitial;
0960   vector<EWAntennaFFres> antVecRes;
0961 
0962   // Trial information.
0963   EWAntenna* antTrial{};
0964   double q2Trial;
0965   bool lastWasSplitSav, lastWasDecSav, lastWasInitialSav, lastWasBelowCut;
0966   int ISav, KSav;
0967 
0968   // Map of branchings.
0969   unordered_map<pair<int, int>, vector<EWBranching> >
0970     *brMapFinal{}, *brMapInitial{}, *brMapResonance{};
0971 
0972   // Cluster maps for spectator selection.
0973   unordered_map<pair<int, int>, vector<pair<int, int> > >
0974     *cluMapFinal{}, *cluMapInitial{};
0975 
0976   // Amplitude calculator.
0977   AmpCalculator* ampCalcPtr{};
0978 
0979   // Flags.
0980   bool isInit, doVetoHardEmissions;
0981   int verbose;
0982   double vetoHardEmissionsDeltaR2;
0983 
0984 };
0985 
0986 //==========================================================================
0987 
0988 // Top-level class for the electroweak shower module.
0989 
0990 class VinciaEW : public VinciaModule {
0991 
0992 public:
0993 
0994   // Constructor.
0995   VinciaEW(): isLoaded{false} {;}
0996 
0997   // Initialize pointers (called at construction time).
0998   void initPtr(Info* infoPtrIn, VinciaCommon* vinComPtrIn) override {
0999     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.
1119   EWParticleData ewData;
1120 
1121   // Amplitude calculator.
1122   AmpCalculator ampCalc;
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