diff --git a/CepGen/Cards/PythonHandler.cpp b/CepGen/Cards/PythonHandler.cpp
index dab6e8d..29d9b9c 100644
--- a/CepGen/Cards/PythonHandler.cpp
+++ b/CepGen/Cards/PythonHandler.cpp
@@ -1,392 +1,393 @@
#include "CepGen/Cards/PythonHandler.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/TamingFunction.h"
#include "CepGen/Core/ParametersList.h"
#include "CepGen/Core/Integrator.h"
#include "CepGen/Processes/ProcessesHandler.h"
#include "CepGen/Hadronisers/HadronisersHandler.h"
#include "CepGen/IO/ExportHandler.h"
#include "CepGen/StructureFunctions/StructureFunctions.h"
#include "CepGen/Physics/GluonGrid.h"
#include "CepGen/Physics/PDG.h"
#include <algorithm>
#if PY_MAJOR_VERSION < 3
# define PYTHON2
#endif
namespace cepgen
{
namespace card
{
//----- specialization for CepGen input cards
PythonHandler::PythonHandler( const char* file )
{
setenv( "PYTHONPATH", ".:Cards:test:../Cards", 1 );
setenv( "PYTHONDONTWRITEBYTECODE", "1", 1 );
CG_DEBUG( "PythonHandler" )
<< "Python PATH: " << getenv( "PYTHONPATH" ) << ".";
std::string filename = pythonPath( 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_DEBUG( "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 ) );
//--- additional particles definition
PyObject* pextp = PyObject_GetAttrString( cfg, PDGLIST_NAME ); // new
if ( pextp ) {
parseExtraParticles( pextp );
Py_CLEAR( pextp );
}
//--- process definition
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 = element( 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 );
params_.setProcess( cepgen::proc::ProcessesHandler::get().build( proc_name, proc_params ) );
//--- process kinematics
PyObject* pin_kinematics = element( process, "inKinematics" ); // borrowed
if ( pin_kinematics )
parseIncomingKinematics( pin_kinematics );
PyObject* pout_kinematics = element( process, "outKinematics" ); // borrowed
if ( pout_kinematics )
parseOutgoingKinematics( pout_kinematics );
//--- taming functions
PyObject* ptam = element( process, "tamingFunctions" ); // borrowed
if ( ptam )
for ( const auto& p : getVector<ParametersList>( ptam ) )
params_.taming_functions->add( p.get<std::string>( "variable" ), p.get<std::string>( "expression" ) );
Py_CLEAR( process );
PyObject* plog = PyObject_GetAttrString( cfg, LOGGER_NAME ); // new
if ( plog ) {
parseLogging( plog );
Py_CLEAR( plog );
}
//--- hadroniser parameters
PyObject* phad = PyObject_GetAttrString( cfg, HADR_NAME ); // new
if ( phad ) {
parseHadroniser( phad );
Py_CLEAR( phad );
}
//--- generation parameters
PyObject* pint = PyObject_GetAttrString( cfg, INTEGRATOR_NAME ); // new
if ( pint ) {
parseIntegrator( pint );
Py_CLEAR( pint );
}
PyObject* pgen = PyObject_GetAttrString( cfg, GENERATOR_NAME ); // new
if ( pgen ) {
parseGenerator( pgen );
Py_CLEAR( pgen );
}
PyObject* pout = PyObject_GetAttrString( cfg, OUTPUT_NAME ); // new
if ( pout ) {
parseOutputModule( pout );
Py_CLEAR( pout );
}
//--- finalisation
Py_CLEAR( cfg );
}
PythonHandler::~PythonHandler()
{
if ( Py_IsInitialized() )
Py_Finalize();
}
void
PythonHandler::parseIncomingKinematics( PyObject* kin )
{
//--- retrieve the beams PDG ids
std::vector<ParametersList> beams_pdg;
fillParameter( kin, "pdgIds", beams_pdg );
if ( !beams_pdg.empty() ) {
if ( beams_pdg.size() != 2 )
throwPythonError( Form( "Invalid list of PDG ids retrieved for incoming beams:\n\t2 PDG ids are expected, %d provided!", beams_pdg.size() ) );
params_.kinematics.incoming_beams. first.pdg = (pdgid_t)beams_pdg.at( 0 ).get<int>( "pdgid" );
params_.kinematics.incoming_beams.second.pdg = (pdgid_t)beams_pdg.at( 1 ).get<int>( "pdgid" );
}
//--- incoming beams kinematics
std::vector<double> beams_pz;
fillParameter( kin, "pz", beams_pz );
if ( !beams_pz.empty() ) {
if ( beams_pz.size() != 2 )
throwPythonError( Form( "Invalid list of pz's retrieved for incoming beams:\n\t2 pz's are expected, %d provided!", beams_pz.size() ) );
params_.kinematics.incoming_beams. first.pz = beams_pz.at( 0 );
params_.kinematics.incoming_beams.second.pz = beams_pz.at( 1 );
}
double sqrt_s = -1.;
fillParameter( kin, "cmEnergy", sqrt_s );
if ( sqrt_s != -1. )
params_.kinematics.setSqrtS( sqrt_s );
//--- structure functions set for incoming beams
PyObject* psf = element( kin, "structureFunctions" ); // borrowed
if ( psf )
params_.kinematics.structure_functions = strfun::StructureFunctionsHandler::get().build( get<ParametersList>( psf ) );
//--- types of parton fluxes for kt-factorisation
std::vector<int> kt_fluxes;
fillParameter( kin, "ktFluxes", kt_fluxes );
if ( !kt_fluxes.empty() ) {
params_.kinematics.incoming_beams.first.kt_flux = (KTFlux)kt_fluxes.at( 0 );
params_.kinematics.incoming_beams.second.kt_flux = ( kt_fluxes.size() > 1 )
? (KTFlux)kt_fluxes.at( 1 )
: (KTFlux)kt_fluxes.at( 0 );
}
//--- specify where to look for the grid path for gluon emission
std::string kmr_grid_path;
fillParameter( kin, "kmrGridPath", kmr_grid_path );
if ( !kmr_grid_path.empty() )
kmr::GluonGrid::get( kmr_grid_path.c_str() );
//--- parse heavy ions beams
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::parseOutgoingKinematics( PyObject* kin )
{
std::vector<int> parts;
fillParameter( kin, "minFinalState", parts );
for ( const auto& pdg : parts )
params_.kinematics.minimum_final_state.emplace_back( (pdgid_t)pdg );
ParametersList part_cuts;
fillParameter( kin, "cuts", part_cuts );
for ( const auto& part : part_cuts.keys() ) {
const auto pdg = (pdgid_t)stoi( part );
const auto& cuts = part_cuts.get<ParametersList>( part );
if ( cuts.has<Limits>( "pt" ) )
params_.kinematics.cuts.central_particles[pdg].pt_single = cuts.get<Limits>( "pt" );
if ( cuts.has<Limits>( "energy" ) )
params_.kinematics.cuts.central_particles[pdg].energy_single = cuts.get<Limits>( "energy" );
if ( cuts.has<Limits>( "eta" ) )
params_.kinematics.cuts.central_particles[pdg].eta_single = cuts.get<Limits>( "eta" );
if ( cuts.has<Limits>( "rapidity" ) )
params_.kinematics.cuts.central_particles[pdg].rapidity_single = cuts.get<Limits>( "rapidity" );
}
// for LPAIR/collinear matrix elements
fillParameter( kin, "q2", params_.kinematics.cuts.initial.q2 );
// for the kT factorised matrix elements
fillParameter( kin, "qt", params_.kinematics.cuts.initial.qt );
fillParameter( kin, "phiqt", params_.kinematics.cuts.initial.phi_qt );
fillParameter( kin, "ptdiff", params_.kinematics.cuts.central.pt_diff );
fillParameter( kin, "phiptdiff", params_.kinematics.cuts.central.phi_pt_diff );
fillParameter( kin, "rapiditydiff", params_.kinematics.cuts.central.rapidity_diff );
// generic phase space limits
fillParameter( kin, "rapidity", params_.kinematics.cuts.central.rapidity_single );
fillParameter( kin, "eta", params_.kinematics.cuts.central.eta_single );
fillParameter( kin, "pt", params_.kinematics.cuts.central.pt_single );
fillParameter( kin, "ptsum", params_.kinematics.cuts.central.pt_sum );
fillParameter( kin, "invmass", params_.kinematics.cuts.central.mass_sum );
fillParameter( kin, "mx", params_.kinematics.cuts.remnants.mass_single );
fillParameter( kin, "yj", params_.kinematics.cuts.remnants.rapidity_single );
Limits lim_xi;
fillParameter( kin, "xi", lim_xi );
if ( lim_xi.valid() )
- params_.kinematics.cuts.remnants.energy_single = ( lim_xi+(-1.) )*( -params_.kinematics.incoming_beams.first.pz );
+ //params_.kinematics.cuts.remnants.energy_single = ( lim_xi+(-1.) )*( -params_.kinematics.incoming_beams.first.pz );
+ params_.kinematics.cuts.remnants.energy_single = -( lim_xi-1. )*params_.kinematics.incoming_beams.first.pz;
}
void
PythonHandler::parseLogging( PyObject* log )
{
int log_level = 0;
fillParameter( log, "level", log_level );
utils::Logger::get().level = (utils::Logger::Level)log_level;
std::vector<std::string> enabled_modules;
fillParameter( log, "enabledModules", enabled_modules );
for ( const auto& mod : enabled_modules )
utils::Logger::get().addExceptionRule( mod );
}
void
PythonHandler::parseIntegrator( PyObject* integr )
{
if ( !PyDict_Check( integr ) )
throwPythonError( "Integrator object should be a dictionary!" );
PyObject* palgo = element( 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_.integration().type = IntegratorType::plain;
else if ( algo == "Vegas" ) {
params_.integration().type = IntegratorType::Vegas;
fillParameter( integr, "alpha", (double&)params_.integration().vegas.alpha );
fillParameter( integr, "iterations", params_.integration().vegas.iterations );
fillParameter( integr, "mode", (int&)params_.integration().vegas.mode );
fillParameter( integr, "verbosity", (int&)params_.integration().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_.integration().vegas.ostream = stderr;
else if ( vegas_logging_output == "cout" )
// redirect all debugging information to the standard stream
params_.integration().vegas.ostream = stdout;
else
params_.integration().vegas.ostream = fopen( vegas_logging_output.c_str(), "w" );
}
else if ( algo == "MISER" ) {
params_.integration().type = IntegratorType::MISER;
fillParameter( integr, "estimateFraction", (double&)params_.integration().miser.estimate_frac );
fillParameter( integr, "minCalls", params_.integration().miser.min_calls );
fillParameter( integr, "minCallsPerBisection", params_.integration().miser.min_calls_per_bisection );
fillParameter( integr, "alpha", (double&)params_.integration().miser.alpha );
fillParameter( integr, "dither", (double&)params_.integration().miser.dither );
}
else
throwPythonError( Form( "Invalid integration() algorithm: %s", algo.c_str() ) );
fillParameter( integr, "numFunctionCalls", params_.integration().ncvg );
fillParameter( integr, "seed", (unsigned long&)params_.integration().rng_seed );
unsigned int rng_engine;
fillParameter( integr, "rngEngine", rng_engine );
switch ( rng_engine ) {
case 0: default: params_.integration().rng_engine = (gsl_rng_type*)gsl_rng_mt19937; break;
case 1: params_.integration().rng_engine = (gsl_rng_type*)gsl_rng_taus2; break;
case 2: params_.integration().rng_engine = (gsl_rng_type*)gsl_rng_gfsr4; break;
case 3: params_.integration().rng_engine = (gsl_rng_type*)gsl_rng_ranlxs0; break;
}
fillParameter( integr, "chiSqCut", params_.integration().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::parseHadroniser( PyObject* hadr )
{
if ( !PyDict_Check( hadr ) )
throwPythonError( "Hadroniser object should be a dictionary!" );
PyObject* pname = element( hadr, MODULE_NAME ); // borrowed
if ( !pname )
throwPythonError( "Hadroniser name is required!" );
std::string hadr_name = get<std::string>( pname );
params_.setHadroniser( cepgen::hadr::HadronisersHandler::get().build( hadr_name, get<ParametersList>( hadr ) ) );
auto h = params_.hadroniser();
h->setParameters( params_ );
{ //--- before calling the init() method
std::vector<std::string> config;
fillParameter( hadr, "preConfiguration", config );
h->readStrings( config );
}
h->init();
{ //--- after init() has been called
std::vector<std::string> config;
fillParameter( hadr, "processConfiguration", config );
for ( const auto& block : config ) {
std::vector<std::string> config_blk;
fillParameter( hadr, block.c_str(), config_blk );
h->readStrings( config_blk );
}
}
}
void
PythonHandler::parseOutputModule( PyObject* pout )
{
if ( !is<ParametersList>( pout ) )
throwPythonError( "Invalid type for output parameters list!" );
PyObject* pname = element( pout, MODULE_NAME ); // borrowed
if ( !pname )
throwPythonError( "Output module name is required!" );
params_.setOutputModule( io::ExportHandler::get().build( get<std::string>( pname ), get<ParametersList>( pout ) ) );
}
void
PythonHandler::parseExtraParticles( PyObject* pparts )
{
if ( !is<ParametersList>( pparts ) )
throwPythonError( "Extra particles definition object should be a parameters list!" );
const auto& parts = get<ParametersList>( pparts );
for ( const auto& k : parts.keys() ) {
const auto& part = parts.get<ParticleProperties>( k );
if ( part.pdgid == 0 || part.mass < 0. )
continue;
CG_DEBUG( "PythonHandler:particles" )
<< "Adding a new particle with name \"" << part.name << "\" to the PDG dictionary.";
PDG::get().define( part );
}
}
}
}
diff --git a/CepGen/Core/Integrand.cpp b/CepGen/Core/Integrand.cpp
index 0fae045..4109c12 100644
--- a/CepGen/Core/Integrand.cpp
+++ b/CepGen/Core/Integrand.cpp
@@ -1,202 +1,202 @@
#include "CepGen/Core/Timer.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/TamingFunction.h"
#include "CepGen/Event/Event.h"
#include "CepGen/Event/Particle.h"
#include "CepGen/Physics/Kinematics.h"
#include "CepGen/Physics/PDG.h"
#include "CepGen/Processes/GenericProcess.h"
#include "CepGen/Hadronisers/GenericHadroniser.h"
#include "CepGen/IO/GenericExportHandler.h"
#include "CepGen/Parameters.h"
#include <cmath>
#include <sstream>
#include <fstream>
namespace cepgen
{
namespace integrand
{
utils::Logger::Level log_level;
utils::Timer tmr;
double
eval( double* x, size_t ndim, void* func_params )
{
log_level = utils::Logger::get().level;
std::shared_ptr<Event> ev;
Parameters* params = nullptr;
if ( !func_params || !( params = static_cast<Parameters*>( func_params ) ) )
throw CG_FATAL( "Integrand" ) << "Failed to retrieve the run parameters!";
proc::GenericProcess* proc = params->process();
if ( !proc )
throw CG_FATAL( "Integrand" ) << "Failed to retrieve the process!";
//================================================================
// start the timer
//================================================================
tmr.reset();
//================================================================
// prepare the event content prior to the process generation
//================================================================
if ( proc->hasEvent() ) // event is not empty
ev = proc->event();
params->prepareRun();
//================================================================
// specify the phase space point to probe
//================================================================
proc->setPoint( ndim, x );
//================================================================
// from this step on, the phase space point is supposed to be set
//================================================================
proc->beforeComputeWeight();
double weight = proc->computeWeight();
//================================================================
// invalidate any unphysical behaviour
//================================================================
if ( weight <= 0. )
return 0.;
//================================================================
// speed up the integration process if no event is to be generated
//================================================================
if ( !ev )
return weight;
if ( !params->storage()
&& !params->taming_functions
&& !params->hadroniser()
&& params->kinematics.cuts.central_particles.empty() )
return weight;
//================================================================
// fill in the process' Event object
//================================================================
proc->fillKinematics();
//================================================================
// once the kinematics variables have been populated, can apply
// the collection of taming functions
//================================================================
if ( params->taming_functions ) {
if ( params->taming_functions->has( "m_central" )
|| params->taming_functions->has( "pt_central" ) ) {
// build the kinematics of the central system
Particle::Momentum central_system;
for ( const auto& part : (*ev)[Particle::CentralSystem] )
central_system += part.momentum();
// tame the cross-section by the reweighting function
if ( params->taming_functions->has( "m_central" ) )
weight *= params->taming_functions->eval( "m_central", central_system.mass() );
if ( params->taming_functions->has( "pt_central" ) )
weight *= params->taming_functions->eval( "pt_central", central_system.pt() );
}
if ( params->taming_functions->has( "q2" ) ) {
weight *= params->taming_functions->eval( "q2", -ev->getOneByRole( Particle::Parton1 ).momentum().mass() );
weight *= params->taming_functions->eval( "q2", -ev->getOneByRole( Particle::Parton2 ).momentum().mass() );
}
}
if ( weight <= 0. )
return 0.;
//================================================================
// set the CepGen part of the event generation
//================================================================
if ( params->storage() )
ev->time_generation = tmr.elapsed();
//================================================================
// event hadronisation and resonances decay
//================================================================
if ( params->hadroniser() ) {
double br = -1.;
if ( !params->hadroniser()->run( *ev, br, params->storage() ) || br == 0. )
return 0.;
weight *= br; // branching fraction for all decays
}
//================================================================
// apply cuts on final state system (after hadronisation!)
// (polish your cuts, as this might be very time-consuming...)
//================================================================
if ( !params->kinematics.cuts.central_particles.empty() ) {
for ( const auto& part : (*ev)[Particle::CentralSystem] ) {
// retrieve all cuts associated to this final state particle
if ( params->kinematics.cuts.central_particles.count( part.pdgId() ) == 0 )
continue;
const auto& cuts_pdgid = params->kinematics.cuts.central_particles.at( part.pdgId() );
// apply these cuts on the given particle
if ( !cuts_pdgid.pt_single.passes( part.momentum().pt() ) )
return 0.;
if ( !cuts_pdgid.energy_single.passes( part.momentum().energy() ) )
return 0.;
if ( !cuts_pdgid.eta_single.passes( part.momentum().eta() ) )
return 0.;
if ( !cuts_pdgid.rapidity_single.passes( part.momentum().rapidity() ) )
return 0.;
}
}
for ( const auto& system : { Particle::OutgoingBeam1, Particle::OutgoingBeam2 } )
for ( const auto& part : (*ev)[system] ) {
if ( part.status() != Particle::Status::FinalState )
continue;
- if ( !params->kinematics.cuts.remnants.energy_single.passes( fabs( part.momentum().energy() ) ) )
+ if ( !params->kinematics.cuts.remnants.energy_single.passes( part.momentum().energy() ) )
return 0.;
if ( !params->kinematics.cuts.remnants.rapidity_single.passes( fabs( part.momentum().rapidity() ) ) )
return 0.;
}
//================================================================
// store the last event in parameters block for a later usage
//================================================================
if ( params->storage() ) {
proc->last_event = ev;
proc->last_event->time_total = tmr.elapsed();
CG_DEBUG( "Integrand" )
<< "[process 0x" << std::hex << proc << std::dec << "] "
<< "Individual time (gen+hadr+cuts): " << proc->last_event->time_total*1.e3 << " ms";
}
//================================================================
// a bit of useful debugging
//================================================================
if ( CG_LOG_MATCH( "Integrand", debugInsideLoop ) ) {
std::ostringstream oss;
for ( unsigned short i = 0; i < ndim; ++i )
oss << Form( "%10.8f ", x[i] );
CG_DEBUG( "Integrand" )
<< "f value for dim-" << ndim << " point ( " << oss.str() << "): "
<< weight;
}
return weight;
}
}
}
diff --git a/CepGen/IO/TextHandler.cpp b/CepGen/IO/TextHandler.cpp
index e4031da..f0a2e6f 100644
--- a/CepGen/IO/TextHandler.cpp
+++ b/CepGen/IO/TextHandler.cpp
@@ -1,330 +1,331 @@
#include "CepGen/IO/ExportHandler.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/ParametersList.h"
#include "CepGen/Core/utils.h"
#include "CepGen/Event/Event.h"
#include "CepGen/Parameters.h"
#include "CepGen/Version.h"
#include <gsl/gsl_histogram.h>
#include <iomanip>
#include <fstream>
#include <regex>
namespace cepgen
{
namespace io
{
/**
* \brief Handler for the generic text file output
* \author Laurent Forthomme <laurent.forthomme@cern.ch>
* \date Jul 2019
*/
class TextHandler : public GenericExportHandler
{
public:
explicit TextHandler( const ParametersList& );
~TextHandler();
void initialise( const Parameters& ) override;
void setCrossSection( double xsec, double ) override { xsec_ = xsec; }
void operator<<( const Event& ) override;
private:
short extractVariableProperties( const std::string& );
std::string textHistogram( const std::string&, const gsl_histogram* ) const;
/// Retrieve a named variable from a particle
double variable( const Particle&, const std::string& ) const;
/// Retrieve a named variable from the whole event
double variable( const Event&, const std::string& ) const;
static const std::regex rgx_select_id_, rgx_select_role_;
static constexpr double INVALID_OUTPUT = -999.;
static constexpr size_t PLOT_WIDTH = 50;
static constexpr char PLOT_CHAR = '#';
std::ofstream file_, hist_file_;
const std::vector<std::string> variables_;
const bool save_banner_, save_variables_;
const bool show_hists_, save_hists_;
const ParametersList hist_variables_;
const std::string separator_;
//--- variables definition
std::unordered_map<short,std::string> variables_name_;
std::unordered_map<short,bool> variable_stored_;
typedef std::pair<unsigned short,std::string> IndexedVariable;
std::unordered_map<short,std::vector<IndexedVariable> > variables_per_id_;
std::unordered_map<Particle::Role,std::vector<IndexedVariable> > variables_per_role_;
std::vector<IndexedVariable> variables_for_event_;
unsigned short num_vars_;
std::ostringstream oss_vars_;
double xsec_;
//--- auxiliary helper maps
const std::unordered_map<std::string,Particle::Role> role_str_ = {
{ "ib1", Particle::Role::IncomingBeam1 }, { "ib2", Particle::Role::IncomingBeam2 },
{ "ob1", Particle::Role::OutgoingBeam1 }, { "ob2", Particle::Role::OutgoingBeam2 },
{ "pa1", Particle::Role::Parton1 }, { "pa2", Particle::Role::Parton2 },
{ "cs", Particle::Role::CentralSystem },
{ "int", Particle::Role::Intermediate }
};
typedef double( Particle::Momentum::*pMethod )(void) const;
/// Mapping of string variables to momentum getter methods
const std::unordered_map<std::string,pMethod> m_mom_str_ = {
{ "px", &Particle::Momentum::px },
{ "py", &Particle::Momentum::py },
{ "pz", &Particle::Momentum::pz },
{ "pt", &Particle::Momentum::pt },
{ "eta", &Particle::Momentum::eta },
{ "phi", &Particle::Momentum::phi },
{ "m", &Particle::Momentum::mass },
{ "e", &Particle::Momentum::energy },
{ "p", &Particle::Momentum::p },
{ "pt2", &Particle::Momentum::pt2 },
{ "th", &Particle::Momentum::theta },
{ "y", &Particle::Momentum::rapidity }
};
//--- kinematic variables
double sqrts_;
unsigned long num_evts_;
struct gsl_histogram_deleter
{
void operator()( gsl_histogram* h ) {
gsl_histogram_free( h );
}
};
std::unordered_map<short,std::unique_ptr<gsl_histogram,gsl_histogram_deleter> > hists_;
};
const std::regex TextHandler::rgx_select_id_( "(\\w+)\\((\\d+)\\)" );
const std::regex TextHandler::rgx_select_role_( "(\\w+)\\(([a-z]+\\d?)\\)" );
TextHandler::TextHandler( const ParametersList& params ) :
GenericExportHandler( "text" ),
file_ ( params.get<std::string>( "filename", "output.txt" ) ),
variables_ ( params.get<std::vector<std::string> >( "variables" ) ),
save_banner_ ( params.get<bool>( "saveBanner", true ) ),
save_variables_( params.get<bool>( "saveVariables", true ) ),
show_hists_ ( params.get<bool>( "showHistograms", true ) ),
save_hists_ ( params.get<bool>( "saveHistograms", false ) ),
hist_variables_( params.get<ParametersList>( "histVariables" ) ),
separator_ ( params.get<std::string>( "separator", "\t" ) ),
num_vars_( 0 ), xsec_( 1. )
{
//--- first extract list of variables to store in output file
oss_vars_.clear();
std::string sep;
for ( const auto& var : variables_ ) {
auto id = extractVariableProperties( var );
if ( id >= 0 ) {
oss_vars_ << sep << var, sep = separator_;
variable_stored_[id] = true;
}
}
//--- then extract list of variables to be plotted in histogram
for ( const auto& var : hist_variables_.keys() ) {
auto id = extractVariableProperties( var );
if ( id < 0 )
continue;
const auto& hvar = hist_variables_.get<ParametersList>( var );
const int nbins = hvar.get<int>( "nbins", 10 );
const double min = hvar.get<double>( "low", 0. ), max = hvar.get<double>( "high", 1. );
hists_[id].reset( gsl_histogram_alloc( nbins ) );
gsl_histogram_set_ranges_uniform( hists_[id].get(), min, max );
CG_INFO( "TextHandler" )
<< "Booking a histogram with " << nbins << " bin" << utils::s( nbins )
<< " between " << min << " and " << max << " for \"" << var << "\".";
}
if ( save_hists_ && !hists_.empty() )
hist_file_.open( "lastrun.hists.txt" );
}
TextHandler::~TextHandler()
{
//--- histograms printout
for ( const auto& var : hist_variables_.keys() ) {
const auto& vn = std::find_if( variables_name_.begin(), variables_name_.end(),
[&var]( auto&& p ) { return p.second == var; } );
if ( vn == variables_name_.end() ) {
CG_WARNING( "TextHandler" )
<< "Failed to retrieve variable \"" << var << "\" for plotting.";
continue;
}
const auto& hist = hists_.at( vn->first ).get();
gsl_histogram_scale( hist, xsec_/( num_evts_+1 ) );
if ( show_hists_ )
CG_INFO( "TextHandler" )
<< textHistogram( var, hist );
if ( save_hists_ )
hist_file_ << "\n" << textHistogram( var, hist ) << "\n";
}
//--- finalisation of the output file
file_.close();
}
void
TextHandler::initialise( const Parameters& params )
{
sqrts_ = params.kinematics.sqrtS();
num_evts_ = 0ul;
if ( save_banner_ )
file_ << banner( params, "#" ) << "\n";
if ( save_variables_ )
file_ << "# " << oss_vars_.str() << "\n";
if ( save_hists_ && !hists_.empty() )
hist_file_ << banner( params, "#" ) << "\n";
}
void
TextHandler::operator<<( const Event& ev )
{
std::vector<double> vars( num_vars_ );
//--- extract and order the variables to be retrieved
//--- particle-level variables (indexed by integer id)
for ( const auto& id_vars : variables_per_id_ ) {
const auto& part = ev[id_vars.first];
//--- loop over the list of variables for this particle
for ( const auto& var : id_vars.second )
vars[var.first] = variable( part, var.second );
}
//--- particle-level variables (indexed by role)
for ( const auto& role_vars : variables_per_role_ ) {
const auto& part = ev[role_vars.first][0];
//--- loop over the list of variables for this particle
for ( const auto& var : role_vars.second )
vars[var.first] = variable( part, var.second );
}
//--- event-level variables
for ( const auto& var : variables_for_event_ )
vars[var.first] = variable( ev, var.second );
//--- write down the variables list in the file
std::string sep;
unsigned short i = 0;
for ( const auto& var : vars ) {
if ( variable_stored_.count( i ) > 0 && variable_stored_.at( i ) )
file_ << sep << var, sep = separator_;
if ( hists_.count( i ) > 0 )
gsl_histogram_increment( hists_.at( i ).get(), var );
++i;
}
file_ << "\n";
++num_evts_;
}
double
TextHandler::variable( const Particle& part, const std::string& var ) const
{
if ( m_mom_str_.count( var ) ) {
auto meth = m_mom_str_.at( var );
return ( part.momentum().*meth )();
}
if ( var == "xi" ) return 1.-part.momentum().energy()*2./sqrts_;
if ( var == "pdg" ) return (double)part.integerPdgId();
if ( var == "charge" ) return part.charge();
if ( var == "status" ) return (double)part.status();
CG_WARNING( "TextHandler" )
<< "Failed to retrieve variable \"" << var << "\".";
return INVALID_OUTPUT;
}
double
TextHandler::variable( const Event& ev, const std::string& var ) const
{
if ( var == "np" )
return (double)ev.size();
if ( var == "nev" )
return (double)num_evts_+1;
if ( var == "nob1" || var == "nob2" ) {
unsigned short out = 0.;
for ( const auto& part : ev[
var == "nob1"
? Particle::Role::OutgoingBeam1
: Particle::Role::OutgoingBeam2
] )
if ( (int)part.status() > 0 )
out++;
return (double)out;
}
if ( var == "tgen" )
return ev.time_generation;
if ( var == "ttot" )
return ev.time_total;
CG_WARNING( "TextHandler" )
<< "Failed to retrieve the event-level variable \"" << var << "\".";
return INVALID_OUTPUT;
}
short
TextHandler::extractVariableProperties( const std::string& var )
{
const auto& vn = std::find_if( variables_name_.begin(), variables_name_.end(),
[&var]( auto&& p ) { return p.second == var; } );
if ( vn != variables_name_.end() )
return vn->first;
std::smatch sm;
if ( std::regex_match( var, sm, rgx_select_id_ ) )
variables_per_id_[std::stod( sm[2].str() )].emplace_back( std::make_pair( num_vars_, sm[1].str() ) );
else if ( std::regex_match( var, sm, rgx_select_role_ ) ) {
const auto& str_role = sm[2].str();
if ( role_str_.count( str_role ) == 0 ) {
CG_WARNING( "TextHandler" )
<< "Invalid particle role retrieved from configuration: \"" << str_role << "\".\n\t"
<< "Skipping the variable \"" << var << "\" in the output module.";
return -1;
}
variables_per_role_[role_str_.at( str_role )].emplace_back( std::make_pair( num_vars_, sm[1].str() ) );
}
else // event-level variables
variables_for_event_.emplace_back( std::make_pair( num_vars_, var ) );
variables_name_[num_vars_] = var;
return num_vars_++;
}
std::string
TextHandler::textHistogram( const std::string& var, const gsl_histogram* hist ) const
{
std::ostringstream os;
const size_t nbins = gsl_histogram_bins( hist );
const double max_bin = gsl_histogram_max_val( hist );
const double inv_max_bin = max_bin > 0. ? 1./max_bin : 0.;
- const std::string sep( 15, ' ' );
+ const std::string sep( 17, ' ' );
os
<< "plot of \"" << var << "\"\n"
- << sep << std::string( PLOT_WIDTH-16-var.size(), ' ' )
+ << sep << std::string( PLOT_WIDTH-15-var.size(), ' ' )
<< "d(sig)/d" << var << " (pb/bin)\n"
<< sep << Form( "%-5.2f", gsl_histogram_min_val( hist ) )
- << std::string( PLOT_WIDTH-9, ' ' )
+ << std::string( PLOT_WIDTH-8, ' ' )
<< Form( "%5.2f", gsl_histogram_max_val( hist ) ) << "\n"
<< sep << std::string( PLOT_WIDTH+2, '.' ); // abscissa axis
for ( size_t i = 0; i < nbins; ++i ) {
double min, max;
gsl_histogram_get_range( hist, i, &min, &max );
const double value = gsl_histogram_get( hist, i );
const int val = value*PLOT_WIDTH*inv_max_bin;
os
- << "\n" << Form( "[%6.2f,%6.2f):", min, max )
+ << "\n" << Form( "[%7.2f,%7.2f):", min, max )
<< std::string( val, PLOT_CHAR ) << std::string( PLOT_WIDTH-val, ' ' )
<< ": " << Form( "%6.2f", value );
}
+ const double bin_width = ( gsl_histogram_max( hist )-gsl_histogram_min( hist ) )/nbins;
os
<< "\n"
<< Form( "%15s", var.c_str() ) << ":" << std::string( PLOT_WIDTH, '.' ) << ":\n" // 2nd abscissa axis
<< "\t("
- << "bin width=" << ( gsl_histogram_max( hist )-gsl_histogram_min( hist ) )/nbins << ", "
+ << "bin width=" << bin_width << " unit" << utils::s( (int)bin_width ) << ", "
<< "mean=" << gsl_histogram_mean( hist ) << ", "
<< "st.dev.=" << gsl_histogram_sigma( hist )
<< ")";
return os.str();
}
}
}
REGISTER_IO_MODULE( text, TextHandler )
diff --git a/CepGen/Physics/Limits.cpp b/CepGen/Physics/Limits.cpp
index 5ea5805..af80887 100644
--- a/CepGen/Physics/Limits.cpp
+++ b/CepGen/Physics/Limits.cpp
@@ -1,132 +1,154 @@
#include "CepGen/Physics/Limits.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/utils.h"
namespace cepgen
{
Limits::Limits( double min, double max ) :
std::pair<double,double>( min, max )
{}
Limits::Limits( const Limits& rhs ) :
std::pair<double,double>( rhs.first, rhs.second )
{}
+ Limits
+ Limits::operator-() const
+ {
+ return Limits( -second, -first );
+ }
+
Limits&
Limits::operator+=( double c )
{
first += c;
second += c;
return *this;
}
Limits&
+ Limits::operator-=( double c )
+ {
+ first -= c;
+ second -= c;
+ return *this;
+ }
+
+ Limits&
Limits::operator*=( double c )
{
first *= c;
second *= c;
if ( c < 0. ) {
const double tmp = first; // will be optimised by compiler anyway...
first = second;
second = tmp;
}
return *this;
}
void
Limits::in( double low, double up )
{
first = low;
second = up;
}
double
Limits::range() const
{
if ( !hasMin() || !hasMax() )
return 0.;
return second-first;
}
bool
Limits::hasMin() const
{
return first != INVALID;
}
bool
Limits::hasMax() const
{
return second != INVALID;
}
bool
Limits::passes( double val ) const
{
if ( hasMin() && val < min() )
return false;
if ( hasMax() && val > max() )
return false;
return true;
}
bool
Limits::valid() const
{
return hasMin() || hasMax();
}
void
Limits::save( bool& on, double& lmin, double& lmax ) const
{
on = false; lmin = lmax = 0.;
if ( !valid() )
return;
on = true;
if ( hasMin() )
lmin = min();
if ( hasMax() )
lmax = max();
if ( lmin == lmax )
on = false;
}
double
Limits::x( double v ) const
{
if ( v < 0. || v > 1. )
throw CG_ERROR( "Limits:shoot" )
<< "x must be comprised between 0 and 1; x value = " << v << ".";
if ( !valid() )
return INVALID;
return first + ( second-first ) * v;
}
std::ostream&
operator<<( std::ostream& os, const Limits& lim )
{
if ( !lim.hasMin() && !lim.hasMax() )
return os << "no cuts";
if ( !lim.hasMin() )
return os << Form( "below %g", lim.max() );
if ( !lim.hasMax() )
return os << Form( "above %g", lim.min() );
return os << Form( "%g to %g", lim.min(), lim.max() );
}
+
Limits
operator+( Limits lim, double c )
{
lim += c;
return lim;
}
Limits
+ operator-( Limits lim, double c )
+ {
+ lim -= c;
+ return lim;
+ }
+
+ Limits
operator*( Limits lim, double c )
{
lim *= c;
return lim;
}
}
diff --git a/CepGen/Physics/Limits.h b/CepGen/Physics/Limits.h
index 7c982b7..1ec246a 100644
--- a/CepGen/Physics/Limits.h
+++ b/CepGen/Physics/Limits.h
@@ -1,56 +1,59 @@
#ifndef CepGen_Physics_Limits_h
#define CepGen_Physics_Limits_h
#include <utility>
#include <ostream>
namespace cepgen
{
/// Validity interval for a variable
class Limits : private std::pair<double,double>
{
public:
/// Define lower and upper limits on a quantity
Limits( double min = INVALID, double max = INVALID );
Limits( const Limits& );
+ Limits operator-() const;
Limits& operator+=( double c );
+ Limits& operator-=( double c );
Limits& operator*=( double c );
friend Limits operator+( Limits lim, double c );
+ friend Limits operator-( Limits lim, double c );
friend Limits operator*( Limits lim, double c );
/// Lower limit to apply on the variable
double min() const { return first; }
/// Lower limit to apply on the variable
double& min() { return first; }
/// Upper limit to apply on the variable
double max() const { return second; }
/// Upper limit to apply on the variable
double& max() { return second; }
/// Export the limits into external variables
void save( bool& on, double& lmin, double& lmax ) const;
/// Find the [0,1] value scaled between minimum and maximum
double x( double v ) const;
/// Specify the lower and upper limits on the variable
void in( double low, double up );
/// Full variable range allowed
double range() const;
/// Have a lower limit?
bool hasMin() const;
/// Have an upper limit?
bool hasMax() const;
/// Check if the value is inside limits' boundaries
bool passes( double val ) const;
/// Is there a lower and upper limit?
bool valid() const;
/// Human-readable expression of the limits
friend std::ostream& operator<<( std::ostream&, const Limits& );
private:
static constexpr double INVALID = -999.999;
};
}
#endif