diff --git a/Cards/Config/Core.py b/Cards/Config/Core.py
index f239c64..82be4de 100644
--- a/Cards/Config/Core.py
+++ b/Cards/Config/Core.py
@@ -1,173 +1,176 @@
#!/usr/bin/env python
class PrintHelper(object):
_indent = 0
_indent_size = 4
def indent(self):
self._indent += self._indent_size
def unindent(self):
self._indent -= self._indent_size
def indentation(self):
return ' '*self._indent
class Parameters(dict):
'''A raw list of steering parameters'''
__getattr__ = dict.get
__setattr__ = dict.__setitem__
__delattr__ = dict.__delitem__
def __init__(self, *args, **kwargs):
self.update(*args, **kwargs)
super(Parameters, self).__init__(*args, **kwargs)
def __deepcopy__(self, memo):
from copy import deepcopy
return Parameters([(deepcopy(k, memo), deepcopy(v, memo)) for k, v in self.items()])
def dump(self, printer=PrintHelper()):
out = self.__class__.__name__+'(\n'
for k, v in self.items():
printer.indent()
out += ('%s%s = ' % (printer.indentation(), k))
if v.__class__.__name__ not in ['Parameters', 'Module']:
out += v.__repr__()
else:
out += v.dump(printer)
out += ',\n'
printer.unindent()
out += printer.indentation()+')'
return out
def __repr__(self):
return self.dump()
def clone(self, *args, **kwargs):
from copy import deepcopy
out = deepcopy(self)
for k in kwargs:
out[k] = kwargs.get(k)
return type(self)(out)
def load(self, mod):
mod = mod.replace('/', '.')
module = __import__(mod)
self.extend(sys.modules[mod])
class Module(Parameters):
'''A named parameters set to steer a generic module'''
def __init__(self, name, *args, **kwargs):
super(Module, self).__init__(*args, **kwargs)
self.mod_name = name
def __len__(self):
return dict.__len__(self)-1 # discard the name key
def dump(self, printer=PrintHelper()):
out = self.__class__.__name__+'('+self.mod_name.__repr__()+'\n'
mod_repr = self.clone('')
mod_repr.pop('mod_name', None)
for k, v in mod_repr.items():
printer.indent()
out += ('%s%s = ' % (printer.indentation(), k))
if v.__class__.__name__ not in ['Parameters', 'Module']:
out += v.__repr__()
else:
out += v.dump(printer)
out += ',\n'
printer.unindent()
out += printer.indentation()+')'
return out
def __repr__(self):
return self.dump()
def clone(self, name, **kwargs):
out = Parameters(self).clone(**kwargs)
out.mod_name = name
return out
class Logging:
'''Logging verbosity'''
Nothing = 0
Error = 1
Warning = 2
Information = 3
Debug = 4
DebugInsideLoop = 5
class StructureFunctions:
'''Types of structure functions supported'''
Electron = Parameters(id=1)
ElasticProton = Parameters(id=2)
SuriYennie = Parameters(id=11)
SzczurekUleshchenko = Parameters(id=12)
BlockDurandHa = Parameters(id=13)
FioreBrasse = Parameters(id=101)
ChristyBosted = Parameters(id=102)
CLAS = Parameters(id=103)
ALLM91 = Parameters(id=201)
ALLM97 = Parameters(id=202)
GD07p = Parameters(id=203)
GD11p = Parameters(id=204)
- MSTWgrid = Parameters(id=205)
+ MSTWgrid = Parameters(
+ id = 205,
+ gridPath = 'External/F2_Luxlike_fit/mstw_f2_scan_nnlo.dat',
+ )
LUXlike = Parameters(
id = 301,
#Q2cut = 10.,
#W2limits = (4.,1.),
continuumSF = GD11p,
resonancesSF = ChristyBosted,
)
GenericLHAPDF = Parameters(
id = 401,
pdfSet = 'LUXqed17_plus_PDF4LHC15_nnlo_100',
numFlavours = 5,
)
class KTFlux:
PhotonElastic = 0
PhotonElasticBudnev = 10
PhotonInelasticBudnev = 1
PhotonInelasticBudnev = 11
PhotonElasticHI = 100
GluonKMR = 20
class ProcessMode:
'''Types of processes supported'''
ElectronProton = 0
ElasticElastic = 1
ElasticInelastic = 2
InelasticElastic = 3
InelasticInelastic = 4
ProtonElectron = 5
ElectronElectron = 6
if __name__ == '__main__':
import unittest
class TestTypes(unittest.TestCase):
def testModules(self):
mod = Module('empty')
self.assertEqual(len(mod), 0)
mod.param1 = 'foo'
self.assertEqual(len(mod), 1)
# playing with modules clones
mod_copy = mod.clone('notEmpty', param1 = 'boo', param2 = 'bar')
self.assertEqual(mod.param1, 'foo')
self.assertEqual(mod_copy.param1, 'boo')
self.assertEqual(mod_copy.param2, 'bar')
self.assertEqual(mod.param1+mod_copy.param2, 'foobar')
def testParameters(self):
params = Parameters(
first = 'foo',
second = 'bar',
third = 42,
fourth = (1, 2),
)
params_copy = params.clone(
second = 'bak',
)
self.assertEqual(len(params), 4)
self.assertEqual(params.first, params['first'])
self.assertEqual(params['second'], 'bar')
self.assertTrue(int(params.third) == params.third)
self.assertEqual(len(params.fourth), 2)
self.assertEqual(params.second, 'bar')
# playing with parameters clones
self.assertEqual(params_copy.second, 'bak')
# check that the clone does not change value if the origin does
# (i.e. we indeed have a deep copy and not a shallow one...)
params.third = 43
self.assertEqual(params.third, 43)
self.assertEqual(params_copy.third, 42)
unittest.main()
diff --git a/CepGen/Cards/PythonHandler.cpp b/CepGen/Cards/PythonHandler.cpp
index 61dc4ad..b2e1c6b 100644
--- a/CepGen/Cards/PythonHandler.cpp
+++ b/CepGen/Cards/PythonHandler.cpp
@@ -1,699 +1,704 @@
#include "PythonHandler.h"
#include "CepGen/Core/Exception.h"
#ifdef PYTHON
#include "CepGen/Core/TamingFunction.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Processes/GamGamLL.h"
#include "CepGen/Processes/PPtoFF.h"
#include "CepGen/Processes/PPtoWW.h"
#include "CepGen/Processes/FortranKTProcess.h"
#include "CepGen/StructureFunctions/StructureFunctionsBuilder.h"
#include "CepGen/StructureFunctions/GenericLHAPDF.h"
+#include "CepGen/StructureFunctions/MSTWGrid.h"
#include "CepGen/StructureFunctions/Schaefer.h"
#include "CepGen/Hadronisers/Pythia8Hadroniser.h"
#include <algorithm>
#include <frameobject.h>
extern "C"
{
extern void nucl_to_ff_( double& );
}
#if PY_MAJOR_VERSION < 3
# define PYTHON2
#endif
namespace CepGen
{
namespace Cards
{
//----- specialization for CepGen input cards
PythonHandler::PythonHandler( const char* file )
{
setenv( "PYTHONPATH", ".:..:Cards", 1 );
std::string filename = getPythonPath( file );
const size_t fn_len = filename.length()+1;
//Py_DebugFlag = 1;
//Py_VerboseFlag = 1;
#ifdef PYTHON2
char* sfilename = new char[fn_len];
snprintf( sfilename, fn_len, "%s", filename.c_str() );
#else
wchar_t* sfilename = new wchar_t[fn_len];
swprintf( sfilename, fn_len, L"%s", filename.c_str() );
#endif
if ( sfilename )
Py_SetProgramName( sfilename );
Py_InitializeEx( 1 );
if ( sfilename )
delete [] sfilename;
if ( !Py_IsInitialized() )
throw CG_FATAL( "PythonHandler" ) << "Failed to initialise the Python cards parser!";
CG_INFO( "PythonHandler" )
<< "Initialised the Python cards parser\n\t"
<< "Python version: " << Py_GetVersion() << "\n\t"
<< "Platform: " << Py_GetPlatform() << ".";
PyObject* cfg = PyImport_ImportModule( filename.c_str() ); // new
if ( !cfg )
throwPythonError( Form( "Failed to parse the configuration card %s", file ) );
PyObject* process = PyObject_GetAttrString( cfg, PROCESS_NAME ); // new
if ( !process )
throwPythonError( Form( "Failed to extract a \"%s\" keyword from the configuration card %s", PROCESS_NAME, file ) );
//--- type of process to consider
PyObject* pproc_name = getElement( process, MODULE_NAME ); // borrowed
if ( !pproc_name )
throwPythonError( Form( "Failed to extract the process name from the configuration card %s", file ) );
const std::string proc_name = decode( pproc_name );
if ( proc_name == "lpair" )
params_.setProcess( new Process::GamGamLL );
else if ( proc_name == "pptoll" || proc_name == "pptoff" )
params_.setProcess( new Process::PPtoFF );
else if ( proc_name == "pptoww" )
params_.setProcess( new Process::PPtoWW );
else if ( proc_name == "patoll" )
params_.setProcess( new Process::FortranKTProcess( "nucltoff", "(p/A)(p/A) ↝ (g/ɣ)ɣ → f⁺f¯", nucl_to_ff_ ) );
else throw CG_FATAL( "PythonHandler" ) << "Unrecognised process: " << proc_name << ".";
//--- process mode
fillParameter( process, "mode", (int&)params_.kinematics.mode );
//--- process kinematics
PyObject* pin_kinematics = getElement( process, "inKinematics" ); // borrowed
if ( pin_kinematics )
parseIncomingKinematics( pin_kinematics );
PyObject* pout_kinematics = getElement( process, "outKinematics" ); // borrowed
if ( pout_kinematics )
parseOutgoingKinematics( pout_kinematics );
Py_CLEAR( process );
PyObject* plog = PyObject_GetAttrString( cfg, "logger" ); // new
if ( plog ) {
parseLogging( plog );
Py_CLEAR( plog );
}
//--- hadroniser parameters
PyObject* phad = PyObject_GetAttrString( cfg, "hadroniser" ); // new
if ( phad ) {
parseHadroniser( phad );
Py_CLEAR( phad );
}
//--- generation parameters
PyObject* pint = PyObject_GetAttrString( cfg, "integrator" ); // new
if ( pint ) {
parseIntegrator( pint );
Py_CLEAR( pint );
}
PyObject* pgen = PyObject_GetAttrString( cfg, "generator" ); // new
if ( pgen ) {
parseGenerator( pgen );
Py_CLEAR( pgen );
}
//--- taming functions
PyObject* ptam = PyObject_GetAttrString( cfg, "tamingFunctions" ); // new
if ( ptam ) {
parseTamingFunctions( ptam );
Py_CLEAR( ptam );
}
PyObject* pkmr_path = PyObject_GetAttrString( cfg, "kmrGridPath" ); // new
if ( pkmr_path ) {
//strncpy( Constants::kmr_grid_path, decode( pkmr_path ).c_str(), 256 );
params_.kinematics.kmr_grid_path = decode( pkmr_path );
Py_CLEAR( pkmr_path );
}
//--- finalisation
Py_CLEAR( cfg );
}
PythonHandler::~PythonHandler()
{
if ( Py_IsInitialized() )
Py_Finalize();
}
void
PythonHandler::parseIncomingKinematics( PyObject* kin )
{
PyObject* ppz = getElement( kin, "pz" ); // borrowed
if ( ppz && PyTuple_Check( ppz ) && PyTuple_Size( ppz ) == 2 ) {
double pz0 = PyFloat_AsDouble( PyTuple_GetItem( ppz, 0 ) );
double pz1 = PyFloat_AsDouble( PyTuple_GetItem( ppz, 1 ) );
params_.kinematics.incoming_beams.first.pz = pz0;
params_.kinematics.incoming_beams.second.pz = pz1;
}
double sqrt_s = -1.;
fillParameter( kin, "cmEnergy", sqrt_s );
if ( sqrt_s != -1. )
params_.kinematics.setSqrtS( sqrt_s );
PyObject* psf = getElement( kin, "structureFunctions" ); // borrowed
if ( psf )
parseStructureFunctions( psf );
std::vector<int> kt_fluxes;
fillParameter( kin, "ktFluxes", kt_fluxes );
if ( kt_fluxes.size() > 0 )
params_.kinematics.incoming_beams.first.kt_flux = kt_fluxes.at( 0 );
if ( kt_fluxes.size() > 1 )
params_.kinematics.incoming_beams.second.kt_flux = kt_fluxes.at( 1 );
std::vector<int> hi_beam1, hi_beam2;
fillParameter( kin, "heavyIonA", hi_beam1 );
if ( hi_beam1.size() == 2 )
params_.kinematics.incoming_beams.first.hi = Kinematics::HeavyIon{ (unsigned short)hi_beam1[0], (unsigned short)hi_beam1[1] };
fillParameter( kin, "heavyIonB", hi_beam2 );
if ( hi_beam2.size() == 2 )
params_.kinematics.incoming_beams.second.hi = Kinematics::HeavyIon{ (unsigned short)hi_beam2[0], (unsigned short)hi_beam2[1] };
}
void
PythonHandler::parseStructureFunctions( PyObject* psf )
{
int str_fun = 0;
fillParameter( psf, "id", str_fun );
params_.kinematics.structure_functions.reset( StructureFunctionsBuilder::get( (SF::Type)str_fun ) );
switch( (SF::Type)str_fun ) {
case SF::Type::GenericLHAPDF: {
auto sf = dynamic_cast<SF::GenericLHAPDF*>( params_.kinematics.structure_functions.get() );
fillParameter( psf, "pdfSet", sf->params.pdf_set );
fillParameter( psf, "numFlavours", (unsigned int&)sf->params.num_flavours );
} break;
+ case SF::Type::MSTWgrid: {
+ auto sf = dynamic_cast<MSTW::Grid*>( params_.kinematics.structure_functions.get() );
+ fillParameter( psf, "gridPath", sf->params.grid_path );
+ } break;
case SF::Type::Schaefer: {
auto sf = dynamic_cast<SF::Schaefer*>( params_.kinematics.structure_functions.get() );
fillParameter( psf, "Q2cut", sf->params.q2_cut );
std::vector<double> w2_lims;
fillParameter( psf, "W2limits", w2_lims );
if ( w2_lims.size() != 0 ) {
if ( w2_lims.size() != 2 )
throwPythonError( Form( "Invalid size for W2limits attribute: %d != 2!", w2_lims.size() ) );
else {
sf->params.w2_lo = *std::min_element( w2_lims.begin(), w2_lims.end() );
sf->params.w2_hi = *std::max_element( w2_lims.begin(), w2_lims.end() );
}
}
PyObject* pcsf = getElement( psf, "continuumSF" ); // borrowed
if ( pcsf )
fillParameter( pcsf, "id", sf->params.cont_model );
PyObject* prsf = getElement( psf, "resonancesSF" ); // borrowed
if ( prsf )
fillParameter( prsf, "id", sf->params.res_model );
fillParameter( psf, "higherTwist", (bool&)sf->params.higher_twist );
} break;
default: break;
}
}
void
PythonHandler::parseOutgoingKinematics( PyObject* kin )
{
PyObject* ppair = getElement( kin, "pair" ); // borrowed
if ( ppair ) {
if ( isInteger( ppair ) ) {
PDG pair = (PDG)asInteger( ppair );
params_.kinematics.central_system = { pair, pair };
}
else if ( PyTuple_Check( ppair ) ) {
if ( PyTuple_Size( ppair ) != 2 )
throw CG_FATAL( "PythonHandler" ) << "Invalid value for in_kinematics.pair!";
PDG pair1 = (PDG)asInteger( PyTuple_GetItem( ppair, 0 ) );
PDG pair2 = (PDG)asInteger( PyTuple_GetItem( ppair, 1 ) );
params_.kinematics.central_system = { pair1, pair2 };
}
}
PyObject* pparts = getElement( kin, "minFinalState" ); // borrowed
if ( pparts && PyTuple_Check( pparts ) )
for ( unsigned short i = 0; i < PyTuple_Size( pparts ); ++i )
params_.kinematics.minimum_final_state.emplace_back( (PDG)asInteger( PyTuple_GetItem( pparts, i ) ) );
PyObject* pcuts = getElement( kin, "cuts" ); // borrowed
if ( pcuts )
parseParticlesCuts( pcuts );
// for LPAIR/collinear matrix elements
fillLimits( kin, "q2", params_.kinematics.cuts.initial.q2 );
// for the kT factorised matrix elements
fillLimits( kin, "qt", params_.kinematics.cuts.initial.qt );
fillLimits( kin, "phiqt", params_.kinematics.cuts.initial.phi_qt );
fillLimits( kin, "ptdiff", params_.kinematics.cuts.central.pt_diff );
fillLimits( kin, "phiptdiff", params_.kinematics.cuts.central.phi_pt_diff );
fillLimits( kin, "rapiditydiff", params_.kinematics.cuts.central.rapidity_diff );
// generic phase space limits
fillLimits( kin, "rapidity", params_.kinematics.cuts.central.rapidity_single );
fillLimits( kin, "eta", params_.kinematics.cuts.central.eta_single );
fillLimits( kin, "pt", params_.kinematics.cuts.central.pt_single );
fillLimits( kin, "ptsum", params_.kinematics.cuts.central.pt_sum );
fillLimits( kin, "invmass", params_.kinematics.cuts.central.mass_sum );
fillLimits( kin, "mx", params_.kinematics.cuts.remnants.mass_single );
}
void
PythonHandler::parseParticlesCuts( PyObject* cuts )
{
if ( !PyDict_Check( cuts ) )
throwPythonError( "Particle cuts object should be a dictionary!" );
PyObject* pkey = nullptr, *pvalue = nullptr;
Py_ssize_t pos = 0;
while ( PyDict_Next( cuts, &pos, &pkey, &pvalue ) ) {
const PDG pdg = (PDG)asInteger( pkey );
fillLimits( pvalue, "pt", params_.kinematics.cuts.central_particles[pdg].pt_single );
fillLimits( pvalue, "energy", params_.kinematics.cuts.central_particles[pdg].energy_single );
fillLimits( pvalue, "eta", params_.kinematics.cuts.central_particles[pdg].eta_single );
fillLimits( pvalue, "rapidity", params_.kinematics.cuts.central_particles[pdg].rapidity_single );
}
}
void
PythonHandler::parseLogging( PyObject* log )
{
fillParameter( log, "level", (int&)Logger::get().level );
std::vector<std::string> enabled_modules;
fillParameter( log, "enabledModules", enabled_modules );
for ( const auto& mod : enabled_modules )
Logger::get().addExceptionRule( mod );
}
void
PythonHandler::parseIntegrator( PyObject* integr )
{
if ( !PyDict_Check( integr ) )
throwPythonError( "Integrator object should be a dictionary!" );
PyObject* palgo = getElement( integr, MODULE_NAME ); // borrowed
if ( !palgo )
throwPythonError( "Failed to retrieve the integration algorithm name!" );
std::string algo = decode( palgo );
if ( algo == "plain" )
params_.integrator.type = Integrator::Type::plain;
else if ( algo == "Vegas" ) {
params_.integrator.type = Integrator::Type::Vegas;
fillParameter( integr, "alpha", (double&)params_.integrator.vegas.alpha );
fillParameter( integr, "iterations", params_.integrator.vegas.iterations );
fillParameter( integr, "mode", (int&)params_.integrator.vegas.mode );
fillParameter( integr, "verbosity", (int&)params_.integrator.vegas.verbose );
std::string vegas_logging_output = "cerr";
fillParameter( integr, "loggingOutput", vegas_logging_output );
if ( vegas_logging_output == "cerr" )
// redirect all debugging information to the error stream
params_.integrator.vegas.ostream = stderr;
else if ( vegas_logging_output == "cout" )
// redirect all debugging information to the standard stream
params_.integrator.vegas.ostream = stdout;
else
params_.integrator.vegas.ostream = fopen( vegas_logging_output.c_str(), "w" );
}
else if ( algo == "MISER" ) {
params_.integrator.type = Integrator::Type::MISER;
fillParameter( integr, "estimateFraction", (double&)params_.integrator.miser.estimate_frac );
fillParameter( integr, "minCalls", params_.integrator.miser.min_calls );
fillParameter( integr, "minCallsPerBisection", params_.integrator.miser.min_calls_per_bisection );
fillParameter( integr, "alpha", (double&)params_.integrator.miser.alpha );
fillParameter( integr, "dither", (double&)params_.integrator.miser.dither );
}
else
throwPythonError( Form( "Invalid integration algorithm: %s", algo.c_str() ) );
fillParameter( integr, "numFunctionCalls", params_.integrator.ncvg );
fillParameter( integr, "seed", (unsigned long&)params_.integrator.rng_seed );
unsigned int rng_engine;
fillParameter( integr, "rngEngine", rng_engine );
switch ( rng_engine ) {
case 0: default: params_.integrator.rng_engine = (gsl_rng_type*)gsl_rng_mt19937; break;
case 1: params_.integrator.rng_engine = (gsl_rng_type*)gsl_rng_taus2; break;
case 2: params_.integrator.rng_engine = (gsl_rng_type*)gsl_rng_gfsr4; break;
case 3: params_.integrator.rng_engine = (gsl_rng_type*)gsl_rng_ranlxs0; break;
}
fillParameter( integr, "chiSqCut", params_.integrator.vegas_chisq_cut );
}
void
PythonHandler::parseGenerator( PyObject* gen )
{
if ( !PyDict_Check( gen ) )
throwPythonError( "Generation information object should be a dictionary!" );
params_.generation.enabled = true;
fillParameter( gen, "treat", params_.generation.treat );
fillParameter( gen, "numEvents", params_.generation.maxgen );
fillParameter( gen, "printEvery", params_.generation.gen_print_every );
fillParameter( gen, "numThreads", params_.generation.num_threads );
fillParameter( gen, "numPoints", params_.generation.num_points );
}
void
PythonHandler::parseTamingFunctions( PyObject* tf )
{
if ( !PyList_Check( tf ) )
throwPythonError( "Taming functions list should be a list!" );
for ( Py_ssize_t i = 0; i < PyList_Size( tf ); ++i ) {
PyObject* pit = PyList_GetItem( tf, i ); // borrowed
if ( !pit )
continue;
if ( !PyDict_Check( pit ) )
throwPythonError( Form( "Item %d has invalid type %s", i, pit->ob_type->tp_name ) );
PyObject* pvar = getElement( pit, "variable" ), *pexpr = getElement( pit, "expression" ); // borrowed
params_.taming_functions->add( decode( pvar ).c_str(), decode( pexpr ).c_str() );
}
}
void
PythonHandler::parseHadroniser( PyObject* hadr )
{
if ( !PyDict_Check( hadr ) )
throwPythonError( "Hadroniser object should be a dictionary!" );
PyObject* pname = getElement( hadr, MODULE_NAME ); // borrowed
if ( !pname )
throwPythonError( "Hadroniser name is required!" );
std::string hadr_name = decode( pname );
fillParameter( hadr, "maxTrials", params_.hadroniser_max_trials );
PyObject* pseed = getElement( hadr, "seed" ); // borrowed
long long seed = -1ll;
if ( pseed && isInteger( pseed ) ) {
seed = PyLong_AsLongLong( pseed );
CG_DEBUG( "PythonHandler:hadroniser" ) << "Hadroniser seed set to " << seed;
}
if ( hadr_name == "pythia8" ) {
params_.setHadroniser( new Hadroniser::Pythia8Hadroniser( params_ ) );
std::vector<std::string> config;
auto pythia8 = dynamic_cast<Hadroniser::Pythia8Hadroniser*>( params_.hadroniser() );
pythia8->setSeed( seed );
fillParameter( hadr, "pythiaPreConfiguration", config );
pythia8->readStrings( config );
pythia8->init();
fillParameter( hadr, "pythiaConfiguration", config );
pythia8->readStrings( config );
fillParameter( hadr, "pythiaProcessConfiguration", config );
pythia8->readStrings( config );
}
}
//------------------------------------------------------------------
// Python API helpers
//------------------------------------------------------------------
std::string
PythonHandler::getPythonPath( const char* file )
{
std::string s_filename = file;
s_filename = s_filename.substr( 0, s_filename.find_last_of( "." ) ); // remove the extension
std::replace( s_filename.begin(), s_filename.end(), '/', '.' ); // replace all '/' by '.'
return s_filename;
}
void
PythonHandler::throwPythonError( const std::string& message )
{
PyObject* ptype = nullptr, *pvalue = nullptr, *ptraceback_obj = nullptr;
// retrieve error indicator and clear it to handle ourself the error
PyErr_Fetch( &ptype, &pvalue, &ptraceback_obj );
PyErr_Clear();
// ensure the objects retrieved are properly normalised and point to compatible objects
PyErr_NormalizeException( &ptype, &pvalue, &ptraceback_obj );
std::ostringstream oss; oss << message;
if ( ptype != nullptr ) { // we can start the traceback
oss << "\n\tError: "
#ifdef PYTHON2
<< PyString_AsString( PyObject_Str( pvalue ) ); // deprecated in python v3+
#else
<< PyUnicode_AsUTF8( PyObject_Str( pvalue ) );
#endif
PyTracebackObject* ptraceback = (PyTracebackObject*)ptraceback_obj;
string tabul = "↪ ";
if ( ptraceback != nullptr ) {
while ( ptraceback->tb_next != nullptr ) {
PyFrameObject* pframe = ptraceback->tb_frame;
if ( pframe != nullptr ) {
int line = PyCode_Addr2Line( pframe->f_code, pframe->f_lasti );
#ifdef PYTHON2
const char* filename = PyString_AsString( pframe->f_code->co_filename );
const char* funcname = PyString_AsString( pframe->f_code->co_name );
#else
const char* filename = PyUnicode_AsUTF8( pframe->f_code->co_filename );
const char* funcname = PyUnicode_AsUTF8( pframe->f_code->co_name );
#endif
oss << Form( "\n\t%s%s on %s (line %d)", tabul.c_str(), boldify( funcname ).c_str(), filename, line );
}
else
oss << Form( "\n\t%s issue in line %d", tabul.c_str(), ptraceback->tb_lineno );
tabul = string( " " )+tabul;
ptraceback = ptraceback->tb_next;
}
}
}
Py_Finalize();
throw CG_FATAL( "PythonHandler:error" ) << oss.str();
}
std::string
PythonHandler::decode( PyObject* obj )
{
std::string out;
#ifdef PYTHON2
out = PyString_AsString( obj ); // deprecated in python v3+
#else
PyObject* pstr = PyUnicode_AsEncodedString( obj, "utf-8", "strict" ); // new
if ( !pstr )
throwPythonError( "Failed to decode a Python object!" );
out = PyBytes_AS_STRING( pstr );
Py_CLEAR( pstr );
#endif
return out;
}
PyObject*
PythonHandler::encode( const char* str )
{
PyObject* obj = PyUnicode_FromString( str ); // new
if ( !obj )
throwPythonError( Form( "Failed to encode the following string:\n\t%s", str ) );
return obj;
}
PyObject*
PythonHandler::getElement( PyObject* obj, const char* key )
{
PyObject* pout = nullptr, *nink = encode( key );
if ( !nink )
return pout;
pout = PyDict_GetItem( obj, nink ); // borrowed
Py_CLEAR( nink );
if ( pout )
CG_DEBUG( "PythonHandler:getElement" )
<< "retrieved " << pout->ob_type->tp_name << " element \"" << key << "\" "
<< "from " << obj->ob_type->tp_name << " object\n\t"
<< "new reference count: " << pout->ob_refcnt;
else
CG_DEBUG( "PythonHandler:getElement" )
<< "did not retrieve a valid element \"" << key << "\"";
return pout;
}
void
PythonHandler::fillLimits( PyObject* obj, const char* key, Limits& lim )
{
PyObject* pobj = getElement( obj, key ); // borrowed
if ( !pobj )
return;
if ( !PyTuple_Check( pobj ) )
throw CG_FATAL( "PythonHandler:fillLimits" ) << "Invalid value retrieved for " << key << ".";
if ( PyTuple_Size( pobj ) < 1 )
throw CG_FATAL( "PythonHandler:fillLimits" ) << "Invalid number of values unpacked for " << key << "!";
double min = PyFloat_AsDouble( PyTuple_GetItem( pobj, 0 ) );
lim.min() = min;
if ( PyTuple_Size( pobj ) > 1 ) {
double max = PyFloat_AsDouble( PyTuple_GetItem( pobj, 1 ) );
if ( max != -1 )
lim.max() = max;
}
}
void
PythonHandler::fillParameter( PyObject* parent, const char* key, bool& out )
{
PyObject* pobj = getElement( parent, key ); // borrowed
if ( !pobj )
return;
if ( !PyBool_Check( pobj ) )
throwPythonError( Form( "Object \"%s\" has invalid type %s", key, pobj->ob_type->tp_name ) );
fillParameter( parent, key, (int&)out );
}
void
PythonHandler::fillParameter( PyObject* parent, const char* key, int& out )
{
PyObject* pobj = getElement( parent, key ); // borrowed
if ( !pobj )
return;
#ifdef PYTHON2
if ( !PyInt_Check( pobj ) && !PyBool_Check( pobj ) )
throwPythonError( Form( "Object \"%s\" has invalid type %s", key, pobj->ob_type->tp_name ) );
out = PyInt_AsLong( pobj );
#else
if ( !PyLong_Check( pobj ) && !PyBool_Check( pobj ) )
throwPythonError( Form( "Object \"%s\" has invalid type %s", key, pobj->ob_type->tp_name ) );
out = PyLong_AsLong( pobj );
#endif
}
void
PythonHandler::fillParameter( PyObject* parent, const char* key, unsigned long& out )
{
PyObject* pobj = getElement( parent, key ); // borrowed
if ( !pobj )
return;
if ( !PyLong_Check( pobj )
#ifdef PYTHON2
&& !PyInt_Check( pobj )
#endif
)
throwPythonError( Form( "Object \"%s\" has invalid type %s", key, pobj->ob_type->tp_name ) );
if ( PyLong_Check( pobj ) )
out = PyLong_AsUnsignedLong( pobj );
#ifdef PYTHON2
else if ( PyInt_Check( pobj ) )
out = PyInt_AsUnsignedLongMask( pobj );
#endif
}
void
PythonHandler::fillParameter( PyObject* parent, const char* key, unsigned int& out )
{
PyObject* pobj = getElement( parent, key ); // borrowed
if ( !pobj )
return;
#ifdef PYTHON2
if ( !PyInt_Check( pobj ) )
throwPythonError( Form( "Object \"%s\" has invalid type %s", key, pobj->ob_type->tp_name ) );
out = PyInt_AsUnsignedLongMask( pobj );
#else
if ( !PyLong_Check( pobj ) )
throwPythonError( Form( "Object \"%s\" has invalid type %s", key, pobj->ob_type->tp_name ) );
out = PyLong_AsUnsignedLong( pobj );
#endif
}
void
PythonHandler::fillParameter( PyObject* parent, const char* key, double& out )
{
PyObject* pobj = getElement( parent, key ); // borrowed
if ( !pobj )
return;
if ( !PyFloat_Check( pobj ) )
throwPythonError( Form( "Object \"%s\" has invalid type %s", key, pobj->ob_type->tp_name ) );
out = PyFloat_AsDouble( pobj );
}
void
PythonHandler::fillParameter( PyObject* parent, const char* key, std::string& out )
{
PyObject* pobj = getElement( parent, key ); // borrowed
if ( !pobj )
return;
if ( !
#ifdef PYTHON2
PyString_Check( pobj )
#else
PyUnicode_Check( pobj )
#endif
)
throwPythonError( Form( "Object \"%s\" has invalid type %s", key, pobj->ob_type->tp_name ) );
out = decode( pobj );
}
void
PythonHandler::fillParameter( PyObject* parent, const char* key, std::vector<double>& out )
{
out.clear();
PyObject* pobj = getElement( parent, key ); // borrowed
if ( !pobj )
return;
if ( !PyTuple_Check( pobj ) )
throwPythonError( Form( "Object \"%s\" has invalid type %s", key, pobj->ob_type->tp_name ) );
for ( Py_ssize_t i = 0; i < PyTuple_Size( pobj ); ++i ) {
PyObject* pit = PyTuple_GetItem( pobj, i ); // borrowed
if ( !PyFloat_Check( pit ) )
continue;
out.emplace_back( PyFloat_AsDouble( pit ) );
}
}
void
PythonHandler::fillParameter( PyObject* parent, const char* key, std::vector<std::string>& out )
{
out.clear();
PyObject* pobj = getElement( parent, key ); // borrowed
if ( !pobj )
return;
if ( !PyTuple_Check( pobj ) )
throwPythonError( Form( "Object \"%s\" has invalid type %s", key, pobj->ob_type->tp_name ) );
for ( Py_ssize_t i = 0; i < PyTuple_Size( pobj ); ++i ) {
PyObject* pit = PyTuple_GetItem( pobj, i ); // borrowed
out.emplace_back( decode( pit ) );
}
}
void
PythonHandler::fillParameter( PyObject* parent, const char* key, std::vector<int>& out )
{
out.clear();
PyObject* pobj = getElement( parent, key ); // borrowed
if ( !pobj )
return;
if ( !PyTuple_Check( pobj ) )
throwPythonError( Form( "Object \"%s\" has invalid type", key ) );
for ( Py_ssize_t i = 0; i < PyTuple_Size( pobj ); ++i ) {
PyObject* pit = PyTuple_GetItem( pobj, i );
if ( !isInteger( pit ) )
throwPythonError( Form( "Object %d has invalid type", i ) );
out.emplace_back( asInteger( pit ) );
}
}
bool
PythonHandler::isInteger( PyObject* obj )
{
#ifdef PYTHON2
return PyInt_Check( obj );
#else
return PyLong_Check( obj );
#endif
}
int
PythonHandler::asInteger( PyObject* obj )
{
#ifdef PYTHON2
return PyInt_AsLong( obj );
#else
return PyLong_AsLong( obj );
#endif
}
}
}
#endif
diff --git a/CepGen/IO/FortranInterface.cpp b/CepGen/IO/FortranInterface.cpp
index f099364..7b27b6c 100644
--- a/CepGen/IO/FortranInterface.cpp
+++ b/CepGen/IO/FortranInterface.cpp
@@ -1,73 +1,67 @@
#include "CepGen/StructureFunctions/StructureFunctionsBuilder.h"
#include "CepGen/StructureFunctions/MSTWGrid.h"
#include "CepGen/Processes/GenericKTProcess.h"
#include "CepGen/Core/Exception.h"
#ifdef __cplusplus
extern "C" {
#endif
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;
- if ( sf_mode == SF::Type::MSTWgrid ) {
- StructureFunctions sf = MSTW::Grid::get( "External/F2_Luxlike_fit/mstw_f2_scan_nnlo.dat" ).eval( q2, xbj );
- f2 = sf.F2;
- fl = sf.FL;
- return;
- }
StructureFunctions* sf = StructureFunctionsBuilder::get( sf_mode );
StructureFunctions val = ( *sf )( q2, xbj );
f2 = val.F2;
fl = val.FL;
delete sf;
}
/*bool lhapdf_init = false;
double
cepgen_coll_flux_( int& fmode, double& q2, double& x )
{
if ( !lhapdf_init ) {
LHAPDF::initPDFSet( set, LHAPDF::LHGRID, 0 );
}
}*/
double
cepgen_kt_flux_( int& fmode, double& kt2, double& x, int& sfmode, double& mx )
{
using namespace CepGen;
using namespace CepGen::Process;
return GenericKTProcess::flux( (GenericKTProcess::Flux)fmode, kt2, x,
*StructureFunctionsBuilder::get( (SF::Type)sfmode ), mx );
}
double
cepgen_kt_flux_hi_( int& fmode, double& kt2, double& x, int& a, int& z )
{
using namespace CepGen;
using namespace CepGen::Process;
return GenericKTProcess::flux( (GenericKTProcess::Flux)fmode, kt2, x,
Kinematics::HeavyIon{ ( unsigned short )a, ( unsigned short )z } );
}
double
cepgen_particle_mass_( int& pdg_id )
{
return CepGen::ParticleProperties::mass( (CepGen::PDG)pdg_id );
}
double
cepgen_particle_charge_( int& pdg_id )
{
return CepGen::ParticleProperties::charge( pdg_id );
}
#ifdef __cplusplus
}
#endif
diff --git a/CepGen/StructureFunctions/MSTWGrid.cpp b/CepGen/StructureFunctions/MSTWGrid.cpp
index 205c48c..e73648f 100644
--- a/CepGen/StructureFunctions/MSTWGrid.cpp
+++ b/CepGen/StructureFunctions/MSTWGrid.cpp
@@ -1,115 +1,118 @@
#include "CepGen/StructureFunctions/MSTWGrid.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 Grid::good_magic = 0x5754534d; // MSTW in ASCII
Grid&
Grid::get( const char* filename )
{
- static Grid instance( filename );
+ Parameterisation p;
+ p.grid_path = filename;
+ static Grid instance( p );
return instance;
}
- Grid::Grid( const char* filename ) :
- CepGen::GridHandler<2>( CepGen::GridType::kLogarithmic )
+ 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( filename, std::ios::binary | std::ios::in );
+ 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 \"" << filename << "%s\"!";
+ 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() ) << "].";
}
- CepGen::StructureFunctions
- Grid::eval( double q2, double xbj ) const
+ Grid&
+ Grid::operator()( double q2, double xbj )
{
const sfval_t val = CepGen::GridHandler<2>::eval( q2, xbj );
- CepGen::StructureFunctions sf;
- sf.F2 = val.value[0];
- sf.FL = val.value[1];
- return sf;
+ 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
index 7f83c33..1180500 100644
--- a/CepGen/StructureFunctions/MSTWGrid.h
+++ b/CepGen/StructureFunctions/MSTWGrid.h
@@ -1,52 +1,58 @@
-#ifndef CepGen_IO_MSTWGridHandler_h
-#define CepGen_IO_MSTWGridHandler_h
+#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::GridHandler<2>
+ 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* filename = "External/F2_Luxlike_fit/mstw_f2_scan_nnlo.dat" );
+ 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$
- CepGen::StructureFunctions eval( double q2, double xbj ) const;
+ 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 char* );
+ explicit Grid( const Parameterisation& );
static const unsigned int good_magic;
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 54a3180..12edcac 100644
--- a/CepGen/StructureFunctions/StructureFunctionsBuilder.cpp
+++ b/CepGen/StructureFunctions/StructureFunctionsBuilder.cpp
@@ -1,38 +1,38 @@
#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/StructureFunctions/MSTWGrid.h"
namespace CepGen
{
StructureFunctions*
StructureFunctionsBuilder::get( const SF::Type& sf_type )
{
switch ( sf_type ) {
case SF::Type::Electron:
- case SF::Type::MSTWgrid: //FIXME
- default:
- case SF::Type::ElasticProton: return new StructureFunctions();
+ case SF::Type::ElasticProton:
+ default: return new StructureFunctions();
case SF::Type::SzczurekUleshchenko: return new SF::SzczurekUleshchenko();
case SF::Type::SuriYennie: return new SF::SuriYennie();
case SF::Type::FioreBrasse: return new SF::FioreBrasse();
case SF::Type::ChristyBosted: return new SF::ChristyBosted();
case SF::Type::CLAS: return new SF::CLAS();
case SF::Type::BlockDurandHa: return new SF::BlockDurandHa();
case SF::Type::ALLM91: return new SF::ALLM( SF::ALLM::Parameterisation::allm91() );
case SF::Type::ALLM97: return new SF::ALLM( SF::ALLM::Parameterisation::allm97() );
case SF::Type::GD07p: return new SF::ALLM( SF::ALLM::Parameterisation::gd07p() );
case SF::Type::GD11p: return new SF::ALLM( SF::ALLM::Parameterisation::gd11p() );
case SF::Type::Schaefer: return new SF::Schaefer();
case SF::Type::GenericLHAPDF: return new SF::GenericLHAPDF();
+ case SF::Type::MSTWgrid: return &MSTW::Grid::get();
}
}
}