Decay/General/GeneralThreeBodyDecayer.h

0001 // -*- C++ -*-
0002 #ifndef HERWIG_GeneralThreeBodyDecayer_H
0003 #define HERWIG_GeneralThreeBodyDecayer_H
0004 //
0005 // This is the declaration of the GeneralThreeBodyDecayer class.
0006 //
0007
0008 #include "Herwig/Decay/DecayIntegrator.h"
0009 #include "Herwig/Models/General/TBDiagram.h"
0010 #include "GeneralThreeBodyDecayer.fh"
0011
0012 namespace Herwig {
0013 using namespace ThePEG;
0014
0015 /**
0016  * Here is the documentation of the GeneralThreeBodyDecayer class.
0017  *
0018  * @see \ref GeneralThreeBodyDecayerInterfaces "The interfaces"
0019  * defined for GeneralThreeBodyDecayer.
0020  */
0021 class GeneralThreeBodyDecayer: public DecayIntegrator {
0022
0023 public:
0024
0025   /** A ParticleData ptr and (possible) mass pair.*/
0026   typedef pair<tcPDPtr, Energy> PMPair;
0027
0028
0029 public:
0030
0031   /**
0032    * The default constructor.
0033    */
0034   GeneralThreeBodyDecayer() : nflow_(999), widthOpt_(1),
0035                   refTag_(), refTagCC_(), iflow_(999),
0036                   intOpt_(0), relerr_(1e-2)
0037   {}
0038
0039   /** @name Virtual functions required by the Decayer class. */
0040   //@{
0041   /**
0042    * For a given decay mode and a given particle instance, perform the
0043    * decay and return the decay products. As this is the base class this
0044    * is not implemented.
0045    * @return The vector of particles produced in the decay.
0046    */
0047   virtual ParticleVector decay(const Particle & parent,
0048                    const tPDVector & children) const;
0049
0050   /**
0051    * Which of the possible decays is required
0052    * @param cc Is this mode the charge conjugate
0053    * @param parent The decaying particle
0054    * @param children The decay products
0055    */
0056   virtual int modeNumber(bool & cc, tcPDPtr parent,const tPDVector & children) const;
0057
0058   /**
0059    * The matrix element to be integrated for the three-body decays as a function
0060    * of the invariant masses of pairs of the outgoing particles.
0061    * @param imode The mode for which the matrix element is needed.
0062    * @param q2 The scale, \e i.e. the mass squared of the decaying particle.
0063    * @param s3 The invariant mass squared of particles 1 and 2, \f$s_3=m^2_{12}\f$.
0064    * @param s2 The invariant mass squared of particles 1 and 3, \f$s_2=m^2_{13}\f$.
0065    * @param s1 The invariant mass squared of particles 2 and 3, \f$s_1=m^2_{23}\f$.
0066    * @param m1 The mass of the first  outgoing particle.
0067    * @param m2 The mass of the second outgoing particle.
0068    * @param m3 The mass of the third  outgoing particle.
0069    * @return The matrix element
0070    */
0071   virtual double threeBodyMatrixElement(const int imode,  const Energy2 q2,
0072                     const Energy2 s3, const Energy2 s2,
0073                     const Energy2 s1, const Energy  m1,
0074                     const Energy  m2, const Energy  m3) const;
0075
0076   /**
0077    * Function to return partial Width
0078    * @param inpart The decaying particle.
0079    * @param outa First  decay product.
0080    * @param outb Second decay product.
0081    * @param outc Third  decay product.
0082    */
0083   virtual Energy partialWidth(PMPair inpart, PMPair outa,
0084                   PMPair outb, PMPair outc) const;
0085
0086   /**
0087    * An overidden member to calculate a branching ratio for a certain
0088    * particle instance.
0089    * @param dm The DecayMode of the particle
0090    * @param p The particle object
0091    * @param oldbrat The branching fraction given in the DecayMode object
0092    */
0093   virtual double brat(const DecayMode & dm, const Particle & p,
0094               double oldbrat) const;
0095
0096   //@}
0097
0098   /**
0099    *  Set the diagrams
0100    */
0101   bool setDecayInfo(PDPtr incoming,vector<PDPtr> outgoing,
0102             const vector<TBDiagram> & process,
0103             double symfac);
0104
0105 public:
0106
0107   /** @name Functions used by the persistent I/O system. */
0108   //@{
0109   /**
0110    * Function used to write out object persistently.
0111    * @param os the persistent output stream written to.
0112    */
0113   void persistentOutput(PersistentOStream & os) const;
0114
0115   /**
0116    * Function used to read in object persistently.
0117    * @param is the persistent input stream read from.
0118    * @param version the version number of the object when written.
0119    */
0120   void persistentInput(PersistentIStream & is, int version);
0121   //@}
0122
0123   /**
0124    * The standard Init function used to initialize the interfaces.
0125    * Called exactly once for each class by the class description system
0126    * before the main function starts or
0127    * when this class is dynamically loaded.
0128    */
0129   static void Init();
0130
0131 protected:
0132
0133   /** @name Standard Interfaced functions. */
0134   //@{
0135   /**
0136    * Initialize this object after the setup phase before saving an
0137    * EventGenerator to disk.
0138    * @throws InitException if object could not be initialized properly.
0139    */
0140   virtual void doinit();
0141
0142   /**
0143    * Initialize this object. Called in the run phase just before
0144    * a run begins.
0145    */
0146   virtual void doinitrun();
0147   //@}
0148
0149   /**
0150    *   Set up the diagrams etc
0151    */
0152   virtual void setupDiagrams(bool checkKinematics);
0153
0154 protected:
0155
0156   /**
0157    * Access the TBDiagrams that store the required information
0158    * to create the diagrams
0159    */
0160   const vector<TBDiagram> & getProcessInfo() const {
0161     return diagrams_;
0162   }
0163
0164   /**
0165    *  Incoming particle
0166    */
0167   PDPtr incoming() const { return incoming_; }
0168
0169   /**
0170    *  Outgoing particles
0171    */
0172   const vector<tPDPtr> & outgoing() const {  return outgoing_; }
0173
0174   /**
0175    *  Number of colour flows
0176    */
0177   unsigned int numberOfFlows() const { return nflow_; }
0178
0179   /**
0180    * Set up the colour factors
0181    */
0182   bool setColourFactors(double symfac);
0183
0184   /**
0185    * Return the matrix of colour factors
0186    */
0187   const vector<DVector> & getColourFactors() const {  return colour_; }
0188
0189   /**
0190    * Return the matrix of colour factors
0191    */
0192   const vector<DVector> & getLargeNcColourFactors() const {
0193     return colourLargeNC_;
0194   }
0195
0196   /**
0197    *  Get the mapping between the phase-space channel and the diagram
0198    */
0199   const vector<unsigned int> & diagramMap() const {
0200     return diagmap_;
0201   }
0202
0203   /**
0204    *  Option for the handling of the widths of the intermediate particles
0205    */
0206   unsigned int widthOption() const { return widthOpt_; }
0207
0208   /**
0209    * Set colour connections
0210    * @param parent Parent particle
0211    * @param out Particle vector containing particles to
0212    * connect colour lines
0213    */
0214   void colourConnections(const Particle & parent,
0215              const ParticleVector & out) const;
0216
0217   /**
0218    *  Method to construct the channels for the integrator to give the partial width
0219    * @param intype  Types of the channels
0220    * @param inmass Mass for the channels
0221    * @param inwidth Width for the channels
0222    * @param inpow Power for the channels
0223    * @param inweights Weights for the channels
0224    */
0225   void constructIntegratorChannels(vector<int> & intype, vector<Energy> & inmass,
0226                    vector<Energy> & inwidth, vector<double> & inpow,
0227                    vector<double> & inweights) const;
0228
0229   /**
0230    *  Set the colour flow
0231    * @param flow The value for the colour flow
0232    */
0233   void colourFlow(unsigned int flow) const { iflow_ = flow; }
0234
0235   /**
0236    *  Set the colour flow
0237    */
0238   unsigned int const & colourFlow() const { return iflow_; }
0239
0240   /**
0241    *  Relative error for GQ integration
0242    */
0243   double relativeError() const {return relerr_;}
0244
0245 private:
0246
0247   /**
0248    * The assignment operator is private and must never be called.
0249    * In fact, it should not even be implemented.
0250    */
0251   GeneralThreeBodyDecayer & operator=(const GeneralThreeBodyDecayer &) = delete;
0252
0253 private:
0254
0255   /**
0256    *  Store the incoming particle
0257    */
0258   PDPtr incoming_;
0259
0260   /**
0261    *  Outgoing particles
0262    */
0263   vector<tPDPtr> outgoing_;
0264
0265   /**
0266    *  Store the diagrams for the decay
0267    */
0268   vector<TBDiagram> diagrams_;
0269
0270   /**
0271    *  Map between the diagrams and the phase-space channels
0272    */
0273   vector<unsigned int> diagmap_;
0274
0275   /**
0276    * Store colour factors for ME calc.
0277    */
0278   vector<DVector> colour_;
0279
0280   /**
0281    *  Store cololur factors for ME calc at large N_c
0282    */
0283   vector<DVector> colourLargeNC_;
0284
0285   /**
0286    * The number of colourflows.
0287    */
0288   unsigned int nflow_;
0289
0290   /**
0291    *  Reference to object to calculate the partial width
0292    */
0293   mutable WidthCalculatorBasePtr widthCalc_;
0294
0295   /**
0296    *  Option for the treatment of the widths
0297    */
0298   unsigned int widthOpt_;
0299
0300   /**
0301    * Store a decay tag for this mode that can be tested when
0302    * trying to determine whether it can be generated by
0303    * this Decayer
0304    */
0305   string refTag_;
0306
0307   /**
0308    * Store a decay tag for the cc-mode that can be tested when
0309    * trying to determine whether it can be generated by
0310    * this Decayer
0311    */
0312   string refTagCC_;
0313
0314   /**
0315    *  The colour flow
0316    */
0317   mutable unsigned int iflow_;
0318
0319   /**
0320    *  Option for the construction of the gaussian integrator
0321    */
0322   unsigned int intOpt_;
0323
0324   /**
0325    *  Relative error for GQ integration of partial width
0326    */
0327   double relerr_;
0328 };
0329
0330 }
0331
0332 #endif /* HERWIG_GeneralThreeBodyDecayer_H */