No OneTemporary
Actions

Size

36 KB

Subscribers

None

View Options

	diff --git a/Cards/pptoww_cfg.py b/Cards/pptoww_cfg.py
	index d9e2a44..55756ca 100644
	--- a/Cards/pptoww_cfg.py
	+++ b/Cards/pptoww_cfg.py
	@@ -1,67 +1,67 @@
	import Config.Core as cepgen
	from Config.Integration.vegas_cff import integrator
	#from Config.logger_cfi import logger
	#logger.enabledModules += ('Hadroniser.configure', 'Generator.*',)
	#--------------------------------------------------------------------
	# Pythia 6 example (with fully leptonic WW decay)
	#--------------------------------------------------------------------
	#from Config.Hadronisation.pythia6_cff import pythia6 as hadroniser
	#from Config.Hadronisation.pythia6Defaults_cfi import WDecayToEMu
	#hadroniser.wDecays = WDecayToEMu
	#hadroniser.processConfiguration += ('wDecays',)
	#hadroniser.remnantsFragmentation = False
	#--------------------------------------------------------------------
	# Pythia 8 example (with fully leptonic WW decay)
	#--------------------------------------------------------------------
	#from Config.Hadronisation.pythia8_cff import pythia8
	#hadroniser = pythia8.clone('pythia8',
	# pythiaConfiguration = (
	# # process-specific
	# '13:onMode = off', # disable muon decays
	# '24:onMode = off', # disable all W decays, but...
	# #'24:onIfAny = 11 13', # enable e-nue + mu-numu final states
	# '24:onPosIfAny = 11', # enable W- -> e- + nu_e decay
	# '24:onNegIfAny = 13', # enable W+ -> mu+ + nu_mu decay
	# ),
	# processConfiguration = pythia8.processConfiguration+('pythiaConfiguration',),
	#)

	import Config.ktProcess_cfi as kt
	process = kt.process.clone('pptoww',
	processParameters = cepgen.Parameters(
	mode = cepgen.ProcessMode.InelasticInelastic,
	polarisationStates = 0, # full
	),
	inKinematics = cepgen.Parameters(
	cmEnergy = 13.e3,
	#structureFunctions = cepgen.StructureFunctions.SzczurekUleshchenko,
	#structureFunctions = cepgen.StructureFunctions.ALLM97,
	structureFunctions = cepgen.StructureFunctions.LUXlike,
	),
	outKinematics = kt.process.outKinematics.clone(
	mx = (1.07, 1000.),
	qt = (0., 1000.),
	#--- extra cuts on the pt(W+) and pt(W-) plane
	ptdiff = (0., 2000.),
	#--- extra cuts on the W+W- system
	invmass = (0.,),
	ptsum = (0.,),
	#--- cuts on single particles' level
	cuts = {
	# cuts on the single W level
	24: cepgen.Parameters(pt = (0.,)), # no pt cut on Ws
	# cuts on the W decay products
	# (mimicking LHC-like experimental cuts)
	- 11: cepgen.Parameters(pt = (20.,), eta = (-2.5, 2.5)),
	- 13: cepgen.Parameters(pt = (20.,), eta = (-2.5, 2.5))
	+ #11: cepgen.Parameters(pt = (20.,), eta = (-2.5, 2.5)),
	+ #13: cepgen.Parameters(pt = (20.,), eta = (-2.5, 2.5))
	},
	#xi = (0.02, 0.15),
	)
	)

	#--- generation parameters
	from Config.generator_cff import generator
	generator = generator.clone(
	numEvents = 10000,
	printEvery = 1000,
	)
	diff --git a/CepGen/Core/Integrator.cpp b/CepGen/Core/Integrator.cpp
	index 711463a..5071230 100644
	--- a/CepGen/Core/Integrator.cpp
	+++ b/CepGen/Core/Integrator.cpp
	@@ -1,469 +1,469 @@
	#include "CepGen/Core/Integrator.h"
	#include "CepGen/Core/GridParameters.h"
	#include "CepGen/Core/utils.h"
	#include "CepGen/Core/Exception.h"

	#include "CepGen/Parameters.h"
	#include "CepGen/Processes/GenericProcess.h"
	#include "CepGen/Hadronisers/GenericHadroniser.h"
	#include "CepGen/Event/Event.h"

	#include <thread>
	#include <math.h>

	#include <gsl/gsl_monte_miser.h>
	#define COORD(s,i,j) ((s)->xi[(i)*(s)->dim + (j)])

	namespace cepgen
	{
	Integrator::Integrator( unsigned int ndim, double integrand( double, size_t, void ), Parameters& params ) :
	ps_bin_( INVALID_BIN ), input_params_( params ),
	function_( new gsl_monte_function{ integrand, ndim, (void*)&input_params_ } ),
	rng_( gsl_rng_alloc( input_params_.integration().rng_engine ) ),
	grid_( new GridParameters( ndim ) )
	{
	//--- initialise the random number generator

	unsigned long seed = ( input_params_.integration().rng_seed > 0 )
	? input_params_.integration().rng_seed
	: time( nullptr ); // seed with time
	gsl_rng_set( rng_.get(), seed );

	//--- a bit of printout for debugging

	CG_DEBUG( "Integrator:build" )
	<< "Number of integration dimensions: " << function_->dim << ",\n\t"
	<< "Number of function calls: " << input_params_.integration().ncvg << ",\n\t"
	<< "Random numbers generator: " << gsl_rng_name( rng_.get() ) << ".";
	switch ( input_params_.integration().type ) {
	case IntegratorType::Vegas:
	CG_DEBUG( "Integrator:build" ) << "Vegas parameters:\n\t"
	<< "Number of iterations in Vegas: " << input_params_.integration().vegas.iterations << ",\n\t"
	<< "α-value: " << input_params_.integration().vegas.alpha << ",\n\t"
	<< "Verbosity: " << input_params_.integration().vegas.verbose << ",\n\t"
	<< "Grid interpolation mode: " << (Integrator::VegasMode)input_params_.integration().vegas.mode << ".";
	break;
	case IntegratorType::MISER:
	CG_DEBUG( "Integrator:build" ) << "MISER parameters:\n\t"
	<< "Number of calls: " << input_params_.integration().miser.min_calls << ", "
	<< "per bisection: " << input_params_.integration().miser.min_calls_per_bisection << ",\n\t"
	<< "Estimate fraction: " << input_params_.integration().miser.estimate_frac << ",\n\t"
	<< "α-value: " << input_params_.integration().miser.alpha << ",\n\t"
	<< "Dither: " << input_params_.integration().miser.dither << ".";
	break;
	case IntegratorType::plain:
	break;
	}
	}

	Integrator::~Integrator()
	{}

	//-----------------------------------------------------------------------------------------------
	// integration part
	//-----------------------------------------------------------------------------------------------

	void
	Integrator::integrate( double& result, double& abserr )
	{
	int res = -1;

	//--- integration bounds
	std::vector<double> x_low( function_->dim, 0. ), x_up( function_->dim, 1. );

	//--- launch integration
	switch ( input_params_.integration().type ) {
	case IntegratorType::plain: {
	std::unique_ptr<gsl_monte_plain_state,void()( gsl_monte_plain_state )>
	pln_state( gsl_monte_plain_alloc( function_->dim ), gsl_monte_plain_free );
	res = gsl_monte_plain_integrate( function_.get(),
	&x_low[0], &x_up[0],
	function_->dim, input_params_.integration().ncvg,
	rng_.get(), pln_state.get(),
	&result, &abserr );
	} break;
	case IntegratorType::Vegas: {
	//----- warmup (prepare the grid)
	warmupVegas( x_low, x_up, 25000 );
	//----- integration
	unsigned short it_chisq = 0;
	do {
	res = gsl_monte_vegas_integrate( function_.get(),
	&x_low[0], &x_up[0],
	function_->dim, 0.2 * input_params_.integration().ncvg,
	rng_.get(), veg_state_.get(),
	&result, &abserr );
	CG_LOG( "Integrator:integrate" )
	<< "\t>> at call " << ( ++it_chisq ) << ": "
	<< Form( "average = %10.6f "
	"sigma = %10.6f chi2 = %4.3f.",
	result, abserr,
	gsl_monte_vegas_chisq( veg_state_.get() ) );
	} while ( fabs( gsl_monte_vegas_chisq( veg_state_.get() )-1. )
	> input_params_.integration().vegas_chisq_cut-1. );
	CG_DEBUG( "Integrator:integrate" )
	<< "Vegas grid information:\n\t"
	<< "ran for " << veg_state_->dim << " dimensions, and generated " << veg_state_->bins_max << " bins.\n\t"
	<< "Integration volume: " << veg_state_->vol << ".";
	grid_->r_boxes = std::pow( veg_state_->bins, function_->dim );
	} break;
	case IntegratorType::MISER: {
	std::unique_ptr<gsl_monte_miser_state,void()( gsl_monte_miser_state )>
	mis_state( gsl_monte_miser_alloc( function_->dim ), gsl_monte_miser_free );
	gsl_monte_miser_params_set( mis_state.get(), &input_params_.integration().miser );
	res = gsl_monte_miser_integrate( function_.get(),
	&x_low[0], &x_up[0],
	function_->dim, input_params_.integration().ncvg,
	rng_.get(), mis_state.get(),
	&result, &abserr );
	} break;
	}

	input_params_.integration().result = result;
	input_params_.integration().err_result = abserr;

	if ( input_params_.hadroniser() )
	input_params_.hadroniser()->setCrossSection( result, abserr );

	if ( res != GSL_SUCCESS )
	throw CG_FATAL( "Integrator:integrate" )
	<< "Error while performing the integration!\n\t"
	<< "GSL error: " << gsl_strerror( res ) << ".";
	}

	void
	Integrator::warmupVegas( std::vector<double>& x_low, std::vector<double>& x_up, unsigned int ncall )
	{
	// start by preparing the grid/state
	veg_state_.reset( gsl_monte_vegas_alloc( function_->dim ) );
	gsl_monte_vegas_params_set( veg_state_.get(), &input_params_.integration().vegas );
	// then perform a first integration with the given calls count
	double result = 0., abserr = 0.;
	const int res = gsl_monte_vegas_integrate( function_.get(),
	&x_low[0], &x_up[0],
	function_->dim, ncall, rng_.get(), veg_state_.get(),
	&result, &abserr );
	// ensure the operation was successful
	if ( res != GSL_SUCCESS )
	throw CG_ERROR( "Integrator:vegas" )
	<< "Failed to warm-up the Vegas grid.\n\t"
	<< "GSL error: " << gsl_strerror( res ) << ".";
	CG_INFO( "Integrator:vegas" )
	<< "Finished the Vegas warm-up.";
	}

	//-----------------------------------------------------------------------------------------------
	// events generation part
	//-----------------------------------------------------------------------------------------------

	void
	Integrator::generateOne( std::function<void( const Event&, unsigned long )> callback )
	{
	if ( !grid_->gen_prepared )
	computeGenerationParameters();

	std::vector<double> xtmp;

	//--- correction cycles

	if ( ps_bin_ != INVALID_BIN ) {
	bool has_correction = false;
	while ( !correctionCycle( xtmp, has_correction ) ) {}
	if ( has_correction ) {
	storeEvent( xtmp, callback );
	return;
	}
	}

	double weight = 0.;

	//--- normal generation cycle

	while ( true ) {
	double y = -1.;
	//----- select a and reject if fmax is too small
	while ( true ) {
	// ...
	ps_bin_ = uniform() * grid_->size();
	y = uniform() * grid_->globalMax();
	grid_->num[ps_bin_] += 1;
	if ( y <= grid_->maxValue( ps_bin_ ) )
	break;
	}
	// shoot a point x in this bin
	grid_->shoot( rng_.get(), ps_bin_, xtmp );
	// get weight for selected x value
	weight = eval( xtmp );
	if ( weight <= 0. )
	continue;
	if ( weight > y )
	break;
	}

	if ( weight <= grid_->maxValue( ps_bin_ ) )
	ps_bin_ = INVALID_BIN;
	else {
	//--- if weight is higher than local or global maximum,
	// init correction cycle
	grid_->f_max_old = grid_->maxValue( ps_bin_ );
	grid_->f_max_diff = weight-grid_->f_max_old;
	grid_->setValue( ps_bin_, weight );
	grid_->correc = ( grid_->num[ps_bin_]-1. ) * grid_->f_max_diff / grid_->globalMax() - 1.;

	CG_DEBUG("Integrator::generateOne")
	<< "Correction " << grid_->correc << " will be applied for phase space bin " << ps_bin_ << ".";
	}

	// return with an accepted event
	if ( weight > 0. )
	storeEvent( xtmp, callback );
	}

	void
	Integrator::generate( unsigned long num_events, std::function<void( const Event&, unsigned long )> callback )
	{
	if ( num_events < 1 )
	num_events = input_params_.generation().maxgen;
	try {
	while ( input_params_.numGeneratedEvents() < num_events )
	generateOne( callback );
	} catch ( const Exception& e ) { return; }
	}

	bool
	Integrator::correctionCycle( std::vector<double>& x, bool& has_correction )
	{
	CG_DEBUG_LOOP( "Integrator:correction" )
	<< "Correction cycles are started.\n\t"
	<< "bin = " << ps_bin_ << "\t"
	<< "correc = " << grid_->correc << "\t"
	<< "corre2 = " << grid_->correc2 << ".";

	if ( grid_->correc >= 1. )
	grid_->correc -= 1.;

	if ( uniform() < grid_->correc ) {
	grid_->correc = -1.;
	std::vector<double> xtmp( function_->dim );
	// Select x values in phase space bin
	grid_->shoot( rng_.get(), ps_bin_, xtmp );
	const double weight = eval( xtmp );
	// Parameter for correction of correction
	if ( weight > grid_->maxValue( ps_bin_ ) ) {
	grid_->f_max2 = std::max( grid_->f_max2, weight );
	grid_->correc += 1.;
	grid_->correc2 -= 1.;
	}
	// Accept event
	if ( weight >= grid_->f_max_diff*uniform() + grid_->f_max_old ) {
	x = xtmp;
	has_correction = true;
	return true;
	}
	return false;
	}
	// Correction if too big weight is found while correction
	// (All your bases are belong to us...)
	if ( grid_->f_max2 > grid_->maxValue( ps_bin_ ) ) {
	grid_->f_max_old = grid_->maxValue( ps_bin_ );
	grid_->f_max_diff = grid_->f_max2-grid_->f_max_old;
	grid_->correc = ( grid_->num[ps_bin_]-1. ) * grid_->f_max_diff / grid_->globalMax();
	if ( grid_->f_max2 >= grid_->globalMax() )
	grid_->correc *= grid_->f_max2 / grid_->globalMax();
	grid_->setValue( ps_bin_, grid_->f_max2 );
	grid_->correc -= grid_->correc2;
	grid_->correc2 = 0.;
	grid_->f_max2 = 0.;
	return false;
	}
	return true;
	}

	bool
	Integrator::storeEvent( const std::vector<double>& x, std::function<void( const Event&, unsigned long )> callback )
	{
	//--- start by computing the matrix element for that point
	const double weight = eval( x );

	//--- reject if unphysical
	if ( weight <= 0. )
	return false;

	{
	if ( input_params_.numGeneratedEvents() % input_params_.generation().gen_print_every == 0 ) {
	CG_INFO( "Integrator:store" )
	<< "Generated events: " << input_params_.numGeneratedEvents();
	input_params_.process()->last_event->dump();
	}
	const Event& last_event = *input_params_.process()->last_event;
	if ( callback )
	callback( last_event, input_params_.numGeneratedEvents() );
	input_params_.addGenerationTime( last_event.time_total );
	}
	return true;
	}

	//-----------------------------------------------------------------------------------------------
	// initial preparation run before the generation of unweighted events
	//-----------------------------------------------------------------------------------------------

	void
	Integrator::computeGenerationParameters()
	{
	input_params_.setStorage( false );

	if ( input_params_.generation().treat
	&& input_params_.integration().type != IntegratorType::Vegas ) {
	CG_INFO( "Integrator:setGen" )
	<< "Treat switched on without a proper Vegas grid; running a warm-up beforehand.";
	std::vector<double> x_low( function_->dim, 0. ), x_up( function_->dim, 1. );
	try {
	warmupVegas( x_low, x_up, 25000 );
	} catch ( const Exception& ) {
	throw CG_FATAL( "Integrator::setGen" )
	<< "Failed to perform a Vegas warm-up.\n\t"
	<< "Try to re-run while disabling integrand treatment...";
	}
	}
	CG_INFO( "Integrator:setGen" )
	<< "Preparing the grid (" << input_params_.generation().num_points << " points/bin) "
	<< "for the generation of unweighted events.";

	const double inv_num_points = 1./input_params_.generation().num_points;
	std::vector<double> x( function_->dim, 0. );
	std::vector<unsigned short> n( function_->dim, 0 );;

	// ...
	double sum = 0., sum2 = 0., sum2p = 0.;

	- utils::ProgressBar prog_bar( grid_->size() );
	+ utils::ProgressBar prog_bar( grid_->size(), 5 );

	//--- main loop
	for ( unsigned int i = 0; i < grid_->size(); ++i ) {
	double fsum = 0., fsum2 = 0.;
	for ( unsigned int j = 0; j < input_params_.generation().num_points; ++j ) {
	grid_->shoot( rng_.get(), i, x );
	const double weight = eval( x );
	grid_->setValue( i, weight );
	fsum += weight;
	fsum2 += weight*weight;
	}
	const double av = fsuminv_num_points, av2 = fsum2inv_num_points, sig2 = av2-av*av;
	sum += av;
	sum2 += av2;
	sum2p += sig2;

	// per-bin debugging loop
	if ( CG_LOG_MATCH( "Integrator:setGen", debugInsideLoop ) ) {
	const double sig = sqrt( sig2 );
	const double eff = ( grid_->maxValue( i ) != 0. )
	? grid_->maxValue( i )/av
	: 1.e4;
	CG_DEBUG_LOOP( "Integrator:setGen" )
	<< "n-vector for bin " << i << ": " << utils::repr( grid_->n( i ) ) << "\n\t"
	<< "av = " << av << "\n\t"
	<< "sig = " << sig << "\n\t"
	<< "fmax = " << grid_->maxValue( i ) << "\n\t"
	<< "eff = " << eff;
	}
	prog_bar.update( i+1 );
	} // end of main loop

	const double inv_max = 1./grid_->size();
	sum *= inv_max;
	sum2 *= inv_max;
	sum2p *= inv_max;

	const double sig = sqrt( sum2-sum*sum ), sigp = sqrt( sum2p );

	double eff1 = 0.;
	for ( unsigned int i = 0; i < grid_->size(); ++i )
	eff1 += sum/grid_->size()*grid_->maxValue( i );
	const double eff2 = sum/grid_->globalMax();

	CG_DEBUG( "Integrator:setGen" )
	<< "Average function value = " << sum << "\n\t"
	<< "Average squared function value = " << sum2 << "\n\t"
	<< "Overall standard deviation = " << sig << "\n\t"
	<< "Average standard deviation = " << sigp << "\n\t"
	<< "Maximum function value = " << grid_->globalMax() << "\n\t"
	<< "Average inefficiency = " << eff1 << "\n\t"
	<< "Overall inefficiency = " << eff2;

	grid_->gen_prepared = true;
	input_params_.setStorage( true );
	CG_INFO( "Integrator:setGen" ) << "Grid prepared! Now launching the production.";
	}

	//------------------------------------------------------------------------------------------------
	// helper / alias methods
	//------------------------------------------------------------------------------------------------

	unsigned short
	Integrator::dimensions() const
	{
	if ( !function_ )
	return 0;
	return function_->dim;
	}

	double
	Integrator::eval( const std::vector<double>& x )
	{
	if ( !input_params_.generation().treat )
	return function_->f( (double)&x[0], function_->dim, (void)&input_params_ );
	//--- treatment of the integration grid
	double w = grid_->r_boxes;
	std::vector<double> x_new( x.size() );
	for ( unsigned short j = 0; j < function_->dim; ++j ) {
	//--- find surrounding coordinates and interpolate
	const double z = x[j]*veg_state_->bins;
	const unsigned int id = z; // coordinate of point before
	const double rel_pos = z-id; // position between coordinates (norm.)
	const double bin_width = ( id == 0 )
	? COORD( veg_state_, 1, j )
	: COORD( veg_state_, id+1, j )-COORD( veg_state_, id, j );
	//--- build new coordinate from linear interpolation
	x_new[j] = COORD( veg_state_, id+1, j )-bin_width*( 1.-rel_pos );
	w *= bin_width;
	}
	return wfunction_->f( (double)&x_new[0], function_->dim, (void*)&input_params_ );
	}

	double
	Integrator::uniform() const
	{
	return gsl_rng_uniform( rng_.get() );
	}

	//------------------------------------------------------------------------------------------------

	std::ostream&
	operator<<( std::ostream& os, const IntegratorType& type )
	{
	switch ( type ) {
	case IntegratorType::plain:
	return os << "plain";
	case IntegratorType::Vegas:
	return os << "Vegas";
	case IntegratorType::MISER:
	return os << "MISER";
	}
	return os;
	}

	std::ostream&
	operator<<( std::ostream& os, const Integrator::VegasMode& mode )
	{
	switch ( mode ) {
	case Integrator::VegasMode::importance:
	return os << "importance";
	case Integrator::VegasMode::importanceOnly:
	return os << "importance-only";
	case Integrator::VegasMode::stratified:
	return os << "stratified";
	}
	return os;
	}
	}

	diff --git a/CepGen/Core/Parameters.cpp b/CepGen/Core/Parameters.cpp
	index a66083a..fc699bf 100644
	--- a/CepGen/Core/Parameters.cpp
	+++ b/CepGen/Core/Parameters.cpp
	@@ -1,307 +1,309 @@
	#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_( 0ul )
	{}

	Parameters::Parameters( Parameters& param ) :
	general( param.general ),
	kinematics( param.kinematics ),
	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_ ),
	integration_( param.integration_ ), generation_( param.generation_ )
	{}

	Parameters::Parameters( const Parameters& param ) :
	general( param.general ),
	kinematics( param.kinematics ),
	taming_functions( param.taming_functions ),
	store_( false ), total_gen_time_( param.total_gen_time_ ), num_gen_events_( param.num_gen_events_ ),
	integration_( param.integration_ ), generation_( param.generation_ )
	{}

	Parameters::~Parameters() // required for unique_ptr initialisation!
	{}

	Parameters&
	Parameters::operator=( Parameters param )
	{
	general = param.general;
	kinematics = param.kinematics;
	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_;
	integration_ = param.integration_;
	generation_ = param.generation_;
	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 )
	return;
	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();
	//--- clear the run statistics
	total_gen_time_ = 0.;
	num_gen_events_ = 0ul;
	process_->first_run = false;
	}

	void
	Parameters::addGenerationTime( double gen_time )
	{
	total_gen_time_ += gen_time;
	num_gen_events_++;
	}

	proc::GenericProcess*
	Parameters::process()
	{
	return process_.get();
	}

	const proc::GenericProcess*
	Parameters::process() const
	{
	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* param )
	{
	const bool pretty = true;

	const int wb = 90, wt = 33;

	os << std::left
	<< "\n"
	<< 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( param->processName().c_str() ) : param->processName() );
	if ( param->process_ ) {
	os << ", " << param->process_->description();
	for ( const auto& par : param->process()->parameters().keys() )
	if ( par != "mode" )
	os << "\n" << std::setw( wt ) << "" << par << ": " << param->process_->parameters().getString( par );
	std::ostringstream proc_mode; proc_mode << param->kinematics.mode;
	if ( param->kinematics.mode != KinematicsMode::invalid )
	os << "\n" << std::setw( wt ) << "Subprocess mode" << ( pretty ? boldify( proc_mode.str().c_str() ) : proc_mode.str() ) << "\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( param->generation_.enabled ) : std::to_string( param->generation_.enabled ) ) << "\n"
	<< std::setw( wt ) << "Number of events to generate"
	<< ( pretty ? boldify( param->generation_.maxgen ) : std::to_string( param->generation_.maxgen ) ) << "\n";
	if ( param->generation_.num_threads > 1 )
	os
	<< std::setw( wt ) << "Number of threads" << param->generation_.num_threads << "\n";
	os
	<< std::setw( wt ) << "Number of points to try per bin" << param->generation_.num_points << "\n"
	<< std::setw( wt ) << "Integrand treatment"
	<< ( pretty ? yesno( param->generation_.treat ) : std::to_string( param->generation_.treat ) ) << "\n"
	<< std::setw( wt ) << "Verbosity level " << utils::Logger::get().level << "\n";
	if ( param->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( param->hadroniser_->name().c_str() ) : param->hadroniser_->name() ) << "\n";
	+ << ( pretty ? boldify( param->hadroniser_->name().c_str() ) : param->hadroniser_->name() ) << "\n"
	+ << std::setw( wt ) << "Remnants fragmentation? "
	+ << ( pretty ? yesno( param->hadroniser_->fragmentRemnants() ) : std::to_string( param->hadroniser_->fragmentRemnants() ) ) << "\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 << param->integration_.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" << param->integration_.ncvg << "\n"
	<< std::setw( wt ) << "Random number generator seed" << param->integration_.rng_seed << "\n";
	if ( param->integration_.rng_engine )
	os
	<< std::setw( wt ) << "Random number generator engine"
	<< param->integration_.rng_engine->name << "\n";
	os
	<< "\n"
	<< std::setfill('_') << std::setw( wb+3 ) << "_/¯ EVENTS KINEMATICS ¯\\_" << std::setfill( ' ' ) << "\n\n"
	<< std::setw( wt ) << "Incoming particles"
	<< param->kinematics.incoming_beams.first << ",\n" << std::setw( wt ) << ""
	<< param->kinematics.incoming_beams.second << "\n"
	<< std::setw( wt ) << "C.m. energy (GeV)" << param->kinematics.sqrtS() << "\n";
	if ( param->kinematics.mode != KinematicsMode::ElasticElastic )
	os << std::setw( wt ) << "Structure functions" << *param->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 : param->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 : param->kinematics.cuts.central.list() )
	if ( lim.second.valid() )
	os << std::setw( wt ) << lim.first << lim.second << "\n";
	if ( param->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 : param->kinematics.cuts.central_particles ) {
	os << " * all single " << std::setw( wt-3 ) << PDG::get().name( 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";
	for ( const auto& lim : param->kinematics.cuts.remnants.list() )
	os << "\n" << std::setw( wt ) << lim.first << lim.second;
	os << "\n";
	return os;
	}

	//-----------------------------------------------------------------------------------------------

	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; // random number of dimensions for VEGAS parameters retrieval
	{
	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 ), gen_print_every( 10000 ),
	num_threads( 2 ), num_points( 100 )
	{}

	Parameters::Generation::Generation( const Generation& rhs ) :
	enabled( rhs.enabled ), maxgen( rhs.maxgen ),
	symmetrise( rhs.symmetrise ), treat( rhs.treat ), gen_print_every( rhs.gen_print_every ),
	num_threads( rhs.num_threads ), num_points( rhs.num_points )
	{}
	}

	diff --git a/CepGen/Hadronisers/GenericHadroniser.h b/CepGen/Hadronisers/GenericHadroniser.h
	index 00e6057..82d7ed3 100644
	--- a/CepGen/Hadronisers/GenericHadroniser.h
	+++ b/CepGen/Hadronisers/GenericHadroniser.h
	@@ -1,74 +1,76 @@
	#ifndef CepGen_Hadronisers_GenericHadroniser_h
	#define CepGen_Hadronisers_GenericHadroniser_h

	#include <vector>
	#include <memory>
	#include <iosfwd>

	namespace cepgen
	{
	class Event;
	class Particle;
	class Parameters;
	class ParametersList;
	/// Location for all hadronisers to be run downstream to the events generation
	namespace hadr
	{
	/**
	* \brief Class template to define any hadroniser as a general object with defined methods
	* \author Laurent Forthomme <laurent.forthomme@cern.ch>
	* \date January 2014
	*/
	class GenericHadroniser
	{
	public:
	/// Write out all hadroniser attributes in output stream
	friend std::ostream& operator<<( std::ostream& os, const GenericHadroniser& hadr );
	/// Write out all hadroniser attributes in output stream
	friend std::ostream& operator<<( std::ostream& os, const GenericHadroniser* hadr );

	/// Default constructor for an undefined hadroniser
	explicit GenericHadroniser( const ParametersList&, const std::string& name = "<invalid hadroniser>" );
	virtual ~GenericHadroniser() {}

	/// Parse a configuration string
	virtual void readString( const char* ) {}
	/// Parse a configuration string
	virtual void readString( const std::string& param ) { readString( param.c_str() ); }
	/// Parse a list of configuration strings
	virtual void readStrings( const std::vector<std::string>& params );
	/// Initialise the event hadroniser before its running
	virtual void init() = 0;
	/** \brief Hadronise a full event
	* \param[inout] ev Event to hadronise
	* \param[inout] weight Event weight after hadronisation
	* \param[in] full Perform the full state hadronisation (incl. remnants fragmentation)
	* \return Boolean stating whether or not the hadronisation occured successfully
	*/
	virtual bool run( Event& ev, double& weight, bool full ) = 0;
	/// Specify the process cross section, in pb
	virtual void setCrossSection( double xsec, double xsec_err ) {}

	virtual void setParameters( const Parameters& params ) { params_ = &params; }
	/// \brief Specify a random numbers generator seed for the hadroniser
	/// \param[in] seed A RNG seed
	void setSeed( long long seed ) { seed_ = seed; }
	+ /// Specify whether the beam remnants are to be fragmented
	+ bool fragmentRemnants() const { return remn_fragm_; }

	/// Return a human-readable name for this hadroniser
	std::string name() const;

	protected:
	/// Name of the hadroniser
	std::string name_;
	/// Random numbers generator seed for the hadroniser
	long long seed_;
	/// Maximal number of trials for the hadronisation of the proton(s) remnants
	unsigned short max_trials_;
	const Parameters* params_; // not owning
	/// Switch on/off the remnants fragmentation where applicable
	const bool remn_fragm_;
	};
	}
	}

	#endif

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 7:53 PM (1 d, 5 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805931
Default Alt Text: (36 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions