Commit 53af62b2 authored by Enrico Bothmann's avatar Enrico Bothmann

Refactor AlphaS handling in reweighting and CF_QCD

parent 995f515e
......@@ -420,7 +420,7 @@ std::vector<Variations::PDFs_And_AlphaS> Variations::PDFsAndAlphaSVector(
Variations::PDFs_And_AlphaS::PDFs_And_AlphaS():
p_alphas(MODEL::as->GetAs(PDF::isr::hard_process))
p_alphas(MODEL::as)
{
// Workaround für C++03 (vector constructor is confused when fed
// with NULL as the initial value for a pointer)
......@@ -456,9 +456,9 @@ Variations::PDFs_And_AlphaS::PDFs_And_AlphaS(std::string pdfname, size_t pdfmemb
// obtain AlphaS based on a loaded PDF or a new one (if none is found)
if (aspdf == NULL) {
p_alphas = new MODEL::One_Running_AlphaS(pdfname, pdfmember);
p_alphas = new MODEL::Running_AlphaS(pdfname, pdfmember);
} else {
p_alphas = new MODEL::One_Running_AlphaS(aspdf);
p_alphas = new MODEL::Running_AlphaS(aspdf);
}
if (p_alphas == NULL) {
THROW(fatal_error, "AlphaS for " + pdfname + " could not be initialised.");
......@@ -576,8 +576,8 @@ std::string Variation_Parameters::GenerateName() const
pdfid = p_pdf1->LHEFNumber();
} else if (p_pdf2 != NULL) {
pdfid = p_pdf2->LHEFNumber();
} else if (p_alphas->PDF() != NULL) {
pdfid = p_alphas->PDF()->LHEFNumber();
} else if (p_alphas->GetAs()->PDF() != NULL) {
pdfid = p_alphas->GetAs()->PDF()->LHEFNumber();
} else {
THROW(fatal_error, "Cannot obtain PDF IDF");
}
......
......@@ -16,7 +16,7 @@
#define ENABLE_REWEIGHTING_FACTORS_HISTOGRAMS 0
namespace PDF { class PDF_Base; }
namespace MODEL { class One_Running_AlphaS; }
namespace MODEL { class Running_AlphaS; }
namespace ATOOLS {
......@@ -70,7 +70,7 @@ namespace ATOOLS {
PDFs_And_AlphaS();
PDFs_And_AlphaS(std::string pdfname, size_t pdfmember);
std::vector<PDF::PDF_Base *> m_pdfs;
MODEL::One_Running_AlphaS *p_alphas;
MODEL::Running_AlphaS *p_alphas;
};
void AddParameterExpandingScaleFactors(std::vector<std::string> scalestringparams,
ScaleFactorExpansions::code,
......@@ -121,7 +121,7 @@ namespace ATOOLS {
double showermuR2fac, double showermuF2fac,
PDF::PDF_Base *pdf1,
PDF::PDF_Base *pdf2,
MODEL::One_Running_AlphaS *alphas,
MODEL::Running_AlphaS *alphas,
bool deletepdfsandalphas):
m_muR2fac(muR2fac), m_muF2fac(muF2fac),
m_showermuR2fac(showermuR2fac), m_showermuF2fac(showermuF2fac),
......@@ -142,7 +142,7 @@ namespace ATOOLS {
//! Pointers to the beam PDFs, which can be NULL for (semi-)leptonic events
PDF::PDF_Base * const p_pdf1;
PDF::PDF_Base * const p_pdf2;
MODEL::One_Running_AlphaS * const p_alphas;
MODEL::Running_AlphaS * const p_alphas;
//! Set whether the pointers to the PDFs and the AlphaS should be deleted after usage
const bool m_deletepdfsandalphas;
const std::string m_name;
......
......@@ -7,6 +7,7 @@
#include "ATOOLS/Org/Exception.H"
#include <algorithm>
#include <cassert>
using namespace ATOOLS;
......@@ -20,25 +21,52 @@ namespace CSSHOWER {
class CF_QCD: public SF_Coupling {
protected:
//! Main underlying coupling (set by SetCoupling)
MODEL::Running_AlphaS *p_cpl;
class QCD_Coupling_Info {
//! (Temporary) alternative coupling (set by SetAlternativeUnderlyingCoupling)
MODEL::One_Running_AlphaS *p_altcpl;
public:
//! Alternative renormalisation scale factor (set by SetAlternativeUnderlyingCoupling)
double m_altrsf;
QCD_Coupling_Info():
p_cpl(NULL), m_rsf(1.0), p_cplmax(NULL)
{};
QCD_Coupling_Info(MODEL::Running_AlphaS * cpl, double rsf, std::vector<double> * cplmax=NULL):
p_cpl(cpl), m_rsf(rsf), p_cplmax(cplmax)
{};
MODEL::Running_AlphaS * const Coupling() const { return p_cpl; }
double RSF() const { return m_rsf; }
std::vector<double> * const MaxCoupling() const { return p_cplmax; }
void SetMaxCoupling(std::vector<double> * const cplmax) { p_cplmax = cplmax; }
bool IsValid() const { return (p_cpl != NULL); }
private:
MODEL::Running_AlphaS * p_cpl;
double m_rsf;
std::vector<double> * p_cplmax;
};
/*!
* Underlying couplings set by SetCoupling and SetAlternativeUnderlyingCoupling, respectively
*
* If the alternative coupling is set, it takes precedence. This can be used for reweighting purposes.
*/
QCD_Coupling_Info m_maincplinfo, m_altcplinfo;
const QCD_Coupling_Info & CurrentCouplingInfo() const
{ return (m_altcplinfo.IsValid() ? m_altcplinfo : m_maincplinfo); }
//! Buffer of max alphas values to avoid re-calculations
std::map<MODEL::One_Running_AlphaS *, double> m_altcplmax;
std::map<MODEL::Running_AlphaS *, std::vector<double> > m_altcplmax;
double m_q, m_rsf, m_k0sq;
double m_q, m_k0sq;
int m_scvmode;
public:
inline CF_QCD(const SF_Key &key):
SF_Coupling(key), p_altcpl(NULL), m_altrsf(1.0), m_altcplmax(), m_q(0.), m_rsf(1.), m_k0sq(0.0), m_scvmode(0)
CF_QCD(const SF_Key &key):
SF_Coupling(key),
m_maincplinfo(QCD_Coupling_Info()), m_altcplinfo(QCD_Coupling_Info()),
m_altcplmax(),
m_q(0.), m_k0sq(0.0), m_scvmode(0)
{
if (key.p_v->in[0].StrongCharge()==8 &&
key.p_v->in[1].StrongCharge()==8 &&
......@@ -65,14 +93,13 @@ namespace CSSHOWER {
bool SetCoupling(MODEL::Model_Base *md,
const double &k0sqi,const double &k0sqf,
const double &isfac,const double &fsfac);
template<class T>
double CplMax(T * as, double rsf) const;
double CplMax(MODEL::Running_AlphaS * as, double rsf) const;
double Coupling(const double &scale,const int pol);
bool AllowSpec(const ATOOLS::Flavour &fl);
double CplFac(const double &scale) const;
bool AllowsAlternativeCouplingUsage() const { return true; }
void SetAlternativeUnderlyingCoupling(void *, double sf=1.0);
void SetAlternativeUnderlyingCoupling(MODEL::Running_AlphaS *, double sf=1.0);
};
}
......@@ -84,35 +111,49 @@ bool CF_QCD::SetCoupling(MODEL::Model_Base *md,
const double &k0sqi,const double &k0sqf,
const double &isfac,const double &fsfac)
{
p_cpl=(MODEL::Running_AlphaS*)md->GetScalarFunction("alpha_S");
SetAlternativeUnderlyingCoupling(NULL);
m_altcplmax.clear(); // buffered values are not valid anymore
m_rsf=ToType<double>(rpa->gen.Variable("RENORMALIZATION_SCALE_FACTOR"))
*ToType<double>(rpa->gen.Variable("CSS_SCALE_FACTOR"));
MODEL::Running_AlphaS * cpl = (MODEL::Running_AlphaS *)(md->GetScalarFunction("alpha_S"));
double rsf = ToType<double>(rpa->gen.Variable("RENORMALIZATION_SCALE_FACTOR"));
rsf *= ToType<double>(rpa->gen.Variable("CSS_SCALE_FACTOR"));
m_scvmode=ToType<int>(rpa->gen.Variable("CSS_SCALE_VARIATION_SCHEME"));
// determine prefactors before calling CplMax below
m_cplfac=((m_type/10==1)?fsfac:isfac);
m_k0sq=(m_type/10==1)?k0sqf:k0sqi;
m_cplmax.push_back(CplMax(p_cpl, m_rsf));
m_maincplinfo = QCD_Coupling_Info(cpl, rsf);
m_cplmax.push_back(CplMax(cpl, rsf));
m_cplmax.push_back(0.0);
m_maincplinfo.SetMaxCoupling(&m_cplmax);
// invalidate alternative value and purge associated cache
m_altcplinfo = QCD_Coupling_Info();
m_altcplmax.clear();
return true;
}
void CF_QCD::SetAlternativeUnderlyingCoupling(void *cpl, double sf)
void CF_QCD::SetAlternativeUnderlyingCoupling(MODEL::Running_AlphaS *cpl, double sf)
{
m_altrsf = sf;
assert(cpl != NULL || sf == 1.0);
if (cpl == NULL) {
p_altcpl = NULL;
// invalidate alternative value
m_altcplinfo = QCD_Coupling_Info();
return;
} else {
p_altcpl = static_cast<MODEL::One_Running_AlphaS *>(cpl);
if (m_altcplmax.find(p_altcpl) == m_altcplmax.end()) {
m_altcplmax[p_altcpl] = CplMax(p_altcpl, m_altrsf);
m_altcplinfo = QCD_Coupling_Info(cpl, sf);
if (m_altcplmax.find(cpl) == m_altcplmax.end()) {
std::vector<double> altcplmax;
altcplmax.push_back(CplMax(cpl, sf));
altcplmax.push_back(0.0);
m_altcplmax[cpl] = altcplmax;
}
m_altcplinfo.SetMaxCoupling(&(m_altcplmax[cpl]));
}
}
template<class T>
double CF_QCD::CplMax(T * as, double rsf) const
double CF_QCD::CplMax(MODEL::Running_AlphaS * as, double rsf) const
{
double minscale = Min(1.0, CplFac(m_k0sq)) * m_k0sq;
double ct(0.);
......@@ -125,19 +166,15 @@ double CF_QCD::Coupling(const double &scale,const int pol)
{
DEBUG_FUNC("pol="<<pol);
if (pol!=0) return 0.0;
// use nominal or alternative coupling and scale factor
One_Running_AlphaS * const as = (p_altcpl) ? p_altcpl : p_cpl->GetAs();
const double rsf = (p_altcpl) ? m_altrsf : m_rsf;
QCD_Coupling_Info cplinfo = CurrentCouplingInfo();
if (scale<0.0) return m_last = (*as)(sqr(rpa->gen.Ecms()))*m_q;
double t(CplFac(scale)*scale), scl(CplFac(scale)*scale*rsf);
double cpl=as->BoundedAlphaS(scl);
double t(CplFac(scale)*scale), scl(CplFac(scale)*scale*cplinfo.RSF());
double cpl=cplinfo.Coupling()->BoundedAlphaS(scl);
msg_Debugging()<<"t="<<t<<", \\mu_R^2="<<scl<<std::endl;
msg_Debugging()<<"as(t)="<<as->BoundedAlphaS(t)<<std::endl;
msg_Debugging()<<"as(t)="<<cplinfo.Coupling()->BoundedAlphaS(t)<<std::endl;
if (!IsEqual(scl,t)) {
msg_Debugging()<<"as(\\mu_R^2)="<<cpl<<std::endl;
std::vector<double> ths(as->Thresholds(t,scl));
std::vector<double> ths(cplinfo.Coupling()->Thresholds(t,scl));
ths.push_back((scl>t)?scl:t);
ths.insert(ths.begin(),(scl>t)?t:scl);
if (t<scl) std::reverse(ths.begin(),ths.end());
......@@ -148,14 +185,14 @@ double CF_QCD::Coupling(const double &scale,const int pol)
case 1:
// replace as(t) -> as(t)*prod[1-as/2pi*beta(nf)*log(th[i]/th[i-1])]
for (size_t i(1);i<ths.size();++i) {
ct=cpl/M_PI*as->Beta0((ths[i]+ths[i-1])/2.0)*log(ths[i]/ths[i-1]);
ct=cpl/M_PI*cplinfo.Coupling()->Beta0((ths[i]+ths[i-1])/2.0)*log(ths[i]/ths[i-1]);
fac*=1.0-ct;
}
break;
case 2:
// replace as(t) -> as(t)*[1-sum as/2pi*beta(nf)*log(th[i]/th[i-1])]
for (size_t i(1);i<ths.size();++i)
ct+=cpl/M_PI*as->Beta0((ths[i]+ths[i-1])/2.0)*log(ths[i]/ths[i-1]);
ct+=cpl/M_PI*cplinfo.Coupling()->Beta0((ths[i]+ths[i-1])/2.0)*log(ths[i]/ths[i-1]);
fac=1.-ct;
break;
default:
......@@ -172,17 +209,16 @@ double CF_QCD::Coupling(const double &scale,const int pol)
msg_Debugging()<<"as(\\mu_R^2)*(1-ct)="<<cpl<<std::endl;
}
cpl*=m_q;
const double cplmax = (p_altcpl) ? m_altcplmax[p_altcpl] : m_cplmax.front();
if (cpl>cplmax) {
if (cpl>cplinfo.MaxCoupling()->front()) {
msg_Error()<<METHOD<<"(): Value exceeds maximum at t = "
<<sqrt(t)<<" -> \\mu_R = "<<sqrt(scl)
<<", qmin = "<<sqrt(as->CutQ2())<<std::endl;
return m_last = cplmax;
<<", qmin = "<<sqrt(cplinfo.Coupling()->CutQ2())<<std::endl;
return m_last = cplinfo.MaxCoupling()->front();
}
#ifdef DEBUG__Trial_Weight
msg_Debugging()<<"as weight kt = "<<sqrt(CplFac(scale))<<" * "
<<sqrt(scale)<<", \\alpha_s("<<sqrt(scl)<<") = "
<<(*as)[scl]<<", m_q = "<<m_q<<"\n";
<<(*cplinfo.Coupling())[scl]<<", m_q = "<<m_q<<"\n";
#endif
return m_last = cpl;
}
......@@ -217,8 +253,8 @@ bool CF_QCD::AllowSpec(const ATOOLS::Flavour &fl)
double CF_QCD::CplFac(const double &scale) const
{
if (m_kfmode==0) return m_cplfac;
One_Running_AlphaS * const as = (p_altcpl) ? p_altcpl : p_cpl->GetAs();
double nf=as->Nf(scale);
QCD_Coupling_Info cplinfo = CurrentCouplingInfo();
double nf=cplinfo.Coupling()->Nf(scale);
double kfac=exp(-(67.0-3.0*sqr(M_PI)-10.0/3.0*nf)/(33.0-2.0*nf));
return m_cplfac*kfac;
}
......
......@@ -9,6 +9,7 @@
#include "ATOOLS/Org/Exception.H"
#include <algorithm>
#include <cassert>
using namespace ATOOLS;
......@@ -22,29 +23,51 @@ namespace MCATNLO {
class CF_QCD: public SF_Coupling {
protected:
//! Main underlying coupling (set by SetCoupling)
MODEL::Running_AlphaS *p_cpl;
class QCD_Coupling_Info {
//! (Temporary) alternative coupling (set by SetAlternativeUnderlyingCoupling)
MODEL::One_Running_AlphaS *p_altcpl;
public:
//! Alternative renormalisation scale factor (set by SetAlternativeUnderlyingCoupling)
double m_altrsf;
QCD_Coupling_Info():
p_cpl(NULL), m_rsf(1.0), p_cplmax(NULL)
{};
QCD_Coupling_Info(MODEL::Running_AlphaS * cpl, double rsf, std::vector<double> * cplmax=NULL):
p_cpl(cpl), m_rsf(rsf), p_cplmax(cplmax)
{};
//! Buffer of max alphas values to avoid re-calculations
std::map<MODEL::One_Running_AlphaS *, double> m_altcplmax;
MODEL::Running_AlphaS * const Coupling() const { return p_cpl; }
double RSF() const { return m_rsf; }
std::vector<double> * const MaxCoupling() const { return p_cplmax; }
void SetMaxCoupling(std::vector<double> * const cplmax) { p_cplmax = cplmax; }
bool IsValid() const { return (p_cpl != NULL); }
double m_q, m_rsf, m_k0sq;
private:
double B0(const double &nf) const
{
return 11.0/6.0*s_CA-2.0/3.0*s_TR*nf;
}
MODEL::Running_AlphaS * p_cpl;
double m_rsf;
std::vector<double> * p_cplmax;
};
/*!
* Underlying couplings set by SetCoupling and SetAlternativeUnderlyingCoupling, respectively
*
* If the alternative coupling is set, it takes precedence. This can be used for reweighting purposes.
*/
QCD_Coupling_Info m_maincplinfo, m_altcplinfo;
const QCD_Coupling_Info & CurrentCouplingInfo() const
{ return (m_altcplinfo.IsValid() ? m_altcplinfo : m_maincplinfo); }
//! Buffer of max alphas values to avoid re-calculations
std::map<MODEL::Running_AlphaS *, std::vector<double> > m_altcplmax;
double m_q, m_k0sq;
public:
inline CF_QCD(const SF_Key &key):
SF_Coupling(key), p_altcpl(NULL), m_altrsf(1.0), m_altcplmax(), m_k0sq(0.0)
CF_QCD(const SF_Key &key):
SF_Coupling(key),
m_maincplinfo(QCD_Coupling_Info()), m_altcplinfo(QCD_Coupling_Info()),
m_altcplmax(),
m_q(0.0), m_k0sq(0.0)
{
if (key.p_v->in[0].StrongCharge()==8 &&
key.p_v->in[1].StrongCharge()==8 &&
......@@ -63,24 +86,25 @@ namespace MCATNLO {
}
}
double B0(const double &nf) const
{
return 11.0/6.0*s_CA-2.0/3.0*s_TR*nf;
}
bool SetCoupling(MODEL::Model_Base *md,
const double &k0sqi,const double &k0sqf,
const double &isfac,const double &fsfac);
template<class T>
double CplMax(T * as) const;
double CplMax(MODEL::Running_AlphaS * as, double rsf) const;
double Coupling(const double &scale,const int pol,
ATOOLS::Cluster_Amplitude *const sub);
bool AllowSpec(const ATOOLS::Flavour &fl);
double CplFac(const double &scale) const;
bool AllowsAlternativeCouplingUsage() const { return true; }
void SetAlternativeUnderlyingCoupling(void *, double sf=1.0);
void SetAlternativeUnderlyingCoupling(MODEL::Running_AlphaS *, double sf=1.0);
void ColorPoint(Parton *const p) const;
double ColorWeight(const Color_Info &ci) const;
};
}
......@@ -92,70 +116,80 @@ bool CF_QCD::SetCoupling(MODEL::Model_Base *md,
const double &k0sqi,const double &k0sqf,
const double &isfac,const double &fsfac)
{
p_cpl=(MODEL::Running_AlphaS*)md->GetScalarFunction("alpha_S");
SetAlternativeUnderlyingCoupling(NULL);
m_altcplmax.clear(); // buffered values are not valid anymore
m_rsf=ToType<double>(rpa->gen.Variable("RENORMALIZATION_SCALE_FACTOR"));
MODEL::Running_AlphaS * cpl = (MODEL::Running_AlphaS *)md->GetScalarFunction("alpha_S");
double rsf = ToType<double>(rpa->gen.Variable("RENORMALIZATION_SCALE_FACTOR"));
// determine prefactors before calling CplMax below
m_cplfac=((m_type/10==1)?fsfac:isfac);
m_k0sq=(m_type/10==1)?k0sqf:k0sqi;
m_cplmax.push_back(CplMax(p_cpl));
m_maincplinfo = QCD_Coupling_Info(cpl, rsf);
m_cplmax.push_back(CplMax(cpl, rsf));
m_cplmax.push_back(0.0);
m_maincplinfo.SetMaxCoupling(&m_cplmax);
// invalidate alternative value and purge associated cache
m_altcplinfo = QCD_Coupling_Info();
m_altcplmax.clear();
return true;
}
void CF_QCD::SetAlternativeUnderlyingCoupling(void *cpl, double sf)
void CF_QCD::SetAlternativeUnderlyingCoupling(MODEL::Running_AlphaS *cpl, double sf)
{
m_altrsf = sf;
assert(cpl != NULL || sf == 1.0);
if (cpl == NULL) {
p_altcpl = NULL;
// invalidate alternative value
m_altcplinfo = QCD_Coupling_Info();
return;
} else {
p_altcpl = static_cast<MODEL::One_Running_AlphaS *>(cpl);
if (m_altcplmax.find(p_altcpl) == m_altcplmax.end()) {
m_altcplmax[p_altcpl] = CplMax(p_altcpl);
m_altcplinfo = QCD_Coupling_Info(cpl, sf);
if (m_altcplmax.find(cpl) == m_altcplmax.end()) {
std::vector<double> altcplmax;
altcplmax.push_back(CplMax(cpl, sf));
altcplmax.push_back(0.0);
m_altcplmax[cpl] = altcplmax;
}
m_altcplinfo.SetMaxCoupling(&(m_altcplmax[cpl]));
}
}
template<class T>
double CF_QCD::CplMax(T * as) const
double CF_QCD::CplMax(MODEL::Running_AlphaS * as, double rsf) const
{
double minscale = Min(1.0, CplFac(m_k0sq)) * m_k0sq;
return (*as)(Max(as->CutQ2(), minscale))*m_q;
return as->BoundedAlphaS(minscale) * m_q;
}
double CF_QCD::Coupling(const double &scale,const int pol,
Cluster_Amplitude *const sub)
{
if (pol!=0) return 0.0; // we do not update m_last when polarized
One_Running_AlphaS * const as = (p_altcpl) ? p_altcpl : p_cpl->GetAs();
const double rsf = (p_altcpl) ? m_altrsf : m_rsf;
double t(CplFac(scale)*scale), scl(sub?sub->MuR2():t*rsf);
if (scl<rsf*as->CutQ2()) return m_last = 0.0;
double cpl=(*as)(scl);
if (pol!=0) return 0.0;
QCD_Coupling_Info cplinfo = CurrentCouplingInfo();
double t(CplFac(scale)*scale), scl(sub?sub->MuR2():t*cplinfo.RSF());
if (scl<cplinfo.RSF()*cplinfo.Coupling()->CutQ2()) return m_last = 0.0;
double cpl=(*cplinfo.Coupling())(scl);
if (sub==NULL && !IsEqual(scl,t)) {
std::vector<double> ths(as->Thresholds(t,scl));
std::vector<double> ths(cplinfo.Coupling()->Thresholds(t,scl));
if (scl>t) std::reverse(ths.begin(),ths.end());
if (ths.empty() || !IsEqual(t,ths.back())) ths.push_back(t);
if (!IsEqual(scl,ths.front())) ths.insert(ths.begin(),scl);
for (size_t i(1);i<ths.size();++i) {
double nf=as->Nf((ths[i]+ths[i-1])/2.0);
double nf=cplinfo.Coupling()->Nf((ths[i]+ths[i-1])/2.0);
double L=log(ths[i]/ths[i-1]), ct=cpl/(2.0*M_PI)*B0(nf)*L;
cpl*=1.0-ct;
}
}
cpl*=m_q*s_qfac;
const double cplmax = (p_altcpl) ? s_qfac*m_altcplmax[p_altcpl] : s_qfac*m_cplmax.front();
if (cpl>cplmax) {
if (cpl>cplinfo.MaxCoupling()->front()) {
msg_Error()<<METHOD<<"(): Value exceeds maximum at k_T = "
<<sqrt(scale)<<" -> q = "<<sqrt(scl)<<"."<<std::endl;
return m_last = cplmax;
return m_last = cplinfo.MaxCoupling()->front();
}
#ifdef DEBUG__Trial_Weight
msg_Debugging()<<"as weight kt = "<<(sub?1.0:sqrt(CplFac(scale)))
<<" * "<<(sub?sqrt(scl):sqrt(scale))<<", \\alpha_s("
<<sqrt(scl)<<") = "<<(*as)[scl]
<<sqrt(scl)<<") = "<<(*cplinfo.Coupling())[scl]
<<", m_q = "<<s_qfac<<" * "<<m_q<<"\n";
#endif
return m_last = cpl;
......@@ -170,8 +204,8 @@ double CF_QCD::CplFac(const double &scale) const
{
if (m_kfmode==-1) return 1.0;
if (m_kfmode==0) return m_cplfac;
One_Running_AlphaS * const as = (p_altcpl) ? p_altcpl : p_cpl->GetAs();
double nf=as->Nf(scale);
QCD_Coupling_Info cplinfo = CurrentCouplingInfo();
double nf=cplinfo.Coupling()->Nf(scale);
double kfac=exp(-(67.0-3.0*sqr(M_PI)-10.0/3.0*nf)/(33.0-2.0*nf));
return m_cplfac*kfac;
}
......
......@@ -31,33 +31,12 @@ namespace MODEL {
}
One_Running_AlphaS::One_Running_AlphaS(const std::string pdfname, const int member,
const double as_MZ, const double m2_MZ,
const int order, const int thmode)
{
InitGenericPDF(pdfname, member);
Init(as_MZ, m2_MZ, order, thmode);
}
One_Running_AlphaS::One_Running_AlphaS(PDF::PDF_Base *const pdf,
const double as_MZ, const double m2_MZ,
const int order, const int thmode):
p_pdf(pdf),
m_pdfowned(false)
m_order(order), m_pdf(0), m_mzset(0), m_as_MZ(as_MZ), m_cutq2(0.0), p_pdf(pdf)
{
Init(as_MZ, m2_MZ, order, thmode);
}
void One_Running_AlphaS::Init(const double as_MZ, const double m2_MZ,
const int order, const int thmode)
{
m_pdf = 0;
m_cutq2 = 0.0;
m_as_MZ = as_MZ;
m_m2_MZ = (m2_MZ != 0.0) ? m2_MZ : Flavour(kf_Z).Mass();
m_order = order;
Data_Reader dataread(" ",";","!","=");
dataread.AddComment("#");
......@@ -133,7 +112,6 @@ void One_Running_AlphaS::Init(const double as_MZ, const double m2_MZ,
for (int i(0);i<m_nth;++i) masses[i]=sortmass[i];
int j = 0;
m_mzset = 0;
for (int i=0; i<m_nth; ++j) {
if ((masses[i]>m_m2_MZ)&&(!m_mzset)) {
//insert Z boson (starting point for any evaluation)
......@@ -191,34 +169,9 @@ void One_Running_AlphaS::Init(const double as_MZ, const double m2_MZ,
}
}
void One_Running_AlphaS::InitGenericPDF(const std::string pdfname, const int member)
{
if (m_pdfowned) delete p_pdf;
p_pdf = CreateGenericPDF(pdfname, member);
m_pdfowned = true;
}
PDF::PDF_Base * One_Running_AlphaS::CreateGenericPDF(const std::string name, const int member)
{
// alphaS should be the same for all hadrons, so we can use a proton (as good as any)
if (s_kftable.find(kf_p_plus)==s_kftable.end()) {
s_kftable[kf_p_plus] = new Particle_Info(kf_p_plus,0.938272,0,3,1,1,1,"P+","P^{+}");
}
Data_Reader dataread(" ",";","!","=");
dataread.AddComment("#");
dataread.AddWordSeparator("\t");
PDF::PDF_Arguments args(Flavour(kf_p_plus), &dataread, 0, name, member);
PDF::PDF_Base * pdf = PDF_Base::PDF_Getter_Function::GetObject(name, args);
pdf->SetBounds();
return pdf;
}
One_Running_AlphaS::~One_Running_AlphaS()
{
if (p_thresh!=0) { delete [] p_thresh; p_thresh = NULL; }
if (m_pdfowned) {
delete p_pdf;
}
}
double One_Running_AlphaS::Beta0(const int nf) {
......@@ -479,9 +432,8 @@ Running_AlphaS::Running_AlphaS(const double as_MZ,const double m2_MZ,
const PDF::ISR_Handler_Map &isr)
{
m_defval = as_MZ;
m_type = std::string("Running Coupling");
m_type = "Running Coupling";
m_name = "Alpha_QCD";
// Read possible override
Data_Reader dataread(" ",";","!","=");
dataread.AddComment("#");
......@@ -490,7 +442,7 @@ Running_AlphaS::Running_AlphaS(const double as_MZ,const double m2_MZ,
std::string name(dataread.GetValue<std::string>("ALPHAS_PDF_SET","CT10nlo"));
int member = dataread.GetValue<int>("ALPHAS_PDF_SET_VERSION",0);
member = dataread.GetValue<int>("ALPHAS_PDF_SET_MEMBER",member);
p_overridingpdf = One_Running_AlphaS::CreateGenericPDF(name, member);
InitOverridingPDF(name, member);
} else {
p_overridingpdf = NULL;
}
......@@ -512,6 +464,49 @@ Running_AlphaS::Running_AlphaS(const double as_MZ,const double m2_MZ,
SetActiveAs(PDF::isr::hard_process);
}
Running_AlphaS::Running_AlphaS(const std::string pdfname, const int member,
const double as_MZ, const double m2_MZ,
const int order, const int thmode)
{
m_type = "Running Coupling";
m_name = "Alpha_QCD";
InitOverridingPDF(pdfname, member);
m_alphas.insert(make_pair(PDF::isr::id::none, new One_Running_AlphaS(p_overridingpdf, as_MZ, m2_MZ, order, thmode)));
SetActiveAs(PDF::isr::id::none);
m_defval = AsMZ();
}
Running_AlphaS::Running_AlphaS(PDF::PDF_Base *const pdf,
const double as_MZ, const double m2_MZ,
const int order, const int thmode):
p_overridingpdf(NULL)
{
m_type = "Running Coupling";
m_name = "Alpha_QCD";
m_alphas.insert(make_pair(PDF::isr::id::none, new One_Running_AlphaS(pdf, as_MZ, m2_MZ, order, thmode)));
SetActiveAs(PDF::isr::id::none);
m_defval = AsMZ();
}
void Running_AlphaS::InitOverridingPDF(const std::string name, const int member)
{
if (p_overridingpdf != NULL) {
delete p_overridingpdf;
p_overridingpdf = NULL;
}
// alphaS should be the same for all hadrons, so we can use a proton (as good as any)
if (s_kftable.find(kf_p_plus)==s_kftable.end()) {
s_kftable[kf_p_plus] = new Particle_Info(kf_p_plus,0.938272,0,3,1,1,1,"P+","P^{+}");
}
Data_Reader dataread(" ",";","!","=");
dataread.AddComment("#");
dataread.AddWordSeparator("\t");
PDF::PDF_Arguments args(Flavour(kf_p_plus), &dataread, 0, name, member);
PDF::PDF_Base * pdf = PDF_Base::PDF_Getter_Function::GetObject(name, args);
pdf->SetBounds();
p_overridingpdf = pdf;
}
Running_AlphaS::~Running_AlphaS()
{
for (AlphasMap::iterator it=m_alphas.begin(); it!=m_alphas.end(); ++it) {
......
......@@ -28,18 +28,12 @@ namespace MODEL {
public:
One_Running_AlphaS(const std::string pdfname, const int member = 0,
const double as_MZ = 0.0,const double m2_MZ = 0.0,
const int order = 0, const int thmode = 1);
One_Running_AlphaS(PDF::PDF_Base *const pdf,
const double as_MZ = 0.0,const double m2_MZ = 0.0,
const int order = 0, const int thmode = 1);
~One_Running_AlphaS();
// construct PDFs that are only used for their AlphaS info and values
static PDF::PDF_Base * CreateGenericPDF(const std::string name, const int member);
int Order() { return m_order; }
double Beta0(const double q2) { return Beta0(Nf(q2)); }
double Beta1(const double q2) { return Beta1(Nf(q2)); }
......@@ -76,13 +70,6 @@ namespace MODEL {
double AlphaSLam(const double,const int);
void ContinueAlphaS(int &);
private:
bool m_pdfowned;
void InitGenericPDF(const std::string pdfname, const int member = 0);
void Init(const double as_MZ, const double m2_MZ,
const int order, const int thmode);