Page MenuHomeHEPForge
Files F8308872

AlphaStrong.cpp
No OneTemporary
Actions

AlphaStrong.cpp
View Options

//==============================================================================
// AlphaStrong.cpp
//
// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre version: 1.1
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// Author: Thomas Ullrich
// Last update:
// $Date: 2013-05-29 21:26:19 +0100 (Wed, 29 May 2013) $
// $Author: thomas.ullrich@bnl.gov $
//==============================================================================
//
// Original code from Pythia8:
// Copyright (C) 2010 Torbjorn Sjostrand.
// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
//
// The code was modified for Sartre and made multithread safe. Several
// methods that are not needed in Sartre where removed.
//==============================================================================
#include "AlphaStrong.h"
#include <iostream>
#include <cmath>
#define ALPHASTRONG_REMEMBERS_LAST_CALL
using namespace std;
const double AlphaStrong::mMassC = 1.5;
const double AlphaStrong::mMassB = 4.8;
const double AlphaStrong::mMassZ = 91.188;
const double AlphaStrong::mSafetyMargin1 = 1.07;
const double AlphaStrong::mSafetyMargin2 = 1.33;
AlphaStrong::AlphaStrong(double as, int order)
{
mIsInit = false;
init(as, order) ;
}
void AlphaStrong::init(double value, int order) {
// Order of alpha_s evaluation
mValueRef = value;
mOrder = order;
if (mOrder > 2 || mOrder < 0) {
cout << "AlphaStrong::init(): error, order must be 0, 1, or 2." << endl;
return;
}
// Fix alpha_s
if (mOrder == 0) {
mLambda3Save = mLambda4Save = mLambda5Save = mScale2Min = 0.;
}
// First order alpha_s: match at flavour thresholds
else if (mOrder == 1) {
mLambda5Save = mMassZ * exp( -6. * M_PI / (23. * mValueRef) );
mLambda4Save = mLambda5Save * pow(mMassB/mLambda5Save, 2./25.);
mLambda3Save = mLambda4Save * pow(mMassC/mLambda4Save, 2./27.);
mScale2Min = pow2(mSafetyMargin1 * mLambda3Save);
}
// Second order alpha_s: iterative match at flavour thresholds
else {
double b15 = 348. / 529.;
double b14 = 462. / 625.;
double b13 = 64. / 81.;
double b25 = 224687. / 242208.;
double b24 = 548575. / 426888.;
double b23 = 938709. / 663552.;
double logScale, loglogScale, correction, valueIter;
// Find Lambda_5 at m_Z
mLambda5Save = mMassZ * exp( -6. * M_PI / (23. * mValueRef) );
for (int iter = 0; iter < mMaxIter; ++iter) {
logScale = 2. * log(mMassZ/mLambda5Save);
loglogScale = log(logScale);
correction = 1. - b15 * loglogScale / logScale
+ pow2(b15 / logScale) * (pow2(loglogScale - 0.5) + b25 - 1.25);
valueIter = mValueRef / correction;
mLambda5Save = mMassZ * exp( -6. * M_PI / (23. * valueIter) );
}
// Find Lambda_4 at m_b
double logScaleB = 2. * log(mMassB/mLambda5Save);
double loglogScaleB = log(logScaleB);
double valueB = 12. * M_PI / (23. * logScaleB)
* (1. - b15 * loglogScaleB / logScaleB
+ pow2(b15 / logScaleB) * (pow2(loglogScaleB - 0.5) + b25- 1.25) );
mLambda4Save = mLambda5Save;
for (int iter = 0; iter < mMaxIter; ++iter) {
logScale = 2. * log(mMassB/mLambda4Save);
loglogScale = log(logScale);
correction = 1. - b14 * loglogScale / logScale
+ pow2(b14 / logScale) * (pow2(loglogScale - 0.5) + b24 - 1.25);
valueIter = valueB / correction;
mLambda4Save = mMassB * exp( -6. * M_PI / (25. * valueIter) );
}
// Find Lambda_3 at m_c
double logScaleC = 2. * log(mMassC/mLambda4Save);
double loglogScaleC = log(logScaleC);
double valueC = 12. * M_PI / (25. * logScaleC)
* (1. - b14 * loglogScaleC / logScaleC
+ pow2(b14 / logScaleC) * (pow2(loglogScaleC - 0.5) + b24 - 1.25) );
mLambda3Save = mLambda4Save;
for (int iter = 0; iter < mMaxIter; ++iter) {
logScale = 2. * log(mMassC/mLambda3Save);
loglogScale = log(logScale);
correction = 1. - b13 * loglogScale / logScale
+ pow2(b13 / logScale) * (pow2(loglogScale - 0.5) + b23 - 1.25);
valueIter = valueC / correction;
mLambda3Save = mMassC * exp( -6. * M_PI / (27. * valueIter) );
}
mScale2Min = pow2(mSafetyMargin2 * mLambda3Save);
}
// Save squares of mass and Lambda values as well
mMassC2 = pow2(mMassC);
mMassB2 = pow2(mMassB);
mLambda3Save2 = pow2(mLambda3Save);
mLambda4Save2 = pow2(mLambda4Save);
mLambda5Save2 = pow2(mLambda5Save);
mValueNow = value;
mScale2Now = mMassZ * mMassZ;
mIsInit = true;
}
double AlphaStrong::alphaS( double scale2) {
// Check for initialization and ensure minimal scale2 value
if (!mIsInit) return 0.;
if (scale2 < mScale2Min) scale2 = mScale2Min;
#if defined(ALPHASTRONG_REMEMBERS_LAST_CALL)
// If equal to old scale then same answer
if (scale2 == mScale2Now) return mValueNow;
mScale2Now = scale2;
#endif
// Fix alpha_s
if (mOrder == 0) {
mValueNow = mValueRef;
}
// First order alpha_s: differs by mass region
else if (mOrder == 1) {
if (scale2 > mMassB2)
mValueNow = 12. * M_PI / (23. * log(scale2/mLambda5Save2));
else if (scale2 > mMassC2)
mValueNow = 12. * M_PI / (25. * log(scale2/mLambda4Save2));
else
mValueNow = 12. * M_PI / (27. * log(scale2/mLambda3Save2));
}
// Second order alpha_s: differs by mass region
else {
double Lambda2, b0, b1, b2;
if (scale2 > mMassB2) {
Lambda2 = mLambda5Save2;
b0 = 23.;
b1 = 348. / 529.;
b2 = 224687. / 242208.;
}
else if (scale2 > mMassC2) {
Lambda2 = mLambda4Save2;
b0 = 25.;
b1 = 462. / 625.;
b2 = 548575. / 426888.;
}
else {
Lambda2 = mLambda3Save2;
b0 = 27.;
b1 = 64. / 81.;
b2 = 938709. / 663552.;
}
double logScale = log(scale2/Lambda2);
double loglogScale = log(logScale);
mValueNow = 12. * M_PI / (b0 * logScale)
* ( 1. - b1 * loglogScale / logScale
+ pow2(b1 / logScale) * (pow2(loglogScale - 0.5) + b2 - 1.25) );
}
return mValueNow;
}

File Metadata

Mime Type
text/x-c
Expires
Sat, Dec 21, 1:26 PM (19 h, 5 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4002677
Default Alt Text
AlphaStrong.cpp (7 KB)

Event Timeline