diff --git a/CepGen/Cards/PythonHandler.cpp b/CepGen/Cards/PythonHandler.cpp
index 7c6a963..1182d1c 100644
--- a/CepGen/Cards/PythonHandler.cpp
+++ b/CepGen/Cards/PythonHandler.cpp
@@ -1,595 +1,596 @@
#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/StructureFunctions/StructureFunctionsBuilder.h"
#include "CepGen/Hadronisers/Pythia8Hadroniser.h"
#include <algorithm>
#include <frameobject.h>
#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 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 );
}
//--- 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.inp = { pz0, pz1 };
}
double sqrt_s = -1.;
fillParameter( kin, "cmEnergy", sqrt_s );
if ( sqrt_s != -1. )
params_.kinematics.setSqrtS( sqrt_s );
int str_fun = 0;
fillParameter( kin, "structureFunctions", str_fun );
- //params_.kinematics.structure_functions = ;
+ params_.kinematics.structure_functions = StructureFunctionsBuilder::get( (SF::Type)str_fun );
}
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<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 ) );
}
}
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 89028c0..38151c6 100644
--- a/CepGen/IO/FortranInterface.cpp
+++ b/CepGen/IO/FortranInterface.cpp
@@ -1,50 +1,22 @@
-#include "CepGen/StructureFunctions/SzczurekUleshchenko.h"
-#include "CepGen/StructureFunctions/SuriYennie.h"
-#include "CepGen/StructureFunctions/FioreBrasse.h"
-#include "CepGen/StructureFunctions/ChristyBosted.h"
-#include "CepGen/StructureFunctions/CLAS.h"
-#include "CepGen/StructureFunctions/BlockDurandHa.h"
-#include "CepGen/StructureFunctions/ALLM.h"
-#include "CepGen/StructureFunctions/Schaefer.h"
-#include "CepGen/StructureFunctions/GenericLHAPDF.h"
+#include "CepGen/StructureFunctions/StructureFunctionsBuilder.h"
#include "CepGen/IO/MSTWGridHandler.h"
extern "C"
{
void
cepgen_structure_functions_( int& sfmode, double& q2, double& xbj, double& f2, double& fl )
{
using namespace CepGen;
- StructureFunctions sf;
- switch ( (SF::Type)sfmode ) {
- case SF::Type::Electron:
- case SF::Type::ElasticProton:
- case SF::Type::Invalid: {
- f2 = fl = 0.;
- return;
- }
- case SF::Type::SzczurekUleshchenko: sf = SF::SzczurekUleshchenko(); break;
- case SF::Type::SuriYennie: sf = SF::SuriYennie(); break;
- case SF::Type::FioreBrasse: sf = SF::FioreBrasse(); break;
- case SF::Type::ChristyBosted: sf = SF::ChristyBosted(); break;
- case SF::Type::CLAS: sf = SF::CLAS(); break;
- case SF::Type::BlockDurandHa: sf = SF::BlockDurandHa(); break;
- case SF::Type::ALLM91: sf = SF::ALLM( SF::ALLM::Parameterisation::allm91() ); break;
- case SF::Type::ALLM97: sf = SF::ALLM( SF::ALLM::Parameterisation::allm97() ); break;
- case SF::Type::GD07p: sf = SF::ALLM( SF::ALLM::Parameterisation::gd07p() ); break;
- case SF::Type::GD11p: sf = SF::ALLM( SF::ALLM::Parameterisation::gd11p() ); break;
- case SF::Type::Schaefer: sf = SF::Schaefer(); break;
- case SF::Type::GenericLHAPDF: sf = SF::GenericLHAPDF(); break;
- case SF::Type::MSTWgrid: {
- sf = MSTW::GridHandler::get().eval( q2, xbj );
- f2 = sf.F2;
- fl = sf.FL;
- return;
- }
+ SF::Type sf_mode = (SF::Type)sfmode;
+ if ( sf_mode == SF::Type::MSTWgrid ) {
+ StructureFunctions sf = MSTW::GridHandler::get().eval( q2, xbj );
+ f2 = sf.F2;
+ fl = sf.FL;
+ return;
}
- sf = sf( q2, xbj );
+ StructureFunctions sf = StructureFunctionsBuilder::get( sf_mode )( q2, xbj );
f2 = sf.F2;
fl = sf.FL;
}
}
diff --git a/CepGen/Physics/Kinematics.cpp b/CepGen/Physics/Kinematics.cpp
index 896b16a..b7d61f0 100644
--- a/CepGen/Physics/Kinematics.cpp
+++ b/CepGen/Physics/Kinematics.cpp
@@ -1,70 +1,75 @@
#include "CepGen/Physics/Kinematics.h"
#include "CepGen/Physics/PDG.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/utils.h"
#include "CepGen/StructureFunctions/SuriYennie.h"
#include "CepGen/Event/Particle.h"
namespace CepGen
{
Kinematics::Kinematics() :
inp( { 6500., 6500. } ), inpdg( { PDG::Proton, PDG::Proton } ),
mode( Mode::ElasticElastic ), structure_functions( SF::SuriYennie() )
{}
+ Kinematics::Kinematics( const Kinematics& kin ) :
+ inp( kin.inp ), inpdg( kin.inpdg ), mode( kin.mode ),
+ structure_functions( kin.structure_functions )
+ {}
+
Kinematics::~Kinematics()
{}
void
Kinematics::setSqrtS( double sqrts )
{
const double pin = 0.5 * sqrts;
inp = { pin, pin };
}
double
Kinematics::sqrtS() const
{
return inp.first + inp.second;
}
//-----------------------------------------------------------------------------------------------
// User-friendly displayers
//-----------------------------------------------------------------------------------------------
std::ostream&
operator<<( std::ostream& os, const Kinematics::Mode& pm )
{
switch ( pm ) {
case Kinematics::Mode::ElectronElectron:
return os << "electron/electron";
case Kinematics::Mode::ElectronProton:
return os << "electron/proton";
case Kinematics::Mode::ProtonElectron:
return os << "proton/electron";
case Kinematics::Mode::ElasticElastic:
return os << "elastic/elastic";
case Kinematics::Mode::InelasticElastic:
return os << "inelastic/elastic";
case Kinematics::Mode::ElasticInelastic:
return os << "elastic/inelastic";
case Kinematics::Mode::InelasticInelastic:
return os << "inelastic/inelastic";
}
return os;
}
//------------------------------------------------------------------------------------------------
// List of kinematics limits
//------------------------------------------------------------------------------------------------
Kinematics::CutsList::CutsList()
{
initial.q2 = { 0., 1.e5 };
central.pt_single.min() = 0.;
remnants.mass_single = { 1.07, 320. };
}
}
diff --git a/CepGen/Physics/Kinematics.h b/CepGen/Physics/Kinematics.h
index 4f40228..6edd99d 100644
--- a/CepGen/Physics/Kinematics.h
+++ b/CepGen/Physics/Kinematics.h
@@ -1,96 +1,97 @@
#ifndef CepGen_Physics_Kinematics_h
#define CepGen_Physics_Kinematics_h
#include <iomanip>
#include <algorithm>
#include "CepGen/Core/Logger.h"
#include "CepGen/StructureFunctions/StructureFunctions.h"
#include "CepGen/Physics/ParticleProperties.h"
#include "CepGen/Physics/Cuts.h"
#include "CepGen/Physics/Limits.h"
#include <vector>
#include <unordered_map>
namespace CepGen
{
enum class PDG;
/// List of kinematic constraints to apply on the process phase space.
class Kinematics
{
public:
Kinematics();
+ Kinematics( const Kinematics& );
~Kinematics();
/// Type of kinematics to consider for the process
enum class Mode {
ElectronProton = 0, ///< electron-proton elastic case
ElasticElastic = 1, ///< proton-proton elastic case
ElasticInelastic = 2, ///< proton-proton single-dissociative (or inelastic-elastic) case
InelasticElastic = 3, ///< proton-proton single-dissociative (or elastic-inelastic) case
InelasticInelastic = 4, ///< proton-proton double-dissociative case
ProtonElectron,
ElectronElectron
};
/// Human-readable format of a process mode (elastic/dissociative parts)
friend std::ostream& operator<<( std::ostream&, const Mode& );
/// Incoming particles' momentum (in \f$\text{GeV}/c\f$)
std::pair<double,double> inp;
/// Set the incoming particles' momenta (if the collision is symmetric)
void setSqrtS( double sqrts );
/// Process centre of mass energy
double sqrtS() const;
/// Beam/primary particle's PDG identifier
std::pair<PDG,PDG> inpdg;
/// PDG id of the outgoing central particles
std::vector<PDG> central_system;
/// Minimum list of central particles required
std::vector<PDG> minimum_final_state;
/// Type of kinematics to consider for the phase space
Mode mode;
/// Type of structure functions to consider
StructureFunctions structure_functions;
struct CutsList
{
CutsList();
template<class T,bool>
struct hasher
{
inline size_t operator()( const T& t ) const {
return std::hash<T>()( t );
}
};
template<class T>
struct hasher<T, true>
{
inline size_t operator() ( const T& t ) {
typedef typename std::underlying_type<T>::type enumType;
return std::hash<enumType>()( static_cast<enumType>( t ) );
}
};
template<class T>
struct EnumHash
{
inline size_t operator()( const T& t ) const {
return hasher<T,std::is_enum<T>::value>()( t );
}
};
/// Cuts on the initial particles kinematics
Cuts initial;
/// Cuts on the central system produced
Cuts central;
std::unordered_map<PDG,Cuts,EnumHash<PDG> > central_particles;
/// Cuts on the beam remnants system
Cuts remnants;
};
CutsList cuts;
};
}
#endif
diff --git a/CepGen/StructureFunctions/ALLM.cpp b/CepGen/StructureFunctions/ALLM.cpp
index 227218e..6866a51 100644
--- a/CepGen/StructureFunctions/ALLM.cpp
+++ b/CepGen/StructureFunctions/ALLM.cpp
@@ -1,167 +1,168 @@
#include "ALLM.h"
#include "CepGen/Physics/ParticleProperties.h"
#include <cmath>
namespace CepGen
{
namespace SF
{
ALLM::Parameterisation
ALLM::Parameterisation::allm91()
{
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.67639, 0.49027, 2.66275 },
{ 0.60408, 0.17353, 1.61812 },
{ 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;
p.type = Type::ALLM91;
return p;
}
ALLM::Parameterisation
ALLM::Parameterisation::allm97()
{
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.80107, 0.97307, 3.4924 },
{ 0.58400, 0.37888, 2.6063 },
{ 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;
p.type = Type::ALLM97;
return p;
}
ALLM::Parameterisation
ALLM::Parameterisation::hht_allm()
{
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( {
{ -1.04, 2.97, 0.163 },
{ 0.706, 0.185, -16.4 },
{ -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;
p.type = Type::Invalid;
return p;
}
ALLM::Parameterisation
ALLM::Parameterisation::hht_allm_ft()
{
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.636, 3.37, -0.660 },
{ 0.882, 0.082, -8.5 },
{ 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;
p.type = Type::Invalid;
return p;
}
ALLM::Parameterisation
ALLM::Parameterisation::gd07p()
{
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( {
{ 0.838, 2.36, 1.77 },
{ 0.374, 0.998, 0.775 },
{ 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;
p.type = Type::GD07p;
return p;
}
ALLM::Parameterisation
ALLM::Parameterisation::gd11p()
{
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.3633, 2.256, 2.209 },
{ 0.3425, 1.0603, 0.5164 },
{ -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;
p.type = Type::GD11p;
return p;
}
ALLM::ALLM( const ALLM::Parameterisation& param ) :
StructureFunctions( param.type ), params_( param )
{}
ALLM
ALLM::operator()( double q2, double xbj, const SigmaRatio& rcomp ) const
{
const double W2_eff = q2*( 1.-xbj )/xbj;
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+params_.q02 )/ params_.lambda2 ), xlog2 = log( params_.q02/params_.lambda2 );
const double t = log( xlog1/xlog2 );
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 = 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 );
ALLM allm;
allm.F2 = q2/( q2+params_.m02 ) * ( F2_Pom + F2_Reg );
allm.computeFL( q2, xbj, rcomp );
return allm;
}
}
}
+
diff --git a/CepGen/StructureFunctions/StructureFunctions.cpp b/CepGen/StructureFunctions/StructureFunctions.cpp
index edfd54f..05ef043 100644
--- a/CepGen/StructureFunctions/StructureFunctions.cpp
+++ b/CepGen/StructureFunctions/StructureFunctions.cpp
@@ -1,74 +1,75 @@
#include "CepGen/StructureFunctions/StructureFunctions.h"
#include "CepGen/Physics/PDG.h"
#include "CepGen/Physics/ParticleProperties.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/utils.h"
#include <iostream>
namespace CepGen
{
const double StructureFunctions::mp_ = ParticleProperties::mass( PDG::Proton );
const double StructureFunctions::mp2_ = StructureFunctions::mp_*StructureFunctions::mp_;
double
StructureFunctions::F1( double q2, double xbj ) const
{
if ( xbj == 0. || q2 == 0. ) {
CG_ERROR( "StructureFunctions:F1" )
<< "Invalid range for Q² = " << q2 << " or xBj = " << xbj << ".";
return 0.;
}
const double F1 = 0.5*( ( 1+4.*xbj*xbj*mp2_/q2 )*F2 - FL )/xbj;
CG_DEBUG_LOOP( "StructureFunctions:F1" )
<< "F1 for Q² = " << q2 << ", xBj = " << xbj << ": " << F1 << "\n\t"
<< "(F2 = " << F2 << ", FL = " << FL << ").";
return F1;
}
void
StructureFunctions::computeFL( double q2, double xbj, const SF::SigmaRatio& ratio )
{
double r_error = 0.;
computeFL( q2, xbj, ratio( q2, xbj, r_error ) );
}
void
StructureFunctions::computeFL( double q2, double xbj, double r )
{
const double tau = 4.*xbj*xbj*mp2_/q2;
FL = F2 * ( 1.+tau ) * ( r/( 1.+r ) );
}
/// Human-readable format of a structure function object
std::ostream&
operator<<( std::ostream& os, const StructureFunctions& sf )
{
return os << "F2 = " << sf.F2 << ", FL = " << sf.FL;
}
/// Human-readable format of a structure function type
std::ostream&
operator<<( std::ostream& os, const SF::Type& sf )
{
switch ( sf ) {
case SF::Type::Invalid: return os << "[INVALID]";
case SF::Type::Electron: return os << "electron";
case SF::Type::ElasticProton: return os << "elastic proton";
case SF::Type::SuriYennie: return os << "Suri-Yennie";
case SF::Type::SzczurekUleshchenko: return os << "Szczurek-Uleshchenko";
case SF::Type::FioreBrasse: return os << "Fiore-Brasse";
case SF::Type::ChristyBosted: return os << "Christy-Bosted";
case SF::Type::CLAS: return os << "CLAS";
case SF::Type::BlockDurandHa: return os << "BDH";
case SF::Type::ALLM91: return os << "ALLM;91";
case SF::Type::ALLM97: return os << "ALLM;97";
case SF::Type::GD07p: return os << "ALLM;GD07p";
case SF::Type::GD11p: return os << "ALLM;GD11p";
case SF::Type::Schaefer: return os << "Schaefer";
case SF::Type::MSTWgrid: return os << "MSTW (grid)";
+ case SF::Type::GenericLHAPDF: return os << "LHAPDF (generic)";
}
return os;
}
}
diff --git a/CepGen/StructureFunctions/StructureFunctions.h b/CepGen/StructureFunctions/StructureFunctions.h
index b314deb..3b62fb3 100644
--- a/CepGen/StructureFunctions/StructureFunctions.h
+++ b/CepGen/StructureFunctions/StructureFunctions.h
@@ -1,58 +1,65 @@
#ifndef CepGen_StructureFunctions_StructureFunctions_h
#define CepGen_StructureFunctions_StructureFunctions_h
+#include "CepGen/StructureFunctions/SigmaRatio.h"
+
#include <iostream>
-#include "SigmaRatio.h"
+#include <map>
namespace CepGen
{
namespace SF
{
/// Proton structure function to be used in the outgoing state description
/// \note Values correspond to the LPAIR legacy steering card values
enum struct Type {
Invalid = 0,
Electron = 1,
ElasticProton = 2,
SuriYennie = 11,
SzczurekUleshchenko = 12,
BlockDurandHa = 13,
FioreBrasse = 101,
ChristyBosted = 102,
CLAS = 103,
ALLM91 = 201,
ALLM97 = 202,
GD07p = 203,
GD11p = 204,
MSTWgrid = 205,
Schaefer = 301,
GenericLHAPDF = 401,
};
}
std::ostream& operator<<( std::ostream&, const SF::Type& );
+ class StructureFunctionsFactory;
class StructureFunctions
{
public:
+ StructureFunctions( const StructureFunctions& sf ) :
+ type( sf.type ), F2( sf.F2 ), FL( sf.FL ), old_vals_( sf.old_vals_ ) {}
StructureFunctions( const SF::Type& type = SF::Type::Invalid, double f2 = 0., double fl = 0. ) :
type( type ), F2( f2 ), FL( fl ), old_vals_({ 0., 0. }) {}
+ static StructureFunctions builder( const SF::Type& );
+
virtual StructureFunctions& operator()( double q2, double xbj ) { return *this; }
void computeFL( double q2, double xbj, const SF::SigmaRatio& ratio = SF::E143Ratio() );
void computeFL( double q2, double xbj, double r );
double F1( double q2, double xbj ) const;
SF::Type type;
double F2, FL;
protected:
static const double mp_, mp2_;
std::pair<double,double> old_vals_;
private:
std::string name_;
};
std::ostream& operator<<( std::ostream&, const StructureFunctions& );
}
#endif
diff --git a/CepGen/StructureFunctions/StructureFunctionsBuilder.cpp b/CepGen/StructureFunctions/StructureFunctionsBuilder.cpp
new file mode 100644
index 0000000..f876f99
--- /dev/null
+++ b/CepGen/StructureFunctions/StructureFunctionsBuilder.cpp
@@ -0,0 +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/IO/MSTWGridHandler.h"
+
+namespace CepGen
+{
+ StructureFunctions
+ StructureFunctionsBuilder::get( const SF::Type& sf_type )
+ {
+ std::cout << sf_type << std::endl;
+ switch ( sf_type ) {
+ case SF::Type::Electron:
+ case SF::Type::MSTWgrid: //FIXME
+ default:
+ case SF::Type::ElasticProton: return StructureFunctions();
+ case SF::Type::SzczurekUleshchenko: return SF::SzczurekUleshchenko();
+ case SF::Type::SuriYennie: return SF::SuriYennie();
+ case SF::Type::FioreBrasse: return SF::FioreBrasse();
+ case SF::Type::ChristyBosted: return SF::ChristyBosted();
+ case SF::Type::CLAS: return SF::CLAS();
+ case SF::Type::BlockDurandHa: return SF::BlockDurandHa();
+ case SF::Type::ALLM91: return SF::ALLM( SF::ALLM::Parameterisation::allm91() );
+ case SF::Type::ALLM97: return SF::ALLM( SF::ALLM::Parameterisation::allm97() );
+ case SF::Type::GD07p: return SF::ALLM( SF::ALLM::Parameterisation::gd07p() );
+ case SF::Type::GD11p: return SF::ALLM( SF::ALLM::Parameterisation::gd11p() );
+ case SF::Type::Schaefer: return SF::Schaefer();
+ }
+ }
+}
diff --git a/CepGen/StructureFunctions/StructureFunctionsBuilder.h b/CepGen/StructureFunctions/StructureFunctionsBuilder.h
new file mode 100644
index 0000000..8a58944
--- /dev/null
+++ b/CepGen/StructureFunctions/StructureFunctionsBuilder.h
@@ -0,0 +1,22 @@
+#ifndef CepGen_StructureFunctions_StructureFunctionsBuilder_h
+#define CepGen_StructureFunctions_StructureFunctionsBuilder_h
+
+#include "CepGen/StructureFunctions/StructureFunctions.h"
+
+namespace CepGen
+{
+ /// Helper class to generate any supported set of structure functions
+ class StructureFunctionsBuilder
+ {
+ public:
+ StructureFunctionsBuilder() {}
+ ~StructureFunctionsBuilder() {}
+
+ /// Build structure functions from the modelling type
+ static StructureFunctions get( const SF::Type& );
+ /// Build structure functions from the modelling name
+ static StructureFunctions get( const char* );
+ };
+}
+
+#endif
diff --git a/CepGen/StructureFunctions/SuriYennie.cpp b/CepGen/StructureFunctions/SuriYennie.cpp
index 850ff99..dc60836 100644
--- a/CepGen/StructureFunctions/SuriYennie.cpp
+++ b/CepGen/StructureFunctions/SuriYennie.cpp
@@ -1,65 +1,67 @@
#include "SuriYennie.h"
#include "CepGen/Physics/ParticleProperties.h"
+#include "CepGen/Core/Exception.h"
#include <math.h>
namespace CepGen
{
namespace SF
{
SuriYennie::Parameterisation
SuriYennie::Parameterisation::standard()
{
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;
}
SuriYennie::Parameterisation
SuriYennie::Parameterisation::alternative()
{
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;
}
SuriYennie::SuriYennie( const SuriYennie::Parameterisation& param ) :
F1( 0. ), FE( 0. ), FM( 0. ), params_( param )
{}
SuriYennie&
SuriYennie::operator()( double q2, double xbj )
{
+ CG_INFO( "SY" );
std::pair<double,double> nv = { q2, xbj };
if ( nv == old_vals_ )
return *this;
old_vals_ = nv;
const double inv_q2 = 1./q2;
const double dm2 = q2*( 1.-xbj )/xbj;
const double mx2 = mp2_ + dm2; // [GeV^2]
const double en = dm2+q2, en2 = en*en; // [GeV^2]
const double nu = 0.5 * en / mp_, Xpr = q2/( q2+mx2 ), tau = 0.25 * q2/mp2_;
const double mq = params_.rho2+q2;
FM = ( params_.C1*dm2*pow( params_.rho2/mq, 2 )
+ ( params_.C2*mp2_*pow( 1.-Xpr, 4 ) ) / ( 1.+Xpr*( Xpr*params_.Cp-2.*params_.Bp ) ) ) * inv_q2;
FE = ( tau*FM + params_.D1*dm2*q2*params_.rho2/mp2_*pow( dm2/mq, 2 )/en2 ) / ( 1.+0.25*en2/mp2_*inv_q2 );
//const double w2 = 2.*mp*FE;
F1 = 0.5*FM*q2/mp_;
F2 = 2.*nu*FE;
return *this;
}
}
}