diff --git a/CepGen/IO/FortranInterface.cpp b/CepGen/IO/FortranInterface.cpp
index 8fd67b2..b71aff5 100644
--- a/CepGen/IO/FortranInterface.cpp
+++ b/CepGen/IO/FortranInterface.cpp
@@ -1,20 +1,20 @@
#include "CepGen/StructureFunctions/StructureFunctionsBuilder.h"
-#include "CepGen/IO/MSTWGridHandler.h"
+#include "CepGen/StructureFunctions/StructureFunctions.h"
#include "CepGen/Core/Exception.h"
extern "C"
{
void
cepgen_structure_functions_( int& sfmode, double& q2, double& xbj, double& f2, double& fl )
{
using namespace CepGen;
SF::Type sf_mode = (SF::Type)sfmode;
CG_DEBUG( "cepgen_structure_functions" ) << sf_mode;
StructureFunctions& val = ( *StructureFunctionsBuilder::get( sf_mode ) )( q2, xbj );
f2 = val.F2;
fl = val.FL;
}
}
diff --git a/CepGen/IO/GridHandler.h b/CepGen/IO/GridHandler.h
new file mode 100644
index 0000000..3f217c5
--- /dev/null
+++ b/CepGen/IO/GridHandler.h
@@ -0,0 +1,182 @@
+#ifndef CepGen_IO_GridHandler_h
+#define CepGen_IO_GridHandler_h
+
+#include <gsl/gsl_version.h>
+#ifdef GSL_MAJOR_VERSION
+# if GSL_MAJOR_VERSION > 2 || ( GSL_MAJOR_VERSION == 2 && GSL_MINOR_VERSION >= 1 )
+# include <gsl/gsl_interp2d.h>
+# include <gsl/gsl_spline2d.h>
+# define GOOD_GSL 1
+# endif
+#endif
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_math.h>
+
+#include "CepGen/Core/Exception.h"
+#include <set>
+
+namespace CepGen
+{
+ class StructureFunctions;
+ enum struct GridType
+ {
+ kLinear = 0,
+ kLogarithmic = 1,
+ kSquare = 2
+ };
+ /// A generic class for 2-dimensional grid interpolation
+ /// \param N Number of values handled per point
+ template <size_t N=1>
+ class GridHandler
+ {
+ public:
+ /// A point and its associated value(s)
+ struct grid_t
+ {
+ /// First coordinate
+ float x;
+ /// Second coordinate
+ float y;
+ /// (list of) value(s) for this point
+ double value[N];
+ };
+
+ public:
+ explicit GridHandler( const GridType& grid_type ) :
+#ifdef GOOD_GSL
+ xacc_( gsl_interp_accel_alloc(), gsl_interp_accel_free ),
+ yacc_( gsl_interp_accel_alloc(), gsl_interp_accel_free ),
+#endif
+ grid_type_( grid_type )
+ {}
+ ~GridHandler() {
+#ifdef GOOD_GSL
+ for ( unsigned short i = 0; i < splines_.size(); ++i ) {
+ gsl_spline2d_free( splines_[i] );
+ if ( values_[i] )
+ delete[] values_[i];
+ }
+#endif
+ }
+
+ /// Interpolate a point to a given coordinate
+ grid_t eval( double x, double y ) const {
+ grid_t ev;
+ if ( grid_type_ == GridType::kLogarithmic ) {
+ x = log10( x );
+ y = log10( y );
+ }
+ else if ( grid_type_ == GridType::kSquare ) {
+ x *= x;
+ y *= y;
+ }
+#ifdef GOOD_GSL
+ unsigned short i = 0;
+ for ( auto& sp : splines_ )
+ if ( gsl_spline2d_eval_e( sp, x, y, xacc_.get(), yacc_.get(), &ev.value[i++] ) != GSL_SUCCESS )
+ CG_WARNING( "GridHandler" )
+ << "Failed to evaluate the grid value "
+ << "for x = " << x << " / y = " << y << ".";
+#else
+ double x_1 = -1., x_2 = -1., y_1 = -1., y_2 = -1.;
+ std::array<double,N> ext_11, ext_12, ext_21, ext_22;
+ for ( unsigned short i = 0; i < x_vals_.size()-1; ++i ) {
+ if ( x >= x_vals_.at( i ) && x < x_vals_.at( i+1 ) ) {
+ x_1 = x_vals_.at( i );
+ x_2 = x_vals_.at( i+1 );
+ break;
+ }
+ }
+ for ( unsigned int i = 0; i < y_vals_.size()-1; ++i ) {
+ if ( y >= y_vals_.at( i ) && y < y_vals_.at( i+1 ) ) {
+ y_1 = y_vals_.at( i );
+ y_2 = y_vals_.at( i+1 );
+ break;
+ }
+ }
+ for ( const auto& val : values_raw_ ) {
+ if ( x_1 == val.x && y_1 == val.y )
+ for ( unsigned short i = 0; i < ext_11.size(); ++i )
+ ext_11[i] = val.value[i];
+ if ( x_2 == val.x && y_1 == val.y )
+ for ( unsigned short i = 0; i < ext_12.size(); ++i )
+ ext_12[i] = val.value[i];
+ if ( x_1 == val.x && y_2 == val.y )
+ for ( unsigned short i = 0; i < ext_21.size(); ++i )
+ ext_21[i] = val.value[i];
+ if ( x_2 == val.x && y_2 == val.y )
+ for ( unsigned short i = 0; i < ext_22.size(); ++i )
+ ext_22[i] = val.value[i];
+ }
+ // now that we have the boundaries, we may interpolate
+ const double x2x1 = x_2-x_1;
+ const double y2y1 = y_2-y_1;
+ const double x2x = x_2-x;
+ const double y2y = y_2-y;
+ const double xx1 = x-x_1;
+ const double yy1 = y-y_1;
+ for ( unsigned short i = 0; i < ext_11.size(); ++i )
+ ev.value[i] = 1.0 / ( x2x1*y2y1 ) * ( ext_11[i]*x2x*y2y + ext_21[i]*xx1*y2y + ext_12[i]*x2x*yy1 + ext_22[i]*xx1*yy1 );
+#endif
+ return ev;
+ }
+
+ /// Return the list of values handled in the grid
+ std::vector<grid_t> values() const { return values_raw_; }
+
+ protected:
+ /// Initialise the grid and all useful interpolators/accelerators
+ void initGSL( const std::set<double>& x_vals, const std::set<double>& y_vals ) {
+#ifdef GOOD_GSL
+ gsl_set_error_handler_off();
+ const gsl_interp2d_type* type = gsl_interp2d_bilinear;
+
+ for ( unsigned short i = 0; i < splines_.size(); ++i ) {
+ values_[i] = new double[x_vals.size() * y_vals.size()];
+ splines_[i] = gsl_spline2d_alloc( type, x_vals.size(), y_vals.size() );
+ }
+
+ // second loop over all points to populate the grid
+ for ( const auto& val : values_raw_ ) {
+ // retrieve the index of the Q2/xbj bin in the set
+ const unsigned short id_x = std::distance( x_vals.begin(), x_vals.lower_bound( log10( val.x ) ) ),
+ id_y = std::distance( y_vals.begin(), y_vals.lower_bound( log10( val.y ) ) );
+ unsigned short i = 0;
+ for ( auto& sp : splines_ ) {
+ gsl_spline2d_set( sp, values_[i], id_x, id_y, val.value[i] );
+ ++i;
+ }
+ }
+
+ // initialise splines objects
+ std::vector<double> q2_vec( x_vals.begin(), x_vals.end() ), xbj_vec( y_vals.begin(), y_vals.end() );
+ for ( unsigned short i = 0; i < splines_.size(); ++i )
+ gsl_spline2d_init( splines_[i], &q2_vec[0], &xbj_vec[0], values_[i], x_vals.size(), y_vals.size() );
+#else
+ CG_WARNING( "GridHandler" )
+ << "GSL version ≥ 2.1 is required for spline bilinear interpolation.\n\t"
+ << "Version " << GSL_VERSION << " is installed on this system!\n\t"
+ << "Will use a linear approximation instead.\n\t"
+ << "You may check the numerical validity of this approach...";
+#endif
+ }
+
+ /// List of coordinates and associated value(s) in the grid
+ std::vector<grid_t> values_raw_;
+#ifdef GOOD_GSL
+ std::array<gsl_spline2d*,N> splines_;
+ std::unique_ptr<gsl_interp_accel,void(*)( gsl_interp_accel* )> xacc_, yacc_;
+ std::array<double*,N> values_;
+#else
+ /// List of coordinates handled in the grid
+ std::vector<double> x_vals_;
+ /// List of coordinates handled in the grid
+ std::vector<double> y_vals_;
+#endif
+ GridType grid_type_;
+
+ };
+}
+
+#endif
+
diff --git a/CepGen/IO/MSTWGridHandler.cpp b/CepGen/IO/MSTWGridHandler.cpp
deleted file mode 100644
index 37be5a1..0000000
--- a/CepGen/IO/MSTWGridHandler.cpp
+++ /dev/null
@@ -1,222 +0,0 @@
-#include "MSTWGridHandler.h"
-#include "CepGen/Core/Exception.h"
-#include "CepGen/Core/utils.h"
-#include "CepGen/StructureFunctions/StructureFunctions.h"
-
-#include <fstream>
-
-#include <gsl/gsl_errno.h>
-#include <gsl/gsl_math.h>
-
-namespace MSTW
-{
- const unsigned int GridHandler::good_magic = 0x5754534d; // MSTW in ASCII
-
- GridHandler&
- GridHandler::get( const char* filename )
- {
- static GridHandler instance( filename );
- return instance;
- }
-
- GridHandler::GridHandler( const char* filename )
-#ifdef GOOD_GSL
- : xacc_( 0 ), yacc_( 0 ), values_( { { 0, 0 } } )
-#endif
- {
- std::set<double> q2_vals, xbj_vals;
-
- { // file readout part
- std::ifstream file( filename, std::ios::binary | std::ios::in );
- if ( !file.is_open() )
- throw CG_FATAL( "GridHandler" ) << "Impossible to load grid file \"" << filename << "%s\"!";
-
- file.read( reinterpret_cast<char*>( &header_ ), sizeof( header_t ) );
-
- // first checks on the file header
-
- if ( header_.magic != good_magic )
- throw CG_FATAL( "GridHandler" ) << "Wrong magic number retrieved: " << header_.magic << ", expecting " << good_magic << ".";
-
- if ( header_.nucleon != header_t::proton )
- throw CG_FATAL( "GridHandler" ) << "Only proton structure function grids can be retrieved for this purpose!";
-
- // retrieve all points and evaluate grid boundaries
-
- sfval_t val;
- while ( file.read( reinterpret_cast<char*>( &val ), sizeof( sfval_t ) ) ) {
- q2_vals.insert( log10( val.q2 ) );
- xbj_vals.insert( log10( val.xbj ) );
-#ifndef GOOD_GSL
- q2_vals_.emplace_back( val.q2 );
- xbj_vals_.emplace_back( val.xbj );
-#endif
- values_raw_.emplace_back( val );
- }
- file.close();
- }
-
- if ( q2_vals.size() < 2 || xbj_vals.size() < 2 )
- throw CG_FATAL( "GridHandler" ) << "Invalid grid retrieved!";
-
- initGSL( q2_vals, xbj_vals );
-
- {
- std::ostringstream ss_cl, ss_ord, ss_nucl;
- ss_cl << header_.cl;
- ss_ord << header_.order;
- ss_nucl << header_.nucleon;
- CG_INFO( "GridHandler" )
- << "MSTW@" << header_.order << " grid evaluator built "
- << "for " << header_.nucleon << " structure functions (" << header_.cl << ")\n\t"
- << " Q² in range [" << pow( 10., *q2_vals.begin() ) << ":" << pow( 10., *q2_vals.rbegin() ) << "]\n\t"
- << "xBj in range [" << pow( 10., *xbj_vals.begin() ) << ":" << pow( 10., *xbj_vals.rbegin() ) << "].";
- }
- }
-
- GridHandler::~GridHandler()
- {
-#ifdef GOOD_GSL
- for ( unsigned short i = 0; i < num_functions_; ++i ) {
- gsl_spline2d_free( splines_[i] );
- if ( values_[i] ) delete[] values_[i];
- }
- if ( xacc_ ) gsl_interp_accel_free( xacc_ );
- if ( yacc_ ) gsl_interp_accel_free( yacc_ );
-#endif
- }
-
- void
- GridHandler::initGSL( const std::set<double>& q2_vals, const std::set<double>& xbj_vals )
- {
-#ifdef GOOD_GSL
- gsl_set_error_handler_off();
- const gsl_interp2d_type* T = gsl_interp2d_bilinear;
-
- for ( unsigned short i = 0; i < num_functions_; ++i ) {
- values_[i] = new double[q2_vals.size() * xbj_vals.size()];
- splines_[i] = gsl_spline2d_alloc( T, q2_vals.size(), xbj_vals.size() );
- }
- xacc_ = gsl_interp_accel_alloc();
- yacc_ = gsl_interp_accel_alloc();
-
- // second loop over all points to populate the grid
- for ( const auto& val : values_raw_ ) {
- // retrieve the index of the Q2/xbj bin in the set
- const unsigned short id_q2 = std::distance( q2_vals.begin(), q2_vals.lower_bound( log10( val.q2 ) ) ),
- id_xbj = std::distance( xbj_vals.begin(), xbj_vals.lower_bound( log10( val.xbj ) ) );
- gsl_spline2d_set( splines_[F2], values_[F2], id_q2, id_xbj, val.f2 );
- gsl_spline2d_set( splines_[FL], values_[FL], id_q2, id_xbj, val.fl );
- }
-
- // initialise splines objects
- std::vector<double> q2_vec( q2_vals.begin(), q2_vals.end() ), xbj_vec( xbj_vals.begin(), xbj_vals.end() );
- for ( unsigned short i = 0; i < num_functions_; ++i ) {
- gsl_spline2d_init( splines_[i], &q2_vec[0], &xbj_vec[0], values_[i], q2_vals.size(), xbj_vals.size() );
- }
-#else
- CG_WARNING( "GridHandler" )
- << "GSL version ≥ 2.1 is required for spline bilinear interpolation.\n\t"
- << "Version " << GSL_VERSION << " is installed on this system!\n\t"
- << "Will use a linear approximation instead.\n\t"
- << "You may check the numerical validity of this approach...";
-#endif
- }
-
- CepGen::StructureFunctions
- GridHandler::eval( double q2, double xbj ) const
- {
- CepGen::StructureFunctions ev;
-#ifdef GOOD_GSL
- if ( gsl_spline2d_eval_e( splines_[F2], log10( q2 ), log10( xbj ), xacc_, yacc_, &ev.F2 ) != GSL_SUCCESS
- || gsl_spline2d_eval_e( splines_[FL], log10( q2 ), log10( xbj ), xacc_, yacc_, &ev.FL ) != GSL_SUCCESS ) {
- CG_WARNING( "GridHandler" )
- << "Failed to evaluate the structure functions "
- << "for Q² = " << q2 << " GeV² / xbj = " << xbj << ".";
- return ev;
- }
-#else
- double q2_1 = -1., q2_2 = -1., xbj_1 = -1., xbj_2 = -1.;
- double f2_11 = 0., f2_12 = 0., f2_21 = 0., f2_22 = 0.;
- double fl_11 = 0., fl_12 = 0., fl_21 = 0., fl_22 = 0.;
- for ( unsigned int i = 0; i < q2_vals_.size()-1; ++i ) {
- if ( q2 >= q2_vals_.at( i ) && q2 < q2_vals_.at( i+1 ) ) {
- q2_1 = q2_vals_.at( i );
- q2_2 = q2_vals_.at( i+1 );
- break;
- }
- }
- for ( unsigned int i = 0; i < xbj_vals_.size()-1; ++i ) {
- if ( xbj >= xbj_vals_.at( i ) && xbj < xbj_vals_.at( i+1 ) ) {
- xbj_1 = xbj_vals_.at( i );
- xbj_2 = xbj_vals_.at( i+1 );
- break;
- }
- }
- for ( const auto& val : values_raw_ ) {
- if ( q2_1 == val.q2 && xbj_1 == val.xbj ) {
- f2_11 = val.f2;
- fl_11 = val.fl;
- }
- if ( q2_2 == val.q2 && xbj_1 == val.xbj ) {
- f2_12 = val.f2;
- fl_12 = val.fl;
- }
- if ( q2_1 == val.q2 && xbj_2 == val.xbj ) {
- f2_21 = val.f2;
- fl_21 = val.fl;
- }
- if ( q2_2 == val.q2 && xbj_2 == val.xbj ) {
- f2_22 = val.f2;
- fl_22 = val.fl;
- }
- }
- const double x2x1 = xbj_2-xbj_1;
- const double y2y1 = q2_2-q2_1;
- const double x2x = xbj_2-xbj;
- const double y2y = q2_2-q2;
- const double yy1 = q2-q2_1;
- const double xx1 = xbj-xbj_1;
- ev.F2 = 1.0 / (x2x1 * y2y1) * ( f2_11*x2x*y2y + f2_21*xx1*y2y + f2_12*x2x*yy1 + f2_22*xx1*yy1 );
- ev.FL = 1.0 / (x2x1 * y2y1) * ( fl_11*x2x*y2y + fl_21*xx1*y2y + fl_12*x2x*yy1 + fl_22*xx1*yy1 );
-#endif
- return ev;
- }
-
- std::ostream&
- operator<<( std::ostream& os, const GridHandler::sfval_t& val )
- {
- return os << Form( "Q² = %.5e GeV²\txbj = %.4f\tF₂ = % .6e\tFL = % .6e", val.q2, val.xbj, val.f2, val.fl );
- }
-
- std::ostream&
- operator<<( std::ostream& os, const GridHandler::header_t::order_t& order )
- {
- switch ( order ) {
- case GridHandler::header_t::lo: return os << "LO";
- case GridHandler::header_t::nlo: return os << "nLO";
- case GridHandler::header_t::nnlo: return os << "nnLO";
- }
- return os;
- }
-
- std::ostream&
- operator<<( std::ostream& os, const GridHandler::header_t::cl_t& cl )
- {
- switch ( cl ) {
- case GridHandler::header_t::cl68: return os << "68% C.L.";
- case GridHandler::header_t::cl95: return os << "95% C.L.";
- }
- return os;
- }
-
- std::ostream&
- operator<<( std::ostream& os, const GridHandler::header_t::nucleon_t& nucl )
- {
- switch ( nucl ) {
- case GridHandler::header_t::proton: return os << "proton";
- case GridHandler::header_t::neutron: return os << "neutron";
- }
- return os;
- }
-}
diff --git a/CepGen/IO/MSTWGridHandler.h b/CepGen/IO/MSTWGridHandler.h
deleted file mode 100644
index 33fcd8c..0000000
--- a/CepGen/IO/MSTWGridHandler.h
+++ /dev/null
@@ -1,75 +0,0 @@
-#ifndef CepGen_IO_MSTWGridHandler_h
-#define CepGen_IO_MSTWGridHandler_h
-
-#include <gsl/gsl_version.h>
-
-#ifdef GSL_MAJOR_VERSION
-#if GSL_MAJOR_VERSION > 2 || ( GSL_MAJOR_VERSION == 2 && GSL_MINOR_VERSION >= 1 )
-#include <gsl/gsl_interp2d.h>
-#include <gsl/gsl_spline2d.h>
-#define GOOD_GSL 1
-#endif
-#endif
-
-#include <array>
-#include <vector>
-#include <set>
-
-namespace CepGen { class StructureFunctions; }
-namespace MSTW
-{
- class GridHandler
- {
- public:
- struct sfval_t {
- float q2, xbj;
- double f2, fl;
- };
- struct header_t {
- enum order_t : unsigned short { lo = 0, nlo = 1, nnlo = 2 };
- enum cl_t : unsigned short { cl68 = 0, cl95 = 1 };
- enum nucleon_t : unsigned short { proton = 1, neutron = 2 };
- unsigned int magic;
- order_t order;
- cl_t cl;
- nucleon_t nucleon;
- };
-
- public:
- static GridHandler& get( const char* filename = "External/F2_Luxlike_fit/mstw_f2_scan_nnlo.dat" );
- ~GridHandler();
-
- CepGen::StructureFunctions eval( double q2, double xbj ) const;
-
- header_t header() const { return header_; }
- std::vector<sfval_t> values() const { return values_raw_; }
-
- private:
- explicit GridHandler( const char* );
- void initGSL( const std::set<double>& q2_vals, const std::set<double>& xbj_vals );
-
- enum spline_type { F2 = 0, FL = 1, num_functions_ };
- static const unsigned int good_magic;
-
- header_t header_;
- std::vector<sfval_t> values_raw_;
-#ifdef GOOD_GSL
- std::array<gsl_spline2d*,2> splines_;
- gsl_interp_accel* xacc_, *yacc_;
- std::array<double*,2> values_;
-#else
- std::vector<double> xbj_vals_, q2_vals_;
-#endif
-
- public:
- GridHandler( const GridHandler& ) = delete;
- void operator=( const GridHandler& ) = delete;
- };
- std::ostream& operator<<( std::ostream&, const GridHandler::sfval_t& );
- std::ostream& operator<<( std::ostream&, const GridHandler::header_t::order_t& );
- std::ostream& operator<<( std::ostream&, const GridHandler::header_t::cl_t& );
- std::ostream& operator<<( std::ostream&, const GridHandler::header_t::nucleon_t& );
-}
-
-#endif
-
diff --git a/CepGen/StructureFunctions/MSTWGrid.cpp b/CepGen/StructureFunctions/MSTWGrid.cpp
new file mode 100644
index 0000000..1d7292a
--- /dev/null
+++ b/CepGen/StructureFunctions/MSTWGrid.cpp
@@ -0,0 +1,119 @@
+#include "CepGen/StructureFunctions/MSTWGrid.h"
+#include "CepGen/Core/Exception.h"
+#include "CepGen/Core/utils.h"
+#include "CepGen/StructureFunctions/StructureFunctions.h"
+
+#include <fstream>
+#include <atomic>
+
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_math.h>
+
+namespace MSTW
+{
+ const unsigned int Grid::good_magic = 0x5754534d; // MSTW in ASCII
+
+ Grid&
+ Grid::get( const char* filename )
+ {
+ Parameterisation p;
+ p.grid_path = filename;
+ static Grid instance( p );
+ return instance;
+ }
+
+ Grid::Grid( const Parameterisation& param ) :
+ CepGen::StructureFunctions( CepGen::SF::Type::MSTWgrid ),
+ CepGen::GridHandler<2>( CepGen::GridType::kLogarithmic ),
+ params( param )
+ {
+ std::set<double> q2_vals, xbj_vals;
+
+ { // file readout part
+ std::ifstream file( params.grid_path, std::ios::binary | std::ios::in );
+ if ( !file.is_open() )
+ throw CG_FATAL( "Grid" ) << "Impossible to load grid file \"" << params.grid_path << "\"!";
+
+ file.read( reinterpret_cast<char*>( &header_ ), sizeof( header_t ) );
+
+ // first checks on the file header
+
+ if ( header_.magic != good_magic )
+ throw CG_FATAL( "Grid" ) << "Wrong magic number retrieved: " << header_.magic << ", expecting " << good_magic << ".";
+
+ if ( header_.nucleon != header_t::proton )
+ throw CG_FATAL( "Grid" ) << "Only proton structure function grids can be retrieved for this purpose!";
+
+ // retrieve all points and evaluate grid boundaries
+
+ sfval_t val;
+ while ( file.read( reinterpret_cast<char*>( &val ), sizeof( sfval_t ) ) ) {
+ q2_vals.insert( log10( val.x ) );
+ xbj_vals.insert( log10( val.y ) );
+#ifndef GOOD_GSL
+ x_vals_.emplace_back( val.x );
+ y_vals_.emplace_back( val.y );
+#endif
+ values_raw_.emplace_back( val );
+ }
+ file.close();
+ }
+
+ if ( q2_vals.size() < 2 || xbj_vals.size() < 2 )
+ throw CG_FATAL( "Grid" ) << "Invalid grid retrieved!";
+
+ initGSL( q2_vals, xbj_vals );
+
+ CG_INFO( "Grid" )
+ << "MSTW@" << header_.order << " grid evaluator built "
+ << "for " << header_.nucleon << " structure functions (" << header_.cl << ")\n\t"
+ << " Q² in range [" << pow( 10., *q2_vals.begin() ) << ":" << pow( 10., *q2_vals.rbegin() ) << "]\n\t"
+ << "xBj in range [" << pow( 10., *xbj_vals.begin() ) << ":" << pow( 10., *xbj_vals.rbegin() ) << "].";
+ }
+
+ Grid&
+ Grid::operator()( double q2, double xbj )
+ {
+ const sfval_t val = CepGen::GridHandler<2>::eval( q2, xbj );
+ F2 = val.value[0];
+ FL = val.value[1];
+ return *this;
+ }
+
+ std::ostream&
+ operator<<( std::ostream& os, const sfval_t& val )
+ {
+ return os << Form( "Q² = %.5e GeV²\txbj = %.4f\tF₂ = % .6e\tFL = % .6e", val.x, val.y, val.value[0], val.value[1] );
+ }
+
+ std::ostream&
+ operator<<( std::ostream& os, const Grid::header_t::order_t& order )
+ {
+ switch ( order ) {
+ case Grid::header_t::lo: return os << "LO";
+ case Grid::header_t::nlo: return os << "nLO";
+ case Grid::header_t::nnlo: return os << "nnLO";
+ }
+ return os;
+ }
+
+ std::ostream&
+ operator<<( std::ostream& os, const Grid::header_t::cl_t& cl )
+ {
+ switch ( cl ) {
+ case Grid::header_t::cl68: return os << "68% C.L.";
+ case Grid::header_t::cl95: return os << "95% C.L.";
+ }
+ return os;
+ }
+
+ std::ostream&
+ operator<<( std::ostream& os, const Grid::header_t::nucleon_t& nucl )
+ {
+ switch ( nucl ) {
+ case Grid::header_t::proton: return os << "proton";
+ case Grid::header_t::neutron: return os << "neutron";
+ }
+ return os;
+ }
+}
diff --git a/CepGen/StructureFunctions/MSTWGrid.h b/CepGen/StructureFunctions/MSTWGrid.h
new file mode 100644
index 0000000..b2f5484
--- /dev/null
+++ b/CepGen/StructureFunctions/MSTWGrid.h
@@ -0,0 +1,61 @@
+#ifndef CepGen_StructureFunctions_MSTWGrid_h
+#define CepGen_StructureFunctions_MSTWGrid_h
+
+#include "CepGen/IO/GridHandler.h"
+#include "CepGen/StructureFunctions/StructureFunctions.h"
+
+/// Martin-Stirling-Thorne-Watt PDFs structure functions
+namespace MSTW
+{
+ /// \note x/y = Q2/xbj given by the parent
+ typedef CepGen::GridHandler<2>::grid_t sfval_t;
+ std::ostream& operator<<( std::ostream&, const sfval_t& );
+
+ /// A \f$F_{2,L}\f$ grid interpolator
+ class Grid : public CepGen::StructureFunctions, private CepGen::GridHandler<2>
+ {
+ public:
+ /// Grid header information as parsed from the file
+ struct header_t {
+ enum order_t : unsigned short { lo = 0, nlo = 1, nnlo = 2 };
+ enum cl_t : unsigned short { cl68 = 0, cl95 = 1 };
+ enum nucleon_t : unsigned short { proton = 1, neutron = 2 };
+ unsigned int magic;
+ order_t order;
+ cl_t cl;
+ nucleon_t nucleon;
+ };
+ struct Parameterisation {
+ Parameterisation() : grid_path( "External/F2_Luxlike_fit/mstw_f2_scan_nnlo.dat" ) {}
+ std::string grid_path;
+ };
+
+ public:
+ /// Retrieve the grid interpolator (singleton)
+ static Grid& get( const char* path = "External/F2_Luxlike_fit/mstw_f2_scan_nnlo.dat" );
+
+ /// Compute the structure functions at a given \f$Q^2/x_{\rm Bj}\f$
+ Grid& operator()( double q2, double xbj ) override;
+ /// Retrieve the grid's header information
+ header_t header() const { return header_; }
+ Parameterisation params;
+
+ public:
+ Grid( const Grid& ) = delete;
+ void operator=( const GridHandler& ) = delete;
+
+ private:
+ explicit Grid( const Parameterisation& = Parameterisation() );
+ static const unsigned int good_magic;
+ static std::shared_ptr<Grid> singl_;
+
+ header_t header_;
+ };
+
+ std::ostream& operator<<( std::ostream&, const Grid::header_t::order_t& );
+ std::ostream& operator<<( std::ostream&, const Grid::header_t::cl_t& );
+ std::ostream& operator<<( std::ostream&, const Grid::header_t::nucleon_t& );
+}
+
+#endif
+
diff --git a/CepGen/StructureFunctions/StructureFunctionsBuilder.cpp b/CepGen/StructureFunctions/StructureFunctionsBuilder.cpp
index b0d2b52..1de7159 100644
--- a/CepGen/StructureFunctions/StructureFunctionsBuilder.cpp
+++ b/CepGen/StructureFunctions/StructureFunctionsBuilder.cpp
@@ -1,40 +1,40 @@
#include "CepGen/StructureFunctions/StructureFunctionsBuilder.h"
#include "CepGen/StructureFunctions/ALLM.h"
#include "CepGen/StructureFunctions/BlockDurandHa.h"
#include "CepGen/StructureFunctions/FioreBrasse.h"
#include "CepGen/StructureFunctions/ChristyBosted.h"
#include "CepGen/StructureFunctions/CLAS.h"
#include "CepGen/StructureFunctions/GenericLHAPDF.h"
#include "CepGen/StructureFunctions/SuriYennie.h"
#include "CepGen/StructureFunctions/SzczurekUleshchenko.h"
#include "CepGen/StructureFunctions/Schaefer.h"
-#include "CepGen/IO/MSTWGridHandler.h"
+#include "CepGen/StructureFunctions/MSTWGrid.h"
namespace CepGen
{
std::shared_ptr<StructureFunctions>
StructureFunctionsBuilder::get( const SF::Type& sf_type )
{
switch ( sf_type ) {
case SF::Type::Electron:
case SF::Type::ElasticProton:
default: return std::make_shared<StructureFunctions>();
case SF::Type::SzczurekUleshchenko: return std::make_shared<SF::SzczurekUleshchenko>();
case SF::Type::SuriYennie: return std::make_shared<SF::SuriYennie>();
case SF::Type::FioreBrasse: return std::make_shared<SF::FioreBrasse>();
case SF::Type::ChristyBosted: return std::make_shared<SF::ChristyBosted>();
case SF::Type::CLAS: return std::make_shared<SF::CLAS>();
case SF::Type::BlockDurandHa: return std::make_shared<SF::BlockDurandHa>();
case SF::Type::ALLM91: return std::make_shared<SF::ALLM>( SF::ALLM::Parameterisation::allm91() );
case SF::Type::ALLM97: return std::make_shared<SF::ALLM>( SF::ALLM::Parameterisation::allm97() );
case SF::Type::GD07p: return std::make_shared<SF::ALLM>( SF::ALLM::Parameterisation::gd07p() );
case SF::Type::GD11p: return std::make_shared<SF::ALLM>( SF::ALLM::Parameterisation::gd11p() );
case SF::Type::Schaefer: return std::make_shared<SF::Schaefer>();
case SF::Type::GenericLHAPDF: return std::make_shared<SF::GenericLHAPDF>();
//--- particular case for the MSTW grid as we are dealing
//--- with a singleton ; hence, no deleter is needed!
case SF::Type::MSTWgrid: return std::shared_ptr<MSTW::Grid>( &MSTW::Grid::get(), [=]( MSTW::Grid* ){} );
}
}
}