Dipole/AlphaS/alpha_s.h

0001 // -*- C++ -*-
0002
0003 // couplings/alpha_s.h is part of matchbox
0004 // (C) 2008 Simon Platzer -- sp@particle.uni-karlsruhe.de
0005
0006 #ifndef matchbox_couplings_alpha_s_h
0007 #define matchbox_couplings_alpha_s_h
0008
0009 #include <string>
0010
0011 #include <array>
0012
0013 #include "ThePEG/Interface/Interfaced.h"
0014 #include "ThePEG/StandardModel/AlphaSBase.h"
0015
0016 #include "gsl.h"
0017
0018 namespace matchbox {
0019
0020   using namespace ThePEG;
0021
0022   template<class AlphaS>
0023   struct solve_lambda_below {
0024
0025     typedef AlphaS alpha_s;
0026
0027     inline solve_lambda_below (alpha_s* a,
0028                    unsigned int n,
0029                    Energy2 lambda2n,
0030                    Energy2 mass2)
0031       : alpha(a), nf_in(n), lambda2_nf_in(lambda2n), threshold(mass2) {}
0032
0033     alpha_s * alpha;
0034     unsigned int nf_in;
0035     Energy2 lambda2_nf_in;
0036     Energy2 threshold;
0037
0038     inline double operator () (double lambda2) {
0039       return ((*alpha)(threshold,lambda2_nf_in,nf_in) -
0040           (*alpha)(threshold,lambda2*MeV2,nf_in-1));
0041     }
0042
0043   };
0044
0045   template<class AlphaS>
0046   struct solve_lambda_above {
0047
0048     typedef AlphaS alpha_s;
0049
0050     inline solve_lambda_above (alpha_s * a,
0051                    unsigned int n,
0052                    Energy2 lambda2n,
0053                    Energy2 mass2)
0054       : alpha(a), nf_in(n), lambda2_nf_in(lambda2n), threshold(mass2) {}
0055
0056     alpha_s * alpha;
0057     unsigned int nf_in;
0058     Energy2 lambda2_nf_in;
0059     Energy2 threshold;
0060
0061     inline double operator () (double lambda2) {
0062       return ((*alpha)(threshold,lambda2_nf_in,nf_in) -
0063           (*alpha)(threshold,lambda2*MeV2,nf_in+1));
0064     }
0065
0066   };
0067
0068   template<class AlphaS>
0069   struct solve_input_lambda {
0070
0071     typedef AlphaS alpha_s;
0072
0073     inline solve_input_lambda (alpha_s * a,
0074                    unsigned int n,
0075                    double inalpha,
0076                    Energy2 inscale)
0077       : alpha(a), nf_in(n), alpha_in(inalpha), scale_in(inscale) {}
0078
0079     alpha_s * alpha;
0080     unsigned int nf_in;
0081     double alpha_in;
0082     Energy2 scale_in;
0083
0084     inline double operator () (double lambda2) {
0085       return ((*alpha)(scale_in,lambda2*MeV2,nf_in) - alpha_in);
0086     }
0087
0088   };
0089
0090   /**
0091    * Base class for the strong coupling.
0092    *
0093    * @see \ref alpha_sInterfaces "The interfaces"
0094    * defined for alpha_s.
0095    */
0096   class alpha_s
0097     : public AlphaSBase {
0098
0099   public:
0100
0101     /**
0102      * The default constructor.
0103      */
0104     alpha_s();
0105
0106   public:
0107
0108     /** @name Virtual functions as required by AlphaSBase. */
0109     //@{
0110     /**
0111      * The \f$\alpha_S\f$. Return the QCD coupling for a given \a scale
0112      * using the given standard model object \a sm.
0113      */
0114     virtual inline double value(Energy2 scale, const StandardModelBase &) const {
0115       return operator() (scale);
0116     }
0117
0118     /**
0119      * Return the flavour thresholds used. The returned vector contains
0120      * (in position <code>i</code>) the scales when the active number of
0121      * flavours changes from <code>i</code> to <code>i+1</code>.
0122      */
0123     virtual inline vector<Energy2> flavourThresholds() const {
0124       assert(!nfvector.empty());
0125       return nfvector;
0126     }
0127
0128     /**
0129      * Return the \f$\Lambda_{QCD}\f$ used for different numbers of
0130      * active flavours.
0131      */
0132     virtual inline vector<Energy> LambdaQCDs() const {
0133       vector<Energy> res;
0134       for (size_t k = 0; k < 7; ++k)
0135     res.push_back(sqrt(lambda_squared_[k]));
0136       return res;
0137     }
0138     //@}
0139
0140   public:
0141
0142     /// return alpha_s as function of scale
0143     inline double operator () (Energy2 scale) const {
0144
0145       if ( fixed_ )
0146     return alpha_s_in_;
0147
0148       assert(matched());
0149       unsigned int active = active_flavours(scale_factor_*scale);
0150       return operator () (scale_factor_*scale,lambda_squared_[active],active);
0151
0152     }
0153
0154     /// return alpha_s as function of scale, QCD scale
0155     /// and number of active flavours
0156     virtual double operator () (Energy2 scale,
0157                 Energy2 lambda2,
0158                 unsigned int nf) const = 0;
0159
0160     /// match thresholds and write alpha_s
0161     /// to specified file; arguments are
0162     /// Q_low/GeV Q_high/GeV n_steps filename
0163     string check (string args);
0164
0165   public:
0166
0167     /// return minimum number of active flavours
0168     inline unsigned int min_active_flavours () const { return min_active_flavours_; }
0169
0170     /// set minimum number of active flavours
0171     inline void min_active_flavours (unsigned int nf) { min_active_flavours_ = nf; }
0172
0173     /// return maximum number of active flavours
0174     inline unsigned int max_active_flavours () const { return max_active_flavours_; }
0175
0176     /// set maximum number of active flavours
0177     inline void max_active_flavours (unsigned int nf) { max_active_flavours_ = nf; }
0178
0179     /// return the number of active flavours at the given scale
0180     inline unsigned int active_flavours (Energy2 scale) const {
0181       unsigned int active = 0;
0182       if (scale > 0.*GeV2) {
0183     while(quark_mass_squared(active) < scale) {
0184       if (++active == max_active_flavours_+1)
0185         break;
0186     }
0187     active -= 1;
0188       } else {
0189     active = 0;
0190       }
0191       return active;
0192     }
0193
0194     /// return the lambda squared for the given number of flavours
0195     inline Energy2 lambda_squared (unsigned int f) const {
0196       assert(f < 7);
0197       return lambda_squared_[f];
0198     }
0199
0200     /// return the mass squared for given flavour
0201     inline Energy2 quark_mass_squared (unsigned int f) const {
0202       assert(f < 7);
0203       return quark_masses_squared_[f];
0204     }
0205
0206     /// set the mass squared for given flavour
0207     inline void quark_mass_squared (unsigned int f, Energy2 m2) {
0208       assert(f < 7);
0209       quark_masses_squared_[f] = m2;
0210       matched_ = false;
0211     }
0212
0213   public:
0214
0215     /// perform the threshold matching
0216     /// given alpha_s value at reference scale
0217     void match_thresholds ();
0218
0219     /// return true, if threshold matching has been
0220     /// performed
0221     inline bool matched () const { return matched_; }
0222
0223   protected:
0224
0225     /** @name Standard Interfaced functions. */
0226     //@{
0227
0228     /**
0229      * Initialize this object after the setup phase before saving an
0230      * EventGenerator to disk.
0231      * @throws InitException if object could not be initialized properly.
0232      */
0233     virtual inline void doinit() {
0234       match_thresholds();
0235       copy(quark_masses_squared_.begin()+1,
0236            quark_masses_squared_.end(),nfvector.begin());
0237       AlphaSBase::doinit();
0238     }
0239
0240     //@}
0241
0242     /// return the scale factor
0243     double scale_factor () const { return scale_factor_; }
0244
0245   public:
0246
0247     /** @name Functions used by the persistent I/O system. */
0248     //@{
0249     /**
0250      * Function used to write out object persistently.
0251      * @name os the persistent output stream written to.
0252      */
0253     void persistentOutput(PersistentOStream & os) const;
0254
0255     /**
0256      * Function used to read in object persistently.
0257      * @name is the persistent input stream read from.
0258      * @name version the version number of the object when written.
0259      */
0260     void persistentInput(PersistentIStream & is, int version);
0261     //@}
0262
0263     /**
0264      * The standard Init function used to initialize the interfaces.
0265      * Called exactly once for each class by the class description system
0266      * before the main function starts or
0267      * when this class is dynamically loaded.
0268      */
0269     static void Init();
0270
0271   private:
0272
0273     /**
0274      * The static object used to initialize the description of this class.
0275      * Indicates that this is an abstract class with persistent data.
0276      */
0277     static AbstractClassDescription<alpha_s> initalpha_s;
0278
0279     /**
0280      * The assignment operator is private and must never be called.
0281      * In fact, it should not even be implemented.
0282      */
0283     alpha_s & operator=(const alpha_s &) = delete;
0284
0285   private:
0286
0287     unsigned int min_active_flavours_;
0288     unsigned int max_active_flavours_;
0289
0290     bool matched_;
0291
0292     double scale_factor_;
0293
0294     std::array<Energy2,7> quark_masses_squared_;
0295     std::array<Energy2,7> lambda_squared_;
0296     vector<Energy2> nfvector=vector<Energy2>(6);
0297
0298     double alpha_s_in_;
0299     Energy scale_in_;
0300
0301     pair<Energy2,Energy2> lambda_range_;
0302
0303     bool fixed_;
0304
0305   };
0306
0307 }
0308
0309 #include "ThePEG/Utilities/ClassTraits.h"
0310
0311 namespace ThePEG {
0312
0313   /** @cond TRAITSPECIALIZATIONS */
0314
0315   /** This template specialization informs ThePEG about the
0316    *  base classes of alpha_s. */
0317   template <>
0318   struct BaseClassTrait<matchbox::alpha_s,1> {
0319     /** Typedef of the first base class of alpha_s. */
0320     typedef AlphaSBase NthBase;
0321   };
0322
0323   /** This template specialization informs ThePEG about the name of
0324    *  the alpha_s class and the shared object where it is defined. */
0325   template <>
0326   struct ClassTraits<matchbox::alpha_s>
0327     : public ClassTraitsBase<matchbox::alpha_s> {
0328     /** Return a platform-independent class name */
0329     static string className() { return "matchbox::alpha_s"; }
0330     /**
0331      * The name of a file containing the dynamic library where the class
0332      * alpha_s is implemented. It may also include several, space-separated,
0333      * libraries if the class alpha_s depends on other classes (base classes
0334      * excepted). In this case the listed libraries will be dynamically
0335      * linked in the order they are specified.
0336      */
0337     static string library() { return "HwDipoleShowerAlphaS.so"; }
0338   };
0339
0340   /** @endcond */
0341
0342 }
0343
0344 #endif /* matchbox_couplings_alpha_s_h */