Commit 679ab589 authored by Frank Siegert's avatar Frank Siegert Committed by Johannes Krause

Introduce CSS_EVOLUTION_SCHEME=2,3 which are equivalent to 0,1 but take masses into account.

parent b7c8d270
......@@ -72,30 +72,43 @@ double Kinematics_FF::GetKT2(const double &Q2,const double &y,const double &z,
const ATOOLS::Flavour &fla,const ATOOLS::Flavour &flc) const
{
double pipj=(Q2-mi2-mj2-mk2)*y;
if (m_evolscheme==0) {
return pipj*z*(1.0-z)-sqr(1.0-z)*mi2-sqr(z)*mj2;
if (m_evolscheme==0 || m_evolscheme==2) {
double kt2=pipj*z*(1.0-z)-sqr(1.0-z)*mi2-sqr(z)*mj2;
if (m_evolscheme==0) return kt2;
if (m_evolscheme==2) return kt2+mi2+mj2;
}
double kt2=pipj*z*(1.0-z);
if (fla.IsFermion()) kt2=pipj*(flc.IsVector()?(1.0-z):z);
else if (flc.IsFermion()) kt2=pipj;
return kt2;
else if (m_evolscheme==1 || m_evolscheme==3) {
double kt2=pipj*z*(1.0-z);
if (fla.IsFermion()) kt2=pipj*(flc.IsVector()?(1.0-z):z);
else if (flc.IsFermion()) kt2=pipj;
if (m_evolscheme==1) return kt2;
if (m_evolscheme==3) return kt2+mi2+mj2;
}
else THROW(fatal_error, "Not implemented");
return 0.0;
}
double Kinematics_FF::GetY(const double &Q2,const double &kt2,const double &z,
double Kinematics_FF::GetY(const double &Q2,const double &_kt2,const double &z,
const double &mi2,const double &mj2,const double &mk2,
const ATOOLS::Flavour &fla,const ATOOLS::Flavour &flc,
const bool force) const
{
if (!force && (z<=0.0 || z>=1.0 || Q2<=mi2+mj2+mk2)) return -1.0;
if (m_evolscheme==0) {
double kt2=_kt2;
if (m_evolscheme==2 || m_evolscheme==3) kt2=kt2-mi2-mj2;
if (m_evolscheme==0 || m_evolscheme==2) {
return (kt2/(z*(1.0-z))+(1.0-z)/z*mi2+z/(1.0-z)*mj2)/(Q2-mi2-mj2-mk2);
}
if (fla.IsFermion()) {
if (flc.IsFermion()) return kt2/z/(Q2-mi2-mj2-mk2);
return kt2/(1.0-z)/(Q2-mi2-mj2-mk2);
else if (m_evolscheme==1 || m_evolscheme==3) {
if (fla.IsFermion()) {
if (flc.IsFermion()) return kt2/z/(Q2-mi2-mj2-mk2);
return kt2/(1.0-z)/(Q2-mi2-mj2-mk2);
}
if (flc.IsFermion()) return kt2/(Q2-mi2-mj2-mk2);
return kt2/(z*(1.0-z))/(Q2-mi2-mj2-mk2);
}
if (flc.IsFermion()) return kt2/(Q2-mi2-mj2-mk2);
return kt2/(z*(1.0-z))/(Q2-mi2-mj2-mk2);
else THROW(fatal_error, "Not implemented");
return 0.0;
}
int Kinematics_FF::MakeKinematics
......@@ -147,30 +160,44 @@ double Kinematics_FI::GetKT2(const double &Q2,const double &y,const double &z,
const ATOOLS::Flavour &fla,const ATOOLS::Flavour &flc) const
{
double pipj=-(Q2-ma2-mi2-mj2)*(1.0-y)/y;
if (m_evolscheme==0) {
return pipj*z*(1.0-z)-sqr(1.0-z)*mi2-sqr(z)*mj2;
if (m_evolscheme==0 || m_evolscheme==2) {
double kt2=pipj*z*(1.0-z)-sqr(1.0-z)*mi2-sqr(z)*mj2;
if (m_evolscheme==0) return kt2;
if (m_evolscheme==2) return kt2+mi2+mj2;
}
else if (m_evolscheme==1 || m_evolscheme==3) {
double kt2=pipj*z*(1.0-z);
if (fla.IsFermion()) kt2=pipj*(flc.IsVector()?(1.0-z):z);
else if (flc.IsFermion()) kt2=pipj;
if (m_evolscheme==1) return kt2;
if (m_evolscheme==3) return kt2+mi2+mj2;
}
double kt2=pipj*z*(1.0-z);
if (fla.IsFermion()) kt2=pipj*(flc.IsVector()?(1.0-z):z);
else if (flc.IsFermion()) kt2=pipj;
return kt2;
else THROW(fatal_error, "Not implemented");
return 0.0;
}
double Kinematics_FI::GetY(const double &Q2,const double &kt2,const double &z,
double Kinematics_FI::GetY(const double &Q2,const double &_kt2,const double &z,
const double &mi2,const double &mj2,const double &ma2,
const ATOOLS::Flavour &fla,const ATOOLS::Flavour &flc,
const bool force) const
{
if (!force && (z<=0.0 || z>=1.0 || Q2>=mi2+mj2+ma2)) return -1.0;
if (m_evolscheme==0) {
double kt2=_kt2;
if (m_evolscheme==2 || m_evolscheme==3) kt2=kt2-mi2-mj2;
if (m_evolscheme==0 || m_evolscheme==2) {
return 1.0/(1.0-(kt2/(z*(1.0-z))+mi2*(1.0-z)/z+mj2*z/(1.0-z))/(Q2-ma2-mi2-mj2));
}
if (fla.IsFermion()) {
if (flc.IsFermion()) return 1.0/(1.0-kt2/z/(Q2-ma2-mi2-mj2));
return 1.0/(1.0-kt2/(1.0-z)/(Q2-ma2-mi2-mj2));
else if (m_evolscheme==1 || m_evolscheme==3) {
if (fla.IsFermion()) {
if (flc.IsFermion()) return 1.0/(1.0-kt2/z/(Q2-ma2-mi2-mj2));
return 1.0/(1.0-kt2/(1.0-z)/(Q2-ma2-mi2-mj2));
}
if (flc.IsFermion()) return 1.0/(1.0-kt2/(Q2-ma2-mi2-mj2));
return 1.0/(1.0-kt2/(z*(1.0-z))/(Q2-ma2-mi2-mj2));
}
if (flc.IsFermion()) return 1.0/(1.0-kt2/(Q2-ma2-mi2-mj2));
return 1.0/(1.0-kt2/(z*(1.0-z))/(Q2-ma2-mi2-mj2));
else THROW(fatal_error, "Not implemented");
return 0.0;
}
int Kinematics_FI::MakeKinematics
......@@ -222,25 +249,38 @@ double Kinematics_IF::GetKT2(const double &Q2,const double &y,const double &z,
const ATOOLS::Flavour &flb,const ATOOLS::Flavour &flc) const
{
double pipj=(Q2-ma2-mi2-mk2)*y/z;
if (m_evolscheme==0) {
return -pipj*(1.0-z)-mi2-sqr(1.0-z)*ma2;
if (m_evolscheme==0 || m_evolscheme==2) {
double kt2=-pipj*(1.0-z)-mi2-sqr(1.0-z)*ma2;
if (m_evolscheme==0) return kt2;
if (m_evolscheme==2) return kt2+mi2+ma2;
}
double kt2=-pipj*(1.0-z);
if (flc.IsFermion()) kt2=-pipj;
return kt2;
else if (m_evolscheme==1 || m_evolscheme==3) {
double kt2=-pipj*(1.0-z);
if (flc.IsFermion()) kt2=-pipj;
if (m_evolscheme==1) return kt2;
if (m_evolscheme==3) return kt2+mi2+ma2;
}
else THROW(fatal_error, "Not implemented");
return 0.0;
}
double Kinematics_IF::GetY(const double &Q2,const double &kt2,const double &z,
double Kinematics_IF::GetY(const double &Q2,const double &_kt2,const double &z,
const double &ma2,const double &mi2,const double &mk2,
const ATOOLS::Flavour &flb,const ATOOLS::Flavour &flc,
const bool force) const
{
if (!force && (z<=0.0 || z>=1.0 || Q2>=ma2+mi2+mk2)) return -1.0;
if (m_evolscheme==0) {
double kt2=_kt2;
if (m_evolscheme==2 || m_evolscheme==3) kt2=kt2-mi2-ma2;
if (m_evolscheme==0 || m_evolscheme==2) {
return -z/(Q2-ma2-mi2-mk2)*((kt2+mi2)/(1.0-z)+(1.0-z)*ma2);
}
if (flc.IsFermion()) return -z/(Q2-ma2-mi2-mk2)*kt2;
return -z/(Q2-ma2-mi2-mk2)*kt2/(1.0-z);
else if (m_evolscheme==1 || m_evolscheme==3) {
if (flc.IsFermion()) return -z/(Q2-ma2-mi2-mk2)*kt2;
return -z/(Q2-ma2-mi2-mk2)*kt2/(1.0-z);
}
else THROW(fatal_error, "Not implemented");
return 0.0;
}
int Kinematics_IF::MakeKinematics
......@@ -300,25 +340,38 @@ double Kinematics_II::GetKT2(const double &Q2,const double &y,const double &z,
const ATOOLS::Flavour &flb,const ATOOLS::Flavour &flc) const
{
double pipj=(Q2-ma2-mi2-mb2)*y/z;
if (m_evolscheme==0) {
return pipj*(1.0-z)-mi2-sqr(1.0-z)*ma2;
if (m_evolscheme==0 || m_evolscheme==2) {
double kt2=pipj*(1.0-z)-mi2-sqr(1.0-z)*ma2;
if (m_evolscheme==0) return kt2;
if (m_evolscheme==2) return kt2+mi2+ma2;
}
double kt2=pipj*(1.0-z);
if (flc.IsFermion()) kt2=pipj;
return kt2;
else if (m_evolscheme==1 || m_evolscheme==3) {
double kt2=pipj*(1.0-z);
if (flc.IsFermion()) kt2=pipj;
if (m_evolscheme==1) return kt2;
if (m_evolscheme==3) return kt2+mi2+ma2;
}
else THROW(fatal_error, "Not implemented");
return 0.0;
}
double Kinematics_II::GetY(const double &Q2,const double &kt2,const double &z,
double Kinematics_II::GetY(const double &Q2,const double &_kt2,const double &z,
const double &ma2,const double &mi2,const double &mb2,
const ATOOLS::Flavour &flb,const ATOOLS::Flavour &flc,
const bool force) const
{
if (!force && (z<=0.0 || z>=1.0 || Q2<=ma2+mi2+mb2)) return -1.0;
if (m_evolscheme==0) {
double kt2=_kt2;
if (m_evolscheme==2 || m_evolscheme==3) kt2=kt2-mi2-ma2;
if (m_evolscheme==0 || m_evolscheme==2) {
return z/(Q2-ma2-mb2-mi2)*((kt2+mi2)/(1.0-z)+(1.0-z)*ma2);
}
if (flc.IsFermion()) return z/(Q2-ma2-mb2-mi2)*kt2;
return z/(Q2-ma2-mb2-mi2)*kt2/(1.0-z);
else if (m_evolscheme==1 || m_evolscheme==3) {
if (flc.IsFermion()) return z/(Q2-ma2-mb2-mi2)*kt2;
return z/(Q2-ma2-mb2-mi2)*kt2/(1.0-z);
}
else THROW(fatal_error, "Not implemented");
return 0.0;
}
int Kinematics_II::MakeKinematics
......
......@@ -15,30 +15,43 @@ double Kinematics_FF::GetKT2(const double &Q2,const double &y,const double &z,
const ATOOLS::Flavour &fla,const ATOOLS::Flavour &flc) const
{
double pipj=(Q2-mi2-mj2-mk2)*y;
if (m_evolscheme==0) {
return pipj*z*(1.0-z)-sqr(1.0-z)*mi2-sqr(z)*mj2;
if (m_evolscheme==0 || m_evolscheme==2) {
double kt2=pipj*z*(1.0-z)-sqr(1.0-z)*mi2-sqr(z)*mj2;
if (m_evolscheme==0) return kt2;
if (m_evolscheme==2) return kt2+mi2+mj2;
}
double kt2=pipj*z*(1.0-z);
if (fla.IsFermion()) kt2=pipj*(flc.IsVector()?(1.0-z):z);
else if (flc.IsFermion()) kt2=pipj;
return kt2;
else if (m_evolscheme==1 || m_evolscheme==3) {
double kt2=pipj*z*(1.0-z);
if (fla.IsFermion()) kt2=pipj*(flc.IsVector()?(1.0-z):z);
else if (flc.IsFermion()) kt2=pipj;
if (m_evolscheme==1) return kt2;
if (m_evolscheme==3) return kt2+mi2+mj2;
}
else THROW(fatal_error, "Not implemented");
return 0.0;
}
double Kinematics_FF::GetY(const double &Q2,const double &kt2,const double &z,
double Kinematics_FF::GetY(const double &Q2,const double &_kt2,const double &z,
const double &mi2,const double &mj2,const double &mk2,
const ATOOLS::Flavour &fla,const ATOOLS::Flavour &flc,
const bool force) const
{
if (!force && (z<=0.0 || z>=1.0 || Q2<=mi2+mj2+mk2)) return -1.0;
if (m_evolscheme==0) {
double kt2=_kt2;
if (m_evolscheme==2 || m_evolscheme==3) kt2=kt2-mi2-mj2;
if (m_evolscheme==0 || m_evolscheme==2) {
return (kt2/(z*(1.0-z))+(1.0-z)/z*mi2+z/(1.0-z)*mj2)/(Q2-mi2-mj2-mk2);
}
if (fla.IsFermion()) {
if (flc.IsFermion()) return kt2/z/(Q2-mi2-mj2-mk2);
return kt2/(1.0-z)/(Q2-mi2-mj2-mk2);
else if (m_evolscheme==1 || m_evolscheme==3) {
if (fla.IsFermion()) {
if (flc.IsFermion()) return kt2/z/(Q2-mi2-mj2-mk2);
return kt2/(1.0-z)/(Q2-mi2-mj2-mk2);
}
if (flc.IsFermion()) return kt2/(Q2-mi2-mj2-mk2);
return kt2/(z*(1.0-z))/(Q2-mi2-mj2-mk2);
}
if (flc.IsFermion()) return kt2/(Q2-mi2-mj2-mk2);
return kt2/(z*(1.0-z))/(Q2-mi2-mj2-mk2);
else THROW(fatal_error, "Not implemented");
return 0.0;
}
int Kinematics_FF::MakeKinematics
......@@ -75,30 +88,44 @@ double Kinematics_FI::GetKT2(const double &Q2,const double &y,const double &z,
const ATOOLS::Flavour &fla,const ATOOLS::Flavour &flc) const
{
double pipj=-(Q2-ma2-mi2-mj2)*(1.0-y)/y;
if (m_evolscheme==0) {
return pipj*z*(1.0-z)-sqr(1.0-z)*mi2-sqr(z)*mj2;
if (m_evolscheme==0 || m_evolscheme==2) {
double kt2=pipj*z*(1.0-z)-sqr(1.0-z)*mi2-sqr(z)*mj2;
if (m_evolscheme==0) return kt2;
if (m_evolscheme==2) return kt2+mi2+mj2;
}
else if (m_evolscheme==1 || m_evolscheme==3) {
double kt2=pipj*z*(1.0-z);
if (fla.IsFermion()) kt2=pipj*(flc.IsVector()?(1.0-z):z);
else if (flc.IsFermion()) kt2=pipj;
if (m_evolscheme==1) return kt2;
if (m_evolscheme==3) return kt2+mi2+mj2;
}
double kt2=pipj*z*(1.0-z);
if (fla.IsFermion()) kt2=pipj*(flc.IsVector()?(1.0-z):z);
else if (flc.IsFermion()) kt2=pipj;
return kt2;
else THROW(fatal_error, "Not implemented");
return 0.0;
}
double Kinematics_FI::GetY(const double &Q2,const double &kt2,const double &z,
double Kinematics_FI::GetY(const double &Q2,const double &_kt2,const double &z,
const double &mi2,const double &mj2,const double &ma2,
const ATOOLS::Flavour &fla,const ATOOLS::Flavour &flc,
const bool force) const
{
if (!force && (z<=0.0 || z>=1.0 || Q2>=mi2+mj2+ma2)) return -1.0;
if (m_evolscheme==0) {
double kt2=_kt2;
if (m_evolscheme==2 || m_evolscheme==3) kt2=kt2-mi2-mj2;
if (m_evolscheme==0 || m_evolscheme==2) {
return 1.0/(1.0-(kt2/(z*(1.0-z))+mi2*(1.0-z)/z+mj2*z/(1.0-z))/(Q2-ma2-mi2-mj2));
}
if (fla.IsFermion()) {
if (flc.IsFermion()) return 1.0/(1.0-kt2/z/(Q2-ma2-mi2-mj2));
return 1.0/(1.0-kt2/(1.0-z)/(Q2-ma2-mi2-mj2));
else if (m_evolscheme==1 || m_evolscheme==3) {
if (fla.IsFermion()) {
if (flc.IsFermion()) return 1.0/(1.0-kt2/z/(Q2-ma2-mi2-mj2));
return 1.0/(1.0-kt2/(1.0-z)/(Q2-ma2-mi2-mj2));
}
if (flc.IsFermion()) return 1.0/(1.0-kt2/(Q2-ma2-mi2-mj2));
return 1.0/(1.0-kt2/(z*(1.0-z))/(Q2-ma2-mi2-mj2));
}
if (flc.IsFermion()) return 1.0/(1.0-kt2/(Q2-ma2-mi2-mj2));
return 1.0/(1.0-kt2/(z*(1.0-z))/(Q2-ma2-mi2-mj2));
else THROW(fatal_error, "Not implemented");
return 0.0;
}
int Kinematics_FI::MakeKinematics
......@@ -135,25 +162,38 @@ double Kinematics_IF::GetKT2(const double &Q2,const double &y,const double &z,
const ATOOLS::Flavour &flb,const ATOOLS::Flavour &flc) const
{
double pipj=(Q2-ma2-mi2-mk2)*y/z;
if (m_evolscheme==0) {
return -pipj*(1.0-z)-mi2-sqr(1.0-z)*ma2;
if (m_evolscheme==0 || m_evolscheme==2) {
double kt2=-pipj*(1.0-z)-mi2-sqr(1.0-z)*ma2;
if (m_evolscheme==0) return kt2;
if (m_evolscheme==2) return kt2+mi2+ma2;
}
double kt2=-pipj*(1.0-z);
if (flc.IsFermion()) kt2=-pipj;
return kt2;
else if (m_evolscheme==1 || m_evolscheme==3) {
double kt2=-pipj*(1.0-z);
if (flc.IsFermion()) kt2=-pipj;
if (m_evolscheme==1) return kt2;
if (m_evolscheme==3) return kt2+mi2+ma2;
}
else THROW(fatal_error, "Not implemented");
return 0.0;
}
double Kinematics_IF::GetY(const double &Q2,const double &kt2,const double &z,
double Kinematics_IF::GetY(const double &Q2,const double &_kt2,const double &z,
const double &ma2,const double &mi2,const double &mk2,
const ATOOLS::Flavour &flb,const ATOOLS::Flavour &flc,
const bool force) const
{
if (!force && (z<=0.0 || z>=1.0 || Q2>=ma2+mi2+mk2)) return -1.0;
if (m_evolscheme==0) {
double kt2=_kt2;
if (m_evolscheme==2 || m_evolscheme==3) kt2=kt2-mi2-ma2;
if (m_evolscheme==0 || m_evolscheme==2) {
return -z/(Q2-ma2-mi2-mk2)*((kt2+mi2)/(1.0-z)+(1.0-z)*ma2);
}
if (flc.IsFermion()) return -z/(Q2-ma2-mi2-mk2)*kt2;
return -z/(Q2-ma2-mi2-mk2)*kt2/(1.0-z);
else if (m_evolscheme==1 || m_evolscheme==3) {
if (flc.IsFermion()) return -z/(Q2-ma2-mi2-mk2)*kt2;
return -z/(Q2-ma2-mi2-mk2)*kt2/(1.0-z);
}
else THROW(fatal_error, "Not implemented");
return 0.0;
}
int Kinematics_IF::MakeKinematics
......@@ -201,25 +241,38 @@ double Kinematics_II::GetKT2(const double &Q2,const double &y,const double &z,
const ATOOLS::Flavour &flb,const ATOOLS::Flavour &flc) const
{
double pipj=(Q2-ma2-mi2-mb2)*y/z;
if (m_evolscheme==0) {
return pipj*(1.0-z)-mi2-sqr(1.0-z)*ma2;
if (m_evolscheme==0 || m_evolscheme==2) {
double kt2=pipj*(1.0-z)-mi2-sqr(1.0-z)*ma2;
if (m_evolscheme==0) return kt2;
if (m_evolscheme==2) return kt2+mi2+ma2;
}
double kt2=pipj*(1.0-z);
if (flc.IsFermion()) kt2=pipj;
return kt2;
else if (m_evolscheme==1 || m_evolscheme==3) {
double kt2=pipj*(1.0-z);
if (flc.IsFermion()) kt2=pipj;
if (m_evolscheme==1) return kt2;
if (m_evolscheme==3) return kt2+mi2+ma2;
}
else THROW(fatal_error, "Not implemented");
return 0.0;
}
double Kinematics_II::GetY(const double &Q2,const double &kt2,const double &z,
double Kinematics_II::GetY(const double &Q2,const double &_kt2,const double &z,
const double &ma2,const double &mi2,const double &mb2,
const ATOOLS::Flavour &flb,const ATOOLS::Flavour &flc,
const bool force) const
{
if (!force && (z<=0.0 || z>=1.0 || Q2<=ma2+mi2+mb2)) return -1.0;
if (m_evolscheme==0) {
double kt2=_kt2;
if (m_evolscheme==2 || m_evolscheme==3) kt2=kt2-mi2-ma2;
if (m_evolscheme==0 || m_evolscheme==2) {
return z/(Q2-ma2-mb2-mi2)*((kt2+mi2)/(1.0-z)+(1.0-z)*ma2);
}
if (flc.IsFermion()) return z/(Q2-ma2-mb2-mi2)*kt2;
return z/(Q2-ma2-mb2-mi2)*kt2/(1.0-z);
else if (m_evolscheme==1 || m_evolscheme==3) {
if (flc.IsFermion()) return z/(Q2-ma2-mb2-mi2)*kt2;
return z/(Q2-ma2-mb2-mi2)*kt2/(1.0-z);
}
else THROW(fatal_error, "Not implemented");
return 0.0;
}
int Kinematics_II::MakeKinematics
......
......@@ -109,7 +109,9 @@ double counting, this has to be disabled as explained in that section.
The evolution variable of the CS shower can be changed using
@option{CSS_EVOLUTION_SCHEME}. Two options are currently implemented,
which correspond to transverse momentum ordering (option 0) and modified
transverse momentum ordering (option 1). The scale at which the strong
transverse momentum ordering (option 1). In addition, modified versions of these
options (option 2 and option 3) are implemented, which take parton masses into
account where applicable. The scale at which the strong
coupling for gluon splitting into quarks is evaluated can be chosen
with @option{CSS_SCALE_SCHEME}, where 0 (default) corresponds to the ordering
parameter and 1 corresponds to invariant mass. Additionally, the CS shower
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment