MatrixElement/Hadron/MEDiffraction.h

0001 // -*- C++ -*-
0002 #ifndef HERWIG_MEDiffraction_H
0003 #define HERWIG_MEDiffraction_H
0004 //
0005 // This is the declaration of the MEDiffraction class.
0006 //
0007
0008 #include "Herwig/MatrixElement/HwMEBase.h"
0009 #include "Herwig/Shower/UEBase.h"
0010
0011 namespace Herwig {
0012
0013 using namespace ThePEG;
0014
0015 /**
0016  * The MEDiffraction class provides a simple colour singlet exchange matrix element
0017  * to be used in the soft component of the multiple scattering model of the
0018  * underlying event
0019  *
0020  * @see \ref MEDiffractionInterfaces "The interfaces"
0021  * defined for MEDiffraction.
0022  */
0023 class MEDiffraction: public HwMEBase {
0024
0025 public:
0026
0027   MEDiffraction();
0028
0029   /** @name Virtual functions required by the MEBase class. */
0030   //@{
0031   /**
0032    * Return the order in \f$\alpha_S\f$ in which this matrix
0033    * element is given.
0034    */
0035   virtual unsigned int orderInAlphaS() const;
0036
0037   /**
0038    * Return the order in \f$\alpha_{EW}\f$ in which this matrix
0039    * element is given.
0040    */
0041   virtual unsigned int orderInAlphaEW() const;
0042
0043   /**
0044    * The matrix element for the kinematical configuration
0045    * previously provided by the last call to setKinematics(), suitably
0046    * scaled by sHat() to give a dimension-less number.
0047    * @return the matrix element scaled with sHat() to give a
0048    * dimensionless number.
0049    */
0050   virtual double me2() const;
0051
0052   /**
0053    * Correction weight to reweight the cross section to the diffractive
0054    * cross section.
0055    */
0056   double correctionweight() const;
0057
0058   /**
0059    * Return the scale associated with the last set phase space point.
0060    */
0061   virtual Energy2 scale() const;
0062
0063   /**
0064    * Set the typed and momenta of the incoming and outgoing partons to
0065    * be used in subsequent calls to me() and colourGeometries()
0066    * according to the associated XComb object. If the function is
0067    * overridden in a sub class the new function must call the base
0068    * class one first.
0069    */
0070   virtual void setKinematics();
0071
0072   /**
0073    * The number of internal degrees of freedom used in the matrix
0074    * element.
0075    */
0076   virtual int nDim() const;
0077
0078   /**
0079    * Generate internal degrees of freedom given nDim() uniform
0080    * random numbers in the interval \f$ ]0,1[ \f$. To help the phase space
0081    * generator, the dSigHatDR should be a smooth function of these
0082    * numbers, although this is not strictly necessary.
0083    * @param r a pointer to the first of nDim() consecutive random numbers.
0084    * @return true if the generation succeeded, otherwise false.
0085    */
0086   virtual bool generateKinematics(const double * r);
0087
0088   /**
0089    * Return the matrix element squared differential in the variables
0090    * given by the last call to generateKinematics().
0091    */
0092   virtual CrossSection dSigHatDR() const;
0093
0094   /**
0095    * Add all possible diagrams with the add() function.
0096    */
0097   virtual void getDiagrams() const;
0098
0099   /**
0100    * Get diagram selector. With the information previously supplied with the
0101    * setKinematics method, a derived class may optionally
0102    * override this method to weight the given diagrams with their
0103    * (although certainly not physical) relative probabilities.
0104    * @param dv the diagrams to be weighted.
0105    * @return a Selector relating the given diagrams to their weights.
0106    */
0107   virtual Selector<DiagramIndex> diagrams(const DiagramVector & dv) const;
0108
0109   /**
0110    * Return a Selector with possible colour geometries for the selected
0111    * diagram weighted by their relative probabilities.
0112    * @param diag the diagram chosen.
0113    * @return the possible colour geometries weighted by their
0114    * relative probabilities.
0115    */
0116   virtual Selector<const ColourLines *>
0117   colourGeometries(tcDiagPtr diag) const;
0118   //@}
0119
0120   /**
0121    * Expect the incoming partons in the laboratory frame
0122    */
0123   /* virtual bool wantCMS() const { return false; } */
0124
0125 public:
0126
0127   /** @name Functions used by the persistent I/O system. */
0128   //@{
0129   /**
0130    * Function used to write out object persistently.
0131    * @param os the persistent output stream written to.
0132    */
0133   void persistentOutput(PersistentOStream & os) const;
0134
0135   /**
0136    * Function used to read in object persistently.
0137    * @param is the persistent input stream read from.
0138    * @param version the version number of the object when written.
0139    */
0140   void persistentInput(PersistentIStream & is, int version);
0141   //@}
0142
0143   /**
0144    * The standard Init function used to initialize the interfaces.
0145    * Called exactly once for each class by the class description system
0146    * before the main function starts or
0147    * when this class is dynamically loaded.
0148    */
0149   static void Init();
0150
0151 protected:
0152
0153   /** @name Standard Interfaced functions. */
0154   //@{
0155   /**
0156    * Initialize this object after the setup phase before saving an
0157    * EventGenerator to disk.
0158    * @throws InitException if object could not be initialized properly.
0159    */
0160   virtual void doinit();
0161
0162   /**
0163    * Initialize this object. Called in the run phase just before a run begins.
0164    */
0165   virtual void doinitrun();
0166   //@}
0167
0168   /** @name Clone Methods. */
0169   //@{
0170   /**
0171    * Make a simple clone of this object.
0172    * @return a pointer to the new object.
0173    */
0174   virtual IBPtr clone() const;
0175
0176   /** Make a clone of this object, possibly modifying the cloned object
0177    * to make it sane.
0178    * @return a pointer to the new object.
0179    */
0180   virtual IBPtr fullclone() const;
0181   //@}
0182
0183 private:
0184
0185   /* The matrix element squared */
0186   double theme2;
0187
0188   /* Use only delta as excited state */
0189   bool deltaOnly;
0190
0191   /* Only produce single-diffractive events */
0192   // Default = off
0193   bool singleDiffractionOnly = 0;
0194
0195   /* Direction of the excited proton */
0196   unsigned int diffDirection;
0197
0198   /* Number of clusters the dissociated proton decays into */
0199   unsigned int dissociationDecay;
0200
0201   /* The mass of the consitutent quark */
0202   Energy mq() const {return Energy(0.325*GeV);}
0203
0204   /* The mass of the constituent diquark */
0205   Energy mqq() const {return Energy(0.650*GeV);}
0206
0207   /* The proton-pomeron slope */
0208   double theprotonPomeronSlope;
0209
0210   /* The soft pomeron intercept */
0211   double thesoftPomeronIntercept;
0212
0213   /* The soft pomeron slope */
0214   double thesoftPomeronSlope;
0215
0216
0217   /**
0218    * Sample the diffractive mass squared M2 and the momentum transfer t
0219    */
0220   pair<pair<Energy2,Energy2>,Energy2> diffractiveMassAndMomentumTransfer() const;
0221
0222
0223   /**
0224    * Random value for the diffractive mass squared M2 according to (M2/s0)^(-intercept)
0225    */
0226   Energy2 randomM2() const;
0227
0228   /**
0229    * Random value for t according to exp(diffSlope*t)
0230    */
0231   Energy2 randomt(Energy2 M2) const;
0232
0233   /**
0234    * Random value for t according to exp(diffSlope*t) for double diffraction
0235    */
0236
0237   Energy2 doublediffrandomt(Energy2 M12, Energy2 M22) const;
0238
0239
0240   /**
0241    * Returns the momenta of the two-body decay of momentum pp
0242    */
0243   pair<Lorentz5Momentum,Lorentz5Momentum> twoBodyDecayMomenta(Lorentz5Momentum pp) const;
0244
0245   /**
0246    * Returns the proton-pomeron slope
0247    */
0248   InvEnergy2 protonPomeronSlope() const;
0249
0250   /**
0251    * Returns the soft pomeron intercept
0252    */
0253   double softPomeronIntercept() const;
0254
0255   //M12 and M22 are masses squared of
0256   //outgoing particles
0257
0258   /**
0259    * Returns the minimal possible value of momentum transfer t given the center
0260    * of mass energy and diffractive masses
0261    */
0262   Energy2 tminfun(Energy2 s, Energy2 M12, Energy2 M22) const;
0263
0264   /**
0265    * Returns the maximal possible value of momentum transfer t given the center
0266    * of mass energy and diffractive masses
0267    */
0268   Energy2 tmaxfun(Energy2 s , Energy2 M12, Energy2 M22) const;
0269
0270   /**
0271    * Returns the minimal possible value of diffractive mass
0272    */
0273   //lowest possible mass given the constituent masses of quark and diquark
0274   Energy2 M2min() const{return sqr(getParticleData(2212)->mass()+mq()+mqq());}
0275
0276   /**
0277    * Returns the maximal possible value of diffractive mass
0278    */
0279   Energy2 M2max() const{
0280     return sqr(generator()->maximumCMEnergy()-getParticleData(2212)->mass());
0281   }//TODO:modify to get proper parameters
0282
0283   InvEnergy2 softPomeronSlope() const;
0284
0285
0286
0287   /* Kallen function */
0288   template<typename A, typename B, typename C>
0289   auto kallen(A a, B b, C c) const -> decltype(a*a)
0290   {
0291     return a*a + b*b + c*c - 2.0*(a*b + b*c + c*a);
0292   }
0293
0294
0295
0296   /**
0297    * The assignment operator is private and must never be called.
0298    * In fact, it should not even be implemented.
0299    */
0300   MEDiffraction & operator=(const MEDiffraction &) = delete;
0301
0302   bool isInRunPhase;
0303
0304
0305   /* The proton mass */
0306   Energy theProtonMass;
0307
0308   /**
0309    * a MPIHandler to administer the creation of several (semihard)
0310    * partonic interactions.
0311    */
0312   UEBasePtr MPIHandler_;
0313
0314 };
0315
0316 }
0317
0318 #endif /* HERWIG_MEDiffraction_H */