diff --git a/CepGen/Core/Generator.cpp b/CepGen/Core/Generator.cpp
index 92c9e56..22adf79 100644
--- a/CepGen/Core/Generator.cpp
+++ b/CepGen/Core/Generator.cpp
@@ -1,185 +1,185 @@
#include "CepGen/Generator.h"
#include "CepGen/Parameters.h"
#include "CepGen/Version.h"
#include "CepGen/Core/Integrator.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/Timer.h"
#include "CepGen/Physics/PDG.h"
#include "CepGen/Processes/ProcessesHandler.h"
#include "CepGen/Event/Event.h"
#include <fstream>
#include <chrono>
#include <atomic>
namespace cepgen
{
namespace utils { std::atomic<int> gSignal; }
Generator::Generator() :
- parameters( std::unique_ptr<Parameters>( new Parameters ) ), result_( -1. ), result_error_( -1. )
+ parameters_( new Parameters ), result_( -1. ), result_error_( -1. )
{
CG_DEBUG( "Generator:init" ) << "Generator initialized";
try {
printHeader();
} catch ( const Exception& e ) {
e.dump();
}
// Random number initialization
std::chrono::system_clock::time_point time = std::chrono::system_clock::now();
srandom( time.time_since_epoch().count() );
}
Generator::Generator( Parameters* ip ) :
- parameters( ip ), result_( -1. ), result_error_( -1. )
+ parameters_( ip ), result_( -1. ), result_error_( -1. )
{}
Generator::~Generator()
{
- if ( parameters->generation.enabled
- && parameters->process()
- && parameters->numGeneratedEvents() > 0 )
+ if ( parameters_->generation.enabled
+ && parameters_->process()
+ && parameters_->numGeneratedEvents() > 0 )
CG_INFO( "Generator" )
<< "Mean generation time / event: "
- << parameters->totalGenerationTime()*1.e3/parameters->numGeneratedEvents()
+ << parameters_->totalGenerationTime()*1.e3/parameters_->numGeneratedEvents()
<< " ms.";
}
size_t
Generator::numDimensions() const
{
- if ( !parameters->process() )
+ if ( !parameters_->process() )
return 0;
- return parameters->process()->numDimensions();
+ return parameters_->process()->numDimensions();
}
void
Generator::clearRun()
{
- integrator_.reset();
- if ( parameters->process() ) {
- parameters->process()->first_run = true;
- parameters->process()->addEventContent();
- parameters->process()->setKinematics( parameters->kinematics );
+ if ( parameters_->process() ) {
+ parameters_->process()->first_run = true;
+ parameters_->process()->addEventContent();
+ parameters_->process()->setKinematics( parameters_->kinematics );
}
result_ = result_error_ = -1.;
{
std::ostringstream os;
for ( const auto& pr : cepgen::proc::ProcessesHandler::get().modules() )
os << " " << pr;
CG_DEBUG( "Generator:clearRun" ) << "Processes handled:" << os.str() << ".";
}
}
void
Generator::setParameters( const Parameters& ip )
{
- parameters.reset( new Parameters( (Parameters&)ip ) ); // copy constructor
+ parameters_.reset( new Parameters( (Parameters&)ip ) ); // copy constructor
}
void
Generator::printHeader()
{
std::string tmp;
std::ostringstream os; os << "version " << version() << std::endl;
std::ifstream hf( "README" );
if ( !hf.good() )
throw CG_WARNING( "Generator" ) << "Failed to open README file.";
while ( true ) {
if ( !hf.good() ) break;
getline( hf, tmp );
os << "\n " << tmp;
}
hf.close();
CG_INFO( "Generator" ) << os.str();
}
double
Generator::computePoint( double* x )
{
- double res = integrand::eval( x, numDimensions(), (void*)parameters.get() );
+ double res = integrand::eval( x, numDimensions(), (void*)parameters_.get() );
std::ostringstream os;
for ( unsigned int i = 0; i < numDimensions(); ++i )
os << x[i] << " ";
CG_DEBUG( "Generator:computePoint" )
<< "Result for x[" << numDimensions() << "] = { " << os.str() << "}:\n\t"
<< res << ".";
return res;
}
void
Generator::computeXsection( double& xsec, double& err )
{
CG_INFO( "Generator" ) << "Starting the computation of the process cross-section.";
- clearRun();
integrate();
xsec = result_;
err = result_error_;
if ( xsec < 1.e-2 )
CG_INFO( "Generator" )
<< "Total cross section: " << xsec*1.e3 << " +/- " << err*1.e3 << " fb.";
else if ( xsec < 5.e2 )
CG_INFO( "Generator" )
<< "Total cross section: " << xsec << " +/- " << err << " pb.";
else if ( xsec < 5.e5 )
CG_INFO( "Generator" )
<< "Total cross section: " << xsec*1.e-3 << " +/- " << err*1.e-3 << " nb.";
else if ( xsec < 5.e8 )
CG_INFO( "Generator" )
<< "Total cross section: " << xsec*1.e-6 << " +/- " << err*1.e-6 << " µb.";
else
CG_INFO( "Generator" )
<< "Total cross section: " << xsec*1.e-9 << " +/- " << err*1.e-9 << " mb.";
}
void
Generator::integrate()
{
+ clearRun();
+
// first destroy and recreate the integrator instance
if ( !integrator_ )
- integrator_ = std::unique_ptr<Integrator>( new Integrator( numDimensions(), integrand::eval, parameters.get() ) );
+ integrator_ = std::unique_ptr<Integrator>( new Integrator( numDimensions(), integrand::eval, parameters_.get() ) );
else if ( integrator_->dimensions() != numDimensions() )
- integrator_.reset( new Integrator( numDimensions(), integrand::eval, parameters.get() ) );
+ integrator_.reset( new Integrator( numDimensions(), integrand::eval, parameters_.get() ) );
CG_DEBUG( "Generator:newInstance" )
<< "New integrator instance created\n\t"
- << "Considered topology: " << parameters->kinematics.mode << " case\n\t"
+ << "Considered topology: " << parameters_->kinematics.mode << " case\n\t"
<< "Will proceed with " << numDimensions() << "-dimensional integration.";
const int res = integrator_->integrate( result_, result_error_ );
if ( res != 0 )
throw CG_FATAL( "Generator" )
<< "Error while computing the cross-section!\n\t"
<< "GSL error: " << gsl_strerror( res ) << ".";
}
std::shared_ptr<Event>
Generator::generateOneEvent()
{
integrator_->generateOne();
- parameters->addGenerationTime( parameters->process()->last_event->time_total );
- return parameters->process()->last_event;
+ parameters_->addGenerationTime( parameters_->process()->last_event->time_total );
+ return parameters_->process()->last_event;
}
void
Generator::generate( std::function<void( const Event&, unsigned long )> callback )
{
const utils::Timer tmr;
CG_INFO( "Generator" )
- << parameters->generation.maxgen << " events will be generated.";
+ << parameters_->generation.maxgen << " events will be generated.";
- integrator_->generate( parameters->generation.maxgen, callback );
+ integrator_->generate( parameters_->generation.maxgen, callback );
const double gen_time_s = tmr.elapsed();
CG_INFO( "Generator" )
- << parameters->generation.ngen << " events generated "
+ << parameters_->generation.ngen << " events generated "
<< "in " << gen_time_s << " s "
- << "(" << gen_time_s/parameters->generation.ngen*1.e3 << " ms/event).";
+ << "(" << gen_time_s/parameters_->generation.ngen*1.e3 << " ms/event).";
}
}
diff --git a/CepGen/Core/Parameters.cpp b/CepGen/Core/Parameters.cpp
index f61b999..88358bb 100644
--- a/CepGen/Core/Parameters.cpp
+++ b/CepGen/Core/Parameters.cpp
@@ -1,281 +1,296 @@
#include "CepGen/Parameters.h"
#include "CepGen/Core/Integrator.h"
#include "CepGen/Core/ParametersList.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Core/TamingFunction.h"
#include "CepGen/Physics/PDG.h"
#include "CepGen/Processes/GenericProcess.h"
#include "CepGen/Hadronisers/GenericHadroniser.h"
#include "CepGen/StructureFunctions/StructureFunctions.h"
#include <iomanip>
namespace cepgen
{
Parameters::Parameters() :
general( new ParametersList ),
taming_functions( new utils::TamingFunctionsCollection ),
store_( false ), total_gen_time_( 0. ), num_gen_events_( 0 )
{}
Parameters::Parameters( Parameters& param ) :
general( param.general ),
kinematics( param.kinematics ), integrator( param.integrator ), generation( param.generation ),
taming_functions( param.taming_functions ),
process_( std::move( param.process_ ) ),
hadroniser_( std::move( param.hadroniser_ ) ),
store_( false ), total_gen_time_( param.total_gen_time_ ), num_gen_events_( param.num_gen_events_ )
{}
Parameters::Parameters( const Parameters& param ) :
general( param.general ),
kinematics( param.kinematics ), integrator( param.integrator ), generation( param.generation ),
taming_functions( param.taming_functions ),
store_( false ), total_gen_time_( param.total_gen_time_ ), num_gen_events_( param.num_gen_events_ )
{}
Parameters::~Parameters() // required for unique_ptr initialisation!
{}
+ Parameters&
+ Parameters::operator=( Parameters param )
+ {
+ general = param.general;
+ kinematics = param.kinematics;
+ integrator = param.integrator;
+ generation = param.generation;
+ taming_functions = param.taming_functions;
+ process_ = std::move( param.process_ );
+ hadroniser_ = std::move( param.hadroniser_ );
+ total_gen_time_ = param.total_gen_time_;
+ num_gen_events_ = param.num_gen_events_;
+ return *this;
+ }
+
void
Parameters::setThetaRange( float thetamin, float thetamax )
{
kinematics.cuts.central.eta_single = {
Particle::thetaToEta( thetamax ),
Particle::thetaToEta( thetamin )
};
CG_DEBUG( "Parameters" )
<< "eta in range: " << kinematics.cuts.central.eta_single
<< " => theta(min) = " << thetamin << ", theta(max) = " << thetamax << ".";
}
void
Parameters::prepareRun()
{
//--- first-run preparation
if ( process_ && process_->first_run ) {
CG_DEBUG( "Parameters" )
<< "Run started for " << process_->name() << " process "
<< "0x" << std::hex << process_.get() << std::dec << ".\n\t"
<< "Process mode considered: " << kinematics.mode << "\n\t"
<< " first beam: " << kinematics.incoming_beams.first << "\n\t"
<< " second beam: " << kinematics.incoming_beams.second << "\n\t"
<< " structure functions: " << kinematics.structure_functions;
if ( process_->hasEvent() )
process_->clearEvent();
process_->first_run = false;
}
//--- clear the run statistics
total_gen_time_ = 0.;
num_gen_events_ = 0;
}
void
Parameters::addGenerationTime( double gen_time )
{
total_gen_time_ += gen_time;
num_gen_events_++;
}
proc::GenericProcess*
Parameters::process()
{
return process_.get();
}
std::string
Parameters::processName() const
{
if ( !process_ )
return "no process";
return process_->name();
}
void
Parameters::setProcess( std::unique_ptr<proc::GenericProcess> proc )
{
process_ = std::move( proc );
}
void
Parameters::setProcess( proc::GenericProcess* proc )
{
if ( !proc )
throw CG_FATAL( "Parameters" )
<< "Trying to clone an invalid process!";
process_.reset( proc );
}
hadr::GenericHadroniser*
Parameters::hadroniser()
{
return hadroniser_.get();
}
std::string
Parameters::hadroniserName() const
{
if ( !hadroniser_ )
return "";
return hadroniser_->name();
}
void
Parameters::setHadroniser( std::unique_ptr<hadr::GenericHadroniser> hadr )
{
hadroniser_ = std::move( hadr );
}
void
Parameters::setHadroniser( hadr::GenericHadroniser* hadr )
{
hadroniser_.reset( hadr );
}
std::ostream&
operator<<( std::ostream& os, const Parameters* p )
{
const bool pretty = true;
const int wb = 90, wt = 33;
os << std::left << "\n";
if ( p->process_ ) {
os
<< std::setfill('_') << std::setw( wb+3 ) << "_/¯ PROCESS INFORMATION ¯\\_" << std::setfill( ' ' ) << "\n"
<< std::right << std::setw( wb ) << std::left << std::endl
<< std::setw( wt ) << "Process to generate"
<< ( pretty ? boldify( p->process_->name().c_str() ) : p->process_->name() ) << "\n"
<< std::setw( wt ) << "" << p->process_->description() << "\n";
for ( const auto& par : p->process_->parameters().keys() )
os << " " << par << ": " << p->process_->parameters().getString( par ) << "\n";
}
os
<< "\n"
<< std::setfill('_') << std::setw( wb+3 ) << "_/¯ RUN INFORMATION ¯\\_" << std::setfill( ' ' ) << "\n"
<< std::right << std::setw( wb ) << std::left << std::endl
<< std::setw( wt ) << "Events generation? "
<< ( pretty ? yesno( p->generation.enabled ) : std::to_string( p->generation.enabled ) ) << "\n"
<< std::setw( wt ) << "Number of events to generate"
<< ( pretty ? boldify( p->generation.maxgen ) : std::to_string( p->generation.maxgen ) ) << "\n";
if ( p->generation.num_threads > 1 )
os
<< std::setw( wt ) << "Number of threads" << p->generation.num_threads << "\n";
os
<< std::setw( wt ) << "Number of points to try per bin" << p->generation.num_points << "\n"
<< std::setw( wt ) << "Integrand treatment"
<< ( pretty ? yesno( p->generation.treat ) : std::to_string( p->generation.treat ) ) << "\n"
<< std::setw( wt ) << "Verbosity level " << utils::Logger::get().level << "\n";
if ( p->hadroniser_ ) {
os
<< "\n"
<< std::setfill( '-' ) << std::setw( wb+6 )
<< ( pretty ? boldify( " Hadronisation algorithm " ) : "Hadronisation algorithm" ) << std::setfill( ' ' ) << "\n\n"
<< std::setw( wt ) << "Name"
<< ( pretty ? boldify( p->hadroniser_->name().c_str() ) : p->hadroniser_->name() ) << "\n";
}
os
<< "\n"
<< std::setfill( '-' ) << std::setw( wb+6 )
<< ( pretty ? boldify( " Integration parameters " ) : "Integration parameters" ) << std::setfill( ' ' ) << "\n\n";
std::ostringstream int_algo; int_algo << p->integrator.type;
os
<< std::setw( wt ) << "Integration algorithm"
<< ( pretty ? boldify( int_algo.str().c_str() ) : int_algo.str() ) << "\n"
<< std::setw( wt ) << "Number of function calls" << p->integrator.ncvg << "\n"
<< std::setw( wt ) << "Random number generator seed" << p->integrator.rng_seed << "\n";
if ( p->integrator.rng_engine )
os
<< std::setw( wt ) << "Random number generator engine"
<< p->integrator.rng_engine->name << "\n";
std::ostringstream proc_mode; proc_mode << p->kinematics.mode;
os
<< "\n"
<< std::setfill('_') << std::setw( wb+3 ) << "_/¯ EVENTS KINEMATICS ¯\\_" << std::setfill( ' ' ) << "\n\n"
<< std::setw( wt ) << "Incoming particles"
<< p->kinematics.incoming_beams.first << ",\n" << std::setw( wt ) << ""
<< p->kinematics.incoming_beams.second << "\n"
<< std::setw( wt ) << "C.m. energy (GeV)" << p->kinematics.sqrtS() << "\n";
if ( p->kinematics.mode != KinematicsMode::invalid )
os << std::setw( wt ) << "Subprocess mode" << ( pretty ? boldify( proc_mode.str().c_str() ) : proc_mode.str() ) << "\n";
if ( p->kinematics.mode != KinematicsMode::ElasticElastic )
os << std::setw( wt ) << "Structure functions" << *p->kinematics.structure_functions << "\n";
os
<< "\n"
<< std::setfill( '-' ) << std::setw( wb+6 ) << ( pretty ? boldify( " Incoming partons " ) : "Incoming partons" ) << std::setfill( ' ' ) << "\n\n";
for ( const auto& lim : p->kinematics.cuts.initial.list() ) // map(particles class, limits)
if ( lim.second.valid() )
os << std::setw( wt ) << lim.first << lim.second << "\n";
os
<< "\n"
<< std::setfill( '-' ) << std::setw( wb+6 ) << ( pretty ? boldify( " Outgoing central system " ) : "Outgoing central system" ) << std::setfill( ' ' ) << "\n\n";
for ( const auto& lim : p->kinematics.cuts.central.list() )
if ( lim.second.valid() )
os << std::setw( wt ) << lim.first << lim.second << "\n";
if ( p->kinematics.cuts.central_particles.size() > 0 ) {
os << std::setw( wt ) << ( pretty ? boldify( ">>> per-particle cuts:" ) : ">>> per-particle cuts:" ) << "\n";
for ( const auto& part_per_lim : p->kinematics.cuts.central_particles ) {
os << " * all single " << std::setw( wt-3 ) << part_per_lim.first << "\n";
for ( const auto& lim : part_per_lim.second.list() )
if ( lim.second.valid() )
os << " - " << std::setw( wt-5 ) << lim.first << lim.second << "\n";
}
}
os << "\n";
os << std::setfill( '-' ) << std::setw( wb+6 ) << ( pretty ? boldify( " Proton / remnants " ) : "Proton / remnants" ) << std::setfill( ' ' ) << "\n\n";
for ( const auto& lim : p->kinematics.cuts.remnants.list() )
os << std::setw( wt ) << lim.first << lim.second << "\n";
return os;
}
std::ostream&
operator<<( std::ostream& os, const Parameters& p )
{
return os << &p;
}
//-----------------------------------------------------------------------------------------------
Parameters::Integration::Integration() :
type( IntegratorType::Vegas ), ncvg( 50000 ),
rng_seed( 0 ), rng_engine( (gsl_rng_type*)gsl_rng_mt19937 ),
vegas_chisq_cut( 1.5 ),
result( -1. ), err_result( -1. )
{
const size_t ndof = 10;
{
std::shared_ptr<gsl_monte_vegas_state> tmp_state( gsl_monte_vegas_alloc( ndof ), gsl_monte_vegas_free );
gsl_monte_vegas_params_get( tmp_state.get(), &vegas );
vegas.iterations = 10;
}
{
std::shared_ptr<gsl_monte_miser_state> tmp_state( gsl_monte_miser_alloc( ndof ), gsl_monte_miser_free );
gsl_monte_miser_params_get( tmp_state.get(), &miser );
}
}
Parameters::Integration::Integration( const Integration& rhs ) :
type( rhs.type ), ncvg( rhs.ncvg ),
rng_seed( rhs.rng_seed ), rng_engine( rhs.rng_engine ),
vegas( rhs.vegas ), vegas_chisq_cut( rhs.vegas_chisq_cut ),
miser( rhs.miser ),
result( -1. ), err_result( -1. )
{}
Parameters::Integration::~Integration()
{
//if ( vegas.ostream && vegas.ostream != stdout && vegas.ostream != stderr )
// fclose( vegas.ostream );
}
//-----------------------------------------------------------------------------------------------
Parameters::Generation::Generation() :
enabled( false ), maxgen( 0 ),
symmetrise( false ), treat( true ), ngen( 0 ), gen_print_every( 10000 ),
num_threads( 2 ), num_points( 100 )
{}
}
diff --git a/CepGen/Generator.h b/CepGen/Generator.h
index a383a58..5118332 100644
--- a/CepGen/Generator.h
+++ b/CepGen/Generator.h
@@ -1,122 +1,126 @@
#ifndef CepGen_Generator_h
#define CepGen_Generator_h
#include <sstream>
#include <memory>
#include <functional>
////////////////////////////////////////////////////////////////////////////////
/**
* \mainpage Foreword
* This Monte Carlo generator was developed as a modern version of the LPAIR code introduced
* in the early 1990s by J. Vermaseren *et al*\cite Baranov:1991yq\cite Vermaseren:1982cz. This latter allows to
* compute the cross-section and to generate events for the \f$\gamma\gamma\to\ell^{+}\ell^{-}\f$
* process in the scope of high energy physics.
*
* Soon after the integration of its matrix element, it was extended as a tool to compute and
* generate events for any generic 2\f$\rightarrow\f$ 3 central exclusive process.
* To do so, the main operation performed here is the integration of the matrix element (given as a
* subset of a GenericProcess object) over the full available phase space.
*
*/
////////////////////////////////////////////////////////////////////////////////
/// Common namespace for this Monte Carlo generator
namespace cepgen
{
namespace integrand
{
/**
* Function to be integrated. It returns the value of the weight for one point
* of the full phase space (or "event"). This weights includes the matrix element
* of the process considered, along with all the kinematic factors, and the cut
* restrictions imposed on this phase space. \f$x\f$ is therefore an array of random
* numbers defined inside its boundaries (as normalised so that \f$\forall i<\mathrm{ndim}\f$,
* \f$0<x_i<1\f$.
*/
double eval( double*, size_t, void* );
}
class Event;
class Integrator;
class Parameters;
////////////////////////////////////////////////////////////////////////////////
/**
* This object represents the core of this Monte Carlo generator, with its
* capability to generate the events (using the embedded Vegas object) and to
* study the phase space in term of the variation of resulting cross section
* while scanning the various parameters (point \f${\bf x}\f$ in the
* multi-dimensional phase space).
*
* The phase space is constrained using the Parameters object given as an
* argument to the constructor, and the differential cross-sections for each
* value of the array \f${\bf x}\f$ are computed in the \a f-function defined
* outside (but populated inside) this object.
*
* This f-function embeds a GenericProcess-inherited object which defines all the
* methods to compute this differential cross-section as well as the in- and outgoing
* kinematics associated to each particle.
*
* \author Laurent Forthomme <laurent.forthomme@cern.ch>
* \date Feb 2013
* \brief Core of the Monte-Carlo generator
*
*/
class Generator {
public:
/// Core of the Monte Carlo integrator and events generator
Generator();
/// Core of the Monte Carlo integrator and events generator
/// \param[in] ip List of input parameters defining the phase space on which to perform the integration
Generator( Parameters *ip );
~Generator();
+
/// Dump this program's header into the standard output stream
void printHeader();
+
+ const Parameters& parameters() const { return *parameters_; }
/// Feed the generator with a Parameters object
void setParameters( const Parameters& ip );
/// Remove all references to a previous generation/run
void clearRun();
/**
* Compute the cross section for the run parameters defined by this object.
* This returns the cross section as well as the absolute error computed along.
* \brief Compute the cross-section for the given process
* \param[out] xsec The computed cross-section, in pb
* \param[out] err The absolute integration error on the computed cross-section, in pb
*/
void computeXsection( double& xsec, double& err );
- /// Integrate the functional over the whole phase space
- void integrate();
/// Last cross section computed by the generator
double crossSection() const { return result_; }
/// Last error on the cross section computed by the generator
double crossSectionError() const { return result_error_; }
//void terminate();
/// Generate one single event given the phase space computed by Vegas in the integration step
/// \return A pointer to the Event object generated in this run
std::shared_ptr<Event> generateOneEvent();
/// Launch the generation of events
void generate( std::function<void( const Event&, unsigned long )> callback = nullptr );
/// Number of dimensions on which the integration is performed
size_t numDimensions() const;
/// Compute one single point from the total phase space
/// \param[in] x the n-dimensional point to compute
/// \return the function value for the given point
double computePoint( double* x );
- /// Physical Parameters used in the events generation and cross-section computation
- std::unique_ptr<Parameters> parameters;
+
private:
+ /// Integrate the functional over the whole phase space
+ void integrate();
+ /// Physical Parameters used in the events generation and cross-section computation
+ std::unique_ptr<Parameters> parameters_;
/// Vegas instance which will integrate the function
std::unique_ptr<Integrator> integrator_;
/// Cross section value computed at the last integration
double result_;
/// Error on the cross section as computed in the last integration
double result_error_;
};
}
#endif
diff --git a/CepGen/Parameters.h b/CepGen/Parameters.h
index 2fe6e5e..0260b0d 100644
--- a/CepGen/Parameters.h
+++ b/CepGen/Parameters.h
@@ -1,149 +1,153 @@
#ifndef CepGen_Parameters_h
#define CepGen_Parameters_h
#include "CepGen/Physics/Kinematics.h"
#include <memory>
#include <gsl/gsl_monte_vegas.h>
#include <gsl/gsl_monte_miser.h>
namespace cepgen
{
class Event;
class ParametersList;
namespace proc { class GenericProcess; }
namespace hadr { class GenericHadroniser; }
namespace utils { class TamingFunctionsCollection; }
enum class IntegratorType;
/// List of parameters used to start and run the simulation job
class Parameters
{
public:
Parameters();
/// Copy constructor (transfers ownership to the process/hadroniser!)
Parameters( Parameters& );
/// Const copy constructor (all but the process and the hadroniser)
Parameters( const Parameters& );
~Parameters(); // required for unique_ptr initialisation!
+
+ /// Assignment operator
+ Parameters& operator=( Parameters );
+
/// Set the polar angle range for the produced leptons
/// \param[in] thetamin The minimal value of \f$\theta\f$ for the outgoing leptons
/// \param[in] thetamax The maximal value of \f$\theta\f$ for the outgoing leptons
void setThetaRange( float thetamin, float thetamax );
/// Dump the input parameters in the terminal
friend std::ostream& operator<<( std::ostream& os, const Parameters* );
friend std::ostream& operator<<( std::ostream& os, const Parameters& );
std::shared_ptr<ParametersList> general;
//----- process to compute
/// Process for which the cross-section will be computed and the events will be generated
proc::GenericProcess* process();
/// Name of the process considered
std::string processName() const;
/// Set the process to study
void setProcess( std::unique_ptr<proc::GenericProcess> proc );
/// Set the process to study
void setProcess( proc::GenericProcess* proc );
//----- events kinematics
/// Events kinematics for phase space definition
Kinematics kinematics;
//----- VEGAS
/// Collection of integrator parameters
struct Integration
{
Integration();
Integration( const Integration& );
~Integration();
IntegratorType type;
/// Number of function calls to be computed for each point
unsigned int ncvg; // ??
/// Random number generator seed
long rng_seed;
/// Random number generator engine
gsl_rng_type* rng_engine;
gsl_monte_vegas_params vegas;
double vegas_chisq_cut;
gsl_monte_miser_params miser;
double result, err_result;
};
/// Integrator parameters
Integration integrator;
//----- events generation
/// Collection of events generation parameters
struct Generation
{
Generation();
/// Are we generating events ? (true) or are we only computing the cross-section ? (false)
bool enabled;
/// Maximal number of events to generate in this run
unsigned int maxgen;
/// Do we want the events to be symmetrised with respect to the \f$z\f$-axis ?
bool symmetrise;
/// Is the integrand to be smoothed for events generation?
bool treat;
/// Number of events already generated in this run
unsigned int ngen;
/// Frequency at which the events are displayed to the end-user
unsigned int gen_print_every;
/// Number of threads to perform the integration
unsigned int num_threads;
/// Number of points to "shoot" in each integration bin by the algorithm
unsigned int num_points;
};
/// Events generation parameters
Generation generation;
/// Specify if the generated events are to be stored
void setStorage( bool store ) { store_ = store; }
/// Are the events generated in this run to be stored in the output file ?
bool storage() const { return store_; }
//----- hadronisation algorithm
/// Hadronisation algorithm to use for the proton(s) fragmentation
hadr::GenericHadroniser* hadroniser();
/// Name of the hadroniser (if applicable)
std::string hadroniserName() const;
/// Set the hadronisation algorithm
void setHadroniser( std::unique_ptr<hadr::GenericHadroniser> hadr );
/// Set the hadronisation algorithm
void setHadroniser( hadr::GenericHadroniser* hadr );
//----- taming functions
/// Functionals to be used to account for rescattering corrections (implemented within the process)
std::shared_ptr<utils::TamingFunctionsCollection> taming_functions;
//----- run operations
/// Reset the total generation time and the number of events generated for this run, prepare kinematics
void prepareRun();
/// Add a new timing into the total generation time
/// \param[in] gen_time Time to add (in seconds)
void addGenerationTime( double gen_time );
/// Return the total generation time for this run (in seconds)
inline double totalGenerationTime() const { return total_gen_time_; }
/// Total number of events already generated in this run
inline unsigned int numGeneratedEvents() const { return num_gen_events_; }
private:
std::unique_ptr<proc::GenericProcess> process_;
std::unique_ptr<hadr::GenericHadroniser> hadroniser_;
bool store_;
/// Total generation time (in seconds)
double total_gen_time_;
/// Number of events already generated
unsigned int num_gen_events_;
};
}
#endif
diff --git a/CepGen/Processes/GenericKTProcess.cpp b/CepGen/Processes/GenericKTProcess.cpp
index f84fc04..da287eb 100644
--- a/CepGen/Processes/GenericKTProcess.cpp
+++ b/CepGen/Processes/GenericKTProcess.cpp
@@ -1,378 +1,380 @@
#include "CepGen/Processes/GenericKTProcess.h"
#include "CepGen/Core/Exception.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 proc
{
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( params, 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 )
{
kin_ = kin;
const KTFlux flux1 = (KTFlux)kin_.incoming_beams.first.kt_flux,
flux2 = (KTFlux)kin_.incoming_beams.second.kt_flux;
if ( kin_.mode == KinematicsMode::invalid ) {
//--- try to extrapolate kinematics mode from unintegrated fluxes
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 )
kin_.mode = KinematicsMode::ElasticElastic;
else if ( el1 )
kin_.mode = KinematicsMode::ElasticInelastic;
else if ( el2 )
kin_.mode = KinematicsMode::InelasticElastic;
else
kin_.mode = KinematicsMode::InelasticInelastic;
}
else {
//--- try to extrapolate unintegrated fluxes from kinematics mode
const HeavyIon hi1( kin_.incoming_beams.first.pdg ), hi2( kin_.incoming_beams.second.pdg );
//==========================================================================================
// ensure the first incoming flux is compatible with the kinematics mode
//==========================================================================================
if ( ( kin_.mode == KinematicsMode::ElasticElastic
|| kin_.mode == KinematicsMode::ElasticInelastic )
&& ( flux1 != KTFlux::P_Photon_Elastic ) ) {
kin_.incoming_beams.first.kt_flux = hi1
? KTFlux::HI_Photon_Elastic
: KTFlux::P_Photon_Elastic;
CG_DEBUG( "GenericKTProcess:kinematics" )
<< "Set the kt flux for first incoming photon to \""
<< kin_.incoming_beams.first.kt_flux << "\".";
}
else if ( flux1 != KTFlux::P_Photon_Inelastic
&& flux1 != KTFlux::P_Photon_Inelastic_Budnev ) {
if ( hi1 )
throw CG_FATAL( "GenericKTProcess:kinematics" )
<< "Inelastic photon emission from HI not yet supported!";
kin_.incoming_beams.first.kt_flux = KTFlux::P_Photon_Inelastic_Budnev;
CG_DEBUG( "GenericKTProcess:kinematics" )
<< "Set the kt flux for first incoming photon to \""
<< kin_.incoming_beams.first.kt_flux << "\".";
}
//==========================================================================================
// ensure the second incoming flux is compatible with the kinematics mode
//==========================================================================================
if ( ( kin_.mode == KinematicsMode::ElasticElastic
|| kin_.mode == KinematicsMode::InelasticElastic )
&& ( flux2 != KTFlux::P_Photon_Elastic ) ) {
kin_.incoming_beams.second.kt_flux = hi2
? KTFlux::HI_Photon_Elastic
: KTFlux::P_Photon_Elastic;
CG_DEBUG( "GenericKTProcess:kinematics" )
<< "Set the kt flux for second incoming photon to \""
<< kin_.incoming_beams.second.kt_flux << "\".";
}
else if ( flux2 != KTFlux::P_Photon_Inelastic
&& flux2 != KTFlux::P_Photon_Inelastic_Budnev ) {
if ( hi2 )
throw CG_FATAL( "GenericKTProcess:kinematics" )
<< "Inelastic photon emission from HI not yet supported!";
kin_.incoming_beams.second.kt_flux = KTFlux::P_Photon_Inelastic_Budnev;
CG_DEBUG( "GenericKTProcess:kinematics" )
<< "Set the kt flux for second incoming photon to \""
<< kin_.incoming_beams.second.kt_flux << "\".";
}
}
+std::cout << kin_.incoming_beams.first << "|" << kin_.incoming_beams.second << std::endl;
+
//============================================================================================
// 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, kin_.cuts.initial.qt, { 1.e-10, 500. }, "First incoming parton virtuality" );
registerVariable( qt2_, Mapping::logarithmic, kin_.cuts.initial.qt, { 1.e-10, 500. }, "Second incoming parton virtuality" );
registerVariable( phi_qt1_, Mapping::linear, kin_.cuts.initial.phi_qt, { 0., 2.*M_PI }, "First incoming parton azimuthal angle" );
registerVariable( phi_qt2_, Mapping::linear, kin_.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_ = event_->getOneByRole( Particle::IncomingBeam1 ).mass();
MY_ = event_->getOneByRole( Particle::IncomingBeam2 ).mass();
if ( kin_.mode == KinematicsMode::InelasticElastic || kin_.mode == KinematicsMode::InelasticInelastic )
registerVariable( MX_, Mapping::square, kin_.cuts.remnants.mass_single, { 1.07, 1000. }, "Positive z proton remnant mass" );
if ( kin_.mode == KinematicsMode::ElasticInelastic || kin_.mode == KinematicsMode::InelasticInelastic )
registerVariable( MY_, Mapping::square, kin_.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
{
const HeavyIon hi1( kin_.incoming_beams.first.pdg ), hi2( kin_.incoming_beams.second.pdg );
//--- compute fluxes according to modelling specified in parameters card
std::pair<double,double> fluxes = {
hi1
? ktFlux( (KTFlux)kin_.incoming_beams.first.kt_flux, x1, q1t2, hi1 )
: ktFlux( (KTFlux)kin_.incoming_beams.first.kt_flux, x1, q1t2, *kin_.structure_functions, MX_ ),
hi2
? ktFlux( (KTFlux)kin_.incoming_beams.second.kt_flux, x2, q2t2, hi2 )
: ktFlux( (KTFlux)kin_.incoming_beams.second.kt_flux, x2, q2t2, *kin_.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 ( kin_.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/test/cepgen.cpp b/test/cepgen.cpp
index 25dc453..c58f1ca 100644
--- a/test/cepgen.cpp
+++ b/test/cepgen.cpp
@@ -1,66 +1,70 @@
//--- steering cards
#include "CepGen/Cards/Handler.h"
#include "CepGen/Generator.h"
#include "CepGen/Core/Exception.h"
//--- necessary include to build the default run
#include "CepGen/Physics/PDG.h"
#include "CepGen/Processes/ProcessesHandler.h"
#include "AbortHandler.h"
#include <iostream>
using namespace std;
/**
* Main caller for this MC generator.
* * loads the configuration files' variables if passed as an argument,
* or a default LPAIR-like configuration,
* * launches the cross-section computation and the events generation.
* \author Laurent Forthomme <laurent.forthomme@cern.ch>
*/
int main( int argc, char* argv[] )
{
//--- first start by defining the generator object
cepgen::Generator gen;
if ( argc < 2 ) {
CG_INFO( "main" ) << "No config file provided. Setting the default parameters.";
//--- default run: LPAIR elastic ɣɣ → µ⁺µ¯ at 13 TeV
cepgen::ParametersList pgen;
pgen.set<int>( "pair", (int)cepgen::PDG::muon );
- gen.parameters->setProcess( cepgen::proc::ProcessesHandler::get().build( "lpair", pgen ) );
- gen.parameters->kinematics.mode = cepgen::KinematicsMode::ElasticElastic;
- gen.parameters->kinematics.cuts.central.pt_single.min() = 15.;
- gen.parameters->kinematics.cuts.central.eta_single = { -2.5, 2.5 };
- gen.parameters->generation.enabled = true;
- gen.parameters->generation.maxgen = 1e3;
+
+ cepgen::Parameters params;
+ params.setProcess( cepgen::proc::ProcessesHandler::get().build( "lpair", pgen ) );
+ params.kinematics.mode = cepgen::KinematicsMode::ElasticElastic;
+ params.kinematics.cuts.central.pt_single.min() = 15.;
+ params.kinematics.cuts.central.eta_single = { -2.5, 2.5 };
+ params.generation.enabled = true;
+ params.generation.maxgen = 1e3;
+ gen.setParameters( params );
}
else
gen.setParameters( cepgen::card::Handler::parse( argv[1] ) );
+
//--- list all parameters
- CG_LOG( "main" ) << gen.parameters.get();
+ CG_LOG( "main" ) << gen.parameters();
cepgen::utils::AbortHandler ctrl_c;
try {
//--- let there be a cross-section...
double xsec = 0., err = 0.;
gen.computeXsection( xsec, err );
- if ( gen.parameters->generation.enabled )
+ if ( gen.parameters().generation.enabled )
//--- events generation starts here
// (one may use a callback function)
gen.generate();
} catch ( const cepgen::utils::RunAbortedException& e ) {
CG_INFO( "main" ) << "Run aborted!";
} catch ( const cepgen::Exception& e ) {
e.dump();
}
return 0;
}