MatrixElement/Powheg/MEPP2GammaGammaPowheg.h

0001 // -*- C++ -*-
0002 #ifndef HERWIG_MEPP2GammaGammaPowheg_H
0003 #define HERWIG_MEPP2GammaGammaPowheg_H
0004 //
0005 // This is the declaration of the MEPP2GammaGammaPowheg class.
0006 //
0007
0008 #include "Herwig/MatrixElement/HwMEBase.h"
0009 #include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"
0010 #include "Herwig/MatrixElement/ProductionMatrixElement.h"
0011 #include "Herwig/Shower/ShowerAlpha.h"
0012
0013 namespace Herwig {
0014
0015 using namespace ThePEG;
0016
0017 /**
0018  * Here is the documentation of the MEPP2GammaGammaPowheg class.
0019  *
0020  * @see \ref MEPP2GammaGammaPowhegInterfaces "The interfaces"
0021  * defined for MEPP2GammaGammaPowheg.
0022  */
0023 class MEPP2GammaGammaPowheg: public HwMEBase {
0024
0025   enum DipoleType {           IIQCD1=2,           IIQCD2=4,
0026             IFQED1=5, FIQED1=6, IFQED2=7, FIQED2=8 };
0027
0028   enum RadiationType {Subtraction,Hard,Shower};
0029
0030 public:
0031
0032   /** @name Standard constructors and destructors. */
0033   //@{
0034   /**
0035    * The default constructor.
0036    */
0037   MEPP2GammaGammaPowheg();
0038   //@}
0039
0040   /** @name Member functions for the generation of hard QCD radiation */
0041   //@{
0042   /**
0043    *  Has a POWHEG style correction
0044    */
0045   virtual POWHEGType hasPOWHEGCorrection() {return ISR;}
0046
0047   /**
0048    *  Apply the POWHEG style correction
0049    */
0050   virtual RealEmissionProcessPtr generateHardest(RealEmissionProcessPtr,
0051                          ShowerInteraction);
0052   //@}
0053
0054 public:
0055
0056   /** @name Virtual functions required by the MEBase class. */
0057   //@{
0058   /**
0059    * Return the order in \f$\alpha_S\f$ in which this matrix
0060    * element is given.
0061    */
0062   virtual unsigned int orderInAlphaS() const;
0063
0064   /**
0065    * Return the order in \f$\alpha_{EW}\f$ in which this matrix
0066    * element is given.
0067    */
0068   virtual unsigned int orderInAlphaEW() const;
0069
0070   /**
0071    * The matrix element for the kinematical configuration
0072    * previously provided by the last call to setKinematics(), suitably
0073    * scaled by sHat() to give a dimension-less number.
0074    * @return the matrix element scaled with sHat() to give a
0075    * dimensionless number.
0076    */
0077   virtual double me2() const;
0078
0079   /**
0080    * Return the scale associated with the last set phase space point.
0081    */
0082   virtual Energy2 scale() const;
0083
0084   /**
0085    * The number of internal degrees of freedom used in the matrix
0086    * element.
0087    */
0088   virtual int nDim() const;
0089
0090   /**
0091    * Generate internal degrees of freedom given nDim() uniform
0092    * random numbers in the interval \f$ ]0,1[ \f$. To help the phase space
0093    * generator, the dSigHatDR should be a smooth function of these
0094    * numbers, although this is not strictly necessary.
0095    * @param r a pointer to the first of nDim() consecutive random numbers.
0096    * @return true if the generation succeeded, otherwise false.
0097    */
0098   virtual bool generateKinematics(const double * r);
0099
0100   /**
0101    * Return the matrix element squared differential in the variables
0102    * given by the last call to generateKinematics().
0103    */
0104   virtual CrossSection dSigHatDR() const;
0105
0106   /**
0107    * Add all possible diagrams with the add() function.
0108    */
0109   virtual void getDiagrams() const;
0110
0111   /**
0112    * Get diagram selector. With the information previously supplied with the
0113    * setKinematics method, a derived class may optionally
0114    * override this method to weight the given diagrams with their
0115    * (although certainly not physical) relative probabilities.
0116    * @param dv the diagrams to be weighted.
0117    * @return a Selector relating the given diagrams to their weights.
0118    */
0119   virtual Selector<DiagramIndex> diagrams(const DiagramVector & dv) const;
0120
0121   /**
0122    * Return a Selector with possible colour geometries for the selected
0123    * diagram weighted by their relative probabilities.
0124    * @param diag the diagram chosen.
0125    * @return the possible colour geometries weighted by their
0126    * relative probabilities.
0127    */
0128   virtual Selector<const ColourLines *>
0129   colourGeometries(tcDiagPtr diag) const;
0130   //@}
0131
0132 public:
0133
0134   /** @name Functions used by the persistent I/O system. */
0135   //@{
0136   /**
0137    * Function used to write out object persistently.
0138    * @param os the persistent output stream written to.
0139    */
0140   void persistentOutput(PersistentOStream & os) const;
0141
0142   /**
0143    * Function used to read in object persistently.
0144    * @param is the persistent input stream read from.
0145    * @param version the version number of the object when written.
0146    */
0147   void persistentInput(PersistentIStream & is, int version);
0148   //@}
0149
0150   /**
0151    * The standard Init function used to initialize the interfaces.
0152    * Called exactly once for each class by the class description system
0153    * before the main function starts or
0154    * when this class is dynamically loaded.
0155    */
0156   static void Init();
0157
0158 protected:
0159
0160   /**
0161    *  Calculate of the full next-to-leading order weight
0162    */
0163   virtual double NLOWeight() const;
0164
0165   /**
0166    *  Leading-order matrix element for \f$q\bar q\to \gamma\gamma\f$
0167    */
0168   double loGammaGammaME(const cPDVector & particles,
0169             const vector<Lorentz5Momentum> & momenta,
0170             bool first=false) const;
0171
0172   /**
0173    *  Leading-order matrix element for \f$qg\to \gamma q\f$
0174    */
0175   double loGammaqME(const cPDVector & particles,
0176         const vector<Lorentz5Momentum> & momenta,
0177         bool first=false) const;
0178
0179   /**
0180    *  Leading-order matrix element for \f$g\bar q\to \gamma \bar q\f$
0181    */
0182   double loGammaqbarME(const cPDVector & particles,
0183                const vector<Lorentz5Momentum> & momenta,
0184                bool first=false) const;
0185
0186   /**
0187    *  Leading-order matrix element for \f$q\bar q\to \gamma g\f$
0188    */
0189   double loGammagME(const cPDVector & particles,
0190             const vector<Lorentz5Momentum> & momenta,
0191             bool first=false) const;
0192
0193   /**
0194    *  Real emission matrix element for \f$q\bar q\to \gamma \gamma g\f$
0195    */
0196   InvEnergy2 realGammaGammagME(const cPDVector & particles,
0197                    const vector<Lorentz5Momentum> & momenta,
0198                    DipoleType dipole, RadiationType rad,
0199                    bool first=false) const;
0200
0201   /**
0202    *  Real emission matrix element for \f$qg\to \gamma \gamma q\f$
0203    */
0204   InvEnergy2 realGammaGammaqME(const cPDVector & particles,
0205                    const vector<Lorentz5Momentum> & momenta,
0206                    DipoleType dipole, RadiationType rad,
0207                    bool first=false) const;
0208
0209   /**
0210    *  Real emission matrix element for \f$g\bar q\to \gamma \gamma \bar q\f$
0211    */
0212   InvEnergy2 realGammaGammaqbarME(const cPDVector & particles,
0213                   const vector<Lorentz5Momentum> & momenta,
0214                   DipoleType dipole, RadiationType rad,
0215                   bool first=false) const;
0216
0217   /**
0218    *  The dipole subtractedvirtual contribution
0219    */
0220   double subtractedVirtual() const;
0221
0222   /**
0223    *  Subtracted real contribution
0224    */
0225   double subtractedReal(pair<double,double> x, double z,
0226             double zJac, double oldqPDF, double newqPDF,
0227             double newgPDF,bool order) const;
0228
0229   /**
0230    *  Calculate of the collinear counterterms
0231    */
0232   //@{
0233   /**
0234    *  Quark collinear counter term
0235    */
0236   double collinearQuark(double x, Energy2 mu2, double jac, double z,
0237             double oldPDF, double newPDF) const;
0238
0239   /**
0240    *  Gluon collinear counter term
0241    */
0242   double collinearGluon(Energy2 mu2, double jac, double z,
0243             double oldPDF, double newPDF) const;
0244   //@}
0245
0246   /**
0247    * The real matrix element divided by \f$2 g_S^2\f$, to be implemented in the
0248    * inheriting classes.
0249    * @param particles The ParticleData objects of the particles
0250    * @param momenta The momenta of the particles
0251    */
0252   double realME(const cPDVector & particles,
0253         const vector<Lorentz5Momentum> & momenta) const;
0254
0255   /**
0256    *  Generate hard QCD emission
0257    */
0258   RealEmissionProcessPtr hardQCDEmission(RealEmissionProcessPtr,
0259                      ParticleVector,
0260                      pair<double,double>);
0261
0262   /**
0263    *  Generate hard QED emission
0264    */
0265   RealEmissionProcessPtr hardQEDEmission(RealEmissionProcessPtr,
0266                      ParticleVector,
0267                      pair<double,double>);
0268
0269   /**
0270    *  The supression function
0271    */
0272   pair<double,double> supressionFunction(Energy pT,Energy scale) const {
0273     if(supressionScale_==0) scale = lambda_;
0274     Energy2 scale2 = sqr(scale), pT2 = sqr(pT);
0275     switch( supressionFunction_ ) {
0276     case 0:
0277       return make_pair(1.,0.);
0278     case 1:
0279       if(pT < scale ) return make_pair(1.,0.);
0280       else            return make_pair(0.,1.);
0281     case 2:
0282       return make_pair(scale2/(pT2+scale2),pT2/(pT2+scale2));
0283     default:
0284       assert(false);
0285       return make_pair(0.,0.);
0286     }
0287   }
0288
0289
0290 protected:
0291
0292   /** @name Clone Methods. */
0293   //@{
0294   /**
0295    * Make a simple clone of this object.
0296    * @return a pointer to the new object.
0297    */
0298   virtual IBPtr clone() const;
0299
0300   /** Make a clone of this object, possibly modifying the cloned object
0301    * to make it sane.
0302    * @return a pointer to the new object.
0303    */
0304   virtual IBPtr fullclone() const;
0305   //@}
0306 protected:
0307
0308   /** @name Standard Interfaced functions. */
0309   //@{
0310   /**
0311    * Initialize this object after the setup phase before saving an
0312    * EventGenerator to disk.
0313    * @throws InitException if object could not be initialized properly.
0314    */
0315   virtual void doinit();
0316   //@}
0317
0318 private:
0319
0320   /**
0321    * The assignment operator is private and must never be called.
0322    * In fact, it should not even be implemented.
0323    */
0324   MEPP2GammaGammaPowheg & operator=(const MEPP2GammaGammaPowheg &) = delete;
0325
0326 private:
0327
0328   /**
0329    *  Vertices
0330    */
0331   //@{
0332   /**
0333    *   FFPVertex
0334    */
0335   AbstractFFVVertexPtr FFPvertex_;
0336
0337   /**
0338    *   FFGVertex
0339    */
0340   AbstractFFVVertexPtr FFGvertex_;
0341   //@}
0342
0343   /**
0344    *  Kinematic variables for the real radiation
0345    */
0346   //@{
0347   /**
0348    *  First  variable
0349    */
0350   mutable double zTilde_;
0351
0352   /**
0353    *  Second variable
0354    */
0355   mutable double vTilde_;
0356
0357   /**
0358    *  Azimuthal angle
0359    */
0360   mutable double phi_;
0361   //@}
0362
0363   /**
0364    *  Whether to generate the positive, negative or leading order contribution
0365    */
0366   unsigned int contrib_;
0367
0368   /**
0369    *  Power for sampling \f$x_p\f$
0370    */
0371   double power_;
0372
0373   /**
0374    *  Pointer to the gluon ParticleData object
0375    */
0376   tcPDPtr gluon_;
0377
0378   /**
0379    *  Processes
0380    */
0381   unsigned int process_;
0382
0383   /**
0384    *  Processes
0385    */
0386   unsigned int threeBodyProcess_;
0387
0388   /**
0389    *  Allowed flavours of the incoming quarks
0390    */
0391   int maxflavour_;
0392
0393   /**
0394    *  Factor for \f$C_F\f$ dependent pieces
0395    */
0396   mutable double CFfact_;
0397
0398   /**
0399    *  Factor for \f$T_R\f$ dependent pieces
0400    */
0401   mutable double TRfact_;
0402
0403   /**
0404    *  Strong coupling
0405    */
0406   mutable double alphaS_;
0407
0408   /**
0409    *  Use a fixed value of \f$\alpha_S\f$
0410    */
0411   bool fixedAlphaS_;
0412
0413   /**
0414    *  Electromagnetic coupling
0415    */
0416   mutable double alphaEM_;
0417
0418   /**
0419    *  Leading-order matrix element
0420    */
0421   mutable double loME_;
0422
0423   /**
0424    *  The matrix element
0425    */
0426   mutable ProductionMatrixElement me_;
0427
0428   /**
0429    *  the selected dipole
0430    */
0431   mutable DipoleType dipole_;
0432
0433   /**
0434    *  Supression Function
0435    */
0436   //@{
0437   /**
0438    *  Choice of the supression function
0439    */
0440   unsigned int supressionFunction_;
0441
0442   /**
0443    *  Choice for the scale in the supression function
0444    */
0445   unsigned int supressionScale_;
0446
0447   /**
0448    *  Scalar for the supression function
0449    */
0450   Energy lambda_;
0451   //@}
0452
0453
0454   /**
0455    *  Hard emission stuff
0456    */
0457   //@{
0458   /**
0459    *  Whether the quark is in the + or - z direction
0460    */
0461   bool quarkplus_;
0462   //@}
0463
0464   /**
0465    *  Properties of the incoming particles
0466    */
0467   //@{
0468   /**
0469    *  Pointers to the BeamParticleData objects
0470    */
0471   vector<tcBeamPtr> beams_;
0472
0473   /**
0474    *  Pointers to the ParticleDataObjects for the partons
0475    */
0476   vector<tcPDPtr> partons_;
0477   //@}
0478
0479   /**
0480    *  Constants for the sampling. The distribution is assumed to have the
0481    *  form \f$\frac{c}{{\rm GeV}}\times\left(\frac{{\rm GeV}}{p_T}\right)^n\f$
0482    */
0483   //@{
0484   /**
0485    *  The prefactor, \f$c\f$ for the \f$q\bar{q}\f$ channel
0486    */
0487   double preQCDqqbarq_;
0488   /**
0489    *  The prefactor, \f$c\f$ for the \f$q\bar{q}\f$ channel
0490    */
0491   double preQCDqqbarqbar_;
0492
0493   /**
0494    *  The prefactor, \f$c\f$ for the \f$qg\f$ channel
0495    */
0496   double preQCDqg_;
0497
0498   /**
0499    *  The prefactor, \f$c\f$ for the \f$g\bar{q}\f$ channel
0500    */
0501   double preQCDgqbar_;
0502
0503   double preQEDqqbarq_;
0504   double preQEDqqbarqbar_;
0505   double preQEDqgq_;
0506   double preQEDgqbarqbar_;
0507
0508   /**
0509    *  The prefactors as a vector for easy use
0510    */
0511   vector<double> prefactor_;
0512   //@}
0513
0514   /**
0515    *  The transverse momentum of the jet
0516    */
0517   Energy minpT_;
0518
0519   /**
0520    *  Pointer to the object calculating the strong coupling
0521    */
0522   ShowerAlphaPtr alphaQCD_;
0523
0524   /**
0525    *  Pointer to the object calculating the EM
0526    */
0527   ShowerAlphaPtr alphaQED_;
0528
0529   /**
0530    *  Scale choice
0531    */
0532   unsigned int scaleChoice_;
0533
0534   /**
0535    *  Scale factor
0536    */
0537   double scalePreFactor_;
0538
0539 };
0540
0541 }
0542
0543 #endif /* HERWIG_MEPP2GammaGammaPowheg_H */