Page MenuHomeHEPForge
Files F10664194

No OneTemporary
Actions

View Options

diff --git a/DipoleShower/AlphaS/alpha_s.cc b/DipoleShower/AlphaS/alpha_s.cc
--- a/DipoleShower/AlphaS/alpha_s.cc
+++ b/DipoleShower/AlphaS/alpha_s.cc
@@ -1,229 +1,233 @@
// -*- C++ -*-
// couplings/alpha_s.cc is part of matchbox
// (C) 2008 Simon Platzer -- sp@particle.uni-karlsruhe.de
#include "alpha_s.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Command.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/PDT/ParticleData.h"
using namespace matchbox;
alpha_s::alpha_s()
: AlphaSBase(), min_active_flavours_(3), max_active_flavours_(6),
matched_(false), scale_factor_(1.), quark_masses_squared_(),
lambda_squared_(), alpha_s_in_(.1176), scale_in_(91.1876*GeV),
lambda_range_(1.*MeV2,1.e6*MeV2), fixed_(false) {
}
alpha_s::~alpha_s() {}
void alpha_s::persistentOutput(PersistentOStream & os) const {
os << min_active_flavours_ << max_active_flavours_ << matched_ << scale_factor_;
for (size_t f = 0; f < 7; ++f)
os << ounit(quark_masses_squared_[f],MeV2)
<< ounit(lambda_squared_[f],MeV2);
os << alpha_s_in_ << ounit(scale_in_,GeV)
<< ounit(lambda_range_.first,MeV2) << ounit(lambda_range_.second,MeV2)
<< fixed_;
}
void alpha_s::persistentInput(PersistentIStream & is, int) {
is >> min_active_flavours_ >> max_active_flavours_ >> matched_ >> scale_factor_;
for (size_t f = 0; f < 7; ++f)
is >> iunit(quark_masses_squared_[f],MeV2)
>> iunit(lambda_squared_[f],MeV2);
is >> alpha_s_in_ >> iunit(scale_in_,GeV)
>> iunit(lambda_range_.first,MeV2) >> iunit(lambda_range_.second,MeV2)
>> fixed_;
}
AbstractClassDescription<alpha_s> alpha_s::initalpha_s;
// Definition of the static class description member.
void alpha_s::Init() {
static ClassDocumentation<alpha_s> documentation
("Base class for strong coupoling as used in matchbox");
static Parameter<alpha_s,unsigned int> interfacemin_active_flavours
("min_active_flavours",
"Minimum number of active flavours",
&alpha_s::min_active_flavours_, 3, 0, 6,
true, false, Interface::limited);
static Parameter<alpha_s,unsigned int> interfacemax_active_flavours
("max_active_flavours",
"Maximum number of active flavours",
&alpha_s::max_active_flavours_, 6, 0, 6,
true, false, Interface::limited);
static Parameter<alpha_s,double> interfaceinput_alpha_s
("input_alpha_s",
"alpha_s value at input scale",
&alpha_s::alpha_s_in_, .1176, 0.0, 1.0,
true, false, Interface::limited);
static Parameter<alpha_s,Energy> interfaceinput_scale
("input_scale",
"Input scale for alpha_s value",
&alpha_s::scale_in_, GeV, 91.1876*GeV, 0.*GeV, 0.*GeV,
true, false, Interface::lowerlim);
static Command<alpha_s> interfacecheck
("check",
"check",
&alpha_s::check, false);
static Parameter<alpha_s,double> interfacescale_factor
("scale_factor",
"scale factor for argument",
&alpha_s::scale_factor_, 1., 0.0, 100.0,
true, false, Interface::limited);
static Switch<alpha_s,bool> interfacefixed
("fixed",
"",
&alpha_s::fixed_, false, false, false);
static SwitchOption interfacefixedOn
(interfacefixed,
"On",
"",
true);
static SwitchOption interfacefixedOff
(interfacefixed,
"Off",
"",
false);
}
string alpha_s::check (string args) {
istringstream argin(args);
double Q_low, Q_high;
long n_steps;
argin >> Q_low >> Q_high >> n_steps;
string fname;
argin >> fname;
cout << "checking alpha_s in range [" << Q_low << "," << Q_high << "] GeV in "
<< n_steps << " steps.\nResults are written to " << fname << "\n";
double step_width = (Q_high-Q_low)/n_steps;
match_thresholds();
cout << "threshold matching results:\n"
<< "(m_Q^2 -> Lambda^2) / GeV^2 for dynamic flavours in range ["
<< min_active_flavours_ << "," << max_active_flavours_ << "]\n";
for (size_t f = 0; f < 7; ++f) {
cout << (quark_masses_squared_[f]/GeV2) << " "
<< (lambda_squared_[f]/GeV2) << "\n";
}
ofstream out (fname.c_str());
for (long k = 0; k <= n_steps; ++k) {
Energy Q = Q_low*GeV + k*step_width*GeV;
out << (Q/GeV) << " " << (operator () (Q*Q)) << "\n";
}
return "alpha_s check finished";
}
void alpha_s::match_thresholds () {
if (matched_)
return;
// get the quark masses
quark_masses_squared_[0] = 0.*MeV2;
for (long f = 1; f < 7; ++f) {
- quark_masses_squared_[static_cast<size_t>(f)]
- = sqr(getParticleData(f)->mass());
+ if ( quarkMasses().empty() )
+ quark_masses_squared_[static_cast<size_t>(f)]
+ = sqr(getParticleData(f)->mass());
+ else
+ quark_masses_squared_[static_cast<size_t>(f)]
+ = sqr(quarkMasses()[static_cast<size_t>(f-1)]);
}
unsigned int active_at_input = active_flavours(sqr(scale_in_));
// solve for input lambda
solve_input_lambda<alpha_s> input_equation (this,active_at_input,alpha_s_in_,sqr(scale_in_));
gsl::bisection_root_solver<solve_input_lambda<alpha_s>,100> input_solver(input_equation);
lambda_squared_[active_at_input] =
MeV2 *
input_solver.solve(make_pair(lambda_range_.first/MeV2,
lambda_range_.second/MeV2));
// get lambdas down to min active flavours
unsigned int below = active_at_input;
while (below > min_active_flavours_) {
solve_lambda_below<alpha_s> match_equation (this,below,
lambda_squared_[below],
quark_masses_squared_[below]);
gsl::bisection_root_solver<solve_lambda_below<alpha_s>,100> match_solver(match_equation);
lambda_squared_[below-1] =
MeV2 *
match_solver.solve(make_pair(lambda_range_.first/MeV2,
lambda_range_.second/MeV2));
--below;
}
// get lambdas up to max active flavours
unsigned int above = active_at_input;
while (above < max_active_flavours_) {
solve_lambda_above<alpha_s> match_equation (this,above,
lambda_squared_[above],
quark_masses_squared_[above+1]);
gsl::bisection_root_solver<solve_lambda_above<alpha_s>,100> match_solver(match_equation);
lambda_squared_[above+1] =
MeV2 *
match_solver.solve(make_pair(lambda_range_.first/MeV2,
lambda_range_.second/MeV2));
++above;
}
if (min_active_flavours_ > 0) {
for (size_t f = 0; f < min_active_flavours_; ++f) {
lambda_squared_[f] = lambda_squared_[min_active_flavours_];
}
}
if (max_active_flavours_ < 6) {
for (size_t f = max_active_flavours_+1; f < 7; ++f) {
lambda_squared_[f] = lambda_squared_[max_active_flavours_];
}
}
matched_ = true;
return;
}
diff --git a/Models/StandardModel/O2AlphaS.cc b/Models/StandardModel/O2AlphaS.cc
--- a/Models/StandardModel/O2AlphaS.cc
+++ b/Models/StandardModel/O2AlphaS.cc
@@ -1,167 +1,170 @@
// -*- C++ -*-
//
// O2AlphaS.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the O2AlphaS class.
//
#include "O2AlphaS.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/PDT/ParticleData.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
void O2AlphaS::persistentOutput(PersistentOStream & os) const {
os << ounit(_lambdaQCD,GeV) << _bcoeff << _ccoeff << ounit(_lambdas,GeV)
<< ounit(_threshold,GeV) << _match << _copt;
}
void O2AlphaS::persistentInput(PersistentIStream & is, int) {
is >> iunit(_lambdaQCD,GeV) >> _bcoeff >> _ccoeff >> iunit(_lambdas,GeV)
>> iunit(_threshold,GeV) >> _match >> _copt;
}
ClassDescription<O2AlphaS> O2AlphaS::initO2AlphaS;
// Definition of the static class description member.
void O2AlphaS::Init() {
static ClassDocumentation<O2AlphaS> documentation
("The O2AlphaS class implements the next-to-leading order alphaS in the same"
" way as in FORTRAN HERWIG");
static Parameter<O2AlphaS,Energy> interfaceLambdaQCD
("LambdaQCD",
"The value of Lambda QCD",
&O2AlphaS::_lambdaQCD, MeV, 180.*MeV, 50.*MeV, 500.0*MeV,
false, false, Interface::limited);
static Switch<O2AlphaS,unsigned int> interfaceLambdaType
("LambdaType",
"Which type of Lambda to use",
&O2AlphaS::_copt, 0, false, false);
static SwitchOption interfaceLambdaTypeMonteCarlo
(interfaceLambdaType,
"MonteCarlo",
"Use a Monte Carlo scheme as in the FORTRAN program",
0);
static SwitchOption interfaceLambdaTypeMSBar
(interfaceLambdaType,
"MSBar",
"Use the MSBar scheme",
1);
}
vector<Energy2> O2AlphaS::flavourThresholds() const {
vector<Energy2> thresholds(_threshold.size());
transform(_threshold.begin(), _threshold.end(),
thresholds.begin(),
sqr<Energy>);
return thresholds;
}
void O2AlphaS::doinit() {
// thresholds
for ( int ix=1; ix<7; ++ix ) {
- tPDPtr p = getParticleData(ix);
- _threshold[ix-1] = p->mass();
+ if ( quarkMasses().empty() ) {
+ tPDPtr p = getParticleData(ix);
+ _threshold[ix-1] = p->mass();
+ } else
+ _threshold[ix-1] = quarkMasses()[ix-1];
}
// d is heavier than u, need to swap
swap(_threshold[0],_threshold[1]);
// beta function coefficients
const double ca = generator()->standardModel()->Nc();
const double cf = (sqr(ca)-1.)/2./ca;
for(unsigned int ix=3;ix<7;++ix)
{
_bcoeff[ix-1]=(11.*ca-2.*ix)/(12.*Constants::pi);
_ccoeff[ix-1]=(17.*sqr(ca)-ix*(5.*ca+3.*cf))/(24.*sqr(Constants::pi))/sqr(_bcoeff[ix-1]);
}
if(_copt==0)
{
double kfac(ca*(67./18.-sqr(Constants::pi)/6.)-25./9.);
_lambdas[5]=_lambdaQCD*exp(kfac/(4.*Constants::pi*_bcoeff[4]))/sqrt(2.);
}
else{_lambdas[5]=_lambdaQCD;}
// calculate the threshold matching
double rho=2.*log(_threshold[5]/_lambdas[5]);
double rat=log(rho)/rho;
_match[5]=(_bcoeff[4]/(1.-_ccoeff[4]*rat)-_bcoeff[5]/(1.-_ccoeff[5]*rat))*rho;
rho=2.*log(_threshold[4]/_lambdas[5]);
rat=log(rho)/rho;
_match[4]=(_bcoeff[4]/(1.-_ccoeff[4]*rat)-_bcoeff[3]/(1.-_ccoeff[3]*rat))*rho;
rho=2.*log(_threshold[3]/_lambdas[5]);
rat=log(rho)/rho;
_match[3]=(_bcoeff[3]/(1.-_ccoeff[3]*rat)-_bcoeff[2]/(1.-_ccoeff[2]*rat))*rho
+_match[4];
// calculate the 4-flavour lambda
_lambdas[4]=_lambdas[5]*pow(_threshold[4]/_lambdas[5],2./25.)*
pow(2.*log(_threshold[4]/_lambdas[5]),963./14375.);
// calculate the 3-flavour lambda
double eps(1.e-6),d35(-1./(_bcoeff[2]*_match[3])),rlf,drh;
unsigned int ix=0;
do
{
rat=log(d35)/d35;
rlf=_bcoeff[2]*d35/(1.-_ccoeff[2]*rat);
drh=_bcoeff[2]*(rlf+_match[3])*sqr(d35)/
((1.-2.*_ccoeff[2]*rat+_ccoeff[2]/d35)*sqr(rlf));
d35=d35-drh;
++ix;
}
while(ix<100&&abs(drh)>eps*d35);
_lambdas[3]=_lambdas[5]*exp(0.5*d35);
AlphaSBase::doinit();
}
vector<Energy> O2AlphaS::LambdaQCDs() const
{
vector<Energy> output(4,_lambdas[3]);
output.push_back(_lambdas[4]);
output.push_back(_lambdas[5]);
output.push_back(_lambdas[5]);
return output;
}
double O2AlphaS::value(Energy2 scale, const StandardModelBase &) const
{
Energy rs=sqrt(scale);
if(scale<sqr(_lambdas[5])) {
cerr << "O2AlphaS called with scale less than Lambda QCD "
<< "scale = " << rs/MeV << " MeV and "
<< "Lambda = " << _lambdas[5]/MeV << " MeV\n";
generator()->log() << "O2AlphaS called with scale less than Lambda QCD "
<< "scale = " << rs/MeV << " MeV and "
<< "Lambda = " << _lambdas[5]/MeV << " MeV\n";
return 0.;
}
double rho=2.*log(rs/_lambdas[5]),rat(log(rho)/rho);
double rlf;
if(rs>_threshold[5]) rlf=_bcoeff[5]*rho/(1.-_ccoeff[5]*rat)+_match[5];
else if(rs>_threshold[4]) rlf=_bcoeff[4]*rho/(1.-_ccoeff[4]*rat);
else if(rs>_threshold[3]) rlf=_bcoeff[3]*rho/(1.-_ccoeff[3]*rat)+_match[4];
else rlf=_bcoeff[2]*rho/(1.-_ccoeff[2]*rat)+_match[3];
// must be possible
if(rlf<=0.) {
generator()->log() << "O2AlphaS coupling less than zero \n";
return 0.;
}
return 1./rlf;
}

File Metadata

Mime Type
text/x-diff
Expires
Thu, Apr 24, 6:32 AM (1 d, 12 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4887745
Default Alt Text
(12 KB)

Event Timeline