BtoSGammaKagan.h

// -*- C++ -*-
//
// BtoSGammaKagan.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_BtoSGammaKagan_H
#define HERWIG_BtoSGammaKagan_H
//
// This is the declaration of the BtoSGammaKagan class.
//
 
#include "Herwig/Utilities/Interpolator.h"
#include "BtoSGammaHadronicMass.h"
#include "ThePEG/Config/Complex.h"
#include "ThePEG/Config/Constants.h"
 
namespace Herwig {
 
using namespace ThePEG;
using namespace ThePEG::Constants;
 
class BtoSGammaKagan: public BtoSGammaHadronicMass {
 
  friend struct KaganIntegrand;
 
public:
 
  BtoSGammaKagan();
 
  virtual Energy hadronicMass(Energy mb,Energy mquark);
 
  virtual void dataBaseOutput(ofstream & os,bool header,bool create) const;
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
public:
 
  double operator ()(double x) const {
    if(_iopt==0) {
      double reg(realG(x/_zratio)),img(imagG(x/_zratio));
      return 16./27.*(1.-x)*(_zratio*_zratio/x/x*(sqr(reg)+sqr(img))+_zratio/x*reg+0.25);
    }
    else {
      return -8./9.*_zratio*(realG(x/_zratio)+0.5*x/_zratio);
    }
  }
  typedef double ValType;
  typedef double ArgType;
 
  InvEnergy smeared(Energy kp) const {
    InvEnergy fermi = exponentialFermiFunction(kp,_fermilambda,_fermia,
                                               _ferminorm,_fermilambda1);
    if(_iopt==1) fermi *=KNLO(_MB*_y/(_mb+kp))*_MB/(_mb+kp);
    return fermi;
  }
 
protected:
 
  virtual IBPtr clone() const {return new_ptr(*this);}
 
  virtual IBPtr fullclone() const {return new_ptr(*this);}
 
protected:
 
  virtual void doinit();
 
  virtual void doinitrun();
 
private:
 
  BtoSGammaKagan & operator=(const BtoSGammaKagan &) = delete;
 
private:
 
  double Delta(double y, double alphaS) const {
    if(y>_ycut) return 0.;
    double ln(log(1.-y));
    return -4./3./pi/(1.-y)*alphaS*(ln+1.75)*exp(-2.*alphaS/3./pi*ln*(ln+3.5));
  }
 
  double semiLeptonicf() const  {
    double z2=sqr(_zratio);
    return 1.-8.*_zratio*(1.-z2)-sqr(z2)-12.*z2*log(_zratio);
  }
 
  double s22(double y) const {return (*_s22inter)(y);}
 
  double s27(double y) const {return (*_s27inter)(y);}
 
  double s77(double y) const  {
    if(y>_ycut) y=_ycut; 
    return 1./3.*(7.+y*(1.-2.*y)-2.*(1.+y)*log(1.-y));
  }
 
  double s78(double y) const  {
    if(y>_ycut) y=_ycut;
    return 8./9.*((1.-y)/y*log(1.-y)+1.+0.25*y*y);
  }
 
  double s88(double y) const  {
    double ratio(_mb/_ms),y2(sqr(y));
    if(y>_ycut) y=_ycut;
    return 1./27.*(2.*(2.-2.*y+y2)/y*(log(1.-y)+2.*log(ratio))-2.*y2-y-8.*(1.-y)/y);
  }
 
  double realG(double t) const  {
    if(t<4.) {
      double at(atan(sqrt(t/(4.-t))));
      return -2.*sqr(at);
    }
    else {
      double ln(log(0.5*(sqrt(t)+sqrt(t-4.))));
      return 2.*(sqr(ln)-0.25*sqr(pi));
    }
  }
 
  double imagG(double t) const {
    if(t<4.) return 0.;
    else     return -2.*pi*log(0.5*(sqrt(t)+sqrt(t-4.)));
  }
 
  double alphaS(Energy Q)  {
    double lo(1.-0.5*_beta0*_alphaSZ/pi*log(_mz/Q));
    return _alphaSZ/lo*(1.-0.25*_beta1/_beta0*_alphaSZ/pi*log(lo)/lo);
  }
 
  void calculateWilsonCoefficients();
 
  double KNLO(double y) const {
    return _delta*Delta(y,_alphaSM)
      +_alphaSM/pi*(s22(y)*sqr(_c20)+s77(y)*sqr(_c70)
                    +s88(y)*sqr(_c80)+s78(y)*_c70*_c80
                    +s27(y)*_c20*(_c70-_c80/3.));
  }
 
private:
 
 
  bool _initialize;
 
  Energy _mt;
 
  Energy _mb;
 
  Energy _mc;
 
  Energy _ms;
 
  double _msovermb;
 
  double _zratio;
 
  Energy2 _lambda2;
 
  Energy _mw;
 
  Energy _mz;
 
  Energy _MB;
 
  double _c20;
 
  double _c70;
 
  double _c80;
 
  double _beta0;
 
  double _beta1;
 
  double _alpha;
 
  double _alphaSZ;
 
  Energy _mub;
 
  double _alphaSM;
 
  double _ckm;
  double _delta;
 
  Interpolator<double,double>::Ptr _s22inter;
 
  Interpolator<double,double>::Ptr _s27inter;
 
  Interpolator<InvEnergy,Energy>::Ptr _pmHinter;
 
  vector<Energy> _mHinter;
 
  vector<InvEnergy> _spectrum;
 
  InvEnergy _spectmax;
 
  unsigned int _maxtry;
 
  Energy _fermilambda;
 
  double _fermia;
 
  InvEnergy _ferminorm;
 
  Energy2 _fermilambda1;
 
  double _ycut;
 
  double _y;
 
  double _deltacut;
 
  unsigned int _nsfunct;
 
  unsigned int _nspect;
 
  unsigned int _iopt;
 
};
 
struct KaganIntegrand {
 
  KaganIntegrand(Ptr<BtoSGammaKagan>::pointer in) : _kagan(in) {};
 
  InvEnergy operator ()(Energy arg) const {return _kagan->smeared(arg);}
  typedef InvEnergy ValType;
  typedef Energy    ArgType;
 
  Ptr<BtoSGammaKagan>::pointer _kagan;
};
}
 
#endif /* HERWIG_BtoSGammaKagan_H */