Commit 0e0f1e61 authored by Jose Manuel Morgado Chávez's avatar Jose Manuel Morgado Chávez
Browse files

Two possibilities: analytical and RT solution

parent df1be303
......@@ -52,8 +52,8 @@ algebraicToyModel::algebraicToyModel(const std::string &className) : PARTONS::GP
algebraicToyModel::algebraicToyModel(const algebraicToyModel& other) : PARTONS::GPDModule(other)
{
m2 = 0.101; // Mass parameter algebraic toy model. Eq (30) Physics Letters B 780 (2018) 287–293.
m2D = 0.099; // D-term t-dependence: Fitting of Phys. Rev. D 97, 014020 (2018) gravitational FFs.
m2 = pow(0.318,2.); // Mass parameter algebraic toy model. Eq (30) Physics Letters B 780 (2018) 287–293.
m2D = pow(0.318,2.); // D-term t-dependence: Fitting of Phys. Rev. D 97, 014020 (2018) gravitational FFs.
// RT.init(); // Initialize Radon transform
}
......@@ -103,64 +103,129 @@ PARTONS::PartonDistribution algebraicToyModel::computeH()
QuarkDistribution quarkDistributionUp(QuarkFlavor::UP);
QuarkDistribution quarkDistributionDown(QuarkFlavor::DOWN);
QuarkDistribution quarkDistributionStrange(QuarkFlavor::STRANGE);
QuarkDistribution quarkDistributionCharm(QuarkFlavor::CHARM);
QuarkDistribution quarkDistributionBottom(QuarkFlavor::BOTTOM);
QuarkDistribution quarkDistributionTop(QuarkFlavor::TOP);
// u-quark
double uVal, uValM;
double uVal1, uValM1;
if ( m_t == 0) //! Zero momentum transfer
{
if ( m_x > m_xi ) // DGLAP>
if ( m_x > m_xi || m_x == m_xi ) // DGLAP>
{
uVal = 30 * pow(1 - m_x, 2.) * ( pow(m_x,2.) - pow(m_xi,2.) ) / pow( 1 - pow(m_xi,2.) , 2.);
uValM = 0.;
} else if ( m_x < -m_xi ) //DGLAP>
uVal1 = uVal;
uValM1 = uValM;
} else if ( m_x < -m_xi || m_x == -m_xi ) //DGLAP>
{
uVal = 0.;
uValM = 30 * pow(1 + m_x, 2.) * ( pow(m_x,2.) - pow(m_xi,2.) ) / pow( 1 - pow(m_xi,2.) , 2.);
uValM = 30 * pow(1 + m_x, 2.) * ( pow(m_x,2.) - pow(m_xi,2.) ) / pow( 1 - pow(m_xi,2.) , 2.);
uVal1 = uVal;
uValM1 = uValM;
} else // ERBL
{
// ============================================================================================
// Compute ERBL GPD (Proper computation: RT)
// ============================================================================================
// Compute double distribution.
if ( DDt0.isZero() )
{
Eigen::VectorXd GPD_DGLAP(RT.x.size());
for ( int i = 0; i < RT.x.size(); i ++ )
{
GPD_DGLAP(i) = 30 * pow(1 - RT.x.at(i), 2.) * ( pow(RT.x.at(i),2.) - pow(RT.xi.at(i),2.) ) / pow( 1 - pow(RT.xi.at(i),2.) , 2.);
}
// // Compute double distribution.
// if ( DDt0.isZero() )
// {
// Eigen::VectorXd GPD_DGLAP(RT.x.size());
// for ( int i = 0; i < RT.x.size(); i ++ )
// {
// GPD_DGLAP(i) = 30 * pow(1 - RT.x.at(i), 2.) * ( pow(RT.x.at(i),2.) - pow(RT.xi.at(i),2.) ) / pow( 1 - pow(RT.xi.at(i),2.) , 2.);
// }
DDt0 = RT.computeDD( GPD_DGLAP );
}
// DDt0 = RT.computeDD( GPD_DGLAP );
// }
// // Compute "gauged" GPD.
// uVal = RT.computeGPD( DDt0, m_x, m_xi );
// uValM = RT.computeGPD( DDt0, -m_x, m_xi );
// // Compute Dterm contribution.
// uVal += RT.computeDterm( DDt0, m_x, m_xi );
// uValM += RT.computeDterm( DDt0, -m_x, m_xi );
// ============================================================================================
// Compute ERBL GPD (Analytic computation)
// ============================================================================================
// Compute "gauged" GPD.
uVal = RT.computeGPD( DDt0, m_x, m_xi );
uValM = RT.computeGPD( DDt0, -m_x, m_xi );
// (Gauged) ERBL GPD t = 0
uVal1 = 7.5 * (1 - m_x) * ( pow(m_xi, 2.) - pow(m_x, 2.) ) * (m_x + 2*m_x*m_xi + pow(m_xi, 2.)) / ( pow(m_xi, 3.)*pow(1 + m_xi, 2.) );
uValM1 = 7.5 * (1 + m_x) * ( pow(m_xi, 2.) - pow(m_x, 2.) ) * (-m_x - 2*m_x*m_xi + pow(m_xi, 2.)) / ( pow(m_xi, 3.)*pow(1 + m_xi, 2.) );
// Compute Dterm contribution.
uVal += RT.computeDterm( DDt0, m_x, m_xi );
uValM += RT.computeDterm( DDt0, -m_x, m_xi );
// D-terms contribution
alpha1 = m_x/m_xi;
dplus1 = 1.125*(1-pow(alpha1,2))*(5*pow(alpha1,2)-1);
dminus1 = -3.75*alpha1*(1-pow(alpha1,2));
// Add D-terms to GPD.
dplus1 /= m_xi;
if ( m_xi >= 0 ) // Conditional expression taking into acount the factor sign(\xi) accompanying dminus.
{
uVal1 += dplus1 + dminus1;
uValM1 += dplus1 - dminus1;
} else
{
uVal1 += dplus1 - dminus1;
uValM1 += dplus1 + dminus1;
}
}
} else //! Non-vanishing momentum transfer. CHECK!!!!!!!!!!!!!!!!!!!
{
c = -m_t*pow(1 - m_x, 2.)/(4*m2*(1 - pow(m_xi,2))); // t-dependence algebraic toy model.
cM = -m_t*pow(1 + m_x, 2.)/(4*m2*(1 - pow(m_xi,2)));
c1 = -m_t/(4*m2);
dt = 1/(1-0.25*m_t/m2D); // D-term t-dependence (monopole parametrization).
if ( m_x > m_xi ) // DGLAP>
if ( m_x > m_xi || m_x == m_xi ) // DGLAP>
{
uVal = 7.5 * pow(1 - m_x, 2.) * ( pow(m_x,2.) - pow(m_xi,2.) ) * (3 + ((1 - 2 * c) * atanh(sqrt(c/(1+c))))/((1 + c) * sqrt(c/(1 + c))) )
/ ( pow( 1 - pow(m_xi,2.) , 2.) * pow(1 + c,2.) );
uValM = 0.;
} else if ( m_x < -m_xi ) // DGLAP>
uVal1 = uVal;
uValM1 = uValM;
if ( m_x == 1 ) // Actually this is the limit x->1 (with \xi<1).
{
uVal = 0.;
uValM = 0.;
uVal1 = uVal;
uValM1 = uValM;
}
} else if ( m_x < -m_xi || m_x == -m_xi ) // DGLAP<
{
uVal = 0.;
uValM = 7.5 * pow(1 + m_x, 2.) * ( pow(m_x,2.) - pow(m_xi,2.) ) * (3 + ((1 - 2 * cM) * atanh(sqrt(cM/(1+cM))))/((1 + cM) * sqrt(cM/(1 + cM))) )
/ ( pow( 1 - pow(m_xi,2.) , 2.) * pow(1 + cM,2.) );
uVal1 = uVal;
uValM1 = uValM;
if ( m_x == -1 ) // Actually this is the limit x->1 (with \xi<1).
{
uVal = 0.;
uValM = 0.;
uVal1 = uVal;
uValM1 = uValM;
}
} else // ERBL
{
// ============================================================================================
......@@ -168,65 +233,116 @@ PARTONS::PartonDistribution algebraicToyModel::computeH()
// ============================================================================================
// Compute double distribution.
if ( DD.isZero() ) // TODO: Map with DDs for different t.
{
Eigen::VectorXd GPD_DGLAP(RT.x.size());
// if ( DD.isZero() ) // TODO: Map with DDs for different t.
// {
// Eigen::VectorXd GPD_DGLAP(RT.x.size());
for ( int i = 0; i < RT.x.size(); i ++ )
{
ca = -m_t*pow(1 - RT.x.at(i), 2.)/(4*m2*(1 - pow(RT.xi.at(i),2)));
// for ( int i = 0; i < RT.x.size(); i ++ )
// {
// ca = -m_t*pow(1 - RT.x.at(i), 2.)/(4*m2*(1 - pow(RT.xi.at(i),2)));
GPD_DGLAP(i) = 7.5 * pow(1 - RT.x.at(i), 2.) * ( pow(RT.x.at(i),2.) - pow(RT.xi.at(i),2.) ) * (3 + ((1 - 2 * ca) * atanh(sqrt(ca/(1+ca))))/((1 + ca) * sqrt(ca/(1 + ca))) )
/ ( pow( 1 - pow(RT.xi.at(i),2.) , 2.) * pow(1 + ca,2.) );
}
// GPD_DGLAP(i) = 7.5 * pow(1 - RT.x.at(i), 2.) * ( pow(RT.x.at(i),2.) - pow(RT.xi.at(i),2.) ) * (3 + ((1 - 2 * ca) * atanh(sqrt(ca/(1+ca))))/((1 + ca) * sqrt(ca/(1 + ca))) )
// / ( pow( 1 - pow(RT.xi.at(i),2.) , 2.) * pow(1 + ca,2.) );
// }
DD = RT.computeDD( GPD_DGLAP );
}
// DD = RT.computeDD( GPD_DGLAP );
// }
// // Compute "gauged" GPD.
// uVal = RT.computeGPD( DD, m_x, m_xi );
// uValM = RT.computeGPD( DD, -m_x, m_xi );
// Compute "gauged" GPD.
uVal = RT.computeGPD( DD, m_x, m_xi );
uValM = RT.computeGPD( DD, -m_x, m_xi );
// // Compute Dterm contribution.
// if ( DDt0.isZero() ) // TODO: Overload computeDterm so that it accepts (DDt0, x, xi) and (x, xi) as arguments. In that way, we can here look for DDt0 and choose the function computeDterm to be called.
// {
// Eigen::VectorXd GPD_DGLAP(RT.x.size());
// for ( int i = 0; i < RT.x.size(); i ++ )
// {
// GPD_DGLAP(i) = 30 * pow(1 - RT.x.at(i), 2.) * ( pow(RT.x.at(i),2.) - pow(RT.xi.at(i),2.) ) / pow( 1 - pow(RT.xi.at(i),2.) , 2.);
// }
// DDt0 = RT.computeDD( GPD_DGLAP );
// }
// uVal += dt*RT.computeDterm( DDt0, m_x, m_xi );
// uValM += dt*RT.computeDterm( DDt0, -m_x, m_xi );
// Compute Dterm contribution.
if ( DDt0.isZero() ) // TODO: Overload computeDterm so that it accepts (DDt0, x, xi) and (x, xi) as arguments. In that way, we can here look for DDt0 and choose the function computeDterm to be called.
// ============================================================================================
// Compute ERBL GPD (Analytic computation)
// ============================================================================================
if ( m_xi == 1 )
{
Eigen::VectorXd GPD_DGLAP(RT.x.size());
uVal1 = -3.75*(1-pow(m_x,2))*( (((1-pow(m_x,2))+c1*(1-4*pow(m_x,3)+3*pow(m_x,4)))/pow(1+c1*(1-pow(m_x,2)),2)) - log(1+c1*(1-pow(m_x,2)))/c1 )
/ (c1*pow(1-m_x,2));
for ( int i = 0; i < RT.x.size(); i ++ )
uValM1 = -3.75*(1-pow(m_x,2))*( (((1-pow(m_x,2))+c1*(1+4*pow(m_x,3)+3*pow(m_x,4)))/pow(1+c1*(1-pow(m_x,2)),2)) - log(1+c1*(1-pow(m_x,2)))/c1 )
/ (c1*pow(1+m_x,2));
if ( m_x == 1 )
{
GPD_DGLAP(i) = 30 * pow(1 - RT.x.at(i), 2.) * ( pow(RT.x.at(i),2.) - pow(RT.xi.at(i),2.) ) / pow( 1 - pow(RT.xi.at(i),2.) , 2.);
uVal1 = 0.;
uValM1 = 0.; // Verify this limit.
}
DDt0 = RT.computeDD( GPD_DGLAP );
} else // The following expression shows a "fake" divergence in the limit \xi->1.
{
// (Gauged) ERBL GPD t < 0
uVal1 = -3.75 * (pow(m_xi,2) - pow(m_x,2)) * (sqrt(c * (c + 1)) * ( -c*m_xi*(1 - pow(m_xi,2))*(1 - m_x)
* (pow(m_xi,4) + 6*m_xi*(1-m_x)*pow(m_x,2)-6*pow(m_xi,3)*(1-m_x)+pow(m_xi,2)*(4-3*(3-m_x)*m_x)+m_x*(4-m_x*(8-5*m_x)))
+ pow(c,2)*pow(1-pow(m_xi,2),2)*(pow(m_xi,3)*(3*m_xi-2)+3*pow(m_x,4)-4*m_xi*pow(m_x,3)-6*(m_xi-1)*m_xi*pow(m_x,2)+2*m_xi*(pow(m_xi,2)-1)*m_x)
+ pow(m_x-1,3)*(pow(m_xi,5)+3*pow(m_xi,4)*(m_x-1)+pow(m_xi,3)*(2-5*m_x)+2*m_xi*m_x))
+ (1-2*c)*pow( c*(1-pow(m_xi,2))*(pow(m_xi,2)-pow(m_x,2))+pow(m_xi,2)*pow((m_x-1),2),2)
*(atanh(sqrt(c/(1+c)))-atanh(sqrt(c/(1+c))*(pow(m_xi,2)-m_x)/(m_xi*(1-m_x)))) )
/ (pow(1+c,2.5)*pow(1-pow(m_xi,2),1.5)*pow(1-m_x,2)*sqrt(c*(1-pow(m_xi,2)))*pow(pow(m_xi,2)+c*(1-pow(m_xi,2))*(pow(m_xi,2)-pow(m_x,2))/pow(1-m_x,2),2));
uValM1 = -3.75 * (pow(m_xi,2) - pow(m_x,2)) * ( sqrt(cM * (cM + 1)) * ( -cM*m_xi*(1 - pow(m_xi,2))*(1 + m_x)
* (pow(m_xi,4) + 6*m_xi*(1+m_x)*pow(m_x,2)-6*pow(m_xi,3)*(1+m_x)+pow(m_xi,2)*(4+3*(3+m_x)*m_x)-m_x*(4+m_x*(8+5*m_x)))
+ pow(cM,2)*pow(1-pow(m_xi,2),2)*(pow(m_xi,3)*(3*m_xi-2)+3*pow(m_x,4)+4*m_xi*pow(m_x,3)-6*(m_xi-1)*m_xi*pow(m_x,2)-2*m_xi*(pow(m_xi,2)-1)*m_x)
- pow(m_x+1,3)*(pow(m_xi,5)-3*pow(m_xi,4)*(m_x+1)+pow(m_xi,3)*(2+5*m_x)-2*m_xi*m_x))
+ (1-2*cM)*pow( cM*(1-pow(m_xi,2))*(pow(m_xi,2)-pow(m_x,2))+pow(m_xi,2)*pow(m_x+1,2),2)
*(atanh(sqrt(cM/(1+cM)))-atanh(sqrt(cM/(1+cM))*(pow(m_xi,2)+m_x)/(m_xi*(1+m_x)))))
/ (pow(1+cM,2.5)*pow(1-pow(m_xi,2),1.5)*pow(1+m_x,2)*sqrt(cM*(1-pow(m_xi,2)))*pow(pow(m_xi,2)+cM*(1-pow(m_xi,2))*(pow(m_xi,2)-pow(m_x,2))/pow(1+m_x,2),2));
}
// D-terms contribution (monopole parametrization)
alpha1 = m_x/m_xi;
dplus1 = 1.125*(1-pow(alpha1,2))*(5*pow(alpha1,2)-1)*dt;
dminus1 = -3.75*alpha1*(1-pow(alpha1,2))*dt;
// Add D-terms to GPD.
dplus1 /= m_xi;
if ( m_xi >= 0 ) // Conditional expression taking into acount the factor sign(\xi) accompanying dminus.
{
uVal1 += dplus1 + dminus1;
uValM1 += dplus1 - dminus1;
} else
{
uVal1 += dplus1 - dminus1;
uValM1 += dplus1 + dminus1;
}
uVal += dt*RT.computeDterm( DDt0, m_x, m_xi );
uValM += dt*RT.computeDterm( DDt0, -m_x, m_xi );
}
}
double uSea = 0.;
double uSeaM = 0.;
quarkDistributionUp.setQuarkDistribution(uVal + uSea);
// Singlet distribution
quarkDistributionUp.setQuarkDistributionPlus(uVal + uSea - uValM - uSeaM);
// Nonsinglet distributiion
quarkDistributionUp.setQuarkDistributionMinus(uVal + uValM);
quarkDistributionUp.setQuarkDistribution(uVal1 + uSea);
quarkDistributionUp.setQuarkDistributionPlus(uVal1 + uSea - uValM1 - uSeaM);
quarkDistributionUp.setQuarkDistributionMinus(uVal1 + uValM1);
// d-quark
double dVal = -uValM;
double dValM = -uVal;
double dVal = -uValM1;
double dValM = -uVal1;
double dSea = 0.;
double dSeaM = 0.;
quarkDistributionDown.setQuarkDistribution(dVal + dSea);
// Singlet distribution
quarkDistributionDown.setQuarkDistributionPlus(dVal + dSea - dValM - dSeaM);
// Nonsinglet distributiion
quarkDistributionDown.setQuarkDistributionMinus(dVal + dValM);
quarkDistributionDown.setQuarkDistributionPlus(dVal + dSea - dValM - dSeaM);
quarkDistributionDown.setQuarkDistributionMinus(dVal + dValM);
// s-quark
double sVal = 0.;
......@@ -238,9 +354,42 @@ PARTONS::PartonDistribution algebraicToyModel::computeH()
quarkDistributionStrange.setQuarkDistributionPlus(sVal + sSea - sValM - sSeaM);
quarkDistributionStrange.setQuarkDistributionMinus(sVal + sValM);
// c-quark
double cVal = 0.;
double cValM = 0.;
double cSea = 0.;
double cSeaM = 0.;
quarkDistributionCharm.setQuarkDistribution(cVal + cSea);
quarkDistributionCharm.setQuarkDistributionPlus(cVal + cSea - cValM - cSeaM);
quarkDistributionCharm.setQuarkDistributionMinus(cVal + cValM);
// b-quark
double bVal = 0.;
double bValM = 0.;
double bSea = 0.;
double bSeaM = 0.;
quarkDistributionBottom.setQuarkDistribution(bVal + bSea);
quarkDistributionBottom.setQuarkDistributionPlus(bVal + bSea - bValM - bSeaM);
quarkDistributionBottom.setQuarkDistributionMinus(bVal + bValM);
// t-quark
double tVal = 0.;
double tValM = 0.;
double tSea = 0.;
double tSeaM = 0.;
quarkDistributionTop.setQuarkDistribution(tVal + tSea);
quarkDistributionTop.setQuarkDistributionPlus(tVal + tSea - tValM - tSeaM);
quarkDistributionTop.setQuarkDistributionMinus(tVal + tValM);
partonDistribution.addQuarkDistribution(quarkDistributionUp);
partonDistribution.addQuarkDistribution(quarkDistributionDown);
partonDistribution.addQuarkDistribution(quarkDistributionStrange);
partonDistribution.addQuarkDistribution(quarkDistributionCharm);
partonDistribution.addQuarkDistribution(quarkDistributionBottom);
partonDistribution.addQuarkDistribution(quarkDistributionTop);
// Gluon distributions
GluonDistribution gluonDistribution(0.);
......
/**
* @file RunningAlphaStrongPI.cpp
* @author Jose Manuel Morgado Chávez (University of Huelva),
* @date January 18th 2021
* @date May, 10th 2021
*
* @brief Evaluation of the process independent strong running coupling constant presented in e.g.: Eur. Phys. J. C (2020) 80:1064
* @version 1.0
......@@ -129,6 +129,9 @@ double RunningAlphaStrongPI::compute() {
Running(m_Mu, Lambda, m_nf);
if (m_alphaS > 1.0)
warn(__func__, ElemUtils::Formatter() << "alpha_s = " << m_alphaS);
return fAlphaS;
}
......
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