alpha_s.h

// -*- C++ -*-
 
// couplings/alpha_s.h is part of matchbox
// (C) 2008 Simon Platzer -- sp@particle.uni-karlsruhe.de
 
#ifndef matchbox_couplings_alpha_s_h
#define matchbox_couplings_alpha_s_h
 
#include <string>
 
#include <array>
 
#include "ThePEG/Interface/Interfaced.h"
#include "ThePEG/StandardModel/AlphaSBase.h"
 
#include "gsl.h"
 
namespace matchbox {
 
  using namespace ThePEG;
 
  template<class AlphaS>
  struct solve_lambda_below {
      
    typedef AlphaS alpha_s;
      
    inline solve_lambda_below (alpha_s* a,
                               unsigned int n,
                               Energy2 lambda2n,
                               Energy2 mass2)
      : alpha(a), nf_in(n), lambda2_nf_in(lambda2n), threshold(mass2) {}
      
    alpha_s * alpha;
    unsigned int nf_in;
    Energy2 lambda2_nf_in;
    Energy2 threshold;
      
    inline double operator () (double lambda2) {
      return ((*alpha)(threshold,lambda2_nf_in,nf_in) -
              (*alpha)(threshold,lambda2*MeV2,nf_in-1));
    }
      
  };
    
  template<class AlphaS>
  struct solve_lambda_above {
      
    typedef AlphaS alpha_s;
      
    inline solve_lambda_above (alpha_s * a,
                               unsigned int n,
                               Energy2 lambda2n,
                               Energy2 mass2)
      : alpha(a), nf_in(n), lambda2_nf_in(lambda2n), threshold(mass2) {}
      
    alpha_s * alpha;
    unsigned int nf_in;
    Energy2 lambda2_nf_in;
    Energy2 threshold;
      
    inline double operator () (double lambda2) {
      return ((*alpha)(threshold,lambda2_nf_in,nf_in) -
              (*alpha)(threshold,lambda2*MeV2,nf_in+1));
    }
      
  };
 
  template<class AlphaS>
  struct solve_input_lambda {
      
    typedef AlphaS alpha_s;
      
    inline solve_input_lambda (alpha_s * a,
                               unsigned int n,
                               double inalpha,
                               Energy2 inscale)
      : alpha(a), nf_in(n), alpha_in(inalpha), scale_in(inscale) {}
      
    alpha_s * alpha;
    unsigned int nf_in;
    double alpha_in;
    Energy2 scale_in;
      
    inline double operator () (double lambda2) {
      return ((*alpha)(scale_in,lambda2*MeV2,nf_in) - alpha_in);
    }
      
  };
 
  class alpha_s
    : public AlphaSBase {
 
  public:
 
    alpha_s();
 
  public:
 
    virtual inline double value(Energy2 scale, const StandardModelBase &) const {
      return operator() (scale);
    }
 
    virtual inline vector<Energy2> flavourThresholds() const {
      assert(!nfvector.empty());
      return nfvector;
    }
 
    virtual inline vector<Energy> LambdaQCDs() const {
      vector<Energy> res;
      for (size_t k = 0; k < 7; ++k)
        res.push_back(sqrt(lambda_squared_[k]));
      return res;
    }
 
  public:
 
    inline double operator () (Energy2 scale) const {
 
      if ( fixed_ )
        return alpha_s_in_;
 
      assert(matched());
      unsigned int active = active_flavours(scale_factor_*scale);
      return operator () (scale_factor_*scale,lambda_squared_[active],active);
 
    }
 
    virtual double operator () (Energy2 scale,
                                Energy2 lambda2,
                                unsigned int nf) const = 0;
 
    string check (string args);
 
  public:
 
    inline unsigned int min_active_flavours () const { return min_active_flavours_; }
 
    inline void min_active_flavours (unsigned int nf) { min_active_flavours_ = nf; }
 
    inline unsigned int max_active_flavours () const { return max_active_flavours_; }
 
    inline void max_active_flavours (unsigned int nf) { max_active_flavours_ = nf; }
 
    inline unsigned int active_flavours (Energy2 scale) const {
      unsigned int active = 0;
      if (scale > 0.*GeV2) {
        while(quark_mass_squared(active) < scale) {
          if (++active == max_active_flavours_+1)
            break;
        }
        active -= 1;
      } else {
        active = 0;
      }
      return active;
    }
 
    inline Energy2 lambda_squared (unsigned int f) const {
      assert(f < 7);
      return lambda_squared_[f];
    }
 
    inline Energy2 quark_mass_squared (unsigned int f) const {
      assert(f < 7);
      return quark_masses_squared_[f];
    }
 
    inline void quark_mass_squared (unsigned int f, Energy2 m2) {
      assert(f < 7);
      quark_masses_squared_[f] = m2;
      matched_ = false;
    }
 
  public:
 
    void match_thresholds ();
 
    inline bool matched () const { return matched_; }
 
  protected:
 
 
    virtual inline void doinit() {
      match_thresholds();
      copy(quark_masses_squared_.begin()+1,
           quark_masses_squared_.end(),nfvector.begin());
      AlphaSBase::doinit();
    }
 
 
    double scale_factor () const { return scale_factor_; }
 
  public:
 
    void persistentOutput(PersistentOStream & os) const;
 
    void persistentInput(PersistentIStream & is, int version);
 
    static void Init();
 
  private:
 
    static AbstractClassDescription<alpha_s> initalpha_s;
 
    alpha_s & operator=(const alpha_s &) = delete;
 
  private:
 
    unsigned int min_active_flavours_;
    unsigned int max_active_flavours_;
 
    bool matched_;
 
    double scale_factor_;
 
    std::array<Energy2,7> quark_masses_squared_;
    std::array<Energy2,7> lambda_squared_;
    vector<Energy2> nfvector=vector<Energy2>(6);
  
    double alpha_s_in_;
    Energy scale_in_;
 
    pair<Energy2,Energy2> lambda_range_;
 
    bool fixed_;
 
  };
 
}
 
#include "ThePEG/Utilities/ClassTraits.h"
 
namespace ThePEG {
 
  template <>
  struct BaseClassTrait<matchbox::alpha_s,1> {
    typedef AlphaSBase NthBase;
  };
 
  template <>
  struct ClassTraits<matchbox::alpha_s>
    : public ClassTraitsBase<matchbox::alpha_s> {
    static string className() { return "matchbox::alpha_s"; }
    static string library() { return "HwDipoleShowerAlphaS.so"; }
  };
 
}
 
#endif /* matchbox_couplings_alpha_s_h */