Decay/Perturbative/SMWDecayer.h

0001 // -*- C++ -*-
0002 //
0003 // SMWDecayer.h is a part of Herwig - A multi-purpose Monte Carlo event generator
0004 // Copyright (C) 2002-2019 The Herwig Collaboration
0005 //
0006 // Herwig is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 #ifndef HERWIG_SMWDecayer_H
0010 #define HERWIG_SMWDecayer_H
0011 //
0012 // This is the declaration of the SMWDecayer class.
0013 //
0014 #include "Herwig/Decay/PerturbativeDecayer.h"
0015 #include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"
0016 #include "ThePEG/Helicity/Vertex/AbstractVVVVertex.h"
0017 #include "Herwig/Decay/PhaseSpaceMode.h"
0018
0019 namespace Herwig {
0020 using namespace ThePEG;
0021 using namespace ThePEG::Helicity;
0022
0023 /** \ingroup Decay
0024  *
0025  *  The <code>SMWDecayer</code> is designed to perform the decay of the
0026  *  W boson to the Standard Model fermions, including the first order
0027  *  electroweak corrections.
0028  *
0029  * @see PerturbativeDecayer
0030  *
0031  */
0032 class SMWDecayer: public PerturbativeDecayer {
0033
0034 public:
0035
0036   /**
0037    * Default constructor.
0038    */
0039   SMWDecayer();
0040
0041 public:
0042
0043   /**
0044    *  Virtual members to be overridden by inheriting classes
0045    *  which implement hard corrections
0046    */
0047   //@{
0048   /**
0049    *  Has an old fashioned ME correction
0050    */
0051   virtual bool hasMECorrection() {return true;}
0052
0053   /**
0054    *  Initialize the ME correction
0055    */
0056   virtual void initializeMECorrection(RealEmissionProcessPtr , double & ,
0057                       double & );
0058
0059   /**
0060    * Apply the soft matrix element correction
0061    * @param parent The initial particle in the current branching
0062    * @param progenitor The progenitor particle of the jet
0063    * @param fs Whether the emission is initial or final-state
0064    * @param highestpT The highest pT so far in the shower
0065    * @param ids ids of the particles produced in the branching
0066    * @param z The momentum fraction of the branching
0067    * @param scale the evolution scale of the branching
0068    * @param pT The transverse momentum of the branching
0069    * @return If true the emission should be vetoed
0070    */
0071   virtual bool softMatrixElementVeto(PPtr parent,
0072                      PPtr progenitor,
0073                      const bool & fs,
0074                      const Energy & highestpT,
0075                      const vector<tcPDPtr> & ids,
0076                      const double & z,
0077                      const Energy & scale,
0078                      const Energy & pT);
0079
0080   /**
0081    *  Has a POWHEG style correction
0082    */
0083   virtual POWHEGType hasPOWHEGCorrection() {return FSR;}
0084
0085 public:
0086
0087   /**
0088    * Which of the possible decays is required
0089    * @param cc Is this mode the charge conjugate
0090    * @param parent The decaying particle
0091    * @param children The decay products
0092    */
0093   virtual int modeNumber(bool & cc, tcPDPtr parent,
0094              const tPDVector & children) const;
0095
0096   /**
0097    * For a given decay mode and a given particle instance, perform the
0098    * decay and return the decay products. As this is the base class this
0099    * is not implemented.
0100    * @return The vector of particles produced in the decay.
0101    */
0102   virtual ParticleVector decay(const Particle & parent,const tPDVector & children) const;
0103
0104   /**
0105    * Return the matrix element squared for a given mode and phase-space channel.
0106    * @param ichan The channel we are calculating the matrix element for.
0107    * @param part The decaying Particle.
0108    * @param outgoing The particles produced in the decay
0109    * @param momenta  The momenta of the particles produced in the decay
0110    * @param meopt Option for the calculation of the matrix element
0111    * @return The matrix element squared for the phase-space configuration.
0112    */
0113   double me2(const int ichan,const Particle & part,
0114          const tPDVector & outgoing,
0115          const vector<Lorentz5Momentum> & momenta,
0116          MEOption meopt) const;
0117
0118   /**
0119    *   Construct the SpinInfos for the particles produced in the decay
0120    */
0121   virtual void constructSpinInfo(const Particle & part,
0122                  ParticleVector outgoing) const;
0123
0124   /**
0125    * Output the setup information for the particle database
0126    * @param os The stream to output the information to
0127    * @param header Whether or not to output the information for MySQL
0128    */
0129   virtual void dataBaseOutput(ofstream & os,bool header) const;
0130
0131 public:
0132
0133   /** @name Functions used by the persistent I/O system. */
0134   //@{
0135   /**
0136    * Function used to write out object persistently.
0137    * @param os the persistent output stream written to.
0138    */
0139   void persistentOutput(PersistentOStream & os) const;
0140
0141   /**
0142    * Function used to read in object persistently.
0143    * @param is the persistent input stream read from.
0144    * @param version the version number of the object when written.
0145    */
0146   void persistentInput(PersistentIStream & is, int version);
0147   //@}
0148
0149   /**
0150    * Standard Init function used to initialize the interfaces.
0151    */
0152   static void Init();
0153
0154 protected:
0155
0156   /** @name Clone Methods. */
0157   //@{
0158   /**
0159    * Make a simple clone of this object.
0160    * @return a pointer to the new object.
0161    */
0162   virtual IBPtr clone() const {return new_ptr(*this);}
0163
0164   /** Make a clone of this object, possibly modifying the cloned object
0165    * to make it sane.
0166    * @return a pointer to the new object.
0167    */
0168   virtual IBPtr fullclone() const {return new_ptr(*this);}
0169   //@}
0170
0171 protected:
0172
0173   /** @name Standard Interfaced functions. */
0174   //@{
0175   /**
0176    * Initialize this object after the setup phase before saving and
0177    * EventGenerator to disk.
0178    * @throws InitException if object could not be initialized properly.
0179    */
0180   virtual void doinit();
0181
0182   /**
0183    * Initialize this object. Called in the run phase just before
0184    * a run begins.
0185    */
0186   virtual void doinitrun();
0187   //@}
0188
0189 protected:
0190
0191   /**
0192    *  Set the \f$\rho\f$ parameter
0193    */
0194   void setRho(double);
0195
0196   /**
0197    *  Set the \f$\tilde{\kappa}\f$ parameters symmetrically
0198    */
0199   void setKtildeSymm();
0200
0201   /**
0202    * Set second \f$\tilde{\kappa}\f$, given the first.
0203    */
0204   void setKtilde2();
0205
0206   /**
0207    *  Translate the variables from \f$x_q,x_{\bar{q}}\f$ to \f$\tilde{\kappa},z\f$
0208    */
0209   //@{
0210   /**
0211    *  Calculate \f$z\f$.
0212    */
0213   double getZfromX(double, double);
0214
0215   /**
0216    *  Calculate \f$\tilde{\kappa}\f$.
0217    */
0218   double getKfromX(double, double);
0219   //@}
0220
0221   /**
0222    * Calculate \f$x_{q},x_{\bar{q}}\f$ from \f$\tilde{\kappa},z\f$.
0223    * @param kt \f$\tilde{\kappa}\f$
0224    * @param z \f$z\f$
0225    * @param x \f$x_{q}\f$
0226    * @param xbar \f$x_{\bar{q}}\f$
0227    */
0228   void getXXbar(double kt, double z, double & x, double & xbar);
0229
0230   /**
0231    *  Soft weight
0232    */
0233   //@{
0234   /**
0235    *  Soft quark weight calculated from \f$x_{q},x_{\bar{q}}\f$
0236    * @param x \f$x_{q}\f$
0237    * @param xbar \f$x_{\bar{q}}\f$
0238    */
0239   double qWeight(double x, double xbar);
0240
0241   /**
0242    *  Soft antiquark weight calculated from \f$x_{q},x_{\bar{q}}\f$
0243    * @param x \f$x_{q}\f$
0244    * @param xbar \f$x_{\bar{q}}\f$
0245    */
0246   double qbarWeight(double x, double xbar);
0247
0248   /**
0249    * Soft quark weight calculated from \f$\tilde{q},z\f$
0250    * @param qtilde  \f$\tilde{q}\f$
0251    * @param z \f$z\f$
0252    */
0253   double qWeightX(Energy qtilde, double z);
0254
0255   /**
0256    * Soft antiquark weight calculated from \f$\tilde{q},z\f$
0257    * @param qtilde  \f$\tilde{q}\f$
0258    * @param z \f$z\f$
0259    */
0260   double qbarWeightX(Energy qtilde, double z);
0261   //@}
0262
0263   /**
0264    * ????
0265    */
0266   double u(double);
0267
0268   /**
0269    *  Vector and axial vector parts of the matrix element
0270    */
0271   //@{
0272   /**
0273    *  Vector part of the matrix element
0274    */
0275   double MEV(double, double);
0276
0277   /**
0278    *  Axial vector part of the matrix element
0279    */
0280   double MEA(double, double);
0281
0282   /**
0283    * The matrix element, given \f$x_1\f$, \f$x_2\f$.
0284    * @param x1 \f$x_1\f$
0285    * @param x2 \f$x_2\f$
0286    */
0287   double PS(double x1, double x2);
0288   //@}
0289
0290 protected:
0291
0292   /**
0293    *  Real emission term, for use in generating the hardest emission
0294    */
0295   double calculateRealEmission(double x1, double x2,
0296                    tcPDVector outgoing,
0297                    vector<Lorentz5Momentum> momenta,
0298                    double phi, double muj, double muk,
0299                    int iemit, bool subtract) const;
0300
0301   /**
0302    *  Calculate the ratio between NLO & LO ME
0303    */
0304   double meRatio(tcPDVector partons,
0305          vector<Lorentz5Momentum> momenta,
0306          unsigned int iemitter,bool subtract) const;
0307
0308   /**
0309    *  Calculate matrix element ratio R/B
0310    */
0311   virtual double matrixElementRatio(const Particle & inpart, const ParticleVector & decay2,
0312                     const ParticleVector & decay3, MEOption meopt,
0313                     ShowerInteraction inter);
0314
0315   /**
0316    *  Calculate the LO ME
0317    */
0318   double loME(const vector<tcPDPtr> & partons,
0319           const vector<Lorentz5Momentum> & momenta) const;
0320
0321   /**
0322    *  Calculate the NLO real emission piece of ME
0323    */
0324   InvEnergy2 realME(const vector<tcPDPtr> & partons,
0325             const vector<Lorentz5Momentum> & momenta,
0326             ShowerInteraction inter) const;
0327
0328 private:
0329
0330   /**
0331    * Private and non-existent assignment operator.
0332    */
0333   SMWDecayer & operator=(const SMWDecayer &) = delete;
0334
0335  private:
0336
0337
0338   /**
0339    *  Pointer to the fermion-antifermion W vertex
0340    */
0341   AbstractFFVVertexPtr FFWVertex_;
0342
0343   /**
0344    *  Pointer to the fermion-antifermion G vertex
0345    */
0346   AbstractFFVVertexPtr FFGVertex_;
0347
0348   /**
0349    *  Pointer to the fermion-antifermion G vertex
0350    */
0351   AbstractFFVVertexPtr FFPVertex_;
0352
0353   /**
0354    *  Pointer to the fermion-antifermion G vertex
0355    */
0356   AbstractVVVVertexPtr WWWVertex_;
0357
0358   /**
0359    * maximum weights for the different integrations
0360    */
0361   //@{
0362   /**
0363    *  Weights for the W to quarks decays.
0364    */
0365   vector<double> quarkWeight_;
0366
0367   /**
0368    *  Weights for the W to leptons decays.
0369    */
0370   vector<double> leptonWeight_;
0371   //@}
0372
0373   /**
0374    *  Spin density matrix for the decay
0375    */
0376   mutable RhoDMatrix rho_;
0377
0378   /**
0379    *  Polarization vectors for the decay
0380    */
0381   mutable vector<VectorWaveFunction> vectors_;
0382
0383   /**
0384    *  Spinors for the decay
0385    */
0386   mutable vector<SpinorWaveFunction> wave_;
0387
0388   /**
0389    *  Barred spinors for the decay
0390    */
0391   mutable vector<SpinorBarWaveFunction> wavebar_;
0392
0393 private:
0394
0395   /**
0396    * CM energy
0397    */
0398   Energy d_Q_;
0399
0400   /**
0401    *  Quark mass
0402    */
0403   Energy d_m_;
0404
0405   /**
0406    * The rho parameter
0407    */
0408   double d_rho_;
0409
0410   /**
0411    * The v parameter
0412    */
0413   double d_v_;
0414
0415   /**
0416    * The initial kappa-tilde values for radiation from the quark
0417    */
0418   double d_kt1_;
0419
0420   /**
0421    * The initial kappa-tilde values for radiation from the antiquark
0422    */
0423   double d_kt2_;
0424
0425   /**
0426    *  Cut-off parameter
0427    */
0428   static const double EPS_;
0429
0430 private:
0431
0432   /**
0433    *  The colour factor
0434    */
0435   double CF_;
0436
0437   /**
0438    *  The W mass
0439    */
0440   mutable Energy mW_;
0441
0442
0443   /**
0444    *  The strong coupling
0445    */
0446   mutable double aS_;
0447
0448   /**
0449    * The scale
0450    */
0451   mutable Energy2 scale_;
0452
0453   /**
0454    *  Stuff for the POWHEG correction
0455    */
0456   //@{
0457   /**
0458    *  ParticleData object for the gluon
0459    */
0460   tcPDPtr gluon_;
0461
0462   /**
0463    *  The ParticleData objects for the fermions
0464    */
0465   vector<tcPDPtr> partons_;
0466
0467   /**
0468    * The fermion momenta
0469    */
0470   vector<Lorentz5Momentum> quark_;
0471
0472   /**
0473    *  The momentum of the radiated gauge boson
0474    */
0475   Lorentz5Momentum gauge_;
0476
0477   /**
0478    *  The W boson
0479    */
0480   PPtr wboson_;
0481
0482   /**
0483    *  W mass squared
0484    */
0485   Energy2 mw2_;
0486   //@}
0487
0488   /**
0489    *  Whether ro return the LO or NLO result
0490    */
0491   bool NLO_;
0492 };
0493
0494 }
0495
0496
0497 #endif /* HERWIG_SMWDecayer_H */