MatrixElement/Hadron/MEQCD2to2Fast.h

0001 // -*- C++ -*-
0002 //
0003 // MEQCD2to2Fast.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_MEQCD2to2Fast_H
0010 #define HERWIG_MEQCD2to2Fast_H
0011 //
0012 // This is the declaration of the MEQCD2to2Fast class.
0013 //
0014
0015 #include "Herwig/MatrixElement/HwMEBase.h"
0016 #include "ThePEG/Repository/UseRandom.h"
0017
0018 namespace Herwig {
0019 using namespace ThePEG;
0020 using namespace ThePEG::Helicity;
0021
0022 /**
0023  * The MEQCD2to2Fast class implements the matrix elements for
0024  * QCD \f$2\to2\f$ scattering processes using hard coded formulae and
0025  * as such can not include spin correlations. It is designed to be a faster
0026  * replacement for MEQCD2to2 for use in the underlying event.
0027  *
0028  * @see \ref MEQCD2to2FastInterfaces "The interfaces"
0029  * defined for MEQCD2to2Fast.
0030  */
0031 class MEQCD2to2Fast: public HwMEBase {
0032
0033 public:
0034
0035   /**
0036    * The default constructor.
0037    */
0038   MEQCD2to2Fast() :_maxflavour(5),_process(0),_strictFlavourScheme(false) {
0039     massOption(vector<unsigned int>(2,0));
0040   }
0041
0042   /** @name Virtual functions required by the MEBase class. */
0043   //@{
0044   /**
0045    * Return the order in \f$\alpha_S\f$ in which this matrix
0046    * element is given.
0047    */
0048   virtual unsigned int orderInAlphaS() const;
0049
0050   /**
0051    * Return the order in \f$\alpha_{EW}\f$ in which this matrix
0052    * element is given.
0053    */
0054   virtual unsigned int orderInAlphaEW() const;
0055
0056   /**
0057    * The matrix element for the kinematical configuration
0058    * previously provided by the last call to setKinematics(), suitably
0059    * scaled by sHat() to give a dimension-less number.
0060    * @return the matrix element scaled with sHat() to give a
0061    * dimensionless number.
0062    */
0063   virtual double me2() const;
0064
0065   /**
0066    * Return the scale associated with the last set phase space point.
0067    */
0068   virtual Energy2 scale() const;
0069
0070   /**
0071    * Add all possible diagrams with the add() function.
0072    */
0073   virtual void getDiagrams() const;
0074
0075   /**
0076    * Get diagram selector. With the information previously supplied with the
0077    * setKinematics method, a derived class may optionally
0078    * override this method to weight the given diagrams with their
0079    * (although certainly not physical) relative probabilities.
0080    * @param dv the diagrams to be weighted.
0081    * @return a Selector relating the given diagrams to their weights.
0082    */
0083   virtual Selector<DiagramIndex> diagrams(const DiagramVector & dv) const;
0084
0085   /**
0086    * Return a Selector with possible colour geometries for the selected
0087    * diagram weighted by their relative probabilities.
0088    * @param diag the diagram chosen.
0089    * @return the possible colour geometries weighted by their
0090    * relative probabilities.
0091    */
0092   virtual Selector<const ColourLines *>
0093   colourGeometries(tcDiagPtr diag) const;
0094   //@}
0095
0096
0097 public:
0098
0099   /** @name Functions used by the persistent I/O system. */
0100   //@{
0101   /**
0102    * Function used to write out object persistently.
0103    * @param os the persistent output stream written to.
0104    */
0105   void persistentOutput(PersistentOStream & os) const;
0106
0107   /**
0108    * Function used to read in object persistently.
0109    * @param is the persistent input stream read from.
0110    * @param version the version number of the object when written.
0111    */
0112   void persistentInput(PersistentIStream & is, int version);
0113   //@}
0114
0115   /**
0116    * The standard Init function used to initialize the interfaces.
0117    * Called exactly once for each class by the class description system
0118    * before the main function starts or
0119    * when this class is dynamically loaded.
0120    */
0121   static void Init();
0122
0123 protected:
0124
0125   /**
0126    *  Members to calculate the matrix elements
0127    */
0128   //@{
0129   /**
0130    * Matrix element for \f$gg\to gg\f$.
0131    */
0132   double gg2ggME() const {
0133     Energy2 u(uHat()),t(tHat()),s(sHat());
0134     double output = 9./4.*(3.-t*u/s/s-s*u/t/t-s*t/u/u);
0135     double flow[3]={(1.-u*t/s/s-s*t/u/u+t*t/s/u),
0136             (1.-t*u/s/s-s*u/t/t+u*u/s/t),
0137             (1.-t*s/u/u-u*s/t/t+s*s/u/t)};
0138     _flow = 1+UseRandom::rnd3(flow[0],flow[1],flow[2]);
0139     double diag[3]={(sqr(u)+sqr(t))/sqr(s),
0140             (sqr(s)+sqr(u))/sqr(t),
0141             (sqr(s)+sqr(t))/sqr(u)};
0142     if(_flow==1)      diag[1]=0;
0143     else if(_flow==2) diag[2]=0;
0144     else if(_flow==3) diag[0]=0;
0145     _diagram=1+UseRandom::rnd3(diag[0],diag[1],diag[2]);
0146     return output;
0147   }
0148
0149   /**
0150    * Matrix element for \f$gg\to q\bar{q}\f$
0151    */
0152   double gg2qqbarME() const {
0153     Energy2 u(uHat()),t(tHat()),s(sHat());
0154     Energy4 u2(sqr(u)),t2(sqr(t)),s2(sqr(s));
0155     double output =(1./6./u/t-3./8./s2)*(t2+u2);
0156     double flow[2]={u2/(u2+t2),t2/(u2+t2)};
0157     _flow = 1+UseRandom::rnd2(flow[0],flow[1]);
0158     _diagram=3+_flow;
0159     return output;
0160   }
0161
0162   /**
0163    * Matrix element for \f$q\bar{q}\to gg\f$
0164    */
0165   double qqbar2ggME() const {
0166     Energy2 u(uHat()),t(tHat()),s(sHat());
0167     Energy4 s2(sqr(s)),u2(sqr(u)),t2(sqr(t));
0168     double output = 0.5*(32./27./u/t-8./3./s2)*(t2+u2);
0169     double flow[2] = {u2/(u2+t2),t2/(t2+u2)};
0170     _flow=1+UseRandom::rnd2(flow[0],flow[1]);
0171     _diagram=6+_flow;
0172     return output;
0173   }
0174
0175   /**
0176    * Matrix element for \f$qg\to qg\f$
0177    */
0178   double qg2qgME() const {
0179     Energy2 u(uHat()),t(tHat()),s(sHat());
0180     Energy4 s2(sqr(s)),u2(sqr(u)),t2(sqr(t));
0181     double output = (-4./9./s/u+1./t2)*(s2+u2);
0182     double flow[2]={u2/(s2+u2),s2/(s2+u2)};
0183     _flow=1+UseRandom::rnd2(flow[0],flow[1]);
0184     _diagram=9+_flow;
0185     return output;
0186   }
0187
0188   /**
0189    * Matrix elements for \f$\bar{q}g\to \bar{q}g\f$.
0190    */
0191   double qbarg2qbargME() const {
0192     // scale
0193     Energy2 u(uHat()),t(tHat()),s(sHat());
0194     Energy4 u2(sqr(u)),s2(sqr(s)); // t2(sqr(t))
0195     double flow[2]={u2/(s2+u2),s2/(s2+u2)};
0196     _flow=1+UseRandom::rnd2(flow[0],flow[1]);
0197     _diagram=12+_flow;
0198     return (-4./9./s/u+1./t/t)*(s*s+u*u);
0199   }
0200
0201   /**
0202    * Matrix element for \f$qq\to qq\f$
0203    */
0204   double qq2qqME() const {
0205     Energy2 u(uHat()),t(tHat());
0206     Energy4 s2(sqr(sHat())),u2(sqr(u)),t2(sqr(t));
0207     double output;
0208     if(mePartonData()[0]->id()==mePartonData()[1]->id()) {
0209       output = 0.5*(4./9.*((s2+u2)/t2+(s2+t2)/u2)
0210             -8./27.*s2/u/t);
0211       double flow[2]={(s2+u2)/t2,(s2+t2)/u2};
0212       _flow=1+UseRandom::rnd2(flow[0],flow[1]);
0213     }
0214     else {
0215       output = 4./9.*(s2+u2)/t2;
0216       _flow=2;
0217     }
0218     _diagram = 15+_flow;
0219     return output;
0220   }
0221
0222   /**
0223    * Matrix element for \f$\bar{q}\bar{q}\to \bar{q}\bar{q}\f$
0224    */
0225   double qbarqbar2qbarqbarME() const {
0226     Energy2 u(uHat()),t(tHat());
0227     Energy4 u2(sqr(u)),t2(sqr(t)),s2(sqr(sHat()));
0228     double output;
0229     if(mePartonData()[0]->id()==mePartonData()[1]->id()) {
0230       output = 0.5*(4./9.*((s2+u2)/t2+(s2+t2)/u2)
0231             -8./27.*s2/u/t);
0232       double flow[2]={(s2+u2)/t2,(s2+t2)/u2};
0233       _flow=1+UseRandom::rnd2(flow[0],flow[1]);
0234     }
0235     else {
0236       output = 4./9.*(s2+u2)/t2;
0237       _flow = 2;
0238     }
0239     _diagram = 17+_flow;
0240     // final part of colour and spin factors
0241     return output;
0242   }
0243
0244   /**
0245    * Matrix element for \f$q\bar{q}\to q\bar{q}\f$
0246    */
0247   double qqbar2qqbarME() const {
0248     // type of process
0249     bool diagon[2]={mePartonData()[0]->id()== -mePartonData()[1]->id(),
0250             mePartonData()[0]->id()==  mePartonData()[2]->id()};
0251     // scale
0252     Energy2 u(uHat()),t(tHat()),s(sHat());
0253     Energy4 s2(sqr(s)),t2(sqr(t)),u2(sqr(u));
0254     double output;
0255     if(diagon[0]&&diagon[1]) {
0256       output= (4./9.*((s2+u2)/t2+(u2+t2)/s2)
0257            -8./27.*u2/s/t);
0258       double flow[2]={(t2+u2)/s2,(s2+u2)/t2};
0259       _flow=1+UseRandom::rnd2(flow[0],flow[1]);
0260     }
0261     else if(diagon[0]) {
0262       output = (4./9.*(t2+u2)/s2);
0263       _flow=1;
0264     }
0265     else {
0266       output = (4./9.*(s2+u2)/t2);
0267       _flow=2;
0268     }
0269     _diagram=19+_flow;
0270     return output;
0271   }
0272   //@}
0273
0274 protected:
0275
0276   /** @name Clone Methods. */
0277   //@{
0278   /**
0279    * Make a simple clone of this object.
0280    * @return a pointer to the new object.
0281    */
0282   virtual IBPtr clone() const;
0283
0284   /** Make a clone of this object, possibly modifying the cloned object
0285    * to make it sane.
0286    * @return a pointer to the new object.
0287    */
0288   virtual IBPtr fullclone() const;
0289   //@}
0290
0291 protected:
0292
0293   /**
0294    * Initialize this object after the setup phase before saving an
0295    * EventGenerator to disk.
0296    * @throws InitException if object could not be initialized properly.
0297    */
0298   virtual void doinit();
0299
0300   //@}
0301
0302 private:
0303
0304   /**
0305    * The assignment operator is private and must never be called.
0306    * In fact, it should not even be implemented.
0307    */
0308   MEQCD2to2Fast & operator=(const MEQCD2to2Fast &) = delete;
0309
0310 private:
0311
0312   /**
0313    *  Maximum numbere of quark flavours to include
0314    */
0315   unsigned int _maxflavour;
0316
0317   /**
0318    *  Processes to include
0319    */
0320   unsigned int _process;
0321
0322   /**
0323    *  Colour flow
0324    */
0325   mutable unsigned int _flow;
0326
0327   /**
0328    *  Diagram
0329    */
0330   mutable unsigned int _diagram;
0331
0332   /**
0333    * Exclude contributions with massive incominbg quarks
0334    */
0335   bool _strictFlavourScheme;
0336
0337 };
0338
0339 }
0340
0341 #endif /* HERWIG_MEQCD2to2Fast_H */