diff --git a/Cards/Config/ktProcess_cfi.py b/Cards/Config/ktProcess_cfi.py
index 4fb9d0e..fd40eb0 100644
--- a/Cards/Config/ktProcess_cfi.py
+++ b/Cards/Config/ktProcess_cfi.py
@@ -1,19 +1,22 @@
import Config.Core as cepgen
from math import pi
class ProtonFlux:
PhotonElastic = 0
PhotonInelastic = 1
PhotonInelasticBudnev = 11
+ GluonKMR = 20
+class HeavyIonFlux:
+ PhotonElastic = 100
process = cepgen.Module('ktProcess',
outKinematics = cepgen.Parameters(
qt = (0., 50.),
phiqt = (0., 2.*pi),
#--- cuts on individual particles defining the central system
rapidity = (-6., 6.),
#--- cuts on the pt(outgoing system) (hyper-)plane
ptdiff = (0., 500.),
phiptdiff = (0., 2.*pi),
),
)
diff --git a/Cards/patoccbar_cfg.py b/Cards/patoccbar_cfg.py
new file mode 100644
index 0000000..c3d6834
--- /dev/null
+++ b/Cards/patoccbar_cfg.py
@@ -0,0 +1,40 @@
+import Config.Core as cepgen
+import Config.ktProcess_cfi as ktfactor
+from Config.integrators_cff import vegas as integrator
+from Config.pdg_cff import PDG
+
+from Config.logger_cfi import logger
+logger.enabledModules += ('GenericKTProcess.registerVariable',)
+
+process = ktfactor.process.clone('patoll',
+ processParameters = cepgen.Parameters(
+ pair = PDG.charm,
+ ),
+ inKinematics = cepgen.Parameters(
+ pz = (6500., 2562.2),
+ structureFunctions = cepgen.StructureFunctions.SuriYennie,
+ #structureFunctions = cepgen.StructureFunctions.FioreBrasse,
+ ktFluxes = (ktfactor.ProtonFlux.GluonKMR, ktfactor.HeavyIonFlux.PhotonElastic),
+ #ktFluxes = (ktfactor.ProtonFlux.PhotonElastic, ktfactor.HeavyIonFlux.PhotonElastic),
+ heavyIonB = (208, 82),
+ kmrGridPath = 'gluon_mmht2014nlo_Watt.dat',
+ ),
+ outKinematics = ktfactor.process.outKinematics.clone(
+ pt = (0.,),
+ energy = (0.,),
+ rapidity = (-7., 9.),
+ #qt = (0.,1000.),
+ #eta = (-2.5, 2.5),
+ mx = (1.07, 1000.),
+ #--- extra cuts on the p1t(l) and p2t(l) plane
+ #ptdiff = (0., 2.5),
+ #--- distance in rapidity between l^+ and l^-
+ #dely = (4., 5.),
+ ),
+)
+
+#--- events generation
+from Config.generator_cff import generator
+generator.numEvents = 10000
+generator.numThreads = 1
+
diff --git a/Cards/patoll_cfg.py b/Cards/patoll_cfg.py
new file mode 100644
index 0000000..25b9830
--- /dev/null
+++ b/Cards/patoll_cfg.py
@@ -0,0 +1,35 @@
+import Config.Core as cepgen
+import Config.ktProcess_cfi as ktfactor
+from Config.integrators_cff import vegas as integrator
+from Config.pdg_cff import PDG
+
+process = ktfactor.process.clone('patoll',
+ processParameters = cepgen.Parameters(
+ pair = PDG.muon,
+ ),
+ inKinematics = cepgen.Parameters(
+ pz = (6500., 2562.2),
+ #structureFunctions = cepgen.StructureFunctions.SuriYennie,
+ #structureFunctions = cepgen.StructureFunctions.FioreBrasse,
+ #structureFunctions = cepgen.StructureFunctions.ALLM91,
+ structureFunctions = cepgen.StructureFunctions.LUXlike,
+ ktFluxes = (ktfactor.ProtonFlux.PhotonInelasticBudnev, ktfactor.HeavyIonFlux.PhotonElastic),
+ heavyIonB = (208, 82),
+ ),
+ outKinematics = ktfactor.process.outKinematics.clone(
+ pt = (4.,),
+ energy = (0.,),
+ rapidity = (-6., 7.),
+ #eta = (-2.5, 2.5),
+ mx = (1.07, 1000.),
+ #--- extra cuts on the p1t(l) and p2t(l) plane
+ #ptdiff = (0., 2.5),
+ #--- distance in rapidity between l^+ and l^-
+ #dely = (4., 5.),
+ ),
+)
+
+#--- events generation
+from Config.generator_cff import generator
+generator.numEvents = 100000
+
diff --git a/CepGen/CMakeLists.txt b/CepGen/CMakeLists.txt
index 55b517c..42483b1 100644
--- a/CepGen/CMakeLists.txt
+++ b/CepGen/CMakeLists.txt
@@ -1,28 +1,34 @@
file(GLOB core_sources Core/*.cpp *.cpp)
file(GLOB phys_sources Physics/*.cpp)
file(GLOB sf_sources StructureFunctions/*.cpp)
include_directories(${PROJECT_SOURCE_DIR})
#----- check the external dependencies for SFs
set(GRV_PATH ${PROJECT_SOURCE_DIR}/External)
file(GLOB grv_sources ${GRV_PATH}/grv_*.f)
if(grv_sources)
message(STATUS "GRV PDFset found in ${grv_sources}!")
add_definitions(-DGRVPDF)
list(APPEND sf_sources ${grv_sources})
else()
message(STATUS "GRV PDFset not found. Will proceed without it")
endif()
+#----- check the external dependencies for physics utilities
+
+if(alphas_sources)
+ list(APPEND phys_sources ${alphas_sources})
+endif()
+
#----- build the object
add_library(CepGenCore SHARED ${core_sources} ${phys_sources} ${sf_sources})
target_link_libraries(CepGenCore ${CEPGEN_EXTERNAL_CORE_REQS})
target_link_libraries(CepGenCore ${CEPGEN_EXTERNAL_STRF_REQS})
#----- installation rules
install(TARGETS CepGenCore DESTINATION lib)
diff --git a/CepGen/Cards/LpairHandler.cpp b/CepGen/Cards/LpairHandler.cpp
index 5c11a03..cfb0b7c 100644
--- a/CepGen/Cards/LpairHandler.cpp
+++ b/CepGen/Cards/LpairHandler.cpp
@@ -1,192 +1,205 @@
#include "CepGen/Cards/LpairHandler.h"
#include "CepGen/Core/ParametersList.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Physics/PDG.h"
#include "CepGen/Processes/GamGamLL.h"
#include "CepGen/Processes/PPtoFF.h"
#include "CepGen/Processes/PPtoWW.h"
#include "CepGen/Hadronisers/Pythia8Hadroniser.h"
#include <fstream>
namespace CepGen
{
namespace Cards
{
const int LpairHandler::kInvalid = 99999;
//----- specialization for LPAIR input cards
LpairHandler::LpairHandler( const char* file ) :
- proc_params_( new ParametersList )
+ proc_params_( new ParametersList ),
+ hi_1_( { 0, 0 } ), hi_2_( { 0, 0 } )
{
std::ifstream f( file, std::fstream::in );
if ( !f.is_open() )
throw CG_FATAL( "LpairHandler" ) << "Failed to parse file \"" << file << "%s\".";
init( ¶ms_ );
std::unordered_map<std::string, std::string> m_params;
std::string key, value;
std::ostringstream os;
while ( f >> key >> value ) {
if ( key[0] == '#' ) continue; // FIXME need to ensure there is no extra space before!
setParameter( key, value );
m_params.insert( { key, value } );
if ( getDescription( key ) != "null" )
os << "\n>> " << key << " = " << std::setw( 15 ) << getParameter( key )
<< " (" << getDescription( key ) << ")";
}
f.close();
if ( proc_name_ == "lpair" )
params_.setProcess( new Process::GamGamLL( *proc_params_ ) );
else if ( proc_name_ == "pptoll" || proc_name_ == "pptoff" )
params_.setProcess( new Process::PPtoFF( *proc_params_ ) );
else if ( proc_name_ == "pptoww" )
params_.setProcess( new Process::PPtoWW( *proc_params_ ) );
else
throw CG_FATAL( "LpairHandler" ) << "Unrecognised process name: " << proc_name_ << "!";
if ( integr_type_ == "plain" )
params_.integrator.type = Integrator::Type::plain;
else if ( integr_type_ == "Vegas" )
params_.integrator.type = Integrator::Type::Vegas;
else if ( integr_type_ == "MISER" )
params_.integrator.type = Integrator::Type::MISER;
else if ( integr_type_ != "" )
throw CG_FATAL( "LpairHandler" ) << "Unrecognized integrator type: " << integr_type_ << "!";
if ( hadr_name_ == "pythia8" )
params_.setHadroniser( new Hadroniser::Pythia8Hadroniser( params_ ) );
if ( m_params.count( "IEND" ) )
setValue<bool>( "IEND", ( std::stoi( m_params["IEND"] ) > 1 ) );
+ HeavyIon hi1{ hi_1_.first, (Element)hi_1_.second }, hi2{ hi_2_.first, (Element)hi_2_.second };
+ if ( hi1 )
+ params_.kinematics.incoming_beams.first.pdg = hi1;
+ if ( hi2 )
+ params_.kinematics.incoming_beams.second.pdg = hi2;
+
CG_INFO( "LpairHandler" ) << "File '" << file << "' succesfully opened!\n\t"
<< "The following parameters are set:" << os.str();
}
void
LpairHandler::init( Parameters* params )
{
//-------------------------------------------------------------------------------------------
// Process/integration/hadronisation parameters
//-------------------------------------------------------------------------------------------
registerParameter<std::string>( "PROC", "Process name to simulate", &proc_name_ );
registerParameter<std::string>( "ITYP", "Integration algorithm", &integr_type_ );
registerParameter<std::string>( "HADR", "Hadronisation algorithm", &hadr_name_ );
+ registerParameter<std::string>( "KMRG", "KMR grid interpolation path", ¶ms_.kinematics.kmr_grid_path );
//-------------------------------------------------------------------------------------------
// General parameters
//-------------------------------------------------------------------------------------------
registerParameter<bool>( "IEND", "Generation type", ¶ms->generation.enabled );
registerParameter<bool>( "NTRT", "Smoothen the integrand", ¶ms->generation.treat );
registerParameter<int>( "DEBG", "Debugging verbosity", (int*)&Logger::get().level );
registerParameter<int>( "NCVG", "Number of function calls", (int*)¶ms->integrator.ncvg );
registerParameter<int>( "ITVG", "Number of integration iterations", (int*)¶ms->integrator.vegas.iterations );
registerParameter<int>( "SEED", "Random generator seed", (int*)¶ms->integrator.rng_seed );
registerParameter<int>( "NTHR", "Number of threads to use for events generation", (int*)¶ms->generation.num_threads );
registerParameter<int>( "MODE", "Subprocess' mode", (int*)¶ms->kinematics.mode );
registerParameter<int>( "NCSG", "Number of points to probe", (int*)¶ms->generation.num_points );
registerParameter<int>( "NGEN", "Number of events to generate", (int*)¶ms->generation.maxgen );
registerParameter<int>( "NPRN", "Number of events before printout", (int*)¶ms->generation.gen_print_every );
//-------------------------------------------------------------------------------------------
// Process-specific parameters
//-------------------------------------------------------------------------------------------
registerParameter<int>( "METH", "Computation method (kT-factorisation)", &proc_params_->operator[]<int>( "method" ) );
registerParameter<int>( "IPOL", "Polarisation states to consider", &proc_params_->operator[]<int>( "polarisationStates" ) );
//-------------------------------------------------------------------------------------------
// Process kinematics parameters
//-------------------------------------------------------------------------------------------
registerParameter<int>( "PMOD", "Outgoing primary particles' mode", (int*)¶ms->kinematics.structure_functions );
registerParameter<int>( "EMOD", "Outgoing primary particles' mode", (int*)¶ms->kinematics.structure_functions );
registerParameter<int>( "PAIR", "Outgoing particles' PDG id", (int*)&proc_params_->operator[]<int>( "pair" ) );
+
+ registerParameter<int>( "INA1", "Heavy ion atomic weight (1st incoming beam)", (int*)&hi_1_.first );
+ registerParameter<int>( "INZ1", "Heavy ion atomic number (1st incoming beam)", (int*)&hi_1_.second );
+ registerParameter<int>( "INA2", "Heavy ion atomic weight (1st incoming beam)", (int*)&hi_2_.first );
+ registerParameter<int>( "INZ2", "Heavy ion atomic number (1st incoming beam)", (int*)&hi_2_.second );
registerParameter<double>( "INP1", "Momentum (1st primary particle)", ¶ms->kinematics.incoming_beams.first.pz );
registerParameter<double>( "INP2", "Momentum (2nd primary particle)", ¶ms->kinematics.incoming_beams.second.pz );
registerParameter<double>( "INPP", "Momentum (1st primary particle)", ¶ms->kinematics.incoming_beams.first.pz );
registerParameter<double>( "INPE", "Momentum (2nd primary particle)", ¶ms->kinematics.incoming_beams.second.pz );
registerParameter<double>( "PTCT", "Minimal transverse momentum (single central outgoing particle)", ¶ms->kinematics.cuts.central.pt_single.min() );
registerParameter<double>( "MSCT", "Minimal central system mass", ¶ms->kinematics.cuts.central.mass_sum.min() );
registerParameter<double>( "ECUT", "Minimal energy (single central outgoing particle)", ¶ms->kinematics.cuts.central.energy_single.min() );
registerParameter<double>( "ETMN", "Minimal pseudo-rapidity (central outgoing particles)", ¶ms->kinematics.cuts.central.eta_single.min() );
registerParameter<double>( "ETMX", "Maximal pseudo-rapidity (central outgoing particles)", ¶ms->kinematics.cuts.central.eta_single.max() );
registerParameter<double>( "YMIN", "Minimal rapidity (central outgoing particles)", ¶ms->kinematics.cuts.central.rapidity_single.min() );
registerParameter<double>( "YMAX", "Maximal rapidity (central outgoing particles)", ¶ms->kinematics.cuts.central.rapidity_single.max() );
registerParameter<double>( "Q2MN", "Minimal Q² = -q² (exchanged parton)", ¶ms->kinematics.cuts.initial.q2.min() );
registerParameter<double>( "Q2MX", "Maximal Q² = -q² (exchanged parton)", ¶ms->kinematics.cuts.initial.q2.max() );
registerParameter<double>( "MXMN", "Minimal invariant mass of proton remnants", ¶ms->kinematics.cuts.remnants.mass_single.min() );
registerParameter<double>( "MXMX", "Maximal invariant mass of proton remnants", ¶ms->kinematics.cuts.remnants.mass_single.max() );
}
void
LpairHandler::store( const char* file )
{
std::ofstream f( file, std::fstream::out | std::fstream::trunc );
if ( !f.is_open() ) {
CG_ERROR( "LpairHandler" ) << "Failed to open file \"" << file << "%s\" for writing.";
}
for ( const auto& it : p_strings_ )
if ( it.second.value )
f << it.first << " = " << *it.second.value << "\n";
for ( const auto& it : p_ints_ )
if ( it.second.value )
f << it.first << " = " << *it.second.value << "\n";
for ( const auto& it : p_doubles_ )
if ( it.second.value )
f << it.first << " = " << *it.second.value << "\n";
for ( const auto& it : p_bools_ )
if ( it.second.value )
f << it.first << " = " << *it.second.value << "\n";
f.close();
}
void
LpairHandler::setParameter( const std::string& key, const std::string& value )
{
try { setValue<double>( key.c_str(), std::stod( value ) ); } catch ( std::invalid_argument& ) {}
try { setValue<int>( key.c_str(), std::stoi( value ) ); } catch ( std::invalid_argument& ) {}
//setValue<bool>( key.c_str(), std::stoi( value ) );
setValue<std::string>( key.c_str(), value );
}
std::string
LpairHandler::getParameter( std::string key ) const
{
double dd = getValue<double>( key.c_str() );
if ( dd != -999. )
return std::to_string( dd );
int ui = getValue<int>( key.c_str() );
if ( ui != 999 )
return std::to_string( ui );
//if ( out = getValue<bool>( key.c_str() ) );
return getValue<std::string>( key.c_str() );
}
std::string
LpairHandler::getDescription( std::string key ) const
{
if ( p_strings_.count( key ) )
return p_strings_.find( key )->second.description;
if ( p_ints_.count( key ) )
return p_ints_.find( key )->second.description;
if ( p_doubles_.count( key ) )
return p_doubles_.find( key )->second.description;
if ( p_bools_.count( key ) )
return p_bools_.find( key )->second.description;
return "null";
}
}
}
diff --git a/CepGen/Cards/LpairHandler.h b/CepGen/Cards/LpairHandler.h
index 79ba803..78cf708 100644
--- a/CepGen/Cards/LpairHandler.h
+++ b/CepGen/Cards/LpairHandler.h
@@ -1,109 +1,110 @@
#ifndef CepGen_Cards_LpairReader_h
#define CepGen_Cards_LpairReader_h
#include "CepGen/Cards/Handler.h"
#include <unordered_map>
using std::string;
namespace CepGen
{
class ParametersList;
namespace Cards
{
/// LPAIR-like steering cards parser and writer
class LpairHandler : public Handler
{
public:
/// Read a LPAIR steering card
explicit LpairHandler( const char* file );
/// Store a configuration into a LPAIR steering card
void store( const char* file );
private:
template<class T> struct Parameter {
Parameter( const char* key, const char* descr, T* value ) : key( key ), description( descr ), value( value ) {}
std::string key, description;
T* value;
};
/// Register a parameter to be steered to a configuration variable
template<class T> void registerParameter( const char* key, const char* description, T* def ) {}
/// Set a parameter value
template<class T> void setValue( const char* key, const T& value ) {}
/// Retrieve a parameter value
template<class T> T getValue( const char* key ) const {}
void setParameter( const std::string& key, const std::string& value );
std::string getParameter( std::string key ) const;
std::string getDescription( std::string key ) const;
static const int kInvalid;
std::unordered_map<std::string, Parameter<std::string> > p_strings_;
std::unordered_map<std::string, Parameter<double> > p_doubles_;
std::unordered_map<std::string, Parameter<int> > p_ints_;
std::unordered_map<std::string, Parameter<bool> > p_bools_;
void init( Parameters* );
std::shared_ptr<ParametersList> proc_params_;
std::string proc_name_, hadr_name_, integr_type_;
+ std::pair<unsigned short,unsigned short> hi_1_, hi_2_;
};
//----- specialised registerers
template<> inline void LpairHandler::registerParameter<std::string>( const char* key, const char* description, std::string* def ) { p_strings_.insert( std::make_pair( key, Parameter<std::string>( key, description, def ) ) ); }
template<> inline void LpairHandler::registerParameter<double>( const char* key, const char* description, double* def ) { p_doubles_.insert( std::make_pair( key, Parameter<double>( key, description, def ) ) ); }
template<> inline void LpairHandler::registerParameter<int>( const char* key, const char* description, int* def ) { p_ints_.insert( std::make_pair( key, Parameter<int>( key, description, def ) ) ); }
template<> inline void LpairHandler::registerParameter<bool>( const char* key, const char* description, bool* def ) { p_bools_.insert( std::make_pair( key, Parameter<bool>( key, description, def ) ) ); }
//----- specialised setters
template<> inline void LpairHandler::setValue<std::string>( const char* key, const std::string& value ) {
auto it = p_strings_.find( key );
if ( it != p_strings_.end() ) *it->second.value = value;
}
template<> inline void LpairHandler::setValue<double>( const char* key, const double& value ) {
auto it = p_doubles_.find( key );
if ( it != p_doubles_.end() ) *it->second.value = value;
}
template<> inline void LpairHandler::setValue<int>( const char* key, const int& value ) {
auto it = p_ints_.find( key );
if ( it != p_ints_.end() ) *it->second.value = value;
}
template<> inline void LpairHandler::setValue<bool>( const char* key, const bool& value ) {
auto it = p_bools_.find( key );
if ( it != p_bools_.end() ) *it->second.value = value;
}
//----- specialised getters
/// Retrieve a string parameter value
template<> inline std::string LpairHandler::getValue( const char* key ) const {
const auto& it = p_strings_.find( key );
if ( it != p_strings_.end() ) return *it->second.value;
return "null";
}
/// Retrieve a floating point parameter value
template<> inline double LpairHandler::getValue( const char* key ) const {
const auto& it = p_doubles_.find( key );
if ( it != p_doubles_.end() ) return *it->second.value;
return -999.;
}
/// Retrieve an integer parameter value
template<> inline int LpairHandler::getValue( const char* key ) const {
const auto& it = p_ints_.find( key );
if ( it != p_ints_.end() ) return *it->second.value;
return 999;
}
/// Retrieve a boolean parameter value
template<> inline bool LpairHandler::getValue( const char* key ) const {
const auto& it = p_bools_.find( key );
if ( it != p_bools_.end() ) return *it->second.value;
return true;
}
}
}
#endif
diff --git a/CepGen/Cards/PythonHandler.cpp b/CepGen/Cards/PythonHandler.cpp
index 9f64ef7..f6949c2 100644
--- a/CepGen/Cards/PythonHandler.cpp
+++ b/CepGen/Cards/PythonHandler.cpp
@@ -1,399 +1,415 @@
#include "CepGen/Cards/PythonHandler.h"
#include "CepGen/Core/Exception.h"
#ifdef PYTHON
#include "CepGen/Core/TamingFunction.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/ParametersList.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/LHAPDF.h"
#include "CepGen/StructureFunctions/MSTWGrid.h"
#include "CepGen/StructureFunctions/Schaefer.h"
#include "CepGen/Hadronisers/Pythia8Hadroniser.h"
#include <algorithm>
+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 ) );
//--- list of process-specific parameters
ParametersList proc_params;
fillParameter( process, "processParameters", proc_params );
//--- 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 = get<std::string>( pproc_name );
//--- process mode
params_.kinematics.mode = (KinematicsMode)proc_params.get<int>( "mode", (int)KinematicsMode::invalid );
if ( proc_name == "lpair" )
params_.setProcess( new Process::GamGamLL( proc_params ) );
else if ( proc_name == "pptoll" || proc_name == "pptoff" )
params_.setProcess( new Process::PPtoFF( proc_params ) );
else if ( proc_name == "pptoww" )
params_.setProcess( new Process::PPtoWW( proc_params ) );
+ else if ( proc_name == "patoll" )
+ params_.setProcess( new Process::FortranKTProcess( proc_params, "nucltoff", "(p/A)(p/A) ↝ (g/ɣ)ɣ → f⁺f¯", nucl_to_ff_ ) );
else throw CG_FATAL( "PythonHandler" ) << "Unrecognised process: " << proc_name << ".";
//--- 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 );
//--- taming functions
PyObject* ptam = getElement( process, "tamingFunctions" ); // borrowed
if ( ptam )
parseTamingFunctions( ptam );
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 );
}
//--- finalisation
Py_CLEAR( cfg );
}
PythonHandler::~PythonHandler()
{
if ( Py_IsInitialized() )
Py_Finalize();
}
void
PythonHandler::parseIncomingKinematics( PyObject* kin )
{
//--- retrieve the beams PDG ids
std::vector<double> beams_pz;
fillParameter( kin, "pz", beams_pz );
if ( beams_pz.size() == 2 ) {
params_.kinematics.incoming_beams.first.pz = beams_pz.at( 0 );
params_.kinematics.incoming_beams.second.pz = beams_pz.at( 1 );
}
//--- retrieve the beams longitudinal momentum
std::vector<int> beams_pdg;
fillParameter( kin, "pdgIds", beams_pdg );
if ( beams_pdg.size() == 2 ) {
params_.kinematics.incoming_beams.first.pdg = (PDG)beams_pdg.at( 0 );
params_.kinematics.incoming_beams.second.pdg = (PDG)beams_pdg.at( 1 );
}
double sqrt_s = -1.;
fillParameter( kin, "cmEnergy", sqrt_s );
+ fillParameter( kin, "kmrGridPath", params_.kinematics.kmr_grid_path );
if ( sqrt_s != -1. )
params_.kinematics.setSqrtS( sqrt_s );
PyObject* psf = getElement( kin, "structureFunctions" ); // borrowed
if ( psf )
parseStructureFunctions( psf, params_.kinematics.structure_functions );
std::vector<int> kt_fluxes;
fillParameter( kin, "ktFluxes", kt_fluxes );
if ( kt_fluxes.size() > 0 )
params_.kinematics.incoming_beams.first.kt_flux = (KTFlux)kt_fluxes.at( 0 );
if ( kt_fluxes.size() > 1 )
params_.kinematics.incoming_beams.second.kt_flux = (KTFlux)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.pdg = HeavyIon{ (unsigned short)hi_beam1[0], (Element)hi_beam1[1] };
+ fillParameter( kin, "heavyIonB", hi_beam2 );
+ if ( hi_beam2.size() == 2 )
+ params_.kinematics.incoming_beams.second.pdg = HeavyIon{ (unsigned short)hi_beam2[0], (Element)hi_beam2[1] };
}
void
PythonHandler::parseStructureFunctions( PyObject* psf, std::shared_ptr<StructureFunctions>& sf_handler )
{
int str_fun = 0;
fillParameter( psf, "id", str_fun );
sf_handler = StructureFunctionsBuilder::get( (SF::Type)str_fun );
switch( (SF::Type)str_fun ) {
case SF::Type::LHAPDF: {
auto sf = std::dynamic_pointer_cast<SF::LHAPDF>( params_.kinematics.structure_functions );
fillParameter( psf, "pdfSet", sf->params.pdf_set );
fillParameter( psf, "numFlavours", (unsigned int&)sf->params.num_flavours );
fillParameter( psf, "pdfMember", (unsigned int&)sf->params.pdf_member );
fillParameter( psf, "mode", (unsigned int&)sf->params.mode );
} break;
case SF::Type::MSTWgrid: {
auto sf = std::dynamic_pointer_cast<mstw::Grid>( params_.kinematics.structure_functions );
fillParameter( psf, "gridPath", sf->params.grid_path );
} break;
case SF::Type::Schaefer: {
auto sf = std::dynamic_pointer_cast<SF::Schaefer>( params_.kinematics.structure_functions );
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 )
parseStructureFunctions( pcsf, sf->params.continuum_model );
PyObject* ppsf = getElement( psf, "perturbativeSF" ); // borrowed
if ( ppsf )
parseStructureFunctions( ppsf, sf->params.perturbative_model );
PyObject* prsf = getElement( psf, "resonancesSF" ); // borrowed
if ( prsf )
parseStructureFunctions( prsf, sf->params.resonances_model );
fillParameter( psf, "higherTwist", (bool&)sf->params.higher_twist );
} break;
default: break;
}
}
void
PythonHandler::parseOutgoingKinematics( PyObject* kin )
{
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)get<int>( 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)get<int>( 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 = get<std::string>( 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( get<std::string>( pvar ).c_str(), get<std::string>( 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 = get<std::string>( pname );
fillParameter( hadr, "maxTrials", params_.hadroniser_max_trials );
PyObject* pseed = getElement( hadr, "seed" ); // borrowed
long long seed = -1ll;
if ( pseed && is<int>( 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 );
}
}
}
}
#endif
diff --git a/CepGen/Event/Event.cpp b/CepGen/Event/Event.cpp
index f27e6f6..eda5a35 100644
--- a/CepGen/Event/Event.cpp
+++ b/CepGen/Event/Event.cpp
@@ -1,337 +1,341 @@
#include "CepGen/Event/Event.h"
#include "CepGen/Physics/PDG.h"
+#include "CepGen/Physics/HeavyIon.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/utils.h"
#include <algorithm>
#include <math.h>
namespace CepGen
{
Event::Event() :
num_hadronisation_trials( 0 ),
time_generation( -1. ), time_total( -1. )
{}
Event::Event( const Event& rhs ) :
num_hadronisation_trials( rhs.num_hadronisation_trials ),
time_generation( rhs.time_generation ), time_total( rhs.time_total ),
particles_( rhs.particles_ ),
evtcontent_( rhs.evtcontent_ )
{}
void
Event::clear()
{
particles_.clear();
time_generation = -1.;
time_total = -1.;
}
void
Event::freeze()
{
//--- store a snapshot of the primordial event block
if ( particles_.count( Particle::CentralSystem ) > 0 )
evtcontent_.cs = particles_[Particle::CentralSystem].size();
if ( particles_.count( Particle::OutgoingBeam1 ) > 0 )
evtcontent_.op1 = particles_[Particle::OutgoingBeam1].size();
if ( particles_.count( Particle::OutgoingBeam2 ) > 0 )
evtcontent_.op2 = particles_[Particle::OutgoingBeam2].size();
}
void
Event::restore()
{
//--- remove all particles after the primordial event block
if ( particles_.count( Particle::CentralSystem ) > 0 )
particles_[Particle::CentralSystem].resize( evtcontent_.cs );
if ( particles_.count( Particle::OutgoingBeam1 ) > 0 )
particles_[Particle::OutgoingBeam1].resize( evtcontent_.op1 );
if ( particles_.count( Particle::OutgoingBeam2 ) > 0 )
particles_[Particle::OutgoingBeam2].resize( evtcontent_.op2 );
}
double
Event::cmEnergy() const
{
return CMEnergy( getOneByRole( Particle::IncomingBeam1 ), getOneByRole( Particle::IncomingBeam2 ) );
}
Particles&
Event::getByRole( Particle::Role role )
{
//--- retrieve all particles with a given role
return particles_[role];
}
const Particles&
Event::getByRole( Particle::Role role ) const
{
if ( particles_.count( role ) == 0 )
throw CG_FATAL( "Event" ) << "Failed to retrieve a particle with " << role << " role.";
//--- retrieve all particles with a given role
return particles_.at( role );
}
ParticlesIds
Event::getIdsByRole( Particle::Role role ) const
{
ParticlesIds out;
//--- retrieve all particles ids with a given role
if ( particles_.count( role ) == 0 )
return out;
for ( const auto& part : particles_.at( role ) )
out.insert( part.id() );
return out;
}
Particle&
Event::getOneByRole( Particle::Role role )
{
//--- retrieve the first particle a the given role
Particles& parts_by_role = getByRole( role );
if ( parts_by_role.size() == 0 )
throw CG_FATAL( "Event" ) << "No particle retrieved with " << role << " role.";
if ( parts_by_role.size() > 1 )
throw CG_FATAL( "Event" ) << "More than one particle with " << role << " role: "
<< parts_by_role.size() << " particles.";
return *parts_by_role.begin();
}
const Particle&
Event::getOneByRole( Particle::Role role ) const
{
if ( particles_.count( role ) == 0 )
throw CG_FATAL( "Event" ) << "Failed to retrieve a particle with " << role << " role.";
//--- retrieve the first particle a the given role
const Particles& parts_by_role = particles_.at( role );
if ( parts_by_role.size() == 0 )
throw CG_FATAL( "Event" ) << "No particle retrieved with " << role << " role.";
if ( parts_by_role.size() > 1 )
throw CG_FATAL( "Event" ) << "More than one particle with " << role << " role: "
<< parts_by_role.size() << " particles";
return *parts_by_role.begin();
}
Particle&
Event::operator[]( int id )
{
for ( auto& role_part : particles_ )
for ( auto& part : role_part.second )
if ( part.id() == id )
return part;
throw CG_FATAL( "Event" ) << "Failed to retrieve the particle with id=" << id << ".";
}
const Particle&
Event::getConstById( int id ) const
{
for ( const auto& role_part : particles_ )
for ( const auto& part : role_part.second )
if ( part.id() == id )
return part;
throw CG_FATAL( "Event" ) << "Failed to retrieve the particle with id=" << id << ".";
}
Particles
Event::getByIds( const ParticlesIds& ids ) const
{
Particles out;
for ( const auto& id : ids )
out.emplace_back( getConstById( id ) );
return out;
}
Particles
Event::mothers( const Particle& part )
{
return getByIds( part.mothers() );
}
Particles
Event::daughters( const Particle& part )
{
return getByIds( part.daughters() );
}
ParticleRoles
Event::roles() const
{
ParticleRoles out;
ParticlesMap::const_iterator it, end = particles_.end();
for ( it = particles_.begin(); it != end; it = particles_.upper_bound( it->first ) ) {
out.emplace_back( it->first );
}
return out;
}
Particle&
Event::addParticle( Particle& part, bool replace )
{
CG_DEBUG_LOOP( "Event" ) << "Particle with PDGid = " << part.integerPdgId() << " has role " << part.role();
if ( part.role() <= 0 )
throw CG_FATAL( "Event" ) << "Trying to add a particle with role=" << (int)part.role() << ".";
//--- retrieve the list of particles with the same role
Particles& part_with_same_role = getByRole( part.role() );
//--- specify the id
if ( part_with_same_role.empty() && part.id() < 0 ) part.setId( numParticles() ); // set the id if previously invalid/inexistent
if ( !part_with_same_role.empty() ) {
if ( replace ) part.setId( part_with_same_role[0].id() ); // set the previous id if replacing a particle
else part.setId( numParticles() );
}
//--- add the particle to the collection
if ( replace ) part_with_same_role = Particles( 1, part ); // generate a vector containing only this particle
else part_with_same_role.emplace_back( part );
return part_with_same_role.back();
}
Particle&
Event::addParticle( Particle::Role role, bool replace )
{
Particle np( role, PDG::invalid );
return addParticle( np, replace );
}
size_t
Event::numParticles() const
{
size_t out = 0;
for ( const auto& role_part : particles_ )
out += role_part.second.size();
return out;
}
const Particles
Event::particles() const
{
Particles out;
for ( const auto& role_part : particles_ )
out.insert( out.end(), role_part.second.begin(), role_part.second.end() );
std::sort( out.begin(), out.end() );
return out;
}
const Particles
Event::stableParticles() const
{
Particles out;
for ( const auto& role_part : particles_ )
for ( const auto& part : role_part.second )
if ( (short)part.status() > 0 )
out.emplace_back( part );
std::sort( out.begin(), out.end() );
return out;
}
void
Event::checkKinematics() const
{
// check the kinematics through parentage
for ( const auto& part : particles() ) {
ParticlesIds daughters = part.daughters();
if ( daughters.empty() )
continue;
Particle::Momentum ptot;
for ( const auto& daugh : daughters ) {
const Particle& d = getConstById( daugh );
const ParticlesIds mothers = d.mothers();
ptot += d.momentum();
if ( mothers.size() < 2 )
continue;
for ( const auto& moth : mothers ) {
if ( moth == part.id() )
continue;
ptot -= getConstById( moth ).momentum();
}
}
const double mass_diff = ( ptot-part.momentum() ).mass();
if ( fabs( mass_diff ) > minimal_precision_ ) {
dump();
throw CG_FATAL( "Event" ) << "Error in momentum balance for particle " << part.id() << ": mdiff = " << mass_diff << ".";
}
}
}
void
Event::dump( std::ostream& out, bool stable ) const
{
const Particles parts = ( stable ) ? stableParticles() : particles();
std::ostringstream os;
Particle::Momentum p_total;
for ( const auto& part : parts ) {
const ParticlesIds mothers = part.mothers();
{
std::ostringstream oss_pdg;
if ( part.pdgId() == PDG::invalid && mothers.size() > 0 ) {
for ( unsigned short i = 0; i < mothers.size(); ++i )
oss_pdg << ( i > 0 ? "/" : "" ) << getConstById( *std::next( mothers.begin(), i ) ).pdgId();
os << Form( "\n %2d\t\t%-10s", part.id(), oss_pdg.str().c_str() );
}
else {
- oss_pdg << part.pdgId();
+ if ( (HeavyIon)part.pdgId() )
+ oss_pdg << (HeavyIon)part.pdgId();
+ else
+ oss_pdg << part.pdgId();
os << Form( "\n %2d\t%-+7d %-10s", part.id(), part.integerPdgId(), oss_pdg.str().c_str() );
}
}
os << "\t";
if ( part.charge() != 999. )
os << Form( "%-.2f\t", part.charge() );
else
os << "\t";
{ std::ostringstream oss; oss << part.role(); os << Form( "%-8s %6d\t", oss.str().c_str(), part.status() ); }
if ( !mothers.empty() ) {
std::ostringstream oss;
unsigned short i = 0;
for ( const auto& moth : mothers ) {
oss << ( i > 0 ? "+" : "" ) << moth;
++i;
}
os << Form( "%6s ", oss.str().c_str() );
}
else os << " ";
const Particle::Momentum mom = part.momentum();
os << Form( "% 9.6e % 9.6e % 9.6e % 9.6e % 12.5f", mom.px(), mom.py(), mom.pz(), part.energy(), part.mass() );
// discard non-primary, decayed particles
if ( part.status() >= Particle::Status::Undefined ) {
const int sign = ( part.status() == Particle::Status::Undefined )
? -1
: +1;
p_total += sign*mom;
}
}
//--- set a threshold to the computation precision
p_total.truncate();
//
CG_INFO( "Event" )
<< Form( "Dump of event content:\n"
" Id\tPDG id\tName\t\tCharge\tRole\t Status\tMother\tpx py pz E \t M \n"
" --\t------\t----\t\t------\t----\t ------\t------\t----GeV/c--- ----GeV/c--- ----GeV/c--- ----GeV/c---\t --GeV/c²--"
"%s\n"
" ----------------------------------------------------------------------------------------------------------------------------------\n"
"\t\t\t\t\t\t\tBalance% 9.6e % 9.6e % 9.6e % 9.6e", os.str().c_str(), p_total.px(), p_total.py(), p_total.pz(), p_total.energy() );
}
//------------------------------------------------------------------------------------------------
Event::NumParticles::NumParticles() :
cs( 0 ), op1( 0 ), op2( 0 )
{}
Event::NumParticles::NumParticles( const NumParticles& np ) :
cs( np.cs ), op1( np.op1 ), op2( np.op2 )
{}
}
diff --git a/CepGen/IO/FortranInterface.cpp b/CepGen/IO/FortranInterface.cpp
index e7c6d0f..19dc902 100644
--- a/CepGen/IO/FortranInterface.cpp
+++ b/CepGen/IO/FortranInterface.cpp
@@ -1,58 +1,67 @@
#include "CepGen/StructureFunctions/StructureFunctionsBuilder.h"
#include "CepGen/StructureFunctions/StructureFunctions.h"
#include "CepGen/Physics/KTFlux.h"
+#include "CepGen/Physics/HeavyIon.h"
#include "CepGen/Physics/ParticleProperties.h"
#include "CepGen/Core/Exception.h"
#ifdef __cplusplus
extern "C" {
#endif
void
cepgen_structure_functions_( int& sfmode, double& xbj, double& q2, double& f2, double& fl )
{
using namespace CepGen;
SF::Type sf_mode = (SF::Type)sfmode;
CG_DEBUG( "cepgen_structure_functions" ) << sf_mode;
static StructureFunctions& val = StructureFunctionsBuilder::get( sf_mode )->operator()( xbj, q2 );
f2 = val.F2;
fl = val.FL;
}
double
cepgen_kt_flux_( int& fmode, double& x, double& kt2, int& sfmode, double& mx )
{
using namespace CepGen;
static auto sf = StructureFunctionsBuilder::get( (SF::Type)sfmode );
return ktFlux(
(KTFlux)fmode, x, kt2, *sf, mx );
}
double
+ cepgen_kt_flux_hi_( int& fmode, double& x, double& kt2, int& a, int& z )
+ {
+ using namespace CepGen;
+ return ktFlux(
+ (KTFlux)fmode, x, kt2, HeavyIon{ (unsigned short)a, (Element)z } );
+ }
+
+ double
cepgen_particle_mass_( int& pdg_id )
{
try {
return CepGen::ParticleProperties::mass( (CepGen::PDG)pdg_id );
} catch ( const CepGen::Exception& e ) {
e.dump();
exit( 0 );
}
}
double
cepgen_particle_charge_( int& pdg_id )
{
try {
return CepGen::ParticleProperties::charge( pdg_id );
} catch ( const CepGen::Exception& e ) {
e.dump();
exit( 0 );
}
}
#ifdef __cplusplus
}
#endif
diff --git a/CepGen/Physics/GluonGrid.cpp b/CepGen/Physics/GluonGrid.cpp
new file mode 100644
index 0000000..bcc25e7
--- /dev/null
+++ b/CepGen/Physics/GluonGrid.cpp
@@ -0,0 +1,55 @@
+#include "CepGen/Physics/GluonGrid.h"
+#include "CepGen/Core/Exception.h"
+
+#include <fstream>
+#include <set>
+
+namespace kmr
+{
+ GluonGrid&
+ GluonGrid::get( const char* filename )
+ {
+ Parameterisation p;
+ p.grid_path = filename;
+ static GluonGrid instance( p );
+ return instance;
+ }
+
+ GluonGrid::GluonGrid( const Parameterisation& param ) :
+ CepGen::GridHandler<3,1>( CepGen::GridType::linear ),
+ params( param )
+ {
+ std::set<double> kt2_vals, x_vals, mu2_vals;
+
+ { // file readout part
+ std::ifstream file( params.grid_path, std::ios::in );
+ if ( !file.is_open() )
+ throw CG_FATAL( "GluonGrid" ) << "Impossible to load grid file \"" << params.grid_path << "\"!";
+
+ std::string x_tmp, kt2_tmp, mu2_tmp, fg_tmp;
+ while ( file >> x_tmp >> kt2_tmp >> mu2_tmp >> fg_tmp ) {
+ const double x = stod( x_tmp ), kt2 = stod( kt2_tmp ), mu2 = stod( mu2_tmp ), fg = stod( fg_tmp );
+ x_vals.insert( x );
+ kt2_vals.insert( kt2 );
+ mu2_vals.insert( mu2 );
+ insert( { x, kt2, mu2 }, { fg } );
+ }
+ file.close();
+ }
+
+ init();
+
+ CG_INFO( "GluonGrid" )
+ << "KMR grid evaluator built!\n\t"
+ << " kt² in range [" << *kt2_vals.begin() << ":" << *kt2_vals.rbegin() << "]\n\t"
+ << " x in range [" << *x_vals.begin() << ":" << *x_vals.rbegin() << "]\n\t"
+ << " µ² in range [" << *mu2_vals.begin() << ":" << *mu2_vals.rbegin() << "].";
+ }
+
+ double
+ GluonGrid::operator()( double x, double kt2, double mu2 ) const
+ {
+ return CepGen::GridHandler<3,1>::eval( { x, kt2, mu2 } ).at( 0 );
+ }
+}
+
diff --git a/CepGen/Physics/GluonGrid.h b/CepGen/Physics/GluonGrid.h
new file mode 100644
index 0000000..5580ae5
--- /dev/null
+++ b/CepGen/Physics/GluonGrid.h
@@ -0,0 +1,40 @@
+#ifndef CepGen_Physics_GluonGrid_h
+#define CepGen_Physics_GluonGrid_h
+
+#include "CepGen/IO/GridHandler.h"
+
+#define DEFAULT_KMR_GRID_PATH "gluon_mmht2014nlo_Watt.dat"
+
+/// Kimber-Martin-Ryskin unintegrated gluon densities
+namespace kmr
+{
+ /// A KMR unintegrated gluon densities grid interpolator
+ class GluonGrid : private CepGen::GridHandler<3,1>
+ {
+ public:
+ struct Parameterisation {
+ Parameterisation() : grid_path( DEFAULT_KMR_GRID_PATH ) {}
+ std::string grid_path;
+ };
+
+ public:
+ /// Retrieve the grid interpolator (singleton)
+ static GluonGrid& get( const char* path = DEFAULT_KMR_GRID_PATH );
+
+ /// Compute the gluon flux
+ double operator()( double x, double kt2, double mu2 ) const;
+ Parameterisation params;
+
+ public:
+ GluonGrid( const GluonGrid& ) = delete;
+ void operator=( const GridHandler& ) = delete;
+
+ private:
+ explicit GluonGrid( const Parameterisation& = Parameterisation() );
+ };
+}
+
+#undef DEFAULT_KMR_GRID_PATH
+
+#endif
+
diff --git a/CepGen/Physics/HeavyIon.cpp b/CepGen/Physics/HeavyIon.cpp
new file mode 100644
index 0000000..2477e74
--- /dev/null
+++ b/CepGen/Physics/HeavyIon.cpp
@@ -0,0 +1,46 @@
+#include "CepGen/Physics/HeavyIon.h"
+
+#include <sstream>
+
+namespace CepGen
+{
+ HeavyIon::HeavyIon( const PDG& pdg ) :
+ A( (unsigned int)pdg % 1000 ),
+ Z( (Element)( (unsigned int)pdg/1000000 == 0 ? 0
+ : ( (unsigned int)pdg/1000 ) % 1000 ) )
+ {}
+
+ HeavyIon::operator PDG() const
+ {
+ // Pythia8 convention/10-1e10+1e6
+ return (PDG)( 1000000+1000*(unsigned short)Z+A );
+ }
+
+ std::ostream&
+ operator<<( std::ostream& os, const HeavyIon& hi )
+ {
+ std::ostringstream oss; oss << hi.Z;
+ if ( oss.str().empty() )
+ return os << "HI{Z=" << (unsigned short)hi.Z << ", A=" << hi.A << "}";
+ return os << hi.A << oss.str();
+ }
+
+ std::ostream&
+ operator<<( std::ostream& os, const Element& elem )
+ {
+ switch ( elem ) {
+ case Element::invalid:
+ return os;
+ case Element::H: return os << "H";
+ case Element::C: return os << "C";
+ case Element::O: return os << "O";
+ case Element::Al: return os << "Al";
+ case Element::Cu: return os << "Cu";
+ case Element::Xe: return os << "Xe";
+ case Element::Au: return os << "Au";
+ case Element::Pb: return os << "Pb";
+ case Element::U: return os << "U";
+ }
+ return os;
+ }
+}
diff --git a/CepGen/Physics/HeavyIon.h b/CepGen/Physics/HeavyIon.h
new file mode 100644
index 0000000..e2b996c
--- /dev/null
+++ b/CepGen/Physics/HeavyIon.h
@@ -0,0 +1,44 @@
+#ifndef CepGen_Physics_HeavyIon_h
+#define CepGen_Physics_HeavyIon_h
+
+#include <ostream>
+
+namespace CepGen
+{
+ enum class PDG;
+ /// Enumeration of chemical elements
+ enum class Element
+ {
+ invalid = 0,
+ H = 1, C = 6, O = 8,
+ Al = 13, Cu = 29,
+ Xe = 54, Au = 79, Pb = 82,
+ U = 92
+ };
+ std::ostream& operator<<( std::ostream& os, const Element& elem );
+
+ /// Heavy ion container (Z+A)
+ struct HeavyIon
+ {
+ /// General constructor from mass and atomic number
+ HeavyIon( unsigned short a, const Element& z ) : A( a ), Z( z ) {}
+ /// Build from a custom PDG id
+ HeavyIon( const PDG& pdg );
+ /// Simple proton
+ static inline HeavyIon proton() { return HeavyIon( 1, Element::H ); }
+ /// Convert the HI into a custom PDG id
+ operator PDG() const;
+ /// Check the validity of the heavy ion
+ inline operator bool() const {
+ return Z > Element::invalid && A > 1; // skip the proton
+ }
+ /// Human-readable expression of the ion
+ friend std::ostream& operator<<( std::ostream& os, const HeavyIon& hi );
+ /// Mass number
+ unsigned short A;
+ /// Atomic number
+ Element Z;
+ };
+}
+
+#endif
diff --git a/CepGen/Physics/KTFlux.cpp b/CepGen/Physics/KTFlux.cpp
index 5702fc7..65064b8 100644
--- a/CepGen/Physics/KTFlux.cpp
+++ b/CepGen/Physics/KTFlux.cpp
@@ -1,68 +1,103 @@
#include "CepGen/Physics/KTFlux.h"
#include "CepGen/Physics/FormFactors.h"
#include "CepGen/Physics/ParticleProperties.h"
+#include "CepGen/Physics/HeavyIon.h"
+#include "CepGen/Physics/GluonGrid.h"
+
#include "CepGen/StructureFunctions/StructureFunctions.h"
#include "CepGen/Core/Exception.h"
namespace CepGen
{
const double KTFluxParameters::kMinKTFlux = 1.e-20;
const double KTFluxParameters::kMP = ParticleProperties::mass( PDG::proton );
const double KTFluxParameters::kMP2 = KTFluxParameters::kMP*KTFluxParameters::kMP;
double
ktFlux( const KTFlux& type, double x, double kt2, StructureFunctions& sf, double mx )
{
double flux = 0.;
const double mp2 = KTFluxParameters::kMP2;
switch ( type ) {
case KTFlux::P_Photon_Elastic: {
const double x2 = x*x;
const double q2min = x2*mp2/( 1.-x ), q2 = q2min + kt2/( 1.-x );
// electromagnetic form factors
const auto& ff = FormFactors::protonElastic( q2 );
const double ela1 = ( 1.-x )*( 1.-q2min/q2 );
//const double ela3 = 1.-( q2-kt2 )/q2;
flux = Constants::alphaEM*M_1_PI*( 1.-x )/q2*( ela1*ff.FE + 0.5*x2*ff.FM );
} break;
case KTFlux::P_Photon_Inelastic_Budnev: {
const double mx2 = mx*mx, x2 = x*x;
const double q2min = ( x*( mx2-mp2 ) + x2*mp2 )/( 1.-x ), q2 = q2min + kt2/( 1.-x );
const double xbj = q2 / ( q2+mx2-mp2 );
// structure functions
auto& str_fun = sf( xbj, q2 );
str_fun.computeFL( xbj, q2 );
const double f_D = str_fun.F2/( q2+mx2-mp2 )* ( 1.-x )*( 1.-q2min/q2 );
const double f_C = str_fun.F1( xbj, q2 ) * 2./q2;
flux = Constants::alphaEM*M_1_PI*( 1.-x )/q2*( f_D+0.5*x2*f_C );
} break;
+ case KTFlux::P_Gluon_KMR: {
+ flux = kmr::GluonGrid::get()( log10( x ), log10( kt2 ), 2.*log10( mx ) );
+ } break;
default:
throw CG_FATAL( "GenericKTProcess:flux" ) << "Invalid flux type: " << type;
}
if ( flux < KTFluxParameters::kMinKTFlux )
return 0.;
return flux;
}
+ double
+ ktFlux( const KTFlux& type, double x, double kt2, const HeavyIon& hi )
+ {
+ double flux = 0.;
+ switch ( type ) {
+ case KTFlux::HI_Photon_Elastic: {
+ const double r_a = 1.1*std::pow( hi.A, 1./3 ), a0 = 0.7, m_a = hi.A*KTFluxParameters::kMP;
+ const double q2_ela = ( kt2+x*x*m_a*m_a )/( 1.-x ), cons = sqrt( q2_ela )/0.1973;
+ const double tau = cons*r_a, tau1 = cons*a0;
+ // "Realistic nuclear form-factor" as used in STARLIGHT
+ const double ff1 = 3.*( sin( tau )-tau*cos( tau ) )/pow( tau+1.e-10, 3 );
+ const double ff2 = 1./( 1.+tau1*tau1 );
+ const double ela1 = pow( kt2/( kt2+x*x*m_a*m_a ), 2 );
+ const double ela2 = pow( ff1*ff2, 2 )/*, ela3 = 1.-( q2_ela-kt2 )/q2_ela*/;
+ const unsigned int z = (unsigned short)hi.Z;
+ flux = Constants::alphaEM*M_1_PI*z*z*ela1*ela2/q2_ela;
+ } break;
+ default:
+ throw CG_FATAL("GenericKTProcess:flux") << "Invalid flux type: " << type;
+ }
+ if ( flux < KTFluxParameters::kMinKTFlux )
+ return 0.;
+ return flux;
+ }
+
std::ostream&
operator<<( std::ostream& os, const KTFlux& type )
{
switch ( type ) {
case KTFlux::P_Photon_Elastic:
return os << "elastic photon from proton";
case KTFlux::P_Photon_Inelastic:
return os << "inelastic photon from proton";
case KTFlux::P_Photon_Inelastic_Budnev:
return os << "inelastic photon from proton (Budnev)";
+ case KTFlux::P_Gluon_KMR:
+ return os << "elastic gluon from proton (KMR)";
+ case KTFlux::HI_Photon_Elastic:
+ return os << "elastic photon from HI";
case KTFlux::invalid: default:
return os << "unrecognized flux (" << (int)type << ")";
}
}
}
diff --git a/CepGen/Physics/KTFlux.h b/CepGen/Physics/KTFlux.h
index dce02c5..db8df10 100644
--- a/CepGen/Physics/KTFlux.h
+++ b/CepGen/Physics/KTFlux.h
@@ -1,31 +1,35 @@
#ifndef CepGen_Physics_KTFlux_h
#define CepGen_Physics_KTFlux_h
#include "CepGen/Physics/PDG.h"
#include <ostream>
namespace CepGen
{
class StructureFunctions;
+ class HeavyIon;
struct KTFluxParameters
{
static const double kMinKTFlux, kMP, kMP2;
};
/// Type of incoming partons fluxes
enum class KTFlux
{
invalid = -1,
P_Photon_Elastic = 0,
P_Photon_Inelastic = 1,
P_Photon_Inelastic_Budnev = 11,
+ P_Gluon_KMR = 20,
+ HI_Photon_Elastic = 100
};
std::ostream& operator<<( std::ostream&, const KTFlux& );
/// Get the flux at a given parton x/kT
/// \param[in] x Parton momentum fraction
/// \param[in] kt2 Transverse 2-momentum \f$\mathbf{q}_{\mathrm{T}}^2\f$ of the incoming parton
/// \param[in] sf Structure functions evaluator
/// \param[in] mx Outgoing diffractive proton mass
double ktFlux( const KTFlux& type, double x, double kt2, StructureFunctions& sf, double mx = KTFluxParameters::kMP );
+ double ktFlux( const KTFlux& type, double x, double kt2, const HeavyIon& hi );
}
#endif
diff --git a/CepGen/Physics/Kinematics.cpp b/CepGen/Physics/Kinematics.cpp
index 90050a5..1d225a2 100644
--- a/CepGen/Physics/Kinematics.cpp
+++ b/CepGen/Physics/Kinematics.cpp
@@ -1,85 +1,89 @@
#include "CepGen/Physics/Kinematics.h"
#include "CepGen/Physics/PDG.h"
#include "CepGen/Physics/KTFlux.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() :
incoming_beams( { { 6500., PDG::proton, KTFlux::invalid }, { 6500., PDG::proton, KTFlux::invalid } } ),
mode( KinematicsMode::invalid ), structure_functions( new SF::SuriYennie )
{}
Kinematics::~Kinematics()
{}
void
Kinematics::setSqrtS( double sqrts )
{
incoming_beams.first.pz = incoming_beams.second.pz = 0.5 * sqrts;
}
double
Kinematics::sqrtS() const
{
return incoming_beams.first.pz + incoming_beams.second.pz;
}
//--------------------------------------------------------------------
// User-friendly display of the kinematics mode
//--------------------------------------------------------------------
std::ostream&
operator<<( std::ostream& os, const KinematicsMode& pm )
{
switch ( pm ) {
case KinematicsMode::invalid:
return os << "invalid";
case KinematicsMode::ElectronElectron:
return os << "electron/electron";
case KinematicsMode::ElectronProton:
return os << "electron/proton";
case KinematicsMode::ProtonElectron:
return os << "proton/electron";
case KinematicsMode::ElasticElastic:
return os << "elastic/elastic";
case KinematicsMode::InelasticElastic:
return os << "inelastic/elastic";
case KinematicsMode::ElasticInelastic:
return os << "elastic/inelastic";
case KinematicsMode::InelasticInelastic:
return os << "inelastic/inelastic";
}
return os;
}
//--------------------------------------------------------------------
// User-friendly display of incoming particles
//--------------------------------------------------------------------
std::ostream&
operator<<( std::ostream& os, const Kinematics::Beam& beam )
{
- os << beam.pdg << " (" << beam.pz << " GeV/c)";
+ if ( (HeavyIon)beam.pdg )
+ os << (HeavyIon)beam.pdg;
+ else
+ os << beam.pdg;
+ os << " (" << beam.pz << " GeV/c)";
if ( beam.kt_flux != KTFlux::invalid )
os << " [unint.flux: " << beam.kt_flux << "]";
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 fadb57d..bafb8b8 100644
--- a/CepGen/Physics/Kinematics.h
+++ b/CepGen/Physics/Kinematics.h
@@ -1,77 +1,79 @@
#ifndef CepGen_Physics_Kinematics_h
#define CepGen_Physics_Kinematics_h
#include "CepGen/Core/Hasher.h"
#include "CepGen/Physics/Cuts.h"
+#include "CepGen/Physics/HeavyIon.h"
#include <ostream>
#include <vector>
#include <unordered_map>
#include <memory>
namespace CepGen
{
enum class PDG;
enum class KTFlux;
class StructureFunctions;
/// Type of kinematics to consider for the process
enum class KinematicsMode
{
invalid = -1,
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)
std::ostream& operator<<( std::ostream&, const KinematicsMode& );
/// List of kinematic constraints to apply on the process phase space.
class Kinematics
{
public:
Kinematics();
~Kinematics();
struct Beam
{
/// Incoming particle's momentum (in \f$\text{GeV}/c\f$)
double pz;
PDG pdg;
KTFlux kt_flux;
};
friend std::ostream& operator<<( std::ostream&, const Beam& );
/// Beam/primary particle's kinematics
std::pair<Beam,Beam> incoming_beams;
/// Set the incoming particles' momenta (if the collision is symmetric)
void setSqrtS( double sqrts );
/// Process centre of mass energy
double sqrtS() const;
/// Minimum list of central particles required
std::vector<PDG> minimum_final_state;
/// Type of kinematics to consider for the phase space
KinematicsMode mode;
/// Type of structure functions to consider
std::shared_ptr<StructureFunctions> structure_functions;
struct CutsList
{
CutsList();
/// 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;
+ std::string kmr_grid_path;
};
}
#endif
diff --git a/CepGen/Processes/CMakeLists.txt b/CepGen/Processes/CMakeLists.txt
index 017c1b7..8050ba9 100644
--- a/CepGen/Processes/CMakeLists.txt
+++ b/CepGen/Processes/CMakeLists.txt
@@ -1,7 +1,7 @@
-file(GLOB sources *.cpp)
+file(GLOB sources *.cpp Fortran/*.f)
include_directories(${PROJECT_SOURCE_DIR})
add_library(CepGenProcesses SHARED ${sources})
install(TARGETS CepGenProcesses DESTINATION lib)
diff --git a/CepGen/Processes/Fortran/cepgen_blocks.inc b/CepGen/Processes/Fortran/cepgen_blocks.inc
new file mode 100644
index 0000000..ff15b10
--- /dev/null
+++ b/CepGen/Processes/Fortran/cepgen_blocks.inc
@@ -0,0 +1,61 @@
+c ================================================================
+c F77 process-to-CepGen helper
+c > this header file defines a set of common blocks useful for
+c > the interfacing of any kt-factorised matrix element computation
+c > to a mother CepGen instance
+c ================================================================
+
+c =================================================================
+c collection of fundamental physics constants
+c =================================================================
+ common/constants/am_p,units,pi,alpha_em
+ double precision am_p,units,pi,alpha_em
+
+c =================================================================
+c information on the full process
+c inp1 = proton energy in lab frame
+c inp2 = nucleus energy **per nucleon** in LAB frame
+c Collision is along z-axis
+c =================================================================
+ common/params/icontri,iflux1,iflux2,imethod,sfmod,pdg_l,
+ & a_nuc1,z_nuc1,a_nuc2,z_nuc2,
+ & inp1,inp2
+ integer icontri,iflux1,iflux2,imethod,sfmod,pdg_l,
+ & a_nuc1,z_nuc1,a_nuc2,z_nuc2
+ double precision inp1,inp2
+
+c =================================================================
+c kt-factorisation kinematics
+c =================================================================
+ common/ktkin/q1t,q2t,phiq1t,phiq2t,y1,y2,ptdiff,phiptdiff,
+ & am_x,am_y
+ double precision q1t,q2t,phiq1t,phiq2t,y1,y2,ptdiff,phiptdiff,
+ & am_x,am_y
+
+c =================================================================
+c phase space cuts
+c =================================================================
+ common/kincuts/ipt,iene,ieta,iinvm,iptsum,idely,
+ & pt_min,pt_max,ene_min,ene_max,eta_min,eta_max,
+ & invm_min,invm_max,ptsum_min,ptsum_max,
+ & dely_min,dely_max
+ integer ipt,iene,ieta,iinvm,iptsum,idely
+ double precision pt_min,pt_max,ene_min,ene_max,eta_min,eta_max,
+ & invm_min,invm_max,ptsum_min,ptsum_max,
+ & dely_min,dely_max
+
+c =================================================================
+c generated event kinematics
+c =================================================================
+ common/evtkin/px,py,nout,idum,ipdg,pc
+ integer nout,idum,ipdg(4)
+ double precision px(4),py(4),pc(4,4)
+
+c =================================================================
+c helpers for the evaluation of unintegrated parton fluxes
+c =================================================================
+ external CepGen_kT_flux,CepGen_kT_flux_HI
+ double precision CepGen_kT_flux,CepGen_kT_flux_HI
+ external CepGen_particle_charge,CepGen_particle_mass
+ double precision CepGen_particle_charge,CepGen_particle_mass
+
diff --git a/CepGen/Processes/Fortran/cepgen_print.f b/CepGen/Processes/Fortran/cepgen_print.f
new file mode 100644
index 0000000..38954f6
--- /dev/null
+++ b/CepGen/Processes/Fortran/cepgen_print.f
@@ -0,0 +1,40 @@
+ subroutine cepgen_print
+ implicit none
+ include 'cepgen_blocks.inc'
+ logical params_shown
+ data params_shown/.false./
+ save params_shown
+
+ if(params_shown) return
+
+ print *,'========================================================'
+ print *,'Parameter value(s)'
+ print *,'--------------------------------------------------------'
+ print 101,'Process mode:',icontri
+ print 101,'Computation method:',imethod
+ print 101,'Structure functions:',sfmod
+ print 101,'Central system PDG:',pdg_l
+ print 103,'Beams momenta:',inp1,inp2
+ print 102,'Fluxes modes:',iflux1,iflux2
+ print 104,'Beams (A,Z):',a_nuc1,z_nuc1,a_nuc2,z_nuc2
+ print *,'========================================================'
+ print *,'Cut enabled minimum maximum'
+ print *,'--------------------------------------------------------'
+ print 100,'pt(single)',ipt,pt_min,pt_max
+ print 100,'energy(single)',iene,ene_min,ene_max
+ print 100,'eta(single)',ieta,eta_min,eta_max
+ print 100,'m(sum)',iinvm,invm_min,invm_max
+ print 100,'pt(sum)',iptsum,ptsum_min,ptsum_max
+ print 100,'delta(y)',idely,dely_min,dely_max
+ print *,'========================================================'
+
+ params_shown=.true.
+
+100 format(A26,' ',L2,f12.4,f12.4)
+101 format(A33,I12)
+102 format(A33,I12,I12)
+103 format(A33,f12.2,f12.2)
+104 format(A33,' (',I3,',',I3,'), (',I3,',',I3,')')
+
+ end
+
diff --git a/CepGen/Processes/Fortran/nucl_to_ff.f b/CepGen/Processes/Fortran/nucl_to_ff.f
new file mode 100644
index 0000000..a699f74
--- /dev/null
+++ b/CepGen/Processes/Fortran/nucl_to_ff.f
@@ -0,0 +1,490 @@
+ subroutine nucl_to_ff(aintegrand)
+ implicit none
+ double precision aintegrand
+c =================================================================
+c CepGen common blocks for kinematics definition
+c =================================================================
+ include 'cepgen_blocks.inc'
+ data iflux1,iflux2,sfmod,pdg_l/10,100,11,13/
+ data a_nuc1,z_nuc1,a_nuc2,z_nuc2/1,1,208,82/
+
+c =================================================================
+c local variables
+c =================================================================
+ double precision s,s12
+ double precision alpha1,alpha2,amt1,amt2
+ double precision pt1,pt2,eta1,eta2,dely
+ double precision pt1x,pt1y,pt2x,pt2y
+ double precision ak10,ak1x,ak1y,ak1z,ak20,ak2x,ak2y,ak2z
+ double precision beta1,beta2,x1,x2
+ double precision z1p,z1m,z2p,z2m
+ double precision q1tx,q1ty,q1z,q10,q1t2
+ double precision q2tx,q2ty,q2z,q20,q2t2
+ double precision ptsum
+ double precision that1,that2,that,uhat1,uhat2,uhat
+ double precision term1,term2,term3,term4,term5,term6,term7
+ double precision term8,term9,term10,amat2
+ double precision ak1_x,ak1_y,ak2_x,ak2_y
+ double precision eps12,eps22
+ double precision aux2_1,aux2_2,f1,f2
+ double precision Phi10,Phi102,Phi11_x,Phi11_y,Phi112
+ double precision Phi20,Phi202,Phi21_x,Phi21_y,Phi212
+ double precision Phi11_dot_e,Phi11_cross_e
+ double precision Phi21_dot_e,Phi21_cross_e
+ double precision aintegral
+ integer imat1,imat2
+
+ double precision px_plus,px_minus,py_plus,py_minus
+ double precision r1,r2
+ double precision ptdiffx,ptsumx,ptdiffy,ptsumy
+ double precision invm,invm2,s1_eff,s2_eff
+ double precision t1abs,t2abs
+ double precision am_l,q_l
+ double precision inp_A,inp_B,am_A,am_B
+ double precision p1_plus, p2_minus
+ double precision amat2_1,amat2_2
+ integer iterm11,iterm22,iterm12,itermtt
+
+ double precision coupling
+
+c =================================================================
+c quarks production
+c =================================================================
+#ifdef ALPHA_S
+ double precision t_max,amu2,alphas
+#endif
+ logical first_init
+ data first_init/.true./
+ save first_init,am_l,q_l
+
+ call CepGen_print
+
+ if(first_init) then
+ am_l = CepGen_particle_mass(pdg_l) ! central particles mass
+ q_l = CepGen_particle_charge(pdg_l) ! central particles charge
+ if(iflux1.ge.20.and.iflux1.lt.40) then
+ if(icontri.eq.3.or.icontri.eq.4) then
+ print *,'Invalid process mode for collinear gluon emission!'
+ stop
+ endif
+#ifdef ALPHA_S
+ print *,'Initialisation of the alpha(S) evolution algorithm..'
+ call initAlphaS(0,1.d0,1.d0,0.5d0,
+ & CepGen_particle_mass(4), ! charm
+ & CepGen_particle_mass(5), ! bottom
+ & CepGen_particle_mass(6)) ! top
+#endif
+ endif
+ first_init = .false.
+ endif
+
+c =================================================================
+c FIXME
+c =================================================================
+
+ aintegrand = 0.d0
+ q10 = 0.d0
+ q1z = 0.d0
+ q20 = 0.d0
+ q2z = 0.d0
+
+c =================================================================
+c go to energy of Nucleus = A*inp in order that generator puts out
+c proper momenta in LAB frame
+c =================================================================
+
+ inp_A = inp1*a_nuc1
+ am_A = am_p*a_nuc1
+ inp_B = inp2*a_nuc2
+ am_B = am_p*a_nuc2
+
+c =================================================================
+c four-momenta for incoming beams in LAB !!!!
+c =================================================================
+
+ r1 = dsqrt(1.d0+am_A**2/inp_A**2)
+ r2 = dsqrt(1.d0+am_B**2/inp_B**2)
+
+ ak10 = inp_A*r1
+ ak1x = 0.d0
+ ak1y = 0.d0
+ ak1z = inp_A
+
+ ak20 = inp_B*r2
+ ak2x = 0.d0
+ ak2y = 0.d0
+ ak2z = -inp_B
+
+ s = 4.*inp_A*inp_B*(1.d0 + r1*r2)/2.d0+(am_A**2+am_B**2)
+ s12 = dsqrt(s)
+
+ p1_plus = (ak10+ak1z)/dsqrt(2.d0)
+ p2_minus = (ak20-ak2z)/dsqrt(2.d0)
+
+c terms in the matrix element
+c
+ iterm11 = 1 ! LL
+ iterm22 = 1 ! TT
+ iterm12 = 1 ! LT
+ itermtt = 1 ! TT'
+
+c =================================================================
+c two terms in the Wolfgang's formula for
+c off-shell gamma gamma --> l^+ l^-
+c =================================================================
+ imat1 = 1
+ imat2 = 1
+c =================================================================
+c Outgoing proton final state's mass
+c =================================================================
+ if((icontri.eq.1).or.(icontri.eq.2)) am_x = am_A
+ if((icontri.eq.1).or.(icontri.eq.3)) am_y = am_B
+
+ q1tx = q1t*cos(phiq1t)
+ q1ty = q1t*sin(phiq1t)
+
+ q2tx = q2t*cos(phiq2t)
+ q2ty = q2t*sin(phiq2t)
+
+ ptsumx = q1tx+q2tx
+ ptsumy = q1ty+q2ty
+
+ ptsum = sqrt(ptsumx**2+ptsumy**2)
+
+c =================================================================
+c a window in final state transverse momentum
+c =================================================================
+
+ if(iptsum.eq.1) then
+ if(ptsum.lt.ptsum_min.or.ptsum.gt.ptsum_max) return
+ endif
+
+c =================================================================
+c compute the individual central particles momentum
+c =================================================================
+
+ ptdiffx = ptdiff*cos(phiptdiff)
+ ptdiffy = ptdiff*sin(phiptdiff)
+
+ pt1x = 0.5*(ptsumx+ptdiffx)
+ pt1y = 0.5*(ptsumy+ptdiffy)
+
+ pt2x = 0.5*(ptsumx-ptdiffx)
+ pt2y = 0.5*(ptsumy-ptdiffy)
+
+ pt1 = sqrt(pt1x**2+pt1y**2)
+ pt2 = sqrt(pt2x**2+pt2y**2)
+
+ if(ipt.eq.1) then
+ if(pt1.lt.pt_min.or.pt2.lt.pt_min) return
+ endif
+
+ amt1 = dsqrt(pt1**2+am_l**2)
+ amt2 = dsqrt(pt2**2+am_l**2)
+
+ invm2 = amt1**2 + amt2**2 + 2.d0*amt1*amt2*dcosh(y1-y2)
+ & -ptsum**2
+
+ invm = dsqrt(invm2)
+
+c =================================================================
+c a window in final state invariant mass
+c =================================================================
+
+ if(iinvm.eq.1) then
+ if(invm.lt.invm_min.or.invm.gt.invm_max) return
+ endif
+
+c =================================================================
+c a window in rapidity distance
+c =================================================================
+
+ dely = dabs(y1-y2)
+ if(idely.eq.1) then
+ if(dely.lt.dely_min.or.dely.gt.dely_max) return
+ endif
+
+c =================================================================
+c auxiliary quantities
+c =================================================================
+
+ alpha1 = amt1/(dsqrt(2.d0)*p1_plus)*dexp( y1)
+ alpha2 = amt2/(dsqrt(2.d0)*p1_plus)*dexp( y2)
+ beta1 = amt1/(dsqrt(2.d0)*p2_minus)*dexp(-y1)
+ beta2 = amt2/(dsqrt(2.d0)*p2_minus)*dexp(-y2)
+
+ q1t2 = q1tx**2 + q1ty**2
+ q2t2 = q2tx**2 + q2ty**2
+
+ x1 = alpha1 + alpha2
+ x2 = beta1 + beta2
+
+ z1p = alpha1/x1
+ z1m = alpha2/x1
+ z2p = beta1/x2
+ z2m = beta2/x2
+
+c -----------------------------------------------------------------
+ if(x1.gt.1.0.or.x2.gt.1.0) return
+c -----------------------------------------------------------------
+
+ s1_eff = x1*s - q1t**2
+ s2_eff = x2*s - q2t**2
+
+c -----------------------------------------------------------------
+c Additional conditions for energy-momentum conservation
+c -----------------------------------------------------------------
+ if(((icontri.eq.2).or.(icontri.eq.4))
+ 1 .and.(dsqrt(s1_eff).le.(am_y+invm))) return
+ if(((icontri.eq.3).or.(icontri.eq.4))
+ 1 .and.(dsqrt(s2_eff).le.(am_x+invm))) return
+
+c =================================================================
+c >>> TO THE OUTPUT COMMON BLOCK
+c =================================================================
+
+c =================================================================
+c four-momenta of the outgoing protons (or remnants)
+c =================================================================
+
+ px_plus = (1.d0-x1) * p1_plus
+ px_minus = ((a_nuc1*am_x)**2 + q1tx**2 + q1ty**2)/2.d0/px_plus
+
+ px(1) = - q1tx
+ px(2) = - q1ty
+ px(3) = (px_plus - px_minus)/dsqrt(2.d0)
+ px(4) = (px_plus + px_minus)/dsqrt(2.d0)
+
+ py_minus = (1.d0-x2) * p2_minus
+ py_plus = ((a_nuc2*am_y)**2 + q2tx**2 + q2ty**2)/2.d0/py_minus
+
+ py(1) = - q2tx
+ py(2) = - q2ty
+ py(3) = (py_plus - py_minus)/dsqrt(2.d0)
+ py(4) = (py_plus + py_minus)/dsqrt(2.d0)
+
+ q1t = dsqrt(q1t2)
+ q2t = dsqrt(q2t2)
+
+c =================================================================
+c four-momenta of the outgoing central particles
+c =================================================================
+
+ nout = 2
+
+ ipdg(1) = pdg_l
+ pc(1,1) = pt1x
+ pc(2,1) = pt1y
+ pc(3,1) = alpha1*ak1z + beta1*ak2z
+ pc(4,1) = alpha1*ak10 + beta1*ak20
+
+ ipdg(2) = -pdg_l
+ pc(1,2) = pt2x
+ pc(2,2) = pt2y
+ pc(3,2) = alpha2*ak1z + beta2*ak2z
+ pc(4,2) = alpha2*ak10 + beta2*ak20
+
+ eta1 = 0.5d0*dlog((dsqrt(amt1**2*(dcosh(y1))**2 - am_l**2) +
+ 2 amt1*dsinh(y1))/(dsqrt(amt1**2*(dcosh(y1))**2 - am_l**2)
+ 3 - amt1*dsinh(y1)))
+
+ eta2 = 0.5d0*dlog((dsqrt(amt2**2*(dcosh(y2))**2 - am_l**2) +
+ 2 amt2*dsinh(y2))/(dsqrt(amt2**2*(dcosh(y2))**2 - am_l**2)
+ 3 - amt2*dsinh(y2)))
+
+ if(ieta.eq.1) then
+ if(eta1.lt.eta_min.or.eta1.gt.eta_max) return
+ if(eta2.lt.eta_min.or.eta2.gt.eta_max) return
+ endif
+
+c =================================================================
+c matrix element squared
+c averaged over initial spin polarizations
+c and summed over final spin polarizations
+c (--> see Wolfgang's notes
+c =================================================================
+
+c =================================================================
+c Mendelstam variables
+c =================================================================
+
+ that1 = (q10-pc(1,4))**2
+ & -(q1tx-pc(1,1))**2-(q1ty-pc(2,1))**2-(q1z-pc(3,1))**2
+ uhat1 = (q10-pc(2,4))**2
+ & -(q1tx-pc(1,2))**2-(q1ty-pc(2,2))**2-(q1z-pc(3,2))**2
+ that2 = (q20-pc(2,4))**2
+ & -(q2tx-pc(1,2))**2-(q2ty-pc(2,2))**2-(q2z-pc(3,2))**2
+ uhat2 = (q20-pc(1,4))**2
+ & -(q2tx-pc(1,1))**2-(q2ty-pc(2,1))**2-(q2z-pc(3,1))**2
+
+ that = (that1+that2)/2.d0
+ uhat = (uhat1+uhat2)/2.d0
+
+c =================================================================
+c matrix elements
+c =================================================================
+c How matrix element is calculated
+c imethod = 0: on-shell formula
+c imethod = 1: off-shell formula
+c =================================================================
+ if(imethod.eq.0) then
+c =================================================================
+c on-shell formula for M^2
+c =================================================================
+ term1 = 6.d0*am_l**8
+ term2 = -3.d0*am_l**4*that**2
+ term3 = -14.d0*am_l**4*that*uhat
+ term4 = -3.d0*am_l**4*uhat**2
+ term5 = am_l**2*that**3
+ term6 = 7.d0*am_l**2*that**2*uhat
+ term7 = 7.d0*am_l**2*that*uhat**2
+ term8 = am_l**2*uhat**3
+ term9 = -that**3*uhat
+ term10 = -that*uhat**3
+
+ amat2 = -2.d0*( term1+term2+term3+term4+term5
+ 2 +term6+term7+term8+term9+term10 )
+ 3 / ( (am_l**2-that)**2 * (am_l**2-uhat)**2 )
+
+ elseif(imethod.eq.1)then
+c =================================================================
+c Wolfgang's formulae
+c =================================================================
+
+ ak1_x = z1m*pt1x-z1p*pt2x
+ ak1_y = z1m*pt1y-z1p*pt2y
+
+ ak2_x = z2m*pt1x-z2p*pt2x
+ ak2_y = z2m*pt1y-z2p*pt2y
+
+ !FIXME FIXME FIXME am_p or am_A/B???
+ t1abs = (q1t2+x1*(am_x**2-am_p**2)+x1**2*am_p**2)/(1.d0-x1)
+ t2abs = (q2t2+x2*(am_y**2-am_p**2)+x2**2*am_p**2)/(1.d0-x2)
+
+ eps12 = am_l**2 + z1p*z1m*t1abs
+ eps22 = am_l**2 + z2p*z2m*t2abs
+
+ Phi10 = 1.d0/((ak1_x+z1p*q2tx)**2+(ak1_y+z1p*q2ty)**2+eps12)
+ 2 - 1.d0/((ak1_x-z1m*q2tx)**2+(ak1_y-z1m*q2ty)**2+eps12)
+ Phi11_x = (ak1_x+z1p*q2tx)/
+ 2 ((ak1_x+z1p*q2tx)**2+(ak1_y+z1p*q2ty)**2+eps12)
+ 3 - (ak1_x-z1m*q2tx)/
+ 4 ((ak1_x-z1m*q2tx)**2+(ak1_y-z1m*q2ty)**2+eps12)
+ Phi11_y = (ak1_y+z1p*q2ty)/
+ 2 ((ak1_x+z1p*q2tx)**2+(ak1_y+z1p*q2ty)**2+eps12)
+ 3 - (ak1_y-z1m*q2ty)/
+ 4 ((ak1_x-z1m*q2tx)**2+(ak1_y-z1m*q2ty)**2+eps12)
+
+ Phi102 = Phi10*Phi10
+ Phi112 = Phi11_x**2+Phi11_y**2
+
+ Phi20 = 1.d0/((ak2_x+z2p*q1tx)**2+(ak2_y+z2p*q1ty)**2+eps22)
+ 2 - 1.d0/((ak2_x-z2m*q1tx)**2+(ak2_y-z2m*q1ty)**2+eps22)
+
+ Phi21_x = (ak2_x+z2p*q1tx)/
+ 2 ((ak2_x+z2p*q1tx)**2+(ak2_y+z2p*q1ty)**2+eps22)
+ 3 - (ak2_x-z2m*q1tx)/
+ 4 ((ak2_x-z2m*q1tx)**2+(ak2_y-z2m*q1ty)**2+eps22)
+ Phi21_y = (ak2_y+z2p*q1ty)/
+ 2 ((ak2_x+z2p*q1tx)**2+(ak2_y+z2p*q1ty)**2+eps22)
+ 3 - (ak2_y-z2m*q1ty)/
+ 4 ((ak2_x-z2m*q1tx)**2+(ak2_y-z2m*q1ty)**2+eps22)
+
+ Phi202 = Phi20*Phi20
+ Phi212 = Phi21_x**2+Phi21_y**2
+
+ Phi11_dot_e = (Phi11_x*q1tx + Phi11_y*q1ty)/dsqrt(q1t2)
+ Phi11_cross_e = (Phi11_x*q1ty - Phi11_y*q1tx)/dsqrt(q1t2)
+
+ Phi21_dot_e = (Phi21_x*q2tx +Phi21_y*q2ty)/dsqrt(q2t2)
+ Phi21_cross_e = (Phi21_x*q2ty -Phi21_y*q2tx)/dsqrt(q2t2)
+
+ aux2_1 = iterm11*(am_l**2+4.d0*z1p**2*z1m**2*t1abs)*Phi102
+ 1 +iterm22*( (z1p**2 + z1m**2)*(Phi11_dot_e**2 +
+ 2 Phi11_cross_e**2) )
+ 3 + itermtt*( Phi11_cross_e**2 - Phi11_dot_e**2)
+ 4 - iterm12*4.d0*z1p*z1m*(z1p-z1m)*Phi10
+ 5 *(q1tx*Phi11_x+q1ty*Phi11_y)
+
+ aux2_2 = iterm11*(am_l**2+4.d0*z2p**2*z2m**2*t2abs)*Phi202
+ 1 +iterm22*( (z2p**2 + z2m**2)*(Phi21_dot_e**2 +
+ 2 Phi21_cross_e**2) )
+ 3 + itermtt*( Phi21_cross_e**2 - Phi21_dot_e**2)
+ 4 - iterm12*4.d0*z2p*z2m*(z2p-z2m)*Phi20
+ 5 *(q2tx*Phi21_x+q2ty*Phi21_y)
+
+c =================================================================
+c convention of matrix element as in our kt-factorization
+c for heavy flavours
+c =================================================================
+
+ amat2_1 = (x1*x2*s)**2 * aux2_1 * 2.*z1p*z1m*q1t2 / (q1t2*q2t2)
+ amat2_2 = (x1*x2*s)**2 * aux2_2 * 2.*z2p*z2m*q2t2 / (q1t2*q2t2)
+
+c =================================================================
+c symmetrization
+c =================================================================
+
+ amat2 = (imat1*amat2_1 + imat2*amat2_2)/2.d0
+
+ endif
+
+ coupling = 1.d0
+c =================================================================
+c first parton coupling
+c =================================================================
+ if(iflux1.ge.20.and.iflux1.lt.40) then ! at least one gluon exchanged
+#ifdef ALPHA_S
+ t_max = max(amt1**2,amt2**2)
+ amu2 = max(eps12,t_max)
+ am_x = dsqrt(amu2)
+ coupling = coupling * 4.d0*pi*alphaS(am_x)/2.d0 ! colour flow
+#else
+ print *,'alphaS not linked to this instance!'
+ stop
+#endif
+ else ! photon exchange
+ coupling = coupling * 4.d0*pi*alpha_em*q_l**2
+ endif
+c =================================================================
+c second parton coupling
+c =================================================================
+ coupling = coupling * 4.d0*pi*alpha_em*q_l**2 ! photon exchange
+ coupling = coupling * 3.d0
+
+c ============================================
+c unintegrated parton distributions
+c ============================================
+
+ if(a_nuc1.le.1) then
+ f1 = CepGen_kT_flux(iflux1,x1,q1t2,sfmod,am_x)
+ else
+ f1 = CepGen_kT_flux_HI(iflux1,x1,q1t2,a_nuc1,z_nuc1)
+ endif
+ if(a_nuc2.le.1) then
+ f2 = CepGen_kT_flux(iflux2,x2,q2t2,sfmod,am_y)
+ else
+ f2 = CepGen_kT_flux_HI(iflux2,x2,q2t2,a_nuc2,z_nuc2)
+ endif
+
+c =================================================================
+c factor 2.*pi below from integration over phi_sum
+c factor 1/4 below from jacobian of transformations
+c factors 1/pi and 1/pi due to integration
+c over d^2 kappa_1 d^2 kappa_2 instead d kappa_1^2 d kappa_2^2
+c =================================================================
+
+ aintegral = (2.d0*pi)*1.d0/(16.d0*pi**2*(x1*x2*s)**2) * amat2
+ & * coupling
+ & * f1/pi * f2/pi * (1.d0/4.d0) * units
+ & * 0.5d0 * 4.0d0 / (4.d0*pi)
+
+c *****************************************************************
+c =================================================================
+ aintegrand = aintegral*q1t*q2t*ptdiff
+c =================================================================
+c *****************************************************************
+c print *,aintegrand,aintegral,coupling
+
+ return
+ end
diff --git a/CepGen/Processes/FortranKTProcess.cpp b/CepGen/Processes/FortranKTProcess.cpp
new file mode 100644
index 0000000..c5c6f01
--- /dev/null
+++ b/CepGen/Processes/FortranKTProcess.cpp
@@ -0,0 +1,187 @@
+#include "CepGen/Processes/FortranKTProcess.h"
+#include "CepGen/Core/ParametersList.h"
+
+#include "CepGen/StructureFunctions/StructureFunctions.h"
+#include "CepGen/Event/Event.h"
+
+#include "CepGen/Physics/KTFlux.h"
+#include "CepGen/Physics/Constants.h"
+#include "CepGen/Physics/PDG.h"
+
+extern "C"
+{
+ struct Constants {
+ double m_p, units, pi, alpha_em;
+ };
+ struct Parameters {
+ int icontri, iflux1, iflux2, imethod, sfmod, pdg_l, a_nuc1, z_nuc1, a_nuc2, z_nuc2;
+ double inp1, inp2;
+ };
+ struct KtKinematics {
+ double q1t, q2t, phiq1t, phiq2t, y1, y2, ptdiff, phiptdiff, m_x, m_y;
+ };
+ struct KinematicsCuts {
+ int ipt, iene, ieta, iinvm, iptsum, idely;
+ double pt_min, pt_max, ene_min, ene_max, eta_min, eta_max;
+ double invm_min, invm_max, ptsum_min, ptsum_max;
+ double dely_min, dely_max;
+ };
+ struct EventKinematics {
+ double px[4], py[4];
+ int nout, idum, pdg[4];
+ double pc[4][4];
+ };
+
+ extern Constants constants_;
+ extern Parameters params_;
+ extern KtKinematics ktkin_;
+ extern KinematicsCuts kincuts_;
+ extern EventKinematics evtkin_;
+}
+
+namespace CepGen
+{
+ namespace Process
+ {
+ FortranKTProcess::FortranKTProcess( const ParametersList& params, const char* name, const char* descr, std::function<void( double& )> func ) :
+ GenericKTProcess( params, name, descr, { { PDG::photon, PDG::photon } }, { PDG::muon, PDG::muon } ),
+ pair_( params.get<int>( "pair", 13 ) ),
+ method_( params.get<int>( "method", 1 ) ),
+ func_( func )
+ {
+ constants_.m_p = ParticleProperties::mass( PDG::proton );
+ constants_.units = Constants::GeV2toBarn;
+ constants_.pi = M_PI;
+ constants_.alpha_em = Constants::alphaEM;
+ }
+
+ void
+ FortranKTProcess::preparePhaseSpace()
+ {
+ mom_ip1_ = event_->getOneByRole( Particle::IncomingBeam1 ).momentum();
+ mom_ip2_ = event_->getOneByRole( Particle::IncomingBeam2 ).momentum();
+
+ registerVariable( y1_, Mapping::linear, cuts_.cuts.central.rapidity_single, { -6., 6. }, "First central particle rapidity" );
+ registerVariable( y2_, Mapping::linear, cuts_.cuts.central.rapidity_single, { -6., 6. }, "Second central particle rapidity" );
+ registerVariable( pt_diff_, Mapping::linear, cuts_.cuts.central.pt_diff, { 0., 50. }, "Transverse momentum difference between central particles" );
+ registerVariable( phi_pt_diff_, Mapping::linear, cuts_.cuts.central.phi_pt_diff, { 0., 2.*M_PI }, "Central particles azimuthal angle difference" );
+
+ //===========================================================================================
+ // feed phase space cuts to the common block
+ //===========================================================================================
+
+ cuts_.cuts.central.pt_single.save( (bool&)kincuts_.ipt, kincuts_.pt_min, kincuts_.pt_max );
+ cuts_.cuts.central.energy_single.save( (bool&)kincuts_.iene, kincuts_.ene_min, kincuts_.ene_max );
+ cuts_.cuts.central.eta_single.save( (bool&)kincuts_.ieta, kincuts_.eta_min, kincuts_.eta_max );
+ cuts_.cuts.central.mass_sum.save( (bool&)kincuts_.iinvm, kincuts_.invm_min, kincuts_.invm_max );
+ cuts_.cuts.central.pt_sum.save( (bool&)kincuts_.iptsum, kincuts_.ptsum_min, kincuts_.ptsum_max );
+ cuts_.cuts.central.rapidity_diff.save( (bool&)kincuts_.idely, kincuts_.dely_min, kincuts_.dely_max );
+
+ //===========================================================================================
+ // feed run parameters to the common block
+ //===========================================================================================
+
+ params_.icontri = (int)cuts_.mode;
+ params_.imethod = method_;
+ params_.sfmod = (int)cuts_.structure_functions->type;
+ params_.pdg_l = pair_;
+
+ //-------------------------------------------------------------------------------------------
+ // incoming beams information
+ //-------------------------------------------------------------------------------------------
+
+ params_.inp1 = cuts_.incoming_beams.first.pz;
+ params_.inp2 = cuts_.incoming_beams.second.pz;
+ const HeavyIon in1 = (HeavyIon)cuts_.incoming_beams.first.pdg;
+ if ( in1 ) {
+ params_.a_nuc1 = in1.A;
+ params_.z_nuc1 = (unsigned short)in1.Z;
+ if ( params_.z_nuc1 > 1 ) {
+ event_->getOneByRole( Particle::IncomingBeam1 ).setPdgId( (PDG)in1 );
+ event_->getOneByRole( Particle::OutgoingBeam1 ).setPdgId( (PDG)in1 );
+ }
+ }
+ else
+ params_.a_nuc1 = params_.z_nuc1 = 1;
+
+ const HeavyIon in2 = (HeavyIon)cuts_.incoming_beams.second.pdg;
+ if ( in2 ) {
+ params_.a_nuc2 = in2.A;
+ params_.z_nuc2 = (unsigned short)in2.Z;
+ if ( params_.z_nuc2 > 1 ) {
+ event_->getOneByRole( Particle::IncomingBeam2 ).setPdgId( (PDG)in2 );
+ event_->getOneByRole( Particle::OutgoingBeam2 ).setPdgId( (PDG)in2 );
+ }
+ }
+ else
+ params_.a_nuc2 = params_.z_nuc2 = 1;
+
+ //-------------------------------------------------------------------------------------------
+ // intermediate partons information
+ //-------------------------------------------------------------------------------------------
+
+ params_.iflux1 = (int)cuts_.incoming_beams.first.kt_flux;
+ params_.iflux2 = (int)cuts_.incoming_beams.second.kt_flux;
+ if ( (KTFlux)params_.iflux1 == KTFlux::P_Gluon_KMR )
+ event_->getOneByRole( Particle::Parton1 ).setPdgId( PDG::gluon );
+ if ( (KTFlux)params_.iflux2 == KTFlux::P_Gluon_KMR )
+ event_->getOneByRole( Particle::Parton2 ).setPdgId( PDG::gluon );
+ }
+
+ double
+ FortranKTProcess::computeKTFactorisedMatrixElement()
+ {
+ ktkin_.q1t = qt1_;
+ ktkin_.q2t = qt2_;
+ ktkin_.phiq1t = phi_qt1_;
+ ktkin_.phiq2t = phi_qt2_;
+ ktkin_.y1 = y1_;
+ ktkin_.y2 = y2_;
+ ktkin_.ptdiff = pt_diff_;
+ ktkin_.phiptdiff = phi_pt_diff_;
+ ktkin_.m_x = MX_;
+ ktkin_.m_y = MY_;
+ double weight = 0.;
+ func_( weight );
+ return weight;
+ }
+
+ void
+ FortranKTProcess::fillCentralParticlesKinematics()
+ {
+ //===========================================================================================
+ // outgoing beam remnants
+ //===========================================================================================
+
+ PX_ = Particle::Momentum( evtkin_.px );
+ PY_ = Particle::Momentum( evtkin_.py );
+ // express these momenta per nucleon
+ PX_ *= 1./params_.a_nuc1;
+ PY_ *= 1./params_.a_nuc2;
+
+//std::cout << PX_ << PY_ << std::endl;
+
+ //===========================================================================================
+ // intermediate partons
+ //===========================================================================================
+
+ const Particle::Momentum mom_par1 = mom_ip1_-PX_, mom_par2 = mom_ip2_-PY_;
+ event_->getOneByRole( Particle::Parton1 ).setMomentum( mom_par1 );
+ event_->getOneByRole( Particle::Parton2 ).setMomentum( mom_par2 );
+ event_->getOneByRole( Particle::Intermediate ).setMomentum( mom_par1+mom_par2 );
+
+ //===========================================================================================
+ // central system
+ //===========================================================================================
+
+ Particles& oc = event_->getByRole( Particle::CentralSystem );
+ for ( int i = 0; i < evtkin_.nout; ++i ) {
+ auto& p = oc[i];
+ p.setPdgId( evtkin_.pdg[i] );
+ p.setStatus( Particle::Status::FinalState );
+ p.setMomentum( Particle::Momentum( evtkin_.pc[i] ) );
+ }
+ }
+ }
+}
+
diff --git a/CepGen/Processes/FortranKTProcess.h b/CepGen/Processes/FortranKTProcess.h
new file mode 100644
index 0000000..af15df0
--- /dev/null
+++ b/CepGen/Processes/FortranKTProcess.h
@@ -0,0 +1,35 @@
+#ifndef CepGen_Processes_FortranKTProcess_h
+#define CepGen_Processes_FortranKTProcess_h
+
+#include "GenericKTProcess.h"
+#include <functional>
+
+namespace CepGen
+{
+ namespace Process
+ {
+ /// Compute the matrix element for a generic \f$k_T\f$-factorised process defined in a Fortran subroutine
+ class FortranKTProcess : public GenericKTProcess
+ {
+ public:
+ FortranKTProcess( const ParametersList& params, const char* name, const char* descr, std::function<void(double&)> func );
+ ProcessPtr clone() const override { return ProcessPtr( new FortranKTProcess( *this ) ); }
+
+ private:
+ void preparePhaseSpace() override;
+ double computeKTFactorisedMatrixElement() override;
+ void fillCentralParticlesKinematics() override;
+
+ int pair_;
+ int method_;
+ /// Subroutine to be called for weight computation
+ std::function<void(double&)> func_;
+ double y1_, y2_, pt_diff_, phi_pt_diff_;
+
+ Particle::Momentum mom_ip1_, mom_ip2_;
+ };
+ }
+}
+
+#endif
+
diff --git a/CepGen/Processes/GenericKTProcess.cpp b/CepGen/Processes/GenericKTProcess.cpp
index b8105a4..df22774 100644
--- a/CepGen/Processes/GenericKTProcess.cpp
+++ b/CepGen/Processes/GenericKTProcess.cpp
@@ -1,356 +1,370 @@
#include "CepGen/Processes/GenericKTProcess.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/ParametersList.h"
#include "CepGen/Event/Event.h"
#include "CepGen/StructureFunctions/StructureFunctions.h"
#include "CepGen/StructureFunctions/SigmaRatio.h"
#include "CepGen/Physics/Constants.h"
#include "CepGen/Physics/KTFlux.h"
#include "CepGen/Physics/FormFactors.h"
#include "CepGen/Physics/PDG.h"
#include <iomanip>
namespace CepGen
{
namespace Process
{
GenericKTProcess::GenericKTProcess( const ParametersList& params,
const std::string& name,
const std::string& description,
const std::array<PDG,2>& partons,
const std::vector<PDG>& central ) :
GenericProcess( name, description+" (kT-factorisation approach)" ),
num_dimensions_( 0 ), kt_jacobian_( 0. ),
qt1_( 0. ), phi_qt1_( 0. ), qt2_( 0. ), phi_qt2_( 0. ),
kIntermediateParts( partons ), kProducedParts( central )
{}
void
GenericKTProcess::addEventContent()
{
GenericProcess::setEventContent(
{ // incoming state
{ Particle::IncomingBeam1, PDG::proton },
{ Particle::IncomingBeam2, PDG::proton },
{ Particle::Parton1, kIntermediateParts[0] },
{ Particle::Parton2, kIntermediateParts[1] }
},
{ // outgoing state
{ Particle::OutgoingBeam1, { PDG::proton } },
{ Particle::OutgoingBeam2, { PDG::proton } },
{ Particle::CentralSystem, kProducedParts }
}
);
setExtraContent();
}
unsigned int
GenericKTProcess::numDimensions() const
{
return num_dimensions_;
}
void
GenericKTProcess::setKinematics( const Kinematics& kin )
{
cuts_ = kin;
const KTFlux flux1 = (KTFlux)cuts_.incoming_beams.first.kt_flux,
flux2 = (KTFlux)cuts_.incoming_beams.second.kt_flux;
if ( cuts_.mode == KinematicsMode::invalid ) {
- bool el1 = ( flux1 == KTFlux::P_Photon_Elastic );
- bool el2 = ( flux2 == KTFlux::P_Photon_Elastic );
+ bool el1 = ( flux1 == KTFlux::P_Photon_Elastic
+ || flux1 == KTFlux::HI_Photon_Elastic
+ || flux1 == KTFlux::P_Gluon_KMR );
+ bool el2 = ( flux2 == KTFlux::P_Photon_Elastic
+ || flux2 == KTFlux::HI_Photon_Elastic
+ || flux2 == KTFlux::P_Gluon_KMR );
if ( el1 && el2 )
cuts_.mode = KinematicsMode::ElasticElastic;
else if ( el1 )
cuts_.mode = KinematicsMode::ElasticInelastic;
else if ( el2 )
cuts_.mode = KinematicsMode::InelasticElastic;
else
cuts_.mode = KinematicsMode::InelasticInelastic;
}
else {
//==========================================================================================
// ensure the first incoming flux is compatible with the kinematics mode
//==========================================================================================
if ( ( cuts_.mode == KinematicsMode::ElasticElastic
|| cuts_.mode == KinematicsMode::ElasticInelastic )
&& ( flux1 != KTFlux::P_Photon_Elastic ) ) {
- cuts_.incoming_beams.first.kt_flux = KTFlux::P_Photon_Elastic;
+ cuts_.incoming_beams.first.kt_flux = ( (HeavyIon)cuts_.incoming_beams.first.pdg )
+ ? KTFlux::HI_Photon_Elastic
+ : KTFlux::P_Photon_Elastic;
CG_DEBUG( "GenericKTProcess:kinematics" )
<< "Set the kt flux for first incoming photon to \""
<< cuts_.incoming_beams.first.kt_flux << "\".";
}
else if ( flux1 != KTFlux::P_Photon_Inelastic
&& flux1 != KTFlux::P_Photon_Inelastic_Budnev ) {
+ if ( (HeavyIon)cuts_.incoming_beams.first.pdg )
+ throw CG_FATAL( "GenericKTProcess:kinematics" )
+ << "Inelastic photon emission from HI not yet supported!";
cuts_.incoming_beams.first.kt_flux = KTFlux::P_Photon_Inelastic_Budnev;
CG_DEBUG( "GenericKTProcess:kinematics" )
<< "Set the kt flux for first incoming photon to \""
<< cuts_.incoming_beams.first.kt_flux << "\".";
}
//==========================================================================================
// ensure the second incoming flux is compatible with the kinematics mode
//==========================================================================================
if ( ( cuts_.mode == KinematicsMode::ElasticElastic
|| cuts_.mode == KinematicsMode::InelasticElastic )
&& ( flux2 != KTFlux::P_Photon_Elastic ) ) {
- cuts_.incoming_beams.second.kt_flux = KTFlux::P_Photon_Elastic;
+ cuts_.incoming_beams.second.kt_flux = ( (HeavyIon)cuts_.incoming_beams.second.pdg )
+ ? KTFlux::HI_Photon_Elastic
+ : KTFlux::P_Photon_Elastic;
CG_DEBUG( "GenericKTProcess:kinematics" )
<< "Set the kt flux for second incoming photon to \""
<< cuts_.incoming_beams.second.kt_flux << "\".";
}
else if ( flux2 != KTFlux::P_Photon_Inelastic
&& flux2 != KTFlux::P_Photon_Inelastic_Budnev ) {
+ if ( (HeavyIon)cuts_.incoming_beams.second.pdg )
+ throw CG_FATAL( "GenericKTProcess:kinematics" )
+ << "Inelastic photon emission from HI not yet supported!";
cuts_.incoming_beams.second.kt_flux = KTFlux::P_Photon_Inelastic_Budnev;
CG_DEBUG( "GenericKTProcess:kinematics" )
<< "Set the kt flux for second incoming photon to \""
<< cuts_.incoming_beams.second.kt_flux << "\".";
}
}
//============================================================================================
// initialise the "constant" (wrt x) part of the Jacobian
//============================================================================================
kt_jacobian_ = 1.;
num_dimensions_ = 0;
mapped_variables_.clear();
//============================================================================================
// register the incoming partons' variables
//============================================================================================
registerVariable( qt1_, Mapping::logarithmic, cuts_.cuts.initial.qt, { 1.e-10, 500. }, "First incoming parton virtuality" );
registerVariable( qt2_, Mapping::logarithmic, cuts_.cuts.initial.qt, { 1.e-10, 500. }, "Second incoming parton virtuality" );
registerVariable( phi_qt1_, Mapping::linear, cuts_.cuts.initial.phi_qt, { 0., 2.*M_PI }, "First incoming parton azimuthal angle" );
registerVariable( phi_qt2_, Mapping::linear, cuts_.cuts.initial.phi_qt, { 0., 2.*M_PI }, "Second incoming parton azimuthal angle" );
//============================================================================================
// register all process-dependent variables
//============================================================================================
preparePhaseSpace();
//============================================================================================
// register the outgoing remnants' variables
//============================================================================================
MX_ = MY_ = event_->getOneByRole( Particle::IncomingBeam1 ).mass();
if ( cuts_.mode == KinematicsMode::InelasticElastic || cuts_.mode == KinematicsMode::InelasticInelastic )
registerVariable( MX_, Mapping::square, cuts_.cuts.remnants.mass_single, { 1.07, 1000. }, "Positive z proton remnant mass" );
if ( cuts_.mode == KinematicsMode::ElasticInelastic || cuts_.mode == KinematicsMode::InelasticInelastic )
registerVariable( MY_, Mapping::square, cuts_.cuts.remnants.mass_single, { 1.07, 1000. }, "Negative z proton remnant mass" );
prepareKinematics();
}
double
GenericKTProcess::computeWeight()
{
if ( mapped_variables_.size() == 0 )
throw CG_FATAL( "GenericKTProcess:weight" )
<< "No variables are mapped with this process!";
if ( kt_jacobian_ == 0. )
throw CG_FATAL( "GenericKTProcess:weight" )
<< "Point-independant component of the Jacobian for this "
<< "kt-factorised process is null.\n\t"
<< "Please check the validity of the phase space!";
//============================================================================================
// generate and initialise all variables, and auxiliary (x-dependent) part of the Jacobian
// for this phase space point.
//============================================================================================
const double aux_jacobian = generateVariables();
if ( aux_jacobian <= 0. )
return 0.;
//============================================================================================
// compute the integrand and combine together into a single weight for the phase space point.
//============================================================================================
const double integrand = computeKTFactorisedMatrixElement();
if ( integrand <= 0. )
return 0.;
const double weight = ( kt_jacobian_*aux_jacobian ) * integrand;
CG_DEBUG_LOOP( "GenericKTProcess:weight" )
<< "Jacobian: " << kt_jacobian_ << " * " << aux_jacobian
<< " = " << ( kt_jacobian_*aux_jacobian ) << ".\n\t"
<< "Integrand = " << integrand << "\n\t"
<< "dW = " << weight << ".";
return weight;
}
std::pair<double,double>
GenericKTProcess::incomingFluxes( double x1, double q1t2, double x2, double q2t2 ) const
{
//--- compute fluxes according to modelling specified in parameters card
std::pair<double,double> fluxes = {
ktFlux( (KTFlux)cuts_.incoming_beams.first.kt_flux, x1, q1t2, *cuts_.structure_functions, MX_ ),
ktFlux( (KTFlux)cuts_.incoming_beams.second.kt_flux, x2, q2t2, *cuts_.structure_functions, MY_ )
};
CG_DEBUG_LOOP( "GenericKTProcess:fluxes" )
<< "KT fluxes: " << fluxes.first << " / " << fluxes.second << ".";
return fluxes;
}
void
GenericKTProcess::registerVariable( double& out, const Mapping& type,
const Limits& in, Limits default_limits,
const char* description )
{
Limits lim = in;
out = 0.; // reset the variable
if ( !in.valid() ) {
CG_DEBUG( "GenericKTProcess:registerVariable" )
<< description << " could not be retrieved from the user configuration!\n\t"
<< "Setting it to the default value: " << default_limits << ".";
lim = default_limits;
}
if ( type == Mapping::logarithmic )
lim = {
std::max( log( lim.min() ), -10. ),
std::min( log( lim.max() ), +10. )
};
mapped_variables_.emplace_back( MappingVariable{ description, lim, out, type, num_dimensions_++ } );
switch ( type ) {
case Mapping::square:
kt_jacobian_ *= 2.*lim.range();
break;
default:
kt_jacobian_ *= lim.range();
break;
}
CG_DEBUG( "GenericKTProcess:registerVariable" )
<< description << " has been mapped to variable " << num_dimensions_ << ".\n\t"
<< "Allowed range for integration: " << lim << ".\n\t"
<< "Variable integration mode: " << type << ".";
}
void
GenericKTProcess::dumpVariables() const
{
std::ostringstream os;
for ( const auto& var : mapped_variables_ )
os << "\n\t(" << var.index << ") " << var.type << " mapping (" << var.description << ") in range " << var.limits;
CG_INFO( "GenericKTProcess:dumpVariables" )
<< "List of variables handled by this kt-factorised process:"
<< os.str();
}
double
GenericKTProcess::generateVariables() const
{
double jacobian = 1.;
for ( const auto& cut : mapped_variables_ ) {
if ( !cut.limits.valid() )
continue;
const double xv = x( cut.index ); // between 0 and 1
switch ( cut.type ) {
case Mapping::linear: {
cut.variable = cut.limits.x( xv );
} break;
case Mapping::logarithmic: {
cut.variable = exp( cut.limits.x( xv ) );
jacobian *= cut.variable;
} break;
case Mapping::square: {
cut.variable = cut.limits.x( xv );
jacobian *= cut.variable;
} break;
}
}
if ( CG_EXCEPT_MATCH( "KtProcess:vars", debugInsideLoop ) ) {
std::ostringstream oss;
for ( const auto& cut : mapped_variables_ ) {
oss << "variable " << cut.index
<< " in range " << std::left << std::setw( 20 ) << cut.limits << std::right
<< " has value " << cut.variable << "\n\t";
}
CG_DEBUG_LOOP( "KtProcess:vars" ) << oss.str();
}
return jacobian;
}
void
GenericKTProcess::fillKinematics( bool )
{
fillCentralParticlesKinematics(); // process-dependent!
fillPrimaryParticlesKinematics();
}
void
GenericKTProcess::fillPrimaryParticlesKinematics()
{
//============================================================================================
// outgoing protons
//============================================================================================
Particle& op1 = event_->getOneByRole( Particle::OutgoingBeam1 ),
&op2 = event_->getOneByRole( Particle::OutgoingBeam2 );
op1.setMomentum( PX_ );
op2.setMomentum( PY_ );
switch ( cuts_.mode ) {
case KinematicsMode::ElasticElastic:
op1.setStatus( Particle::Status::FinalState );
op2.setStatus( Particle::Status::FinalState );
break;
case KinematicsMode::ElasticInelastic:
op1.setStatus( Particle::Status::FinalState );
op2.setStatus( Particle::Status::Unfragmented ); op2.setMass( MY_ );
break;
case KinematicsMode::InelasticElastic:
op1.setStatus( Particle::Status::Unfragmented ); op1.setMass( MX_ );
op2.setStatus( Particle::Status::FinalState );
break;
case KinematicsMode::InelasticInelastic:
op1.setStatus( Particle::Status::Unfragmented ); op1.setMass( MX_ );
op2.setStatus( Particle::Status::Unfragmented ); op2.setMass( MY_ );
break;
default: {
throw CG_FATAL( "GenericKTProcess" )
<< "This kT factorisation process is intended for p-on-p collisions! Aborting.";
} break;
}
//============================================================================================
// incoming partons (photons, pomerons, ...)
//============================================================================================
Particle& g1 = event_->getOneByRole( Particle::Parton1 );
g1.setMomentum( event_->getOneByRole( Particle::IncomingBeam1 ).momentum()-PX_, true );
Particle& g2 = event_->getOneByRole( Particle::Parton2 );
g2.setMomentum( event_->getOneByRole( Particle::IncomingBeam2 ).momentum()-PY_, true );
//============================================================================================
// two-parton system
//============================================================================================
event_->getOneByRole( Particle::Intermediate ).setMomentum( g1.momentum()+g2.momentum() );
}
std::ostream&
operator<<( std::ostream& os, const GenericKTProcess::Mapping& type )
{
switch ( type ) {
case GenericKTProcess::Mapping::linear: return os << "linear";
case GenericKTProcess::Mapping::logarithmic: return os << "logarithmic";
case GenericKTProcess::Mapping::square: return os << "squared";
}
return os;
}
}
}
diff --git a/External/alphaS.f b/External/alphaS.f
new file mode 100755
index 0000000..2128285
--- /dev/null
+++ b/External/alphaS.f
@@ -0,0 +1,703 @@
+C----------------------------------------------------------------------
+C-- Stand-alone code for alpha_s cannibalised (with permission)
+C-- from Andreas Vogt's QCD-PEGASUS package (hep-ph/0408244).
+C-- The running coupling alpha_s is obtained at N^mLO (m = 0,1,2,3)
+C-- by solving the renormalisation group equation in the MSbar scheme
+C-- by a fourth-order Runge-Kutta integration. Transitions from
+C-- n_f to n_f+1 flavours are made when the factorisation scale
+C-- mu_f equals the pole masses m_h (h = c,b,t). At exactly
+C-- the thresholds m_{c,b,t}, the number of flavours n_f = {3,4,5}.
+C-- The top quark mass should be set to be very large to evolve with
+C-- a maximum of five flavours. The factorisation scale mu_f may be
+C-- a constant multiple of the renormalisation scale mu_r. The input
+C-- factorisation scale mu_(f,0) should be less than or equal to
+C-- the charm quark mass. However, if it is greater than the
+C-- charm quark mass, the value of alpha_s at mu_(f,0) = 1 GeV will
+C-- be found using a root-finding algorithm.
+C--
+C-- Example of usage.
+C-- First call the initialisation routine (only needed once):
+C--
+C-- IORD = 2 ! perturbative order (N^mLO,m=0,1,2,3)
+C-- FR2 = 1.D0 ! ratio of mu_f^2 to mu_r^2
+C-- MUR = 1.D0 ! input mu_r in GeV
+C-- ASMUR = 0.5D0 ! input value of alpha_s at mu_r
+C-- MC = 1.4D0 ! charm quark mass
+C-- MB = 4.75D0 ! bottom quark mass
+C-- MT = 1.D10 ! top quark mass
+C-- CALL INITALPHAS(IORD, FR2, MUR, ASMUR, MC, MB, MT)
+C--
+C-- Then get alpha_s at a renormalisation scale mu_r with:
+C--
+C-- MUR = 100.D0 ! renormalisation scale in GeV
+C-- ALFAS = ALPHAS(MUR)
+C--
+C----------------------------------------------------------------------
+C-- Comments to Graeme Watt <watt(at)hep.ucl.ac.uk>
+C----------------------------------------------------------------------
+
+ SUBROUTINE INITALPHAS(IORD, FR2, MUR, ASMUR, MC, MB, MT)
+C-- IORD = 0 (LO), 1 (NLO), 2 (NNLO), 3 (NNNLO).
+C-- FR2 = ratio of mu_f^2 to mu_r^2 (must be a fixed value).
+C-- MUR = input renormalisation scale (in GeV) for alpha_s.
+C-- ASMUR = input value of alpha_s at the renormalisation scale MUR.
+C-- MC,MB,MT = heavy quark masses in GeV.
+ IMPLICIT NONE
+ INTEGER IORD,IORDc,MAXF,MODE
+ DOUBLE PRECISION FR2,MUR,ASMUR,MC,MB,MT,EPS,A,B,DZEROX,
+ & R0c,FR2c,MURc,ASMURc,MCc,MBc,MTc,FINDALPHASR0,R0,ASI
+ COMMON / DZEROXcommon / FR2c,MURc,ASMURc,MCc,MBc,MTc,R0c,IORDc
+ PARAMETER(EPS=1.D-10,MAXF=10000,MODE=1)
+ EXTERNAL FINDALPHASR0
+
+ IF (MUR*sqrt(FR2).LE.MC) THEN ! Check that MUF <= MC.
+ R0 = MUR
+ ASI = ASMUR
+ ELSE ! Solve for alpha_s at R0 = 1 GeV.
+C-- Copy variables to common block.
+ R0c = 1.D0/sqrt(FR2)
+ IORDc = IORD
+ FR2c = FR2
+ MURc = MUR
+ ASMURc = ASMUR
+ MCc = MC
+ MBc = MB
+ MTc = MT
+C-- Now get alpha_s(R0) corresponding to alpha_s(MUR).
+ A = 0.02D0 ! lower bound for alpha_s(R0)
+ B = 2.00D0 ! upper bound for alpha_s(R0)
+ R0 = R0c
+ ASI = DZEROX(A,B,EPS,MAXF,FINDALPHASR0,MODE)
+ END IF
+
+ CALL INITALPHASR0(IORD, FR2, R0, ASI, MC, MB, MT)
+
+ RETURN
+ END
+
+C----------------------------------------------------------------------
+
+C-- Find the zero of this function using DZEROX.
+ DOUBLE PRECISION FUNCTION FINDALPHASR0(ASI)
+ IMPLICIT NONE
+ INTEGER IORD
+ DOUBLE PRECISION FR2, R0, ASI, MC, MB, MT, MUR, ASMUR, ALPHAS
+ COMMON / DZEROXcommon / FR2, MUR, ASMUR, MC, MB, MT, R0, IORD
+
+ CALL INITALPHASR0(IORD, FR2, R0, ASI, MC, MB, MT)
+ FINDALPHASR0 = ALPHAS(MUR) - ASMUR ! solve equal to zero
+
+ RETURN
+ END
+
+C----------------------------------------------------------------------
+
+ SUBROUTINE INITALPHASR0(IORD, FR2, R0, ASI, MC, MB, MT)
+C-- IORD = 0 (LO), 1 (NLO), 2 (NNLO), 3 (NNNLO).
+C-- FR2 = ratio of mu_f^2 to mu_r^2 (must be a fixed value).
+C-- R0 = input renormalisation scale (in GeV) for alphas_s.
+C-- ASI = input value of alpha_s at the renormalisation scale R0.
+C-- MC,MB,MT = heavy quark masses in GeV.
+C-- Must have R0*sqrt(FR2) <= MC to call this subroutine.
+ IMPLICIT NONE
+ INTEGER IORD,NAORD,NASTPS,IVFNS,NFF
+ DOUBLE PRECISION FR2,R0,ASI,MC,MB,MT,LOGFR,R20,
+ & PI,ZETA,CF,CA,TR,AS0,M20,MC2,MB2,MT2
+ PARAMETER(PI = 3.1415 92653 58979 D0)
+
+ COMMON / RZETA / ZETA(6)
+ COMMON / COLOUR / CF, CA, TR
+ COMMON / ASINP / AS0, M20
+ COMMON / ASPAR / NAORD, NASTPS
+ COMMON / VARFLV / IVFNS
+ COMMON / NFFIX / NFF
+ COMMON / FRRAT / LOGFR
+
+c
+c ..QCD colour factors
+c
+ CA = 3.D0
+ CF = 4./3.D0
+ TR = 0.5 D0
+c
+c ..The lowest integer values of the Zeta function
+c
+ ZETA(1) = 0.5772 15664 90153 D0
+ ZETA(2) = 1.64493 40668 48226 D0
+ ZETA(3) = 1.20205 69031 59594 D0
+ ZETA(4) = 1.08232 32337 11138 D0
+ ZETA(5) = 1.03692 77551 43370 D0
+ ZETA(6) = 1.01734 30619 84449 D0
+
+ IVFNS = 1 ! variable flavour-number scheme (VFNS)
+c IVFNS = 0 ! fixed flavour-number scheme (FFNS)
+ NFF = 4 ! number of flavours for FFNS
+ NAORD = IORD ! perturbative order of alpha_s
+ NASTPS = 20 ! num. steps in Runge-Kutta integration
+ R20 = R0**2 ! input renormalisation scale
+ MC2 = MC**2 ! mu_f^2 for charm threshold
+ MB2 = MB**2 ! mu_f^2 for bottom threshold
+ MT2 = MT**2 ! mu_f^2 for top threshold
+ LOGFR = LOG(FR2) ! log of ratio of mu_f^2 to mu_r^2
+ M20 = R20 * FR2 ! input factorisation scale
+
+c
+c ..Stop some nonsense
+c
+ IF ( (IVFNS .EQ. 0) .AND. (NFF .LT. 3) ) THEN
+ WRITE (6,*) 'Wrong flavour number for FFNS evolution. STOP'
+ STOP
+ END IF
+ IF ( (IVFNS .EQ. 0) .AND. (NFF .GT. 5) ) THEN
+ WRITE (6,*) 'Wrong flavour number for FFNS evolution. STOP'
+ STOP
+ END IF
+c
+ IF ( NAORD .GT. 3 ) THEN
+ WRITE (6,*) 'Specified order in a_s too high. STOP'
+ STOP
+ END IF
+c
+ IF ( (IVFNS .NE. 0) .AND. (FR2 .GT. 4.001D0) ) THEN
+ WRITE (6,*) 'Too low mu_r for VFNS evolution. STOP'
+ STOP
+ END IF
+c
+ IF ( (IVFNS .EQ. 1) .AND. (M20 .GT. MC2) ) THEN
+ WRITE (6,*) 'Too high mu_0 for VFNS evolution. STOP'
+ STOP
+ END IF
+c
+ IF ( (ASI .GT. 2.D0) .OR. (ASI .LT. 2.D-2) ) THEN
+ WRITE (6,*) 'alpha_s out of range. STOP'
+ STOP
+ END IF
+c
+ IF ( (IVFNS .EQ. 1) .AND. (MC2 .GT. MB2) ) THEN
+ WRITE (6,*) 'Wrong charm-bottom mass hierarchy. STOP'
+ STOP
+ END IF
+ IF ( (IVFNS .EQ. 1) .AND. (MB2 .GT. MT2) ) THEN
+ WRITE (6,*) 'Wrong bottom-top mass hierarchy. STOP'
+ STOP
+ END IF
+c
+
+C-- Store the beta function coefficients in a COMMON block.
+ CALL BETAFCT
+
+C-- Store a_s = alpha_s(mu_r^2)/(4 pi) at the input scale R0.
+ AS0 = ASI / (4.D0* PI)
+
+C-- Store alpha_s at the heavy flavour thresholds in a COMMON block.
+ IF (IVFNS .NE. 0) THEN
+ CALL EVNFTHR (MC2, MB2, MT2)
+ END IF
+
+ RETURN
+ END
+
+C----------------------------------------------------------------------
+
+ DOUBLE PRECISION FUNCTION ALPHAS(MUR)
+ IMPLICIT NONE
+ INTEGER NFF,IVFNS,NF
+ DOUBLE PRECISION PI,LOGFR,AS0,M20,ASC,M2C,ASB,M2B,AST,M2T,M2,MUR,
+ & R2,ASI,ASF,R20,R2T,R2B,R2C,AS
+ PARAMETER ( PI = 3.1415 92653 58979 D0 )
+c
+c ..Input common blocks
+c
+ COMMON / NFFIX / NFF
+ COMMON / VARFLV / IVFNS
+ COMMON / FRRAT / LOGFR
+ COMMON / ASINP / AS0, M20
+ COMMON / ASFTHR / ASC, M2C, ASB, M2B, AST, M2T
+
+ R2 = MUR**2
+ M2 = R2 * EXP(+LOGFR)
+ IF (IVFNS .EQ. 0) THEN
+c
+c Fixed number of flavours
+c
+ NF = NFF
+ R20 = M20 * R2/M2
+ ASI = AS0
+ ASF = AS (R2, R20, AS0, NF)
+c
+ ELSE
+c
+c ..Variable number of flavours
+c
+ IF (M2 .GT. M2T) THEN
+ NF = 6
+ R2T = M2T * R2/M2
+ ASI = AST
+ ASF = AS (R2, R2T, AST, NF)
+c
+ ELSE IF (M2 .GT. M2B) THEN
+ NF = 5
+ R2B = M2B * R2/M2
+ ASI = ASB
+ ASF = AS (R2, R2B, ASB, NF)
+c
+ ELSE IF (M2 .GT. M2C) THEN
+ NF = 4
+ R2C = M2C * R2/M2
+ ASI = ASC
+ ASF = AS (R2, R2C, ASC, NF)
+c
+ ELSE
+ NF = 3
+ R20 = M20 * R2/M2
+ ASI = AS0
+ ASF = AS (R2, R20, AS0, NF)
+c
+ END IF
+c
+ END IF
+c
+c ..Final value of alpha_s
+c
+ ALPHAS = 4.D0*PI*ASF
+
+ RETURN
+ END
+c
+c =================================================================av==
+
+
+c =====================================================================
+c
+c ..The threshold matching of the QCD coupling in the MS(bar) scheme,
+c a_s = alpha_s(mu_r^2)/(4 pi), for NF -> NF + 1 active flavours
+c up to order a_s^4 (NNNLO).
+c
+c ..The value ASNF of a_s for NF flavours at the matching scale, the
+c logarithm LOGRH = ln (mu_r^2/m_H^2) -- where m_H is the pole mass
+c of the heavy quark -- and NF are passed as arguments to the
+c function ASNF1. The order of the expansion NAORD (defined as
+c the 'n' in N^nLO) is provided by the common-block ASPAR.
+c
+c ..The matching coefficients are inverted from Chetyrkin, Kniehl and
+c Steinhauser, Phys. Rev. Lett. 79 (1997) 2184. The QCD colour
+c factors have been hard-wired in these results. The lowest integer
+c values of the Zeta function are given by the common-block RZETA.
+c
+c =====================================================================
+c
+
+ DOUBLE PRECISION FUNCTION ASNF1 (ASNF, LOGRH, NF)
+c
+ IMPLICIT NONE
+ INTEGER NF, NAORD, NASTPS, PRVCLL, K1, K2
+ DOUBLE PRECISION ASNF,LOGRH,ZETA,CMC,CMCI30,CMCF30,CMCF31,
+ & CMCI31,ASP,LRHP
+
+ DIMENSION CMC(3,0:3)
+c
+c ---------------------------------------------------------------------
+c
+c ..Input common-blocks
+c
+ COMMON / ASPAR / NAORD, NASTPS
+ COMMON / RZETA / ZETA(6)
+c
+c ..Variables to be saved for the next call
+c
+ SAVE CMC, CMCI30, CMCF30, CMCF31, CMCI31, PRVCLL
+c
+c ---------------------------------------------------------------------
+c
+c ..The coupling-constant matching coefficients (CMC's) up to NNNLO
+c (calculated and saved in the first call of this routine)
+c
+ IF (PRVCLL .NE. 1) THEN
+c
+ CMC(1,0) = 0.D0
+ CMC(1,1) = 2./3.D0
+c
+ CMC(2,0) = 14./3.D0
+ CMC(2,1) = 38./3.D0
+ CMC(2,2) = 4./9.D0
+c
+ CMCI30 = + 80507./432.D0 * ZETA(3) + 58933./1944.D0
+ 1 + 128./3.D0 * ZETA(2) * (1.+ DLOG(2.D0)/3.D0)
+ CMCF30 = - 64./9.D0 * (ZETA(2) + 2479./3456.D0)
+ CMCI31 = 8941./27.D0
+ CMCF31 = - 409./27.D0
+ CMC(3,2) = 511./9.D0
+ CMC(3,3) = 8./27.D0
+c
+ PRVCLL = 1
+c
+ END IF
+c
+c ---------------------------------------------------------------------
+c
+c ..The N_f dependent CMC's, and the alpha_s matching at order NAORD
+c
+ CMC(3,0) = CMCI30 + NF * CMCF30
+ CMC(3,1) = CMCI31 + NF * CMCF31
+c
+ ASNF1 = ASNF
+ IF (NAORD .EQ. 0) GO TO 1
+ ASP = ASNF
+c
+ DO 11 K1 = 1, NAORD
+ ASP = ASP * ASNF
+ LRHP = 1.D0
+c
+ DO 12 K2 = 0, K1
+ ASNF1 = ASNF1 + ASP * CMC(K1,K2) * LRHP
+ LRHP = LRHP * LOGRH
+c
+ 12 CONTINUE
+ 11 CONTINUE
+c
+c ---------------------------------------------------------------------
+c
+ 1 RETURN
+ END
+
+c
+c =================================================================av==
+c
+c ..The subroutine EVNFTHR for the evolution of a_s = alpha_s/(4 pi)
+c from a three-flavour initial scale to the four- to six-flavour
+c thresholds (identified with the squares of the corresponding quark
+c masses). The results are written to the common-block ASFTHR.
+c
+c ..The input scale M20 = mu_(f,0)^2 and the corresponding value
+c AS0 of a_s are provided by ASINP. The fixed scale logarithm
+c LOGFR = ln (mu_f^2/mu_r^2) is specified in FRRAT. The alpha_s
+c matching is done by the function ASNF1.
+c
+c =====================================================================
+c
+c
+ SUBROUTINE EVNFTHR (MC2, MB2, MT2)
+c
+ IMPLICIT NONE
+ DOUBLE PRECISION MC2, MB2, MT2, M20, M2C, M2B, M2T, R20, R2C,
+ 1 R2B, R2T, AS, ASNF1, AS0, ASC, ASB, AST,
+ 2 ASC3, ASB4, AST5, LOGFR, SC, SB, ST
+c
+c ---------------------------------------------------------------------
+c
+c ..Input common blocks
+c
+ COMMON / ASINP / AS0, M20
+ COMMON / FRRAT / LOGFR
+c
+c ..Output common blocks
+c
+ COMMON / ASFTHR / ASC, M2C, ASB, M2B, AST, M2T
+
+c ---------------------------------------------------------------------
+c
+c ..Coupling constants at and evolution distances to/between thresholds
+c
+ R20 = M20 * EXP(-LOGFR)
+c
+c ..Charm
+c
+ M2C = MC2
+ R2C = M2C * R20/M20
+ ASC3 = AS (R2C, R20, AS0, 3)
+ SC = LOG (AS0 / ASC3)
+ ASC = ASNF1 (ASC3, -LOGFR, 3)
+c
+c ..Bottom
+c
+ M2B = MB2
+ R2B = M2B * R20/M20
+ ASB4 = AS (R2B, R2C, ASC, 4)
+ SB = LOG (ASC / ASB4)
+ ASB = ASNF1 (ASB4, -LOGFR, 4)
+c
+c ..Top
+c
+ M2T = MT2
+ R2T = M2T * R20/M20
+ AST5 = AS (R2T, R2B, ASB, 5)
+ ST = LOG (ASB / AST5)
+ AST = ASNF1 (AST5, -LOGFR, 5)
+
+ RETURN
+ END
+
+c
+c =================================================================av==
+c
+c ..The running coupling of QCD,
+c
+c AS = a_s = alpha_s(mu_r^2)/(4 pi),
+c
+c obtained by integrating the evolution equation for a fixed number
+c of massless flavours NF. Except at leading order (LO), AS is
+c obtained using a fourth-order Runge-Kutta integration.
+c
+c ..The initial and final scales R20 and R2, the value AS0 at
+c R20, and NF are passed as function arguments. The coefficients
+c of the beta function up to a_s^5 (N^3LO) are provided by the
+c common-block BETACOM. The order of the expansion NAORD (defined
+c as the 'n' in N^nLO) and the number of steps NASTPS for the
+c integration beyond LO are given by the common-block ASPAR.
+c
+c =====================================================================
+c
+c
+ DOUBLE PRECISION FUNCTION AS (R2, R20, AS0, NF)
+c
+ IMPLICIT NONE
+ INTEGER NFMIN, NFMAX, NF, NAORD, NASTPS, K1
+ DOUBLE PRECISION R2, R20, AS0, SXTH, BETA0, BETA1, BETA2, BETA3,
+ & FBETA1,FBETA2,FBETA3,A,LRRAT,DLR,XK0,XK1,XK2,XK3
+ PARAMETER (NFMIN = 3, NFMAX = 6)
+ PARAMETER ( SXTH = 0.16666 66666 66666 D0 )
+c
+c ---------------------------------------------------------------------
+c
+c ..Input common-blocks
+c
+ COMMON / ASPAR / NAORD, NASTPS
+ COMMON / BETACOM / BETA0 (NFMIN:NFMAX), BETA1 (NFMIN:NFMAX),
+ , BETA2 (NFMIN:NFMAX), BETA3 (NFMIN:NFMAX)
+c
+c ..The beta functions FBETAn at N^nLO for n = 1, 2, and 3
+c
+ FBETA1(A) = - A**2 * ( BETA0(NF) + A * BETA1(NF) )
+ FBETA2(A) = - A**2 * ( BETA0(NF) + A * ( BETA1(NF)
+ , + A * BETA2(NF) ) )
+ FBETA3(A) = - A**2 * ( BETA0(NF) + A * ( BETA1(NF)
+ , + A * (BETA2(NF) + A * BETA3(NF)) ) )
+c
+c ---------------------------------------------------------------------
+c
+c ..Initial value, evolution distance and step size
+c
+ AS = AS0
+ LRRAT = LOG (R2/R20)
+ DLR = LRRAT / NASTPS
+c
+c ..Solution of the evolution equation depending on NAORD
+c (fourth-order Runge-Kutta beyond the leading order)
+c
+ IF (NAORD .EQ. 0) THEN
+c
+ AS = AS0 / (1.+ BETA0(NF) * AS0 * LRRAT)
+c
+ ELSE IF (NAORD .EQ. 1) THEN
+c
+ DO 2 K1 = 1, NASTPS
+ XK0 = DLR * FBETA1 (AS)
+ XK1 = DLR * FBETA1 (AS + 0.5 * XK0)
+ XK2 = DLR * FBETA1 (AS + 0.5 * XK1)
+ XK3 = DLR * FBETA1 (AS + XK2)
+ AS = AS + SXTH * (XK0 + 2.* XK1 + 2.* XK2 + XK3)
+ 2 CONTINUE
+c
+ ELSE IF (NAORD .EQ. 2) THEN
+c
+ DO 3 K1 = 1, NASTPS
+ XK0 = DLR * FBETA2 (AS)
+ XK1 = DLR * FBETA2 (AS + 0.5 * XK0)
+ XK2 = DLR * FBETA2 (AS + 0.5 * XK1)
+ XK3 = DLR * FBETA2 (AS + XK2)
+ AS = AS + SXTH * (XK0 + 2.* XK1 + 2.* XK2 + XK3)
+ 3 CONTINUE
+c
+ ELSE IF (NAORD .EQ. 3) THEN
+c
+ DO 4 K1 = 1, NASTPS
+ XK0 = DLR * FBETA3 (AS)
+ XK1 = DLR * FBETA3 (AS + 0.5 * XK0)
+ XK2 = DLR * FBETA3 (AS + 0.5 * XK1)
+ XK3 = DLR * FBETA3 (AS + XK2)
+ AS = AS + SXTH * (XK0 + 2.* XK1 + 2.* XK2 + XK3)
+ 4 CONTINUE
+ END IF
+c
+c ---------------------------------------------------------------------
+c
+ RETURN
+ END
+
+c
+c =================================================================av==
+c
+c ..The subroutine BETAFCT for the coefficients BETA0...BETA3 of the
+c beta function of QCD up to order alpha_s^5 (N^3LO), normalized by
+c
+c
+c d a_s / d ln mu_r^2 = - BETA0 a_s^2 - BETA1 a_s^3 - ...
+c
+c with a_s = alpha_s/(4*pi).
+c
+c ..The MSbar coefficients are written to the common-block BETACOM for
+c
+c NF = 3...6 (parameters NFMIN, NFMAX) quark flavours.
+c
+c ..The factors CF, CA and TF are taken from the common-block COLOUR.
+c Beyond NLO the QCD colour factors are hard-wired in this routine,
+c and the numerical coefficients are truncated to six digits.
+c
+c =====================================================================
+c
+c
+ SUBROUTINE BETAFCT
+c
+ IMPLICIT DOUBLE PRECISION (A - Z)
+ INTEGER NFMIN, NFMAX, NF
+ PARAMETER (NFMIN = 3, NFMAX = 6)
+c
+c ---------------------------------------------------------------------
+c
+c ..Input common-block
+c
+ COMMON / COLOUR / CF, CA, TR
+c
+c ..Output common-block
+c
+ COMMON / BETACOM / BETA0 (NFMIN:NFMAX), BETA1 (NFMIN:NFMAX),
+ 1 BETA2 (NFMIN:NFMAX), BETA3 (NFMIN:NFMAX)
+
+c
+c ---------------------------------------------------------------------
+c
+c ..The full LO and NLO coefficients
+c
+ B00 = 11./3.D0 * CA
+ B01 = -4./3.D0 * TR
+ B10 = 34./3.D0 * CA**2
+ B11 = -20./3.D0 * CA*TR - 4.* CF*TR
+c
+c ..Flavour-number loop and output to the array
+c
+ DO 1 NF = NFMIN, NFMAX
+c
+ BETA0(NF) = B00 + B01 * NF
+ BETA1(NF) = B10 + B11 * NF
+c
+ BETA2(NF) = 1428.50 - 279.611 * NF + 6.01852 * NF**2
+ BETA3(NF) = 29243.0 - 6946.30 * NF + 405.089 * NF**2
+ 1 + 1.49931 * NF**3
+c
+c ---------------------------------------------------------------------
+c
+ 1 CONTINUE
+c
+ RETURN
+ END
+c
+c =================================================================av==
+
+
+C-- G.W. DZEROX taken from CERNLIB to find the zero of a function.
+ DOUBLE PRECISION FUNCTION DZEROX(A0,B0,EPS,MAXF,F,MODE)
+ IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+C Based on
+C
+C J.C.P. Bus and T.J. Dekker, Two Efficient Algorithms with
+C Guaranteed Convergence for Finding a Zero of a Function,
+C ACM Trans. Math. Software 1 (1975) 330-345.
+C
+C (MODE = 1: Algorithm M; MODE = 2: Algorithm R)
+ CHARACTER*80 ERRTXT
+ LOGICAL LMT
+ DIMENSION IM1(2),IM2(2),LMT(2)
+ PARAMETER (Z1 = 1, HALF = Z1/2)
+ DATA IM1 /2,3/, IM2 /-1,3/
+ DZEROX = 0.D0 ! G.W. to prevent compiler warning
+ IF(MODE .NE. 1 .AND. MODE .NE. 2) THEN
+ C=0
+ WRITE(ERRTXT,101) MODE
+ WRITE(6,*) ERRTXT
+ GO TO 99
+ ENDIF
+ FA=F(B0)
+ FB=F(A0)
+ IF(FA*FB .GT. 0) THEN
+ C=0
+ WRITE(ERRTXT,102) A0,B0
+ WRITE(6,*) ERRTXT
+ GO TO 99
+ ENDIF
+ ATL=ABS(EPS)
+ B=A0
+ A=B0
+ LMT(2)=.TRUE.
+ MF=2
+ 1 C=A
+ FC=FA
+ 2 IE=0
+ 3 IF(ABS(FC) .LT. ABS(FB)) THEN
+ IF(C .NE. A) THEN
+ D=A
+ FD=FA
+ END IF
+ A=B
+ B=C
+ C=A
+ FA=FB
+ FB=FC
+ FC=FA
+ END IF
+ TOL=ATL*(1+ABS(C))
+ H=HALF*(C+B)
+ HB=H-B
+ IF(ABS(HB) .GT. TOL) THEN
+ IF(IE .GT. IM1(MODE)) THEN
+ W=HB
+ ELSE
+ TOL=TOL*SIGN(Z1,HB)
+ P=(B-A)*FB
+ LMT(1)=IE .LE. 1
+ IF(LMT(MODE)) THEN
+ Q=FA-FB
+ LMT(2)=.FALSE.
+ ELSE
+ FDB=(FD-FB)/(D-B)
+ FDA=(FD-FA)/(D-A)
+ P=FDA*P
+ Q=FDB*FA-FDA*FB
+ END IF
+ IF(P .LT. 0) THEN
+ P=-P
+ Q=-Q
+ END IF
+ IF(IE .EQ. IM2(MODE)) P=P+P
+ IF(P .EQ. 0 .OR. P .LE. Q*TOL) THEN
+ W=TOL
+ ELSEIF(P .LT. HB*Q) THEN
+ W=P/Q
+ ELSE
+ W=HB
+ END IF
+ END IF
+ D=A
+ A=B
+ FD=FA
+ FA=FB
+ B=B+W
+ MF=MF+1
+ IF(MF .GT. MAXF) THEN
+ WRITE(6,*) "Error in DZEROX: TOO MANY FUNCTION CALLS"
+ GO TO 99
+ ENDIF
+ FB=F(B)
+ IF(FB .EQ. 0 .OR. SIGN(Z1,FC) .EQ. SIGN(Z1,FB)) GO TO 1
+ IF(W .EQ. HB) GO TO 2
+ IE=IE+1
+ GO TO 3
+ END IF
+ DZEROX=C
+ 99 CONTINUE
+ RETURN
+ 101 FORMAT('Error in DZEROX: MODE = ',I3,' ILLEGAL')
+ 102 FORMAT('Error in DZEROX: F(A) AND F(B) HAVE THE SAME SIGN, A = ',
+ 1 1P,D15.8,', B = ',D15.8)
+ END
+
+C ---------------------------------------------------------------------
diff --git a/cmake/UseEnvironment.cmake b/cmake/UseEnvironment.cmake
index fde651c..5fa5738 100644
--- a/cmake/UseEnvironment.cmake
+++ b/cmake/UseEnvironment.cmake
@@ -1,97 +1,104 @@
set(LXPLUS_GCC_ENV "source /afs/cern.ch/sw/lcg/external/gcc/4.9.3/x86_64-slc6-gcc49-opt/setup.sh")
+set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -Wall -cpp")
if(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wno-long-long -pedantic-errors -g")
else()
if(CMAKE_CXX_COMPILER_VERSION VERSION_GREATER 4.7)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wno-long-long -pedantic-errors -g")
else()
message(STATUS "gcc version >= 4.7 is required")
if($ENV{HOSTNAME} MATCHES "^lxplus[0-9]+.cern.ch")
message(STATUS "Compiling on LXPLUS. Did you properly source the environment variables?")
message(STATUS "e.g. `${LXPLUS_GCC_ENV}`")
endif()
message(FATAL_ERROR "Please clean up this build environment and try again...")
endif()
endif()
if(NOT CMAKE_VERSION VERSION_LESS 3.1)
set(CMAKE_CXX_STANDARD 11)
else()
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
endif()
if($ENV{HOSTNAME} MATCHES "^lxplus[0-9]+.cern.ch")
set(BASE_DIR "/afs/cern.ch/sw/lcg/external")
list(APPEND CMAKE_PREFIX_PATH "${BASE_DIR}/CMake/2.8.9/Linux-i386/share/cmake-2.8/Modules")
set(GSL_DIR "${BASE_DIR}/GSL/2.2.1/x86_64-slc6-gcc48-opt")
set(HEPMC_DIR "${BASE_DIR}/HepMC/2.06.08/x86_64-slc6-gcc48-opt")
#set(LHAPDF_DIR "${BASE_DIR}/MCGenerators/lhapdf/5.8.9/x86_64-slc6-gcc46-opt")
set(LHAPDF_DIR "${BASE_DIR}/MCGenerators_lcgcmt67c/lhapdf/6.2.0/x86_64-slc6-gcc48-opt")
set(PYTHIA8_DIR "${BASE_DIR}/MCGenerators_lcgcmt67c/pythia8/230/x86_64-slc6-gcc48-opt")
set(PYTHON_DIR "${BASE_DIR}/Python/2.7.4/x86_64-slc6-gcc48-opt")
set(PYTHON_LIBRARY "${PYTHON_DIR}/lib/libpython2.7.so")
set(PYTHON_EXECUTABLE "${PYTHON_DIR}/bin/python")
set(PYTHON_INCLUDE_DIR "${PYTHON_DIR}/include/python2.7")
set(ROOTSYS "${BASE_DIR}/../app/releases/ROOT/6.06.08/x86_64-slc6-gcc49-opt/root")
message(STATUS "Compiling on LXPLUS. Do not forget to source the environment variables!")
message(STATUS "e.g. `${LXPLUS_GCC_ENV}`")
#--- searching for GSL
find_library(GSL_LIB gsl HINTS ${GSL_DIR} PATH_SUFFIXES lib REQUIRED)
find_library(GSL_CBLAS_LIB gslcblas HINTS ${GSL_DIR} PATH_SUFFIXES lib)
#--- searching for LHAPDF
find_library(LHAPDF LHAPDF HINTS ${LHAPDF_DIR} PATH_SUFFIXES lib)
find_path(LHAPDF_INCLUDE LHAPDF HINTS ${LHAPDF_DIR} PATH_SUFFIXES include)
#--- searching for HepMC
find_library(HEPMC_LIB HepMC HINTS ${HEPMC_DIR} PATH_SUFFIXES lib)
find_path(HEPMC_INCLUDE HepMC HINTS ${HEPMC_DIR} PATH_SUFFIXES include)
else()
find_library(GSL_LIB gsl REQUIRED)
find_library(GSL_CBLAS_LIB gslcblas)
find_library(LHAPDF LHAPDF)
find_path(LHAPDF_INCLUDE LHAPDF)
find_library(HEPMC_LIB HepMC)
find_path(HEPMC_INCLUDE HepMC)
endif()
#--- searching for Pythia 8
set(PYTHIA8_DIRS $ENV{PYTHIA8_DIR} ${PYTHIA8_DIR} /usr /usr/local /opt/pythia8)
find_library(PYTHIA8 pythia8 HINTS ${PYTHIA8_DIRS} PATH_SUFFIXES lib)
find_path(PYTHIA8_INCLUDE Pythia8 HINTS ${PYTHIA8_DIRS} PATH_SUFFIXES include include/Pythia8 include/pythia8)
message(STATUS "GSL found in ${GSL_LIB}")
list(APPEND CEPGEN_EXTERNAL_IO_REQS ${GSL_LIB} ${GSL_CBLAS_LIB})
list(APPEND CEPGEN_EXTERNAL_CORE_REQS ${GSL_LIB} ${GSL_CBLAS_LIB})
list(APPEND CEPGEN_EXTERNAL_STRF_REQS ${GSL_LIB} ${GSL_CBLAS_LIB})
if(LHAPDF)
message(STATUS "LHAPDF found in ${LHAPDF}")
list(APPEND CEPGEN_EXTERNAL_STRF_REQS ${LHAPDF})
add_definitions(-DLIBLHAPDF)
include_directories(${LHAPDF_INCLUDE})
endif()
find_package(PythonLibs 2.7)
if(PYTHONLIBS_FOUND)
list(APPEND CEPGEN_EXTERNAL_CARDS_REQS ${PYTHON_LIBRARIES})
add_definitions(-DPYTHON)
message(STATUS "Python v${PYTHONLIBS_VERSION_STRING} found")
include_directories(${PYTHON_INCLUDE_DIRS})
endif()
if(PYTHIA8)
message(STATUS "Pythia8 found in ${PYTHIA8}")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-misleading-indentation")
list(APPEND CEPGEN_EXTERNAL_HADR_REQS ${PYTHIA8} dl)
add_definitions(-DPYTHIA8)
include_directories(${PYTHIA8_INCLUDE})
endif()
if(HEPMC_LIB)
message(STATUS "HepMC found in ${HEPMC_INCLUDE}")
list(APPEND CEPGEN_EXTERNAL_IO_REQS ${HEPMC_LIB})
add_definitions(-DLIBHEPMC)
include_directories(${HEPMC_INCLUDE})
endif()
find_library(MUPARSER muparser)
if(MUPARSER)
message(STATUS "muParser found in ${MUPARSER}")
list(APPEND CEPGEN_EXTERNAL_CORE_REQS ${MUPARSER})
add_definitions(-DMUPARSER)
endif()
+set(ALPHAS_PATH ${PROJECT_SOURCE_DIR}/External)
+file(GLOB alphas_sources ${ALPHAS_PATH}/alphaS.f)
+if(alphas_sources)
+ message(STATUS "alpha(s) evolution found in ${alphas_sources}")
+ add_definitions(-DALPHA_S)
+endif()