Herwig/Decay/PhaseSpaceChannel.h

0001 // -*- C++ -*-
0002 //
0003 // PhaseSpaceChannel.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_PhaseSpaceChannel_H
0010 #define Herwig_PhaseSpaceChannel_H
0011 //
0012 // This is the declaration of the PhaseSpaceChannel class.
0013 //
0014
0015 #include "ThePEG/Config/ThePEG.h"
0016 #include "ThePEG/PDT/ParticleData.h"
0017 #include "ThePEG/Persistency/PersistentOStream.h"
0018 #include "ThePEG/Persistency/PersistentIStream.h"
0019 #include "ThePEG/Utilities/EnumIO.h"
0020 #include "PhaseSpaceMode.fh"
0021 #include "ThePEG/MatrixElement/Tree2toNDiagram.h"
0022
0023 namespace Herwig {
0024
0025 using namespace ThePEG;
0026
0027   /** \ingroup Decay
0028    *
0029    * This class is designed to store the information needed for a given
0030    * phase-space channel for use by the multi-channel phase space decayer
0031    * and perform the generation of the phase space for that channel.
0032    *
0033    * The decay channel is specified as a series of \f$1\to2\f$ decays to either
0034    * external particles or other intermediates. For each intermediate
0035    * the Jacobian to be used can be either a Breit-Wigner(0) or a power-law
0036    * (1).
0037    *
0038    * The class is then capable of generating a phase-space point using this
0039    * channel and computing the weight of a given point for use in a multi-channel
0040    * phase space integration using the <code>PhaseSpaceMode</code> class.
0041    *
0042    * The class is designed so that the phase-space channels can either by specified
0043    * using the addIntermediate method directly or via the repository.
0044    * (In practice at the moment all the channels are constructed by the relevant decayers
0045    *  using the former method at the moment.)
0046    *
0047    * @see PhaseSpaceMode
0048    * @see DecayIntegrator
0049    *
0050    * @author  Peter Richardson
0051    */
0052 class PhaseSpaceChannel  {
0053
0054 public:
0055
0056    friend PersistentOStream;
0057    friend PersistentIStream;
0058
0059 /**
0060  *  Struct for the intermediates in the phase-space channel
0061  */
0062 struct PhaseSpaceResonance {
0063
0064   /**
0065    *   Enum for the jacobians
0066    */
0067   enum Jacobian {BreitWigner,Power,OnShell};
0068
0069    /**
0070     *   The constructor
0071     */
0072 PhaseSpaceResonance() {};
0073    /**
0074     *   The constructor
0075     */
0076 PhaseSpaceResonance(cPDPtr part) :particle(part), mass2(sqr(part->mass())), mWidth(part->mass()*part->width()),
0077                   jacobian(BreitWigner), power(0.), children(make_pair(0,0))
0078 {};
0079   /**
0080    *  The particle
0081    */
0082   cPDPtr particle;
0083
0084   /**
0085    *  Mass squared
0086    */
0087   Energy2 mass2;
0088
0089   /**
0090    *  Mass times width
0091    */
0092   Energy2 mWidth;
0093
0094   /**
0095    *   Type of jacobian
0096    */
0097   Jacobian jacobian;
0098
0099   /**
0100    *   The power for a power law
0101    */
0102   double power;
0103
0104   /**
0105    *  The children
0106    */
0107   pair<int,int> children;
0108
0109   /**
0110    *   The final descendents
0111    */
0112   vector<int> descendents;
0113
0114 };
0115
0116 public:
0117
0118   /**
0119    *  Default constructior
0120    */
0121   PhaseSpaceChannel() : weight_(1.), initialized_(false), skipFirst_(false)  {};
0122
0123   /**
0124    *  Constructor with incoming particles
0125    */
0126   PhaseSpaceChannel(tPhaseSpaceModePtr inm, bool skip=false);
0127
0128   /**
0129    * If less than zero indicate that this channel is competed. Otherwise
0130    * signal the parent of the next added parton.
0131    */
0132   PhaseSpaceChannel & operator , (tPDPtr res) {
0133     if(intermediates_.size()==1&&skipFirst_) {
0134       skipFirst_=false;
0135     }
0136     else
0137       intermediates_.push_back(PhaseSpaceResonance(res));
0138     if(iAdd_<0) return *this;
0139     if(intermediates_[iAdd_].children.first==0)
0140       intermediates_[iAdd_].children.first  = 1-int(intermediates_.size());
0141     else
0142       intermediates_[iAdd_].children.second = 1-int(intermediates_.size());
0143     iAdd_=-1;
0144     return *this;
0145   }
0146
0147   /**
0148    * If less than zero indicate that this channel is competed. Otherwise
0149    * signal the parent of the next added parton.
0150    */
0151   PhaseSpaceChannel & operator , (int o) {
0152     if(iAdd_<0) iAdd_ = o;
0153     else if(o>=0) {
0154       if(intermediates_[iAdd_].children.first==0)
0155     intermediates_[iAdd_].children.first  = o;
0156       else
0157     intermediates_[iAdd_].children.second = o;
0158       iAdd_=-1;
0159     }
0160     else if(o<0) {
0161       assert(false);
0162     }
0163     return *this;
0164   }
0165
0166   /**
0167    *  Set the jacobian for a given resonance
0168    */
0169   void setJacobian(unsigned int ires, PhaseSpaceResonance::Jacobian jac, double power) {
0170     intermediates_[ires].jacobian = jac;
0171     intermediates_[ires].power    = power;
0172   }
0173
0174 public:
0175
0176   /**
0177    *  Initialize the channel
0178    */
0179   void init(tPhaseSpaceModePtr mode);
0180
0181   /**
0182    *  Initialize the channel
0183    */
0184   void initrun(tPhaseSpaceModePtr mode);
0185
0186   /**
0187    *  Check the kinematics
0188    */
0189   bool checkKinematics();
0190
0191   /**
0192    *  The weight
0193    */
0194   const double & weight() const {return weight_;}
0195
0196   /**
0197    *  Set the weight
0198    */
0199   void  weight(double in) {weight_=in;}
0200
0201   /**
0202    * Reset the properties of an intermediate particle. This member is used
0203    * when a Decayer is used a different value for either the mass or width
0204    * of a resonace to that in the ParticleData object. This improves the
0205    * efficiency of the integration.
0206    * @param part A pointer to the particle data object for the intermediate.
0207    * @param mass The new mass of the intermediate
0208    * @param width The new width of the intermediate.
0209    */
0210   void resetIntermediate(tcPDPtr part,Energy mass,Energy width) {
0211     if(!part) return;
0212     for(PhaseSpaceResonance & res : intermediates_) {
0213       if(res.particle!=part) continue;
0214       res.mass2 = sqr(mass);
0215       res.mWidth = mass*width;
0216     }
0217   }
0218
0219   /**
0220    * Generate the momenta of the external particles. This member uses the masses
0221    * of the external particles generated by the PhaseMode class and the
0222    * intermediates for the channel to generate the momenta of the decay products.
0223    * @param pin The momenta of the decay products. This is outputed by the member.
0224    * @param massext The masses of the particles. This is to allow inclusion of
0225    * off-shell effects for the external particles.
0226    */
0227   vector<Lorentz5Momentum> generateMomenta(const Lorentz5Momentum & pin,
0228                        const vector<Energy> & massext) const;
0229
0230
0231   /**
0232    * Generate the weight for this channel given a phase space configuration.
0233    * This member generates the weight for a given phase space configuration
0234    * and is used by the PhaseSpaceMode class to compute the denominator
0235    * of the weight in the multi-channel integration.
0236    * @param output The momenta of the outgoing particles.
0237    */
0238   double generateWeight(const vector<Lorentz5Momentum> & output) const;
0239
0240   /**
0241    * Generate the final-state particles including the intermediate resonances.
0242    * This method takes the outgoing particles and adds the intermediate particles
0243    * specified by this phase-space channel. This is to allow a given set of
0244    * intermediates to be specified even when there is interference between different
0245    * intermediate states.
0246    * @param cc Whether the particles are the mode specified or its charge conjugate.
0247    * @param in The incoming particles.
0248    * @param out The outgoing particles.
0249    *
0250    */
0251   void generateIntermediates(bool cc,const Particle & in, ParticleVector & out);
0252
0253   /**
0254    *   Create a \f$2\to n\f$ diagrams for the channel
0255    */
0256   ThePEG::Ptr<ThePEG::Tree2toNDiagram>::pointer createDiagram() const;
0257
0258 public:
0259
0260   /**
0261    * Output operator to allow the structure to be persistently written
0262    * @param os The output stream
0263    * @param x The intermediate
0264    */
0265   inline friend PersistentOStream & operator<<(PersistentOStream & os,
0266                            const PhaseSpaceChannel  & x) {
0267     os << x.weight_ << x.initialized_ << x.intermediates_;
0268     return os;
0269   }
0270
0271   /**
0272    * Input operator to allow persistently written data to be read in
0273    * @param is The input stream
0274    * @param x The NBVertex
0275    */
0276   inline friend PersistentIStream & operator>>(PersistentIStream & is,
0277                            PhaseSpaceChannel & x) {
0278     is >> x.weight_ >> x.initialized_ >> x.intermediates_;
0279     return is;
0280   }
0281
0282   /**
0283    *  A friend output operator to allow the channel to be outputted for
0284    * debugging purposes
0285    */
0286   friend ostream & operator<<(ostream & os, const PhaseSpaceChannel & channel);
0287
0288 private:
0289
0290   /**
0291    *  Find the external particles which are the children of a given resonance
0292    */
0293   void findChildren(const PhaseSpaceResonance & res,
0294             vector<int> & children) {
0295     if(res.children.first>0)
0296       children.push_back(res.children.first);
0297     else
0298       findChildren(intermediates_[abs(res.children.first)],children);
0299     if(!res.particle) return;
0300     if(res.children.second>0)
0301       children.push_back(res.children.second);
0302     else
0303       findChildren(intermediates_[abs(res.children.second)],children);
0304   }
0305
0306   /**
0307    * Calculate the momenta for a two body decay
0308    * The return value indicates success or failure.
0309    * @param p The momentum of the decaying particle
0310    * @param m1 The mass of the first decay product
0311    * @param m2 The mass of the second decay product
0312    * @param p1 The momentum of the first decay product
0313    * @param p2 The momentum of the second decay product
0314    */
0315   void twoBodyDecay(const Lorentz5Momentum & p,
0316             const Energy m1, const Energy m2,
0317             Lorentz5Momentum & p1, Lorentz5Momentum & p2) const;
0318
0319   /** @name Mass Generation Members */
0320   //@{
0321   /**
0322    * Generate the mass of a resonance.
0323    * @param ires The resonance to be generated.
0324    * @param lower The lower limit on the particle's mass.
0325    * @param upper The upper limit on the particle's mass.
0326    */
0327   Energy generateMass(const PhaseSpaceResonance & res,
0328               Energy lower,Energy upper,
0329               const double & rnd) const;
0330
0331   /**
0332    * Return the weight for a given resonance.
0333    * @param ires The resonance to be generated.
0334    * @param moff The mass of the resonance.
0335    * @param lower The lower limit on the particle's mass.
0336    * @param upper The upper limit on the particle's mass.
0337    */
0338   InvEnergy2 massWeight(const PhaseSpaceResonance & res,
0339             Energy moff,Energy lower,Energy upper) const;
0340   //@}
0341
0342   /**
0343    * Helper function for the weight calculation.
0344    * @param ires The resonance to be generated.
0345    * @param limit The limit on the particle's mass.
0346    */
0347   double atanhelper(const PhaseSpaceResonance & res, Energy limit) const;
0348
0349 private:
0350
0351   /**
0352    *  Pointer to the phase-space mode
0353    */
0354   tPhaseSpaceModePtr mode_;
0355
0356   /**
0357    *  The intermediates
0358    */
0359   vector<PhaseSpaceResonance> intermediates_;
0360
0361   /**
0362    *   Integer to keep track of what we are adding
0363    */
0364   int iAdd_ = -1;
0365
0366   /**
0367    *  The weight
0368    */
0369   double weight_;
0370
0371   /**
0372    *  Whether or not its been initialized
0373    */
0374   bool initialized_;
0375
0376   /**
0377    *  Wheter or not to skiip the first resonance
0378    */
0379   bool skipFirst_;
0380
0381 };
0382
0383
0384 /**
0385  * Output operator to allow the structure to be persistently written
0386  * @param os The output stream
0387  * @param x The intermediate
0388  */
0389 inline PersistentOStream & operator<<(PersistentOStream & os,
0390                       const PhaseSpaceChannel::PhaseSpaceResonance  & x) {
0391   os << x.particle << ounit(x.mass2,GeV2) << ounit(x.mWidth,GeV2)
0392      << oenum(x.jacobian) << x.power << x.children << x.descendents;
0393   return os;
0394 }
0395
0396 /**
0397  * Input operator to allow persistently written data to be read in
0398  * @param is The input stream
0399  * @param x The NBVertex
0400  */
0401 inline PersistentIStream & operator>>(PersistentIStream & is,
0402                       PhaseSpaceChannel::PhaseSpaceResonance & x) {
0403   is >> x.particle >> iunit(x.mass2,GeV2) >> iunit(x.mWidth,GeV2)
0404      >> ienum(x.jacobian) >> x.power >> x.children >> x.descendents;
0405   return is;
0406 }
0407
0408 /**
0409  *  A friend output operator to allow the channel to be outputted for
0410  * debugging purposes
0411  */
0412 ostream & operator<<(ostream & os, const PhaseSpaceChannel & channel);
0413
0414 /**
0415  * exception for this class and those inheriting from it
0416  */
0417 class PhaseSpaceError: public Exception {};
0418
0419 }
0420
0421 #endif /* Herwig_PhaseSpaceChannel_H */