diff --git a/CepGen/Physics/FormFactors.cpp b/CepGen/Physics/FormFactors.cpp
index 17c6d4d..e27b879 100644
--- a/CepGen/Physics/FormFactors.cpp
+++ b/CepGen/Physics/FormFactors.cpp
@@ -1,71 +1,74 @@
#include "FormFactors.h"
namespace CepGen
{
FormFactors
FormFactors::Trivial()
{
return FormFactors( 1.0, 1.0 );
}
FormFactors
FormFactors::ProtonElastic( double q2 )
{
const double mp = Particle::massFromPDGId( Particle::Proton ), mp2 = mp*mp;
const double GE = pow( 1.+q2/0.71, -2. ), GE2 = GE*GE;
const double GM = 2.79*GE, GM2 = GM*GM;
return FormFactors( ( 4.*mp2*GE2 + q2*GM2 ) / ( 4.*mp2 + q2 ), GM2 );
}
FormFactors
FormFactors::ProtonInelastic( const StructureFunctionsType& sf, double q2, double mi2, double mf2 )
{
switch ( sf ) {
case StructureFunctionsType::ElasticProton:
InWarning( "Elastic proton form factors requested! Check your process definition!" );
return FormFactors::ProtonElastic( q2 );
case StructureFunctionsType::SuriYennie:
return FormFactors::SuriYennie( q2, mi2, mf2 );
case StructureFunctionsType::SzczurekUleshchenko:
return FormFactors::SzczurekUleshchenko( q2, mi2, mf2 );
case StructureFunctionsType::Fiore:
case StructureFunctionsType::FioreSea:
case StructureFunctionsType::FioreVal:
return FormFactors::FioreBrasse( q2, mi2, mf2 );
default: throw Exception( __PRETTY_FUNCTION__, "Invalid structure functions required!", FatalError );
}
}
FormFactors
FormFactors::SuriYennie( double q2, double mi2, double mf2 )
{
const double x = q2 / ( q2 + mf2 - mi2 );
- const StructureFunctions sy = SF::SuriYennie( q2, x );
+ SF::SuriYennie suriyennie;
+ const StructureFunctions sy = suriyennie( q2, x );
//std::cout << "---> " << sy.FM << "\t" << sy.F2*x/q2 << "\t" << sy.F2*x*sqrt(mi2)/q2 << std::endl;
return FormFactors( sy.F2 * x * sqrt( mi2 ) / q2, sy.FM ); //FIXME
}
FormFactors
FormFactors::FioreBrasse( double q2, double mi2, double mf2 )
{
const double x = q2 / ( q2 + mf2 - mi2 );
- StructureFunctions sf = SF::FioreBrasse( q2, x );
+ SF::FioreBrasse fb;
+ StructureFunctions sf = fb( q2, x );
return FormFactors( sf.F2 * x / q2, -2.*sf.F1 / q2 );
}
FormFactors
FormFactors::SzczurekUleshchenko( double q2, double mi2, double mf2 )
{
const double x = q2 / ( q2 + mf2 - mi2 );
- StructureFunctions sf = SF::SzczurekUleshchenko( q2, x );
+ SF::SzczurekUleshchenko su;
+ StructureFunctions sf = su( q2, x );
return FormFactors( sf.F2 * x / q2, -2.*sf.F1 / q2 );
}
std::ostream&
operator<<( std::ostream& os, const FormFactors& ff )
{
os << Form( "Form factors: electric: Fe = %.3e ; magnetic: Fm = %.3e", ff.FE, ff.FM ).c_str();
return os;
}
}
diff --git a/CepGen/Processes/GenericKTProcess.cpp b/CepGen/Processes/GenericKTProcess.cpp
index c504d47..dabab26 100644
--- a/CepGen/Processes/GenericKTProcess.cpp
+++ b/CepGen/Processes/GenericKTProcess.cpp
@@ -1,274 +1,279 @@
#include "GenericKTProcess.h"
namespace CepGen
{
namespace Process
{
GenericKTProcess::GenericKTProcess( const std::string& name,
const std::string& description,
const unsigned int& num_user_dimensions,
const std::array<Particle::ParticleCode,2>& partons,
const std::vector<Particle::ParticleCode>& central ) :
GenericProcess( name, description+" (kT-factorisation approach)" ),
kNumUserDimensions( num_user_dimensions ),
kIntermediateParts( partons ), kProducedParts( central )
{}
GenericKTProcess::~GenericKTProcess()
{}
void
GenericKTProcess::addEventContent()
{
GenericProcess::setEventContent(
{ // incoming state
{ Particle::IncomingBeam1, Particle::Proton },
{ Particle::IncomingBeam2, Particle::Proton },
{ Particle::Parton1, kIntermediateParts[0] },
{ Particle::Parton2, kIntermediateParts[1] }
},
{ // outgoing state
{ Particle::OutgoingBeam1, { Particle::Proton } },
{ Particle::OutgoingBeam2, { Particle::Proton } },
{ Particle::CentralSystem, kProducedParts }
}
);
setExtraContent();
}
unsigned int
GenericKTProcess::numDimensions( const Kinematics::ProcessMode& process_mode ) const
{
switch ( process_mode ) {
default:
case Kinematics::ElasticElastic: return kNumRequiredDimensions+kNumUserDimensions;
case Kinematics::ElasticInelastic:
case Kinematics::InelasticElastic: return kNumRequiredDimensions+kNumUserDimensions+1;
case Kinematics::InelasticInelastic: return kNumRequiredDimensions+kNumUserDimensions+2;
}
}
void
GenericKTProcess::addPartonContent()
{
// Incoming partons
qt1_ = exp( log_qmin_+( log_qmax_-log_qmin_ )*x( 0 ) );
qt2_ = exp( log_qmin_+( log_qmax_-log_qmin_ )*x( 1 ) );
phi_qt1_ = 2.*M_PI*x( 2 );
phi_qt2_ = 2.*M_PI*x( 3 );
DebuggingInsideLoop( Form( "photons transverse virtualities (qt):\n\t"
" mag = %f / %f (%.2f < log(qt) < %.2f)\n\t"
" phi = %f / %f",
qt1_, qt2_, log_qmin_, log_qmax_, phi_qt1_, phi_qt2_ ) );
}
void
GenericKTProcess::setKinematics( const Kinematics& kin )
{
cuts_ = kin;
log_qmin_ = -10.; // FIXME //log_qmin_ = std::log( std::sqrt( cuts_.q2min ) );
log_qmax_ = log( cuts_.initial_cuts[Cuts::qt].max() );
}
double
GenericKTProcess::computeWeight()
{
addPartonContent();
prepareKTKinematics();
computeOutgoingPrimaryParticlesMasses();
const double jac = computeJacobian(),
integrand = computeKTFactorisedMatrixElement(),
weight = jac*integrand;
DebuggingInsideLoop( Form( "Jacobian = %f\n\tIntegrand = %f\n\tdW = %f", jac, integrand, weight ) );
return weight;
}
void
GenericKTProcess::computeOutgoingPrimaryParticlesMasses()
{
const unsigned int op_index = kNumRequiredDimensions+kNumUserDimensions;
const Kinematics::Limits remn_mx_cuts = cuts_.remnant_cuts[Cuts::mass];
switch ( cuts_.mode ) {
case Kinematics::ElectronProton: default: {
InError( "This kT factorisation process is intended for p-on-p collisions! Aborting!" );
exit( 0 ); } break;
case Kinematics::ElasticElastic: {
MX_ = event_->getOneByRole( Particle::IncomingBeam1 ).mass();
MY_ = event_->getOneByRole( Particle::IncomingBeam2 ).mass();
} break;
case Kinematics::ElasticInelastic: {
const double mx_min = remn_mx_cuts.min(), mx_range = remn_mx_cuts.range();
MX_ = event_->getOneByRole( Particle::IncomingBeam1 ).mass();
MY_ = mx_min + mx_range*x( op_index );
} break;
case Kinematics::InelasticElastic: {
const double mx_min = remn_mx_cuts.min(), mx_range = remn_mx_cuts.range();
MX_ = mx_min + mx_range*x( op_index );
MY_ = event_->getOneByRole( Particle::IncomingBeam2 ).mass();
} break;
case Kinematics::InelasticInelastic: {
const double mx_min = remn_mx_cuts.min(), mx_range = remn_mx_cuts.range();
MX_ = mx_min + mx_range*x( op_index );
MY_ = mx_min + mx_range*x( op_index+1 );
} break;
}
DebuggingInsideLoop( Form( "outgoing remnants invariant mass: %f / %f (%.2f < M(X/Y) < %.2f)", MX_, MY_, remn_mx_cuts.min(), remn_mx_cuts.max() ) );
}
void
GenericKTProcess::computeIncomingFluxes( double x1, double q1t2, double x2, double q2t2 )
{
flux1_ = flux2_ = 0.;
switch ( cuts_.mode ) {
case Kinematics::ElasticElastic:
flux1_ = elasticFlux( x1, q1t2 );
flux2_ = elasticFlux( x2, q2t2 );
break;
case Kinematics::ElasticInelastic:
flux1_ = elasticFlux( x1, q1t2 );
flux2_ = inelasticFlux( x2, q2t2, MY_ );
break;
case Kinematics::InelasticElastic:
flux1_ = inelasticFlux( x1, q1t2, MX_ );
flux2_ = elasticFlux( x2, q2t2 );
break;
case Kinematics::InelasticInelastic:
flux1_ = inelasticFlux( x1, q1t2, MX_ );
flux2_ = inelasticFlux( x2, q2t2, MY_ );
break;
default: return;
}
flux1_ = std::max( flux1_, 1.e-20 );
flux2_ = std::max( flux2_, 1.e-20 );
DebuggingInsideLoop( Form( "Form factors: %e / %e", flux1_, flux2_ ) );
}
void
GenericKTProcess::fillKinematics( bool )
{
fillPrimaryParticlesKinematics();
fillCentralParticlesKinematics();
}
void
GenericKTProcess::fillPrimaryParticlesKinematics()
{
//=================================================================
// outgoing protons
//=================================================================
Particle& op1 = event_->getOneByRole( Particle::OutgoingBeam1 ),
&op2 = event_->getOneByRole( Particle::OutgoingBeam2 );
op1.setMomentum( PX_ );
op2.setMomentum( PY_ );
switch ( cuts_.mode ) {
case Kinematics::ElasticElastic:
op1.setStatus( Particle::FinalState );
op2.setStatus( Particle::FinalState );
break;
case Kinematics::ElasticInelastic:
op1.setStatus( Particle::FinalState );
op2.setStatus( Particle::Undecayed ); op2.setMass( MY_ );
break;
case Kinematics::InelasticElastic:
op1.setStatus( Particle::Undecayed ); op1.setMass( MX_ );
op2.setStatus( Particle::FinalState );
break;
case Kinematics::InelasticInelastic:
op1.setStatus( Particle::Undecayed ); op1.setMass( MX_ );
op2.setStatus( Particle::Undecayed ); op2.setMass( MY_ );
break;
default: { FatalError( "This kT factorisation process is intended for p-on-p collisions! Aborting!" ); } break;
}
//=================================================================
// incoming partons (photons, pomerons, ...)
//=================================================================
//FIXME ensure the validity of this approach
Particle& g1 = event_->getOneByRole( Particle::Parton1 ),
&g2 = event_->getOneByRole( Particle::Parton2 );
g1.setMomentum( event_->getOneByRole( Particle::IncomingBeam1 ).momentum()-PX_, true );
g1.setStatus( Particle::Incoming );
g2.setMomentum( event_->getOneByRole( Particle::IncomingBeam2 ).momentum()-PY_, true );
g2.setStatus( Particle::Incoming );
}
double
GenericKTProcess::minimalJacobian() const
{
double jac = 1.;
jac *= ( log_qmax_-log_qmin_ )*qt1_; // d(q1t) . q1t
jac *= ( log_qmax_-log_qmin_ )*qt2_; // d(q2t) . q2t
jac *= 2.*M_PI; // d(phi1)
jac *= 2.*M_PI; // d(phi2)
const double mx_range = cuts_.remnant_cuts.at( Cuts::mass ).range();
switch ( cuts_.mode ) {
case Kinematics::ElasticElastic: default: break;
case Kinematics::ElasticInelastic: jac *= 2.* mx_range * MY_; break;
case Kinematics::InelasticElastic: jac *= 2.* mx_range * MX_; break;
case Kinematics::InelasticInelastic: jac *= 2.* mx_range * MX_;
jac *= 2.* mx_range * MY_; break;
} // d(mx/y**2)
return jac;
}
double
GenericKTProcess::elasticFlux( double x, double kt2 )
{
const double mp = Particle::massFromPDGId( Particle::Proton ), mp2 = mp*mp;
const double Q2_min = x*x*mp2/( 1.-x ), Q2_ela = Q2_min + kt2/( 1.-x );
const FormFactors ela = FormFactors::ProtonElastic( Q2_ela );
const double ela1 = ( 1.-x )*( 1.-Q2_min/Q2_ela );
const double ela2 = ela.FE;
//const double ela3 = 1.-( Q2_ela-kt2 )/Q2_ela;
const double f_ela = Constants::alphaEM/M_PI/kt2*( ela1*ela2 + 0.5*x*x*ela.FM );
return f_ela * ( 1.-x ) * kt2 / Q2_ela;
}
double
GenericKTProcess::inelasticFlux( double x, double kt2, double mx )
{
const double mx2 = mx*mx, mp = Particle::massFromPDGId( Particle::Proton ), mp2 = mp*mp;
// F2 structure function
const double Q2min = 1. / ( 1.-x )*( x*( mx2-mp2 ) + x*x*mp2 ),
Q2 = Q2min + kt2 / ( 1.-x );
float xbj = Q2 / ( Q2 + mx2 - mp2 );
StructureFunctions sf;
switch ( cuts_.structure_functions ) {
- case StructureFunctionsType::SzczurekUleshchenko:
- sf = SF::SzczurekUleshchenko( Q2, xbj ); break;
- case StructureFunctionsType::SuriYennie:
- sf = SF::SuriYennie( Q2, xbj ); break;
- case StructureFunctionsType::Fiore:
- sf = SF::FioreBrasse( Q2, xbj ); break;
- case StructureFunctionsType::ALLM91:
- sf = SF::ALLM( Q2, xbj, SF::ALLMParameterisation::allm91() ); break;
- case StructureFunctionsType::ALLM97:
- sf = SF::ALLM( Q2, xbj, SF::ALLMParameterisation::allm97() ); break;
+ case StructureFunctionsType::SzczurekUleshchenko: {
+ SF::SzczurekUleshchenko su;
+ sf = su( Q2, xbj ); } break;
+ case StructureFunctionsType::SuriYennie: {
+ SF::SuriYennie sy;
+ sf = sy( Q2, xbj ); } break;
+ case StructureFunctionsType::Fiore: {
+ SF::FioreBrasse fb;
+ sf = fb( Q2, xbj ); } break;
+ case StructureFunctionsType::ALLM91: {
+ SF::ALLM allm91( SF::ALLM::Parameterisation::allm91() );
+ sf = allm91( Q2, xbj ); } break;
+ case StructureFunctionsType::ALLM97: {
+ SF::ALLM allm97( SF::ALLM::Parameterisation::allm97() );
+ sf = allm97( Q2, xbj ); } break;
default: break; //FIXME
}
// Longitudinal/transverse virtual photon cross section R
// from Sibirtsev and Blunden (Phys Rev C 88,065202 (2013))
const double R = 0.014 * Q2 * ( exp( -0.07*Q2 ) + 41.*exp( -0.8*Q2 ) );
const double F1 = sf.F2 * ( 1.+4*xbj*xbj*mp2/Q2 )/( 1.+R )/( 2.*xbj );
const double ine1 = ( 1.-x )*( 1.-Q2min / Q2 );
const double f_D = sf.F2/( Q2 + mx2 - mp2 ) * ine1, f_C= 2.*F1/Q2;
const double f_ine = Constants::alphaEM/M_PI/kt2*( f_D + 0.5*x*x*f_C );
return f_ine * ( 1.-x ) * kt2 / Q2;
}
}
}
diff --git a/CepGen/StructureFunctions/ALLM.cpp b/CepGen/StructureFunctions/ALLM.cpp
index 4ff74f0..32c7f54 100644
--- a/CepGen/StructureFunctions/ALLM.cpp
+++ b/CepGen/StructureFunctions/ALLM.cpp
@@ -1,154 +1,154 @@
#include "ALLM.h"
namespace CepGen
{
namespace SF
{
- ALLMParameterisation
- ALLMParameterisation::allm91()
+ ALLM::Parameterisation
+ ALLM::Parameterisation::allm91()
{
- ALLMParameterisation p;
- p.pomeron = Parameters(
+ Parameterisation p;
+ p.pomeron = Parameters( {
{ 0.26550, 0.04856, 1.04682 },
{ -0.04503, -0.36407, 8.17091 },
- { 0.49222, 0.52116, 3.5515 } );
- p.reggeon = Parameters(
+ { 0.49222, 0.52116, 3.5515 } } );
+ p.reggeon = Parameters( {
{ 0.67639, 0.49027, 2.66275 },
{ 0.60408, 0.17353, 1.61812 },
- { 1.26066, 1.83624, 0.81141 } );
+ { 1.26066, 1.83624, 0.81141 } } );
p.m02 = 0.30508;
p.mp2 = 10.676;
p.mr2 = 0.20623;
p.q02 = 0.27799;
p.lambda2 = 0.06527;
return p;
}
- ALLMParameterisation
- ALLMParameterisation::allm97()
+ ALLM::Parameterisation
+ ALLM::Parameterisation::allm97()
{
- ALLMParameterisation p;
- p.pomeron = Parameters(
+ Parameterisation p;
+ p.pomeron = Parameters( {
{ 0.28067, 0.22291, 2.1979 },
{ -0.0808, -0.44812, 1.1709 },
- { 0.36292, 1.8917, 1.8439 } );
- p.reggeon = Parameters(
+ { 0.36292, 1.8917, 1.8439 } } );
+ p.reggeon = Parameters( {
{ 0.80107, 0.97307, 3.4924 },
{ 0.58400, 0.37888, 2.6063 },
- { 0.01147, 3.7582, 0.49338 } );
+ { 0.01147, 3.7582, 0.49338 } } );
p.m02 = 0.31985;
p.mp2 = 49.457;
p.mr2 = 0.15052;
p.q02 = 0.52544;
p.lambda2 = 0.06526;
return p;
}
- ALLMParameterisation
- ALLMParameterisation::hht_allm()
+ ALLM::Parameterisation
+ ALLM::Parameterisation::hht_allm()
{
- ALLMParameterisation p;
- p.pomeron = Parameters(
+ Parameterisation p;
+ p.pomeron = Parameters( {
{ 0.412, 0.164, 17.7 },
{ -0.835, -0.446, 10.6 },
- { -45.8, 55.7, -0.031 } );
- p.reggeon = Parameters(
+ { -45.8, 55.7, -0.031 } } );
+ p.reggeon = Parameters( {
{ -1.04, 2.97, 0.163 },
{ 0.706, 0.185, -16.4 },
- { -1.29, 4.51, 1.16 } );
+ { -1.29, 4.51, 1.16 } } );
p.m02 = 0.446;
p.mp2 = 74.2;
p.mr2 = 29.3;
p.q02 = 4.74e-5;
p.lambda2 = 2.2e-8;
return p;
}
- ALLMParameterisation
- ALLMParameterisation::hht_allm_ft()
+ ALLM::Parameterisation
+ ALLM::Parameterisation::hht_allm_ft()
{
- ALLMParameterisation p;
- p.pomeron = Parameters(
+ Parameterisation p;
+ p.pomeron = Parameters( {
{ 0.356, 0.171, 18.6 },
{ -0.075, -0.470, 9.2 },
- { -0.477, 54.0, 0.073 } );
- p.reggeon = Parameters(
+ { -0.477, 54.0, 0.073 } } );
+ p.reggeon = Parameters( {
{ -0.636, 3.37, -0.660 },
{ 0.882, 0.082, -8.5 },
- { 0.339, 3.38, 1.07 } );
+ { 0.339, 3.38, 1.07 } } );
p.m02 = 0.388;
p.mp2 = 50.8;
p.mr2 = 0.838;
p.q02 = 1.87e-5;
p.lambda2 = 4.4e-9;
return p;
}
- ALLMParameterisation
- ALLMParameterisation::gd07p()
+ ALLM::Parameterisation
+ ALLM::Parameterisation::gd07p()
{
- ALLMParameterisation p;
- p.pomeron = Parameters(
+ Parameterisation p;
+ p.pomeron = Parameters( {
{ 0.339, 0.127, 1.16 },
{ -0.105, -0.495, 1.29 },
- { -1.42, 4.51, 0.551 } );
- p.reggeon = Parameters(
+ { -1.42, 4.51, 0.551 } } );
+ p.reggeon = Parameters( {
{ 0.838, 2.36, 1.77 },
{ 0.374, 0.998, 0.775 },
- { 2.71, 1.83, 1.26 } );
+ { 2.71, 1.83, 1.26 } } );
p.m02 = 0.454;
p.mp2 = 30.7;
p.mr2 = 0.117;
p.q02 = 1.15;
p.lambda2 = 0.06527;
return p;
}
- ALLMParameterisation
- ALLMParameterisation::gd11p()
+ ALLM::Parameterisation
+ ALLM::Parameterisation::gd11p()
{
- ALLMParameterisation p;
- p.pomeron = Parameters(
+ Parameterisation p;
+ p.pomeron = Parameters( {
{ 0.3638, 0.1211, 1.166 }, // c
{ -0.11895, -0.4783, 1.353 }, // a
- { 1.0833, 2.656, 1.771 } ); // b
- p.reggeon = Parameters(
+ { 1.0833, 2.656, 1.771 } } ); // b
+ p.reggeon = Parameters( {
{ 1.3633, 2.256, 2.209 },
{ 0.3425, 1.0603, 0.5164 },
- { -10.408, 14.857, 0.07739 } );
+ { -10.408, 14.857, 0.07739 } } );
p.m02 = 0.5063;
p.mp2 = 34.75;
p.mr2 = 0.03190;
p.q02 = 1.374;
p.lambda2 = 0.06527;
return p;
}
StructureFunctions
- ALLM( double q2, double xbj, const ALLMParameterisation& param )
+ ALLM::operator()( double q2, double xbj ) const
{
- const double factor = q2/( q2+param.m02 );
+ const double factor = q2/( q2+params_.m02 );
const double W2_eff = q2*( 1.-xbj )/xbj;
- const double xp = ( q2+param.mp2 )/( q2+W2_eff+param.mp2 ),
- xr = ( q2+param.mr2 )/( q2+W2_eff+param.mr2 );
+ const double xp = ( q2+params_.mp2 )/( q2+W2_eff+params_.mp2 ),
+ xr = ( q2+params_.mr2 )/( q2+W2_eff+params_.mr2 );
- const double xlog1 = log( ( q2+param.q02 )/ param.lambda2 ), xlog2 = log( param.q02/param.lambda2 );
+ const double xlog1 = log( ( q2+params_.q02 )/ params_.lambda2 ), xlog2 = log( params_.q02/params_.lambda2 );
const double t = log( xlog1/xlog2 );
- const double apom = param.pomeron.a[0] + ( param.pomeron.a[0]-param.pomeron.a[1] )*( 1./( 1.+pow( t, param.pomeron.a[2] ) ) - 1. );
- const double bpom = param.pomeron.b[0] + param.pomeron.b[1]*pow( t, param.pomeron.b[2] );
- const double cpom = param.pomeron.c[0] + ( param.pomeron.c[0]-param.pomeron.c[1] )*( 1./( 1.+pow( t, param.pomeron.c[2] ) ) - 1. );
+ const double apom = params_.pomeron.a[0] + ( params_.pomeron.a[0]-params_.pomeron.a[1] )*( 1./( 1.+pow( t, params_.pomeron.a[2] ) ) - 1. );
+ const double bpom = params_.pomeron.b[0] + params_.pomeron.b[1]*pow( t, params_.pomeron.b[2] );
+ const double cpom = params_.pomeron.c[0] + ( params_.pomeron.c[0]-params_.pomeron.c[1] )*( 1./( 1.+pow( t, params_.pomeron.c[2] ) ) - 1. );
- const double areg = param.reggeon.a[0] + param.reggeon.a[1]*pow( t, param.reggeon.a[2] );
- const double breg = param.reggeon.b[0] + param.reggeon.b[1]*pow( t, param.reggeon.b[2] );
- const double creg = param.reggeon.c[0] + param.reggeon.c[1]*pow( t, param.reggeon.c[2] );
+ const double areg = params_.reggeon.a[0] + params_.reggeon.a[1]*pow( t, params_.reggeon.a[2] );
+ const double breg = params_.reggeon.b[0] + params_.reggeon.b[1]*pow( t, params_.reggeon.b[2] );
+ const double creg = params_.reggeon.c[0] + params_.reggeon.c[1]*pow( t, params_.reggeon.c[2] );
const double F2_Pom = cpom*pow( xp, apom )*pow( 1.-xbj, bpom ),
F2_Reg = creg*pow( xr, areg )*pow( 1.-xbj, breg );
StructureFunctions allm;
allm.F2 = factor * ( F2_Pom + F2_Reg );
return allm;
}
}
}
diff --git a/CepGen/StructureFunctions/ALLM.h b/CepGen/StructureFunctions/ALLM.h
index 6884037..9c20cb0 100644
--- a/CepGen/StructureFunctions/ALLM.h
+++ b/CepGen/StructureFunctions/ALLM.h
@@ -1,47 +1,55 @@
#ifndef CepGen_StructureFunctions_ALLM_h
#define CepGen_StructureFunctions_ALLM_h
#include "StructureFunctions.h"
namespace CepGen
{
namespace SF
{
- class ALLMParameterisation
+ class ALLM
{
- private:
- struct Parameters {
- Parameters() :
- a( { 0., 0., 0. } ), b( { 0., 0., 0. } ), c( { 0., 0., 0. } ) {}
- Parameters( std::vector<double> c, std::vector<double> a, std::vector<double> b ) :
- a( a ), b( b ), c( c ) {}
- std::vector<double> a, b, c;
+ public:
+ class Parameterisation
+ {
+ private:
+ struct Parameters {
+ Parameters() :
+ a( { 0., 0., 0. } ), b( { 0., 0., 0. } ), c( { 0., 0., 0. } ) {}
+ Parameters( std::vector<double> c, std::vector<double> a, std::vector<double> b ) :
+ a( a ), b( b ), c( c ) {}
+ std::vector<double> a, b, c;
+ };
+
+ public:
+ /// Pre-HERA data fit (694 data points)
+ static Parameterisation allm91();
+ /// Fixed target and HERA photoproduction total cross sections (1356 points)
+ static Parameterisation allm97();
+ static Parameterisation hht_allm();
+ static Parameterisation hht_allm_ft();
+ static Parameterisation gd07p();
+ static Parameterisation gd11p();
+
+ Parameters pomeron, reggeon;
+ /// Effective photon squared mass
+ double m02;
+ /// Effective pomeron squared mass
+ double mp2;
+ /// Effective reggeon squared mass
+ double mr2;
+ double q02;
+ /// Squared QCD scale
+ double lambda2;
};
- public:
- /// Pre-HERA data fit (694 data points)
- static ALLMParameterisation allm91();
- /// Fixed target and HERA photoproduction total cross sections (1356 points)
- static ALLMParameterisation allm97();
- static ALLMParameterisation hht_allm();
- static ALLMParameterisation hht_allm_ft();
- static ALLMParameterisation gd07p();
- static ALLMParameterisation gd11p();
+ ALLM( const ALLM::Parameterisation& param = ALLM::Parameterisation::allm97() ) : params_( param ) {}
+ StructureFunctions operator()( double q2, double xbj ) const;
- Parameters pomeron, reggeon;
- /// Effective photon squared mass
- double m02;
- /// Effective pomeron squared mass
- double mp2;
- /// Effective reggeon squared mass
- double mr2;
- double q02;
- /// Squared QCD scale
- double lambda2;
+ private:
+ Parameterisation params_;
};
-
- StructureFunctions ALLM( double q2, double xbj, const ALLMParameterisation& param=ALLMParameterisation::allm97() );
}
}
#endif
diff --git a/CepGen/StructureFunctions/BlockDurandHa.cpp b/CepGen/StructureFunctions/BlockDurandHa.cpp
index 149768a..2c3dce1 100644
--- a/CepGen/StructureFunctions/BlockDurandHa.cpp
+++ b/CepGen/StructureFunctions/BlockDurandHa.cpp
@@ -1,40 +1,41 @@
#include "BlockDurandHa.h"
namespace CepGen
{
namespace SF
{
- BlockDurandHaParameters
- BlockDurandHaParameters::standard()
+ BlockDurandHa::Parameterisation
+ BlockDurandHa::Parameterisation::standard()
{
- BlockDurandHaParameters p;
- p.a = { 8.205e-4, -5.148e-2, -4.725e-3 };
- p.b = { 2.217e-3, 1.244e-2, 5.958e-4 };
- p.c = { 0.255e0, 1.475e-1 };
+ Parameterisation p;
+ p.a = { { 8.205e-4, -5.148e-2, -4.725e-3 } };
+ p.b = { { 2.217e-3, 1.244e-2, 5.958e-4 } };
+ p.c = { { 0.255e0, 1.475e-1 } };
p.n = 11.49;
p.lambda = 2.430;
p.mu2 = 2.82;
p.m2 = 0.753;
return p;
}
- StructureFunctions BlockDurandHa( double q2, double xbj, const BlockDurandHaParameters& params )
+ StructureFunctions
+ BlockDurandHa::operator()( double q2, double xbj ) const
{
StructureFunctions bdh;
if ( q2 <= 0. ) return bdh;
- const double tau = q2 / ( q2 + params.mu2 );
- const double xl = log( 1. + q2 / params.mu2 );
+ const double tau = q2 / ( q2 + params_.mu2 );
+ const double xl = log( 1. + q2 / params_.mu2 );
const double xlx = log( tau/xbj );
- const double A = params.a[0] + params.a[1]*xl + params.a[2]*xl*xl;
- const double B = params.b[0] + params.b[1]*xl + params.b[2]*xl*xl;
- const double C = params.c[0] + params.c[1]*xl;
- const double D = q2*( q2+params.lambda*params.m2 ) / pow( q2+params.m2, 2 );
+ const double A = params_.a[0] + params_.a[1]*xl + params_.a[2]*xl*xl;
+ const double B = params_.b[0] + params_.b[1]*xl + params_.b[2]*xl*xl;
+ const double C = params_.c[0] + params_.c[1]*xl;
+ const double D = q2*( q2+params_.lambda*params_.m2 ) / pow( q2+params_.m2, 2 );
- bdh.F2 = D*pow( 1.-xbj, params.n ) * ( C + A*xlx + B*xlx*xlx );
+ bdh.F2 = D*pow( 1.-xbj, params_.n ) * ( C + A*xlx + B*xlx*xlx );
return bdh;
}
}
}
diff --git a/CepGen/StructureFunctions/BlockDurandHa.h b/CepGen/StructureFunctions/BlockDurandHa.h
index 44feb75..0833d68 100644
--- a/CepGen/StructureFunctions/BlockDurandHa.h
+++ b/CepGen/StructureFunctions/BlockDurandHa.h
@@ -1,28 +1,36 @@
#ifndef CepGen_StructureFunctions_BlockDurandHa_h
#define CepGen_StructureFunctions_BlockDurandHa_h
#include "StructureFunctions.h"
#include <array>
namespace CepGen
{
namespace SF
{
- struct BlockDurandHaParameters
+ class BlockDurandHa
{
- std::array<double,3> a, b;
- std::array<double,2> c;
- double n;
- /// Effective mass spread parameter
- double lambda;
- /// Asymptotic log-behaviour transition scale factor
- double mu2;
- /// Squared effective mass (~VM mass)
- double m2;
- static BlockDurandHaParameters standard();
+ public:
+ struct Parameterisation
+ {
+ std::array<double,3> a, b;
+ std::array<double,2> c;
+ double n;
+ /// Effective mass spread parameter
+ double lambda;
+ /// Asymptotic log-behaviour transition scale factor
+ double mu2;
+ /// Squared effective mass (~VM mass)
+ double m2;
+ static Parameterisation standard();
+ };
+ BlockDurandHa( const BlockDurandHa::Parameterisation params = BlockDurandHa::Parameterisation::standard() ) : params_( params ) {}
+ StructureFunctions operator()( double q2, double xbj ) const;
+
+ private:
+ Parameterisation params_;
};
- StructureFunctions BlockDurandHa( double q2, double xbj, const BlockDurandHaParameters& params=BlockDurandHaParameters::standard() );
}
}
#endif
diff --git a/CepGen/StructureFunctions/ChristyBosted.cpp b/CepGen/StructureFunctions/ChristyBosted.cpp
index b1a3649..fec09d9 100644
--- a/CepGen/StructureFunctions/ChristyBosted.cpp
+++ b/CepGen/StructureFunctions/ChristyBosted.cpp
@@ -1,260 +1,256 @@
#include "ChristyBosted.h"
namespace CepGen
{
namespace SF
{
- ChristyBosted::ChristyBosted( const ChristyBostedParameterisation& params ) :
- params_( params )
- {}
-
double
ChristyBosted::resmod507( char sf, double w2, double q2 ) const
{
const double mp = Constants::mp, mp2 = mp*mp,
mpi = Constants::mpi, mpi2 = mpi*mpi,
meta = Particle::massFromPDGId( Particle::Eta ), meta2 = meta*meta;
const double w = sqrt( w2 );
const double xb = q2/( q2+w2-mp2 );
double m0 = 0., q20 = 0.;
if ( sf == 'T' ) { // transverse
m0 = 0.125;
q20 = 0.05;
}
else if ( sf == 'L' ) {
m0 = params_.m0;
q20 = 0.125;
}
else {
InError( "Invalid direction retrieved! Aborting." )
return 0.;
}
const double norm_q2 = 1./0.330/0.330;
const double t = log( log( ( q2+m0 )*norm_q2 )/log( m0*norm_q2 ) );
//--- calculate kinematics needed for threshold relativistic B-W
// equivalent photon energies
const double k = 0.5 * ( w2-mp2 )/mp;
const double kcm = 0.5 * ( w2-mp2 )/w;
const double epicm = 0.5 * ( w2+mpi2-mp2 )/w, ppicm = sqrt( std::max( 0., epicm*epicm-mpi2 ) );
const double epi2cm = 0.5 * ( w2 + 4.*mpi2 - mp2 )/w, ppi2cm = sqrt( std::max( 0., epi2cm*epi2cm-4*mpi2 ) );
const double eetacm = 0.5 * ( w2 + meta2 - mp2 )/w, petacm = sqrt( std::max( 0., eetacm*eetacm-meta2 ) );
std::array<double,7> width, height, pgam;
for ( unsigned short i = 0; i < 7; ++i ) {
- const ChristyBostedParameterisation::ResonanceParameters& res = params_.resonances[i];
+ const Parameterisation::ResonanceParameters& res = params_.resonances[i];
width[i] = res.width;
//--- calculate partial widths
//----- 1-pion decay mode
const double x02 = res.x0*res.x0;
const double partial_width_singlepi = pow( ppicm /res.pcmr( mpi2 ), 2.*res.angular_momentum+1. )
* pow( ( res.pcmr( mpi2 )*res.pcmr( mpi2 )+x02 )/( ppicm *ppicm +x02 ), res.angular_momentum );
//----- 2-pion decay mode
const double partial_width_doublepi = pow( ppi2cm/res.pcmr( 4.*mpi2 ), 2.*( res.angular_momentum+2. ) )
* pow( ( res.pcmr( 4.*mpi2 )*res.pcmr( 4.*mpi2 )+x02 )/( ppi2cm*ppi2cm+x02 ), res.angular_momentum+2 )
* w / res.mass;
//----- eta decay mode (only for S11's)
const double partial_width_eta = ( res.br.eta == 0. ) ? 0. :
pow( petacm/res.pcmr( meta2 ), 2.*res.angular_momentum+1. )
* pow( ( res.pcmr( meta2 )*res.pcmr( meta2 )+x02 )/( petacm*petacm+x02 ), res.angular_momentum );
// virtual photon width
pgam[i] = res.width * pow( kcm/res.kcmr(), 2 ) * ( res.kcmr()*res.kcmr()+x02 )/( kcm*kcm+x02 );
width[i] = ( partial_width_singlepi * res.br.singlepi
+ partial_width_doublepi * res.br.doublepi
+ partial_width_eta * res.br.eta ) * res.width;
//std::cout << i << "\t" << width[i] << "\t" << partial_width_doublepi << "\t" << pow( ppi2cm/res.pcmr( 4.*mpi2 ), 2.*( res.angular_momentum+2. ) ) << std::endl;
//--- resonance Q^2 dependence calculations
if ( sf == 'T' ) height[i] = res.A0_T*( 1.+res.fit_parameters[0]*q2/( 1.+res.fit_parameters[1]*q2 ) )/pow( 1.+q2/0.91, res.fit_parameters[2] );
else if ( sf == 'L' ) height[i] = res.A0_L/( 1.+res.fit_parameters[3]*q2 )*q2*exp( -q2*res.fit_parameters[4] );
//std::cout << sf << "//////" << height[i] << "\t" << res.fit_parameters[4] << std::endl;
height[i] = height[i]*height[i];
}
//--- calculate Breit-Wigners for all resonances
double sig_res = 0.;
for ( unsigned short i = 0; i < 7; ++i ) {
- const ChristyBostedParameterisation::ResonanceParameters res = params_.resonances[i];
+ const Parameterisation::ResonanceParameters res = params_.resonances[i];
const double mass2 = res.mass*res.mass, width2 = width[i]*width[i];
const double sigr = height[i]*res.kr()/k*res.kcmr()/kcm/res.width * ( width[i]*pgam[i] / ( pow( w2-mass2, 2 ) + mass2*width2 ) );
sig_res += sigr;
}
sig_res *= w;
//--- finish resonances / start non-res background calculation
const double xpr = 1./( 1.+( w2-pow( mp+mpi, 2 ) )/( q2+q20 ) );
double sig_nr = 0.;
if ( sf == 'T' ) { // transverse
const double wdif = w - ( mp + mpi );
for ( unsigned short i = 0; i < 2; ++i ) {
const double expo = params_.continuum.transverse[i].fit_parameters[1]
+ params_.continuum.transverse[i].fit_parameters[2]*q2
+ params_.continuum.transverse[i].fit_parameters[3]*q2*q2;
sig_nr += params_.continuum.transverse[i].sig0 / pow( q2+params_.continuum.transverse[i].fit_parameters[0], expo ) * pow( wdif, i+1.5 );
}
sig_nr *= xpr;
}
else if ( sf == 'L' ) { // longitudinal
for ( unsigned short i = 0; i < 1; ++i ) {
const double expo = params_.continuum.longitudinal[i].fit_parameters[0]
+ params_.continuum.longitudinal[i].fit_parameters[1];
sig_nr += params_.continuum.longitudinal[i].sig0
* pow( 1.-xpr, expo )/( 1.-xb )
* pow( q2/( q2+q20 ), params_.continuum.longitudinal[i].fit_parameters[2] )/( q2+q20 )
* pow( xpr, params_.continuum.longitudinal[i].fit_parameters[3]+params_.continuum.longitudinal[i].fit_parameters[4]*t );
}
}
return sig_res + sig_nr;
}
- ChristyBosted::ChristyBostedParameterisation
- ChristyBosted::ChristyBostedParameterisation::standard()
+ ChristyBosted::Parameterisation
+ ChristyBosted::Parameterisation::standard()
{
- ChristyBostedParameterisation params;
+ Parameterisation params;
params.m0 = 4.2802;
params.continuum.transverse = {
ContinuumParameters::DirectionParameters( 246.06, { 0.067469, 1.3501, 0.12054, -0.0038495 } ),
ContinuumParameters::DirectionParameters( -89.360, { 0.20977, 1.5715, 0.090736, 0.010362 } )
};
params.continuum.longitudinal = {
ContinuumParameters::DirectionParameters( 86.746, { 0., 4.0294, 3.1285, 0.33403, 4.9623 } )
};
//--- P33(1232)
ResonanceParameters p33;
p33.br = ResonanceParameters::BranchingRatios( 1., 0., 0. );
p33.angular_momentum = 1.;
//p33.x0 = 0.15;
p33.x0 = 0.14462;
p33.mass = 1.2298;
p33.width = 0.13573;
p33.fit_parameters = { 4.2291, 1.2598, 2.1242, 19.910, 0.22587 };
p33.A0_T = 7.7805;
p33.A0_L = 29.414;
params.resonances.emplace_back( p33 );
//--- S11(1535)
ResonanceParameters s11_1535;
s11_1535.br = ResonanceParameters::BranchingRatios( 0.45, 0.1, 0.45 );
s11_1535.angular_momentum = 0.;
s11_1535.x0 = 0.215;
s11_1535.mass = 1.5304;
s11_1535.width = 0.220;
s11_1535.fit_parameters = { 6823.2, 33521., 2.5686, 0., 0. };
s11_1535.A0_T = 6.3351;
s11_1535.A0_L = 0.;
params.resonances.emplace_back( s11_1535 );
//--- D13(1520)
ResonanceParameters d13;
d13.br = ResonanceParameters::BranchingRatios( 0.65, 0.35, 0. );
d13.angular_momentum = 2.;
d13.x0 = 0.215;
d13.mass = 1.5057;
d13.width = 0.082956;
d13.fit_parameters = { 21.240, 0.055746, 2.4886, 97.046, 0.31042 };
d13.A0_T = 0.60347;
d13.A0_L = 157.92;
params.resonances.emplace_back( d13 );
//--- F15(1680)
ResonanceParameters f15;
f15.br = ResonanceParameters::BranchingRatios( 0.65, 0.35, 0. );
f15.angular_momentum = 3.;
f15.x0 = 0.215;
f15.mass = 1.6980;
f15.width = 0.095782;
f15.fit_parameters = { -0.28789, 0.18607, 0.063534, 0.038200, 1.2182 };
f15.A0_T = 2.3305;
f15.A0_L = 4.2160;
params.resonances.emplace_back( f15 );
//--- S11(1650)
ResonanceParameters s11_1650;
s11_1650.br = ResonanceParameters::BranchingRatios( 0.4, 0.5, 0.1 );
s11_1650.angular_momentum = 0.;
s11_1650.x0 = 0.215;
s11_1650.mass = 1.6650;
s11_1650.width = 0.10936;
s11_1650.fit_parameters = { -0.56175, 0.38964, 0.54883, 0.31393, 2.9997 };
s11_1650.A0_T = 1.9790;
s11_1650.A0_L = 13.764;
params.resonances.emplace_back( s11_1650 );
//--- P11(1440) roper
ResonanceParameters p11;
p11.br = ResonanceParameters::BranchingRatios( 0.65, 0.35, 0. );
p11.angular_momentum = 1.;
p11.x0 = 0.215;
p11.mass = 1.4333;
p11.width = 0.37944;
p11.fit_parameters = { 46.213, 0.19221, 1.9141, 0.053743, 1.3091 };
p11.A0_T = 0.022506;
p11.A0_L = 5.5124;
params.resonances.emplace_back( p11 );
//--- F37(1950)
ResonanceParameters f37;
f37.br = ResonanceParameters::BranchingRatios( 0.5, 0.5, 0. );
f37.angular_momentum = 3.;
f37.x0 = 0.215;
f37.mass = 1.9341;
f37.width = 0.380;
f37.fit_parameters = { 0., 0., 1., 1.8951, 0.51376 };
f37.A0_T = 3.4187;
f37.A0_L = 1.8951;
params.resonances.emplace_back( f37 );
return params;
}
StructureFunctions
ChristyBosted::operator()( double q2, double xbj ) const
{
const double mp2 = Constants::mp*Constants::mp;
const double w2 = mp2 + q2*( 1.-xbj )/xbj;
const double w_min = Constants::mp+Constants::mpi;
StructureFunctions cb;
if ( sqrt( w2 ) < w_min ) return cb;
//-----------------------------
// modification of Christy-Bosted at large q2 as described in the LUXqed paper
//-----------------------------
const double q21 = 30., q20 = 8.;
const double delq2 = q2 - q20;
const double qq = q21 - q20;
const double prefac = 1./( 4.*M_PI*M_PI*Constants::alphaEM ) * ( 1.-xbj );
//------------------------------
double q2_eff = q2, w2_eff = w2;
if ( q2 > q20 ) {
q2_eff = q20 + delq2/( 1.+delq2/qq );
w2_eff = mp2 + q2_eff*( 1.-xbj )/xbj;
}
const double tau = 4.*xbj*xbj*mp2/q2_eff;
const double sigT = resmod507( 'T', w2_eff, q2_eff ), sigL = resmod507( 'L', w2_eff, q2_eff );
cb.F2 = prefac * q2_eff / ( 1+tau ) * ( sigT+sigL ) / Constants::GeV2toBarn * 1.e6;
if ( q2 > q20 ) cb.F2 *= q21/( q21 + delq2 );
const double R = sigL/sigT;
cb.FL = cb.F2*( 1+tau )*R/( 1.+R );
return cb;
}
}
}
diff --git a/CepGen/StructureFunctions/ChristyBosted.h b/CepGen/StructureFunctions/ChristyBosted.h
index 3a1f1cc..cd03870 100644
--- a/CepGen/StructureFunctions/ChristyBosted.h
+++ b/CepGen/StructureFunctions/ChristyBosted.h
@@ -1,79 +1,79 @@
#ifndef CepGen_StructureFunctions_ChristyBosted_h
#define CepGen_StructureFunctions_ChristyBosted_h
#include "StructureFunctions.h"
#include <array>
#include <vector>
namespace CepGen
{
namespace SF
{
class ChristyBosted
{
public:
- struct ChristyBostedParameterisation
+ struct Parameterisation
{
- static ChristyBostedParameterisation standard();
+ static Parameterisation standard();
struct ResonanceParameters
{
struct BranchingRatios
{
BranchingRatios() : singlepi( 0. ), doublepi( 0. ), eta( 0. ) {}
BranchingRatios( double singlepi, double doublepi, double eta ) : singlepi( singlepi ), doublepi( doublepi ), eta( eta ) {}
bool valid() const { return ( singlepi+doublepi+eta == 1. ); }
/// single pion branching ratio
double singlepi;
/// double pion branching ratio
double doublepi;
/// eta meson branching ratio
double eta;
};
ResonanceParameters() : angular_momentum( 0. ), x0( 0. ), mass( 0. ), width( 0. ), A0_T( 0. ), A0_L( 0. ) {}
double kr() const { return 0.5 * ( mass*mass-Constants::mp*Constants::mp )/Constants::mp; }
double ecmr( double m2 ) const { return ( mass == 0. ) ? 0. : 0.5 * ( mass*mass+m2-Constants::mp*Constants::mp ) / mass; }
double kcmr() const { return ecmr( 0. ); }
double pcmr( double m2 ) const { return sqrt( std::max( 0., ecmr( m2 )*ecmr( m2 )-m2 ) ); }
BranchingRatios br;
/// meson angular momentum
double angular_momentum;
/// damping parameter
double x0;
/// mass, in GeV/c2
double mass;
/// full width, in GeV
double width;
double A0_T;
double A0_L;
std::array<double,5> fit_parameters;
};
struct ContinuumParameters
{
struct DirectionParameters
{
DirectionParameters() : sig0( 0. ) {}
DirectionParameters( double sig0, const std::vector<double>& params ) : sig0( sig0 ), fit_parameters( params ) {}
double sig0;
std::vector<double> fit_parameters;
};
std::array<DirectionParameters,2> transverse;
std::array<DirectionParameters,1> longitudinal;
};
double m0;
std::vector<ResonanceParameters> resonances;
ContinuumParameters continuum;
};
- ChristyBosted( const ChristyBostedParameterisation& params = ChristyBostedParameterisation::standard() );
+ ChristyBosted( const ChristyBosted::Parameterisation& params = ChristyBosted::Parameterisation::standard() ) : params_( params ) {}
StructureFunctions operator()( double q2, double xbj ) const;
private:
double resmod507( char sf, double w2, double q2 ) const;
- ChristyBostedParameterisation params_;
+ Parameterisation params_;
};
}
}
#endif
diff --git a/CepGen/StructureFunctions/FioreBrasse.cpp b/CepGen/StructureFunctions/FioreBrasse.cpp
index bc6f117..b3991cd 100644
--- a/CepGen/StructureFunctions/FioreBrasse.cpp
+++ b/CepGen/StructureFunctions/FioreBrasse.cpp
@@ -1,153 +1,152 @@
#include "FioreBrasse.h"
namespace CepGen
{
namespace SF
{
- FioreBrasseParameterisation
- FioreBrasseParameterisation::standard()
+ FioreBrasse::Parameterisation
+ FioreBrasse::Parameterisation::standard()
{
- FioreBrasseParameterisation p;
+ Parameterisation p;
p.s0 = 1.14;
p.norm = 0.021;
p.resonances.emplace_back( -0.8377, 0.95, 0.1473, 1.0, 2.4617, 3./2. );
p.resonances.emplace_back( -0.37, 0.95, 0.1471, 0.5399, 2.4617, 5./2. );
p.resonances.emplace_back( 0.0038, 0.85, 0.1969, 4.2225, 1.5722, 3./2. );
p.resonances.emplace_back( 0.5645, 0.1126, 1.3086, 19.2694, 4.5259, 1. );
return p;
}
- FioreBrasseParameterisation
- FioreBrasseParameterisation::alternative()
+ FioreBrasse::Parameterisation
+ FioreBrasse::Parameterisation::alternative()
{
- FioreBrasseParameterisation p;
+ Parameterisation p;
p.s0 = 1.2871;
p.norm = 0.0207;
p.resonances.emplace_back( -0.8070, 0.9632, 0.1387, 1.0, 2.6066, 3./2. );
p.resonances.emplace_back( -0.3640, 0.9531, 0.1239, 0.6086, 2.6066, 5./2. );
p.resonances.emplace_back( -0.0065, 0.8355, 0.2320, 4.7279, 1.4828, 3./2. );
p.resonances.emplace_back( 0.5484, 0.1373, 1.3139, 14.7267, 4.6041, 1. );
return p;
}
StructureFunctions
- FioreBrasse( double q2, double xbj )
+ FioreBrasse::operator()( double q2, double xbj ) const
{
const double mp = Particle::massFromPDGId( Particle::Proton ), mp2 = mp*mp;
const double akin = 1. + 4.*mp2 * xbj*xbj/q2;
const double prefactor = q2*( 1.-xbj ) / ( 4.*M_PI*Constants::alphaEM*akin );
const double s = q2*( 1.-xbj )/xbj + mp2;
- FioreBrasseParameterisation p = FioreBrasseParameterisation::standard();
double ampli_res = 0., ampli_tot = 0.;
for ( unsigned short i = 0; i < 3; ++i ) { //FIXME 4??
- const FioreBrasseParameterisation::ResonanceParameters res = p.resonances[i];
+ const Parameterisation::ResonanceParameters res = params_.resonances[i];
if ( !res.enabled ) continue;
- const double sqrts0 = sqrt( p.s0 );
+ const double sqrts0 = sqrt( params_.s0 );
std::complex<double> alpha;
- if ( s > p.s0 )
- alpha = std::complex<double>( res.alpha0 + res.alpha2*sqrts0 + res.alpha1*s, res.alpha2*sqrt( s-p.s0 ) );
+ if ( s > params_.s0 )
+ alpha = std::complex<double>( res.alpha0 + res.alpha2*sqrts0 + res.alpha1*s, res.alpha2*sqrt( s-params_.s0 ) );
else
- alpha = std::complex<double>( res.alpha0 + res.alpha1*s + res.alpha2*( sqrts0 - sqrt( p.s0 - s ) ), 0. );
+ alpha = std::complex<double>( res.alpha0 + res.alpha1*s + res.alpha2*( sqrts0 - sqrt( params_.s0 - s ) ), 0. );
double formfactor = 1./pow( 1. + q2/res.q02, 2 );
double denom = pow( res.spin-std::real( alpha ), 2 ) + pow( std::imag( alpha ), 2 );
double ampli_imag = res.a*formfactor*formfactor*std::imag( alpha )/denom;
ampli_res += ampli_imag;
}
{
- const FioreBrasseParameterisation::ResonanceParameters res = p.resonances[3];
+ const Parameterisation::ResonanceParameters res = params_.resonances[3];
double sE = res.alpha2, sqrtsE = sqrt( sE );
std::complex<double> alpha;
if ( s > sE )
alpha = std::complex<double>( res.alpha0 + res.alpha1*sqrtsE, res.alpha1*sqrt( s-sE ) );
else
alpha = std::complex<double>( res.alpha0 + res.alpha1*( sqrtsE - sqrt( sE-s ) ), 0. );
double formfactor = 1./pow( 1. + q2/res.q02, 2 );
double sp = 1.5*res.spin;
double denom = pow( sp-std::real( alpha ), 2 ) + pow( std::imag( alpha ), 2 );
double ampli_bg = res.a*formfactor*formfactor*std::imag( alpha )/denom;
ampli_res += ampli_bg;
}
- ampli_tot = p.norm*ampli_res;
+ ampli_tot = params_.norm*ampli_res;
StructureFunctions fb;
fb.F2 = prefactor*ampli_tot;
return fb;
}
StructureFunctions
- FioreBrasseOld( double q2, double xbj )
+ FioreBrasse::operator()( double q2, double xbj, bool ) const
{
const double mp = Particle::massFromPDGId( Particle::Proton ), mp2 = mp*mp;
//const double m_min = Particle::massFromPDGId(Particle::Proton)+0.135;
const double m_min = mp+Particle::massFromPDGId( Particle::PiZero );
const double mx2 = mp2 + q2*( 1.-xbj )/xbj, mx = sqrt( mx2 );
if ( mx < m_min || mx > 1.99 ) {
InWarning( Form( "Fiore-Brasse form factors to be retrieved for an invalid MX value:\n\t"
"%.2e GeV, while allowed range is [1.07, 1.99] GeV", mx ) );
return StructureFunctions();
}
int n_bin;
double x_bin, dx;
if ( mx < 1.11 ) {
n_bin = 0;
x_bin = mx-m_min;
dx = 1.11-m_min; // Delta w bin sizes
}
else if ( mx < 1.77 ) { // w in [1.11, 1.77[
dx = 0.015; // Delta w bin sizes
n_bin = ( mx-1.11 )/dx + 1;
x_bin = fmod( mx-1.11, dx );
}
else { // w in [1.77, 1.99[
dx = 0.02; // Delta w bin sizes
n_bin = ( mx-1.77 )/dx + 45;
x_bin = fmod( mx-1.77, dx );
}
// values of a, b, c provided from the fits on ep data and retrieved from
// http://dx.doi.org/10.1016/0550-3213(76)90231-5 with 1.110 <= w2 <=1.990
const double a[56] = { 5.045, 5.126, 5.390,5.621, 5.913, 5.955,6.139,6.178,6.125, 5.999,
5.769, 5.622, 5.431,5.288, 5.175, 5.131,5.003,5.065,5.045, 5.078,
5.145, 5.156, 5.234,5.298, 5.371, 5.457,5.543,5.519,5.465, 5.384,
5.341, 5.320, 5.275,5.290, 5.330, 5.375,5.428,5.478,5.443, 5.390,
5.333, 5.296, 5.223,5.159, 5.146, 5.143,5.125,5.158,5.159, 5.178,
5.182, 5.195, 5.160,5.195, 5.163, 5.172 },
b[56] = { 0.798, 1.052, 1.213,1.334,1.397,1.727,1.750,1.878,1.887,1.927,
2.041, 2.089, 2.148,2.205,2.344,2.324,2.535,2.464,2.564,2.610,
2.609, 2.678, 2.771,2.890,2.982,3.157,3.183,3.315,3.375,3.450,
3.477, 3.471, 3.554,3.633,3.695,3.804,3.900,4.047,4.290,4.519,
4.709, 4.757, 4.840,5.017,5.015,5.129,5.285,5.322,5.545,5.623,
5.775, 5.894, 6.138,6.151,6.301,6.542 },
c[56] = { 0.043, 0.024, 0.000,-0.013,-0.023,-0.069,-0.060,-0.080,-0.065,-0.056,
-0.065,-0.056,-0.043,-0.034,-0.054,-0.018,-0.046,-0.015,-0.029,-0.048,
-0.032,-0.045,-0.084,-0.115,-0.105,-0.159,-0.164,-0.181,-0.203,-0.223,
-0.245,-0.254,-0.239,-0.302,-0.299,-0.318,-0.383,-0.393,-0.466,-0.588,
-0.622,-0.568,-0.574,-0.727,-0.665,-0.704,-0.856,-0.798,-1.048,-0.980,
-1.021,-1.092,-1.313,-1.341,-1.266,-1.473 };
const double d = 3.0;
const double nu = 0.5 * ( q2 + mx2 - mp2 ) / mp, nu2 = nu*nu,
logqq0 = 0.5 * log( ( nu2+q2 ) / pow( ( mx2-mp2 ) / ( 2.*mp ), 2 ) );
const double gd2 = pow( 1. / ( 1+q2 / .71 ), 4 ); // dipole form factor of the proton
const double sigLow = ( n_bin == 0 ) ? 0. :
gd2 * exp( a[n_bin-1] + b[n_bin-1]*logqq0 + c[n_bin-1]*pow( fabs( logqq0 ), d ) );
const double sigHigh =
gd2 * exp( a[n_bin] + b[n_bin] *logqq0 + c[n_bin] *pow( fabs( logqq0 ), d ) );
const double sigma_t = sigLow + x_bin*( sigHigh-sigLow )/dx;
const double w1 = ( mx2-mp2 )/( 8.*M_PI*M_PI*mp*Constants::alphaEM )/Constants::GeV2toBarn*1.e6 * sigma_t;
const double w2 = w1 * q2 / ( q2+nu2 );
StructureFunctions fb_old;
fb_old.F1 = w1; //FIXME
fb_old.F2 = w2; //FIXME
return fb_old;
}
}
}
diff --git a/CepGen/StructureFunctions/FioreBrasse.h b/CepGen/StructureFunctions/FioreBrasse.h
index a3bd100..7c6e58e 100644
--- a/CepGen/StructureFunctions/FioreBrasse.h
+++ b/CepGen/StructureFunctions/FioreBrasse.h
@@ -1,37 +1,45 @@
#ifndef CepGen_StructureFunctions_FioreBrasse_h
#define CepGen_StructureFunctions_FioreBrasse_h
#include "StructureFunctions.h"
#include <complex>
namespace CepGen
{
namespace SF
{
- struct FioreBrasseParameterisation
+ class FioreBrasse
{
- static FioreBrasseParameterisation standard();
- static FioreBrasseParameterisation alternative();
+ public:
+ struct Parameterisation
+ {
+ static Parameterisation standard();
+ static Parameterisation alternative();
- struct ResonanceParameters {
- ResonanceParameters( double a0, double a1, double a2, double a, double q02, float spin, bool on=true ) :
- alpha0( a0 ), alpha1( a1 ), alpha2( a2 ), a( a ), q02( q02 ), spin( spin ), enabled( on ) {}
- double alpha0, alpha1, alpha2, a, q02;
- float spin;
- bool enabled;
- };
+ struct ResonanceParameters {
+ ResonanceParameters( double a0, double a1, double a2, double a, double q02, float spin, bool on=true ) :
+ alpha0( a0 ), alpha1( a1 ), alpha2( a2 ), a( a ), q02( q02 ), spin( spin ), enabled( on ) {}
+ double alpha0, alpha1, alpha2, a, q02;
+ float spin;
+ bool enabled;
+ };
- std::vector<ResonanceParameters> resonances;
- double s0, norm;
+ std::vector<ResonanceParameters> resonances;
+ double s0, norm;
+ };
+ /// Fiore-Brasse proton structure functions (F.W Brasse et al., DESY 76/11 (1976),
+ /// http://dx.doi.org/10.1016/0550-3213(76)90231-5)
+ FioreBrasse( const FioreBrasse::Parameterisation& params = FioreBrasse::Parameterisation::standard() ) : params_( params ) {}
+ /// \param[in] q2 Squared 4-momentum transfer
+ /// \param[in] xbj Bjorken's x
+ /// \cite Brasse1976413
+ StructureFunctions operator()( double q2, double xbj ) const;
+ StructureFunctions operator()( double q2, double xbj, bool old ) const;
+
+ private:
+ Parameterisation params_;
};
- /// Fiore-Brasse proton structure functions (F.W Brasse et al., DESY 76/11 (1976),
- /// http://dx.doi.org/10.1016/0550-3213(76)90231-5)
- /// \param[in] q2 Squared 4-momentum transfer
- /// \param[in] xbj Bjorken's x
- /// \cite Brasse1976413
- StructureFunctions FioreBrasse( double q2, double xbj );
- StructureFunctions FioreBrasseOld( double q2, double xbj );
}
}
#endif
diff --git a/CepGen/StructureFunctions/SuriYennie.cpp b/CepGen/StructureFunctions/SuriYennie.cpp
index 254e9a9..8b4318a 100644
--- a/CepGen/StructureFunctions/SuriYennie.cpp
+++ b/CepGen/StructureFunctions/SuriYennie.cpp
@@ -1,51 +1,51 @@
#include "SuriYennie.h"
namespace CepGen
{
namespace SF
{
- SuriYennieParameterisation
- SuriYennieParameterisation::standard()
+ SuriYennie::Parameterisation
+ SuriYennie::Parameterisation::standard()
{
- SuriYennieParameterisation p;
+ Parameterisation p;
p.C1 = 0.86926;
p.C2 = 2.23422;
p.D1 = 0.12549;
p.rho2 = 0.585;
p.Cp = 0.96;
p.Bp = 0.63;
return p;
}
- SuriYennieParameterisation
- SuriYennieParameterisation::alternative()
+ SuriYennie::Parameterisation
+ SuriYennie::Parameterisation::alternative()
{
- SuriYennieParameterisation p;
+ Parameterisation p;
p.C1 = 0.6303;
p.C2 = 2.3049;
p.D1 = 0.04681;
p.rho2 = 1.05;
p.Cp = 1.23;
p.Bp = 0.61;
return p;
}
StructureFunctions
- SuriYennie( double q2, double xbj, const SuriYennieParameterisation& param )
+ SuriYennie::operator()( double q2, double xbj ) const
{
const double mp = Particle::massFromPDGId( Particle::Proton ), mp2 = mp*mp;
const double mx2 = q2 * ( 1.-xbj )/xbj + mp2, // [GeV^2]
nu = 0.5 * ( q2 + mx2 - mp2 ) / mp; // [GeV]
- const double dm2 = mx2-mp2, Xpr = q2/( q2+mx2 ), En = dm2+q2, Tau = 0.25 * q2/mp2, MQ = param.rho2+q2;
+ const double dm2 = mx2-mp2, Xpr = q2/( q2+mx2 ), En = dm2+q2, Tau = 0.25 * q2/mp2, MQ = params_.rho2+q2;
StructureFunctions sy;
- sy.FM = ( param.C1*dm2*pow( param.rho2/MQ, 2 ) + ( param.C2*mp2*pow( 1.-Xpr, 4 ) ) / ( 1.+Xpr*( Xpr*param.Cp-2.*param.Bp ) ) )/q2;
- const double FE = ( Tau*sy.FM + param.D1*dm2*q2*param.rho2/mp2*pow( dm2/MQ/En, 2 ) ) / ( 1.+0.25*En*En/mp2/q2 );
+ sy.FM = ( params_.C1*dm2*pow( params_.rho2/MQ, 2 ) + ( params_.C2*mp2*pow( 1.-Xpr, 4 ) ) / ( 1.+Xpr*( Xpr*params_.Cp-2.*params_.Bp ) ) )/q2;
+ const double FE = ( Tau*sy.FM + params_.D1*dm2*q2*params_.rho2/mp2*pow( dm2/MQ/En, 2 ) ) / ( 1.+0.25*En*En/mp2/q2 );
const double w2 = 2.*mp*FE, w1 = 0.5 * sy.FM*q2/mp;
sy.F2 = nu/mp*w2;
return sy;
}
}
}
diff --git a/CepGen/StructureFunctions/SuriYennie.h b/CepGen/StructureFunctions/SuriYennie.h
index 6c00cef..a6503b6 100644
--- a/CepGen/StructureFunctions/SuriYennie.h
+++ b/CepGen/StructureFunctions/SuriYennie.h
@@ -1,22 +1,30 @@
#ifndef CepGen_StructureFunctions_SuriYennie_h
#define CepGen_StructureFunctions_SuriYennie_h
#include "StructureFunctions.h"
namespace CepGen
{
namespace SF
{
- struct SuriYennieParameterisation {
- // values extracted from experimental fits
- static SuriYennieParameterisation standard();
- static SuriYennieParameterisation alternative();
+ class SuriYennie
+ {
+ public:
+ struct Parameterisation {
+ // values extracted from experimental fits
+ static Parameterisation standard();
+ static Parameterisation alternative();
- double C1, C2, D1, rho2, Cp, Bp;
- };
+ double C1, C2, D1, rho2, Cp, Bp;
+ };
+
+ SuriYennie( const SuriYennie::Parameterisation& param = SuriYennie::Parameterisation::standard() ) : params_( param ) {}
+ StructureFunctions operator()( double q2, double xbj ) const;
- StructureFunctions SuriYennie( double q2, double xbj, const SuriYennieParameterisation& param=SuriYennieParameterisation::standard() );
+ private:
+ Parameterisation params_;
+ };
}
}
#endif
diff --git a/CepGen/StructureFunctions/SzczurekUleshchenko.cpp b/CepGen/StructureFunctions/SzczurekUleshchenko.cpp
index feee450..b471f43 100644
--- a/CepGen/StructureFunctions/SzczurekUleshchenko.cpp
+++ b/CepGen/StructureFunctions/SzczurekUleshchenko.cpp
@@ -1,39 +1,39 @@
#include "SzczurekUleshchenko.h"
namespace CepGen
{
namespace SF
{
StructureFunctions
- SzczurekUleshchenko( double q2, double xbj )
+ SzczurekUleshchenko::operator()( double q2, double xbj ) const
{
#ifndef GRVPDF
FatalError( "Szczurek-Uleshchenko structure functions cannot be computed"
" as GRV PDF set is not linked to this instance!" );
#else
const float q02 = 0.8;
float amu2 = q2+q02; // shift the overall scale
float xuv, xdv, xus, xds, xss, xg;
float xbj_arg = xbj;
grv95lo_( xbj_arg, amu2, xuv, xdv, xus, xds, xss, xg );
DebuggingInsideLoop( Form( "Form factor content at xB = %e (scale = %f GeV^2):\n\t"
" valence quarks: u / d = %e / %e\n\t"
" sea quarks: u / d / s = %e / %e / %e\n\t"
" gluons: = %e",
xbj, amu2, xuv, xdv, xus, xds, xss, xg ) );
// standard partonic structure function
const double F2_aux = 4./9.*( xuv + 2.*xus )
+ 1./9.*( xdv + 2.*xds )
+ 1./9.*( 2.*xss );
StructureFunctions su;
su.F2 = F2_aux * q2 / amu2; // F2 corrected for low Q^2 behaviour
return su;
#endif
}
}
}
diff --git a/CepGen/StructureFunctions/SzczurekUleshchenko.h b/CepGen/StructureFunctions/SzczurekUleshchenko.h
index 2cc58be..869edef 100644
--- a/CepGen/StructureFunctions/SzczurekUleshchenko.h
+++ b/CepGen/StructureFunctions/SzczurekUleshchenko.h
@@ -1,20 +1,25 @@
#ifndef CepGen_StructureFunctions_SzczurekUleshchenko_h
#define CepGen_StructureFunctions_SzczurekUleshchenko_h
#include "StructureFunctions.h"
extern "C"
{
extern void grv95lo_( float&, float&, float&, float&, float&, float&, float&, float& );
//extern void grv95lo_( double&, double&, double&, double&, double&, double&, double&, double& );
}
namespace CepGen
{
namespace SF
{
- StructureFunctions SzczurekUleshchenko( double q2, double xbj );
+ class SzczurekUleshchenko
+ {
+ public:
+ SzczurekUleshchenko() {}
+ StructureFunctions operator()( double q2, double xbj ) const;
+ };
}
}
#endif
diff --git a/test/utils/structure_functions.cxx b/test/utils/structure_functions.cxx
index ecb7b70..778c458 100644
--- a/test/utils/structure_functions.cxx
+++ b/test/utils/structure_functions.cxx
@@ -1,143 +1,147 @@
#include "CepGen/Physics/FormFactors.h"
#include "CepGen/Event/Particle.h"
#include "CepGen/StructureFunctions/ChristyBosted.h"
#include "test/Canvas.h"
#include "TGraph.h"
#include "TMultiGraph.h"
#include <iostream>
using namespace std;
int
main( int argc, char* argv[] )
{
const float min_xbj = 1.e-5, max_xbj = 0.99, q2 = ( argc>1 ) ? atof( argv[1] ) : 2.5;
const char* q2_str = ( argc>2 ) ? argv[2] : std::to_string( q2 ).c_str();
const unsigned int npoints = 5000;
TGraph g_sy_f2, g_fb_f2, g_su_f2, g_bdh_f2, g_cteq_f2, g_mrst_f2;
TGraph g_allm97_f2, g_allm_hht_f2, g_allm_hht_ft_f2;
TGraph g_lux_f2;
TGraph g_cb_f2;
const bool use_logarithmic_x = ( argc>3 ) ? atoi( argv[3] ) : false;
/*CepGen::SF::GenericLHAPDF cteq( "cteq6l1" );
CepGen::SF::GenericLHAPDF mrst( "MRST2004qed_proton" );
CepGen::SF::GenericLHAPDF lux( "LUXqed17_plus_PDF4LHC15_nnlo_100" );*/
+ CepGen::SF::SuriYennie sy;
+ CepGen::SF::SzczurekUleshchenko su;
CepGen::SF::ChristyBosted cb;
+ CepGen::SF::ALLM allm97( CepGen::SF::ALLM::Parameterisation::allm97() );
+ CepGen::SF::FioreBrasse fb;
for ( unsigned int i=0; i<npoints; i++ ) {
float xbj;
if ( use_logarithmic_x ) {
const float min_lxbj = log10( min_xbj ), max_lxbj = log10( max_xbj );
xbj = pow( 10, min_lxbj + i*( max_lxbj-min_lxbj )/( npoints-1 ) );
}
else xbj = min_xbj + i*( max_xbj-min_xbj )/( npoints-1 );
- auto sf_sy = CepGen::SF::SuriYennie( q2, xbj ),
- sf_fb = CepGen::SF::FioreBrasse( q2, xbj ),
- sf_su = CepGen::SF::SzczurekUleshchenko( q2, xbj ),
- sf_allm97 = CepGen::SF::ALLM( q2, xbj, CepGen::SF::ALLMParameterisation::allm97() ),
+ auto sf_sy = sy( q2, xbj ),
+ sf_fb = fb( q2, xbj ),
+ sf_su = su( q2, xbj ),
+ sf_allm97 = allm97( q2, xbj ),
//sf_allm_hht = CepGen::SF::ALLM( q2, xbj, CepGen::SF::ALLMParameterisation::hht_allm() ),
//sf_allm_hht_ft = CepGen::SF::ALLM( q2, xbj, CepGen::SF::ALLMParameterisation::hht_allm_ft() ),
//sf_bdh = CepGen::SF::BlockDurandHa( q2, xbj ),
//sf_cteq = cteq( q2, xbj ),
//sf_mrst = mrst( q2, xbj ),
//sf_lux = lux( q2, xbj ),
sf_cb = cb( q2, xbj );
g_sy_f2.SetPoint( i, xbj, sf_sy.F2 );
g_fb_f2.SetPoint( i, xbj, sf_fb.F2 );
g_su_f2.SetPoint( i, xbj, sf_su.F2 );
//g_bdh_f2.SetPoint( i, xbj, sf_bdh.F2 );
//g_cteq_f2.SetPoint( i, xbj, sf_cteq.F2 );
//g_mrst_f2.SetPoint( i, xbj, sf_mrst.F2 );
//g_lux_f2.SetPoint( i, xbj, sf_lux.F2 );
g_cb_f2.SetPoint( i, xbj, sf_cb.F2 );
g_allm97_f2.SetPoint( i, xbj, sf_allm97.F2 );
//g_allm_hht_f2.SetPoint( i, xbj, sf_allm_hht.F2 );
//g_allm_hht_ft_f2.SetPoint( i, xbj, sf_allm_hht_ft.F2 );
}
CepGen::Canvas c( "test", Form( "CepGen proton structure functions, Q^{2} = %s GeV^{2}", q2_str ) );
c.SetLegendX1( 0.4 );
TMultiGraph mg;
/*g_sy_f2.SetLineStyle( 2 );
g_sy_f2.SetLineWidth( 3 );
mg.Add( &g_sy_f2, "l" );
- c.AddLegendEntry( &g_sy_f2, "Suri-Yennie", "l" );
+ c.AddLegendEntry( &g_sy_f2, "Suri-Yennie", "l" );*/
//g_fb_f2.SetLineStyle( 2 );
g_fb_f2.SetLineColor( kRed+1 );
g_fb_f2.SetLineWidth( 3 );
mg.Add( &g_fb_f2, "l" );
c.AddLegendEntry( &g_fb_f2, "Fiore-Brasse", "l" );
//g_su_f2.SetLineStyle( 2 );
- g_su_f2.SetLineColor( kGreen+2 );
+ /*g_su_f2.SetLineColor( kGreen+2 );
g_su_f2.SetLineWidth( 3 );
mg.Add( &g_su_f2, "l" );
c.AddLegendEntry( &g_su_f2, "Szczurek-Uleshchenko", "l" );*/
g_allm97_f2.SetLineColor( kBlue+1 );
g_allm97_f2.SetLineWidth( 3 );
mg.Add( &g_allm97_f2, "l" );
c.AddLegendEntry( &g_allm97_f2, "Abramowicz et al. 97", "l" );
/*g_allm_hht_f2.SetLineColor( kBlue+1 );
g_allm_hht_f2.SetLineWidth( 3 );
g_allm_hht_f2.SetLineStyle( 2 );
mg.Add( &g_allm_hht_f2, "l" );
c.AddLegendEntry( &g_allm_hht_f2, "Abramowicz et al. HHT", "l" );
g_allm_hht_ft_f2.SetLineColor( kBlue+1 );
g_allm_hht_ft_f2.SetLineWidth( 3 );
g_allm_hht_ft_f2.SetLineStyle( 3 );
mg.Add( &g_allm_hht_ft_f2, "l" );
c.AddLegendEntry( &g_allm_hht_ft_f2, "Abramowicz et al. HHT-FT", "l" );*/
g_cb_f2.SetLineColor( kMagenta );
g_cb_f2.SetLineWidth( 3 );
mg.Add( &g_cb_f2, "l" );
c.AddLegendEntry( &g_cb_f2, "Christy-Bosted", "l" );
/*g_bdh_f2.SetLineColor( kOrange );
g_bdh_f2.SetLineWidth( 3 );
mg.Add( &g_bdh_f2, "l" );
c.AddLegendEntry( &g_bdh_f2, "Block-Durand-Ha", "l" );
g_cteq_f2.SetLineColor( kMagenta );
g_cteq_f2.SetLineWidth( 3 );
mg.Add( &g_cteq_f2, "l" );
c.AddLegendEntry( &g_cteq_f2, "CTEQ6", "l" );
g_mrst_f2.SetLineColor( kMagenta+1 );
g_mrst_f2.SetLineWidth( 3 );
mg.Add( &g_mrst_f2, "l" );
c.AddLegendEntry( &g_mrst_f2, "MRST2004 (QED/proton)", "l" );*/
/*g_lux_f2.SetLineColor( kOrange );
g_lux_f2.SetLineWidth( 3 );
mg.Add( &g_lux_f2, "l" );
c.AddLegendEntry( &g_lux_f2, "LUXqed", "l" );*/
mg.Draw( "alpr" );
mg.SetTitle( "x_{Bj}\\Proton form factor F_{2}" );
c.Prettify( mg.GetHistogram() );
//mg.GetYaxis()->SetRangeUser( 1.e-10, 10. );
mg.GetYaxis()->SetRangeUser( 0., 0.8 );
mg.GetXaxis()->SetLimits( min_xbj, max_xbj );
if ( use_logarithmic_x ) c.SetLogx();
//c.SetLogy();
c.Save( "pdf" );
return 0;
}