Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 10:06:34

0001 // VinciaEW.h is a part of the PYTHIA event generator.
0002 // Copyright (C) 2024 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{}, *beamBPtr{};
0948   Info* infoPtr{};
0949   PartonSystems* partonSystemsPtr{};
0950   Rndm* rndmPtr{};
0951   Settings* settingsPtr{};
0952   Logger* loggerPtr{};
0953   VinciaCommon* vinComPtr{};
0954   AlphaEM* al{};
0955 
0956   // Antennae.
0957   vector<EWAntennaFF> antVecFinal;
0958   vector<EWAntennaII> antVecInitial;
0959   vector<EWAntennaFFres> antVecRes;
0960 
0961   // Trial information.
0962   EWAntenna* antTrial{};
0963   double q2Trial;
0964   bool lastWasSplitSav, lastWasDecSav, lastWasInitialSav, lastWasBelowCut;
0965   int ISav, KSav;
0966 
0967   // Map of branchings.
0968   unordered_map<pair<int, int>, vector<EWBranching> >
0969     *brMapFinal{}, *brMapInitial{}, *brMapResonance{};
0970 
0971   // Cluster maps for spectator selection.
0972   unordered_map<pair<int, int>, vector<pair<int, int> > >
0973     *cluMapFinal{}, *cluMapInitial{};
0974 
0975   // Amplitude calculator.
0976   AmpCalculator* ampCalcPtr{};
0977 
0978   // Flags.
0979   bool isInit, doVetoHardEmissions;
0980   int verbose;
0981   double vetoHardEmissionsDeltaR2;
0982 
0983 };
0984 
0985 //==========================================================================
0986 
0987 // Top-level class for the electroweak shower module.
0988 
0989 class VinciaEW : public VinciaModule {
0990 
0991 public:
0992 
0993   // Constructor.
0994   VinciaEW(): isLoaded{false} {;}
0995 
0996   // Initialize pointers (called at construction time).
0997   void initPtr(Info* infoPtrIn, VinciaCommon* vinComPtrIn) override {
0998     infoPtr = infoPtrIn;
0999     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.
1118   EWParticleData ewData;
1119 
1120   // Amplitude calculator.
1121   AmpCalculator ampCalc;
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