diff --git a/CepGen/Core/Generator.cpp b/CepGen/Core/Generator.cpp
index 7d97c71..b56e13b 100644
--- a/CepGen/Core/Generator.cpp
+++ b/CepGen/Core/Generator.cpp
@@ -1,139 +1,140 @@
#include "CepGen/Generator.h"
#include "CepGen/Parameters.h"
#include "CepGen/Core/Integrator.h"
#include "CepGen/Core/Exception.h"
+#include "CepGen/Processes/GenericProcess.h"
#include "CepGen/Event/Event.h"
#include <fstream>
namespace CepGen
{
Generator::Generator() :
parameters( std::unique_ptr<Parameters>( new Parameters ) ),
cross_section_( -1. ), cross_section_error_( -1. ), has_cross_section_( false )
{
if ( Logger::get().level > Logger::Nothing ) {
Debugging( "Generator initialized" );
try { printHeader(); } catch ( Exception& e ) { e.dump(); }
}
srand( time( 0 ) ); // Random number initialization
}
Generator::Generator( Parameters* ip ) :
parameters( ip )
{}
Generator::~Generator()
{
if ( parameters->generation.enabled && parameters->process() && parameters->process()->numGeneratedEvents()>0 ) {
Information( Form( "Mean generation time / event: %.3f ms", parameters->process()->totalGenerationTime()*1.e3/parameters->process()->numGeneratedEvents() ) );
}
}
size_t
Generator::numDimensions() const
{
if ( !parameters->process() ) return 0;
return parameters->process()->numDimensions( parameters->kinematics.mode );
}
void
Generator::clearRun()
{
integrator_.reset();
parameters->integrator.first_run = true;
has_cross_section_ = false; // force the recreation of the integrator instance
cross_section_ = cross_section_error_ = -1.;
}
void
Generator::setParameters( Parameters& ip )
{
parameters = std::unique_ptr<Parameters>( new Parameters( ip ) ); // copy constructor
}
void
Generator::printHeader()
{
std::string tmp;
std::ostringstream os; os << std::endl;
std::ifstream hf( "README" );
if ( !hf.good() ) throw Exception( __PRETTY_FUNCTION__, "Failed to open README file", JustWarning );
while ( true ) {
if ( !hf.good() ) break;
getline( hf, tmp );
os << "\n " << tmp;
}
hf.close();
Information( os.str().c_str() );
}
double
Generator::computePoint( double* x )
{
prepareFunction();
double res = f( x, numDimensions(), (void*)parameters.get() );
std::ostringstream os;
for ( unsigned int i=0; i<numDimensions(); i++ ) { os << x[i] << " "; }
Debugging( Form( "Result for x[%zu] = ( %s):\n\t%10.6f", numDimensions(), os.str().c_str(), res ) );
return res;
}
void
Generator::computeXsection( double& xsec, double& err )
{
// first destroy and recreate the integrator instance
if ( !integrator_ ) {
integrator_ = std::unique_ptr<Integrator>( new Integrator( numDimensions(), f, parameters.get() ) );
}
else if ( integrator_->dimensions() != numDimensions() ) {
integrator_.reset( new Integrator( numDimensions(), f, parameters.get() ) );
}
if ( Logger::get().level>=Logger::Debug ) {
std::ostringstream topo; topo << parameters->kinematics.mode;
Debugging( Form( "New integrator instance created\n\t"
"Considered topology: %s case\n\t"
"Will proceed with %d-dimensional integration", topo.str().c_str(), numDimensions() ) );
}
Information( "Starting the computation of the process cross-section" );
try { prepareFunction(); } catch ( Exception& e ) { e.dump(); }
has_cross_section_ = ( integrator_->integrate( cross_section_, cross_section_error_ ) == 0 );
xsec = cross_section_;
err = cross_section_error_;
Information( Form( "Total cross section: %f +/- %f pb", xsec, err ) );
}
Event*
Generator::generateOneEvent()
{
bool good = false;
if ( !has_cross_section_ ) {
computeXsection( cross_section_, cross_section_error_ );
}
while ( !good ) { good = integrator_->generateOneEvent(); }
last_event = this->parameters->generation.last_event;
return last_event.get();
}
void
Generator::prepareFunction()
{
if ( !parameters->process() ) {
throw Exception( __PRETTY_FUNCTION__, "No process defined!", FatalError );
}
Kinematics kin = parameters->kinematics;
parameters->process()->addEventContent();
parameters->process()->setKinematics( kin );
Debugging( "Function prepared to be integrated!" );
}
}
diff --git a/CepGen/Core/Integrand.cpp b/CepGen/Core/Integrand.cpp
index 05af43c..e618bc7 100644
--- a/CepGen/Core/Integrand.cpp
+++ b/CepGen/Core/Integrand.cpp
@@ -1,136 +1,137 @@
#include "CepGen/Core/Timer.h"
#include "CepGen/Core/Exception.h"
#include "CepGen/Event/Event.h"
#include "CepGen/Event/Particle.h"
#include "CepGen/Physics/Kinematics.h"
+#include "CepGen/Processes/GenericProcess.h"
#include "CepGen/Parameters.h"
#include <sstream>
namespace CepGen
{
double
f( double* x, size_t ndim, void* params )
{
std::ostringstream os;
Parameters* p = static_cast<Parameters*>( params );
std::shared_ptr<Event> ev = p->process()->event();
if ( p->process()->hasEvent() ) {
p->process()->clearEvent();
const Particle::Momentum p1( 0., 0., p->kinematics.inp.first ), p2( 0., 0., -p->kinematics.inp.second );
p->process()->setIncomingKinematics( p1, p2 ); // at some point introduce non head-on colliding beams?
DebuggingInsideLoop( Form( "Function f called -- some parameters:\n\t"
" pz(p1) = %5.2f pz(p2) = %5.2f\n\t"
" remnant mode: %d",
p->kinematics.inp.first, p->kinematics.inp.second, p->kinematics.structure_functions ) );
if ( p->integrator.first_run ) {
if ( Logger::get().level >= Logger::Debug ) {
std::ostringstream oss; oss << p->kinematics.mode;
Debugging( Form( "Process mode considered: %s", oss.str().c_str() ) );
}
Particle& op1 = ev->getOneByRole( Particle::OutgoingBeam1 ),
&op2 = ev->getOneByRole( Particle::OutgoingBeam2 );
//--- add outgoing protons or remnants
switch ( p->kinematics.mode ) {
case Kinematics::ElectronProton: default: { InError( Form( "Process mode %d not yet handled!", (int)p->kinematics.mode ) ); }
case Kinematics::ElasticElastic: break; // nothing to change in the event
case Kinematics::ElasticInelastic:
op2.setPdgId( Particle::uQuark, 1 ); break;
case Kinematics::InelasticElastic: // set one of the outgoing protons to be fragmented
op1.setPdgId( Particle::uQuark, 1 ); break;
case Kinematics::InelasticInelastic: // set both the outgoing protons to be fragmented
op1.setPdgId( Particle::uQuark, 1 );
op2.setPdgId( Particle::uQuark, 1 );
break;
}
//--- prepare the function to be integrated
p->process()->prepareKinematics();
//--- add central system
Particles& central_system = ev->getByRole( Particle::CentralSystem );
if ( central_system.size() == p->kinematics.central_system.size() ) {
unsigned short i = 0;
for ( Particles::iterator part = central_system.begin(); part != central_system.end(); ++part ) {
if ( p->kinematics.central_system[i] == Particle::invalidParticle ) continue;
part->setPdgId( p->kinematics.central_system[i] );
part->computeMass();
i++;
}
}
p->process()->clearRun();
p->integrator.first_run = false;
}
} // event is not empty
p->process()->setPoint( ndim, x );
if ( Logger::get().level >= Logger::DebugInsideLoop ) {
std::ostringstream oss; for ( unsigned int i=0; i<ndim; i++ ) { oss << x[i] << " "; }
DebuggingInsideLoop( Form( "Computing dim-%d point ( %s)", ndim, oss.str().c_str() ) );
}
//--- from this step on, the phase space point is supposed to be set by 'GenericProcess::setPoint()'
p->process()->beforeComputeWeight();
Timer tmr; // start the timer
double integrand = p->process()->computeWeight();
if ( integrand < 0. ) return 0.;
//--- only fill in the process' Event object if storage is requested
// or if taming functions are to be applied
if ( !p->taming_functions.empty() || p->storage() ) p->process()->fillKinematics();
//--- once the kinematics variables have been populated,
// can apply the collection of taming functions
double taming = 1.0;
if ( p->taming_functions.has( "m_central" ) || p->taming_functions.has( "pt_central" ) ) {
const Particle::Momentum central_system( ev->getByRole( Particle::CentralSystem )[0].momentum() + ev->getByRole( Particle::CentralSystem )[1].momentum() );
taming *= p->taming_functions.eval( "m_central", central_system.mass() );
taming *= p->taming_functions.eval( "pt_central", central_system.pt() );
}
if ( p->taming_functions.has( "q2" ) ) {
taming *= p->taming_functions.eval( "q2", ev->getOneByRole( Particle::Parton1 ).momentum().mass() );
taming *= p->taming_functions.eval( "q2", ev->getOneByRole( Particle::Parton2 ).momentum().mass() );
}
integrand *= taming;
//--- full event content (+ hadronisation) if generating events
if ( p->storage() ) {
ev->time_generation = tmr.elapsed();
ev->time_total = tmr.elapsed();
p->process()->addGenerationTime( ev->time_total );
Debugging( Form( "Generation time: %5.6f sec\n\t"
"Total time (gen+hadr): %5.6f sec",
ev->time_generation,
ev->time_total ) );
p->generation.last_event = ev;
} // generating events
if ( Logger::get().level>=Logger::DebugInsideLoop ) {
os.str( "" ); for ( unsigned int i=0; i<ndim; i++ ) { os << Form( "%10.8f ", x[i] ); }
Debugging( Form( "f value for dim-%d point ( %s): %4.4e", ndim, os.str().c_str(), integrand ) );
}
return integrand;
}
}
diff --git a/CepGen/Core/Integrator.cpp b/CepGen/Core/Integrator.cpp
index 6605f5d..61fc5aa 100644
--- a/CepGen/Core/Integrator.cpp
+++ b/CepGen/Core/Integrator.cpp
@@ -1,329 +1,330 @@
#include "Integrator.h"
#include "CepGen/Parameters.h"
#include "CepGen/Core/Exception.h"
+#include "CepGen/Event/Event.h"
#include <fstream>
namespace CepGen
{
Integrator::Integrator( const unsigned int dim, double f_( double*, size_t, void* ), Parameters* param ) :
ps_bin_( 0 ), input_params_( param ), function_( std::unique_ptr<gsl_monte_function>( new gsl_monte_function ) )
{
//--- function to be integrated
function_->f = f_;
function_->dim = dim;
function_->params = (void*)param;
//--- initialise the random number generator
gsl_rng_env_setup();
rng_ = gsl_rng_alloc( gsl_rng_default );
unsigned long seed = ( param->integrator.seed > 0 )
? param->integrator.seed
: time( nullptr ); // seed with time
gsl_rng_set( rng_, seed );
Debugging( Form( "Number of integration dimensions: %d\n\t"
"Number of iterations: %d\n\t"
"Number of function calls: %d", dim, input_params_->integrator.itvg, input_params_->integrator.ncvg ) );
}
Integrator::~Integrator()
{
if ( rng_ ) gsl_rng_free( rng_ );
}
int
Integrator::integrate( double& result, double& abserr )
{
int res = -1;
gsl_monte_vegas_state* veg_state;
gsl_monte_miser_state* mis_state;
const Integrator::Type algorithm = input_params_->integrator.type;
//--- integration bounds
std::vector<double> x_low( function_->dim, 0. ), x_up( function_->dim, 1. );
//--- prepare integrator
if ( algorithm == Vegas ) veg_state = gsl_monte_vegas_alloc( function_->dim );
else if ( algorithm == MISER ) mis_state = gsl_monte_miser_alloc( function_->dim );
if ( algorithm == Vegas ) {
//----- Vegas warmup (prepare the grid)
if ( !grid_.grid_prepared ) {
res = gsl_monte_vegas_integrate( function_.get(), &x_low[0], &x_up[0], function_->dim, 10000, rng_, veg_state, &result, &abserr );
grid_.grid_prepared = true;
}
//----- integration
for ( unsigned short i = 0; i < input_params_->integrator.itvg; i++ ) {
res = gsl_monte_vegas_integrate( function_.get(), &x_low[0], &x_up[0], function_->dim, 0.2 * input_params_->integrator.ncvg, rng_, veg_state, &result, &abserr );
PrintMessage( Form( ">> Iteration %2d: average = %10.6f sigma = %10.6f chi2 = %4.3f", i+1, result, abserr, gsl_monte_vegas_chisq( veg_state ) ) );
}
}
//----- integration
else if ( algorithm == MISER )
res = gsl_monte_miser_integrate( function_.get(), &x_low[0], &x_up[0], function_->dim, input_params_->integrator.ncvg, rng_, mis_state, &result, &abserr );
//--- clean integrator
if ( algorithm == Vegas ) gsl_monte_vegas_free( veg_state );
else if ( algorithm == MISER ) gsl_monte_miser_free( mis_state );
return res;
}
void
Integrator::generate()
{
std::ofstream of;
std::string fn;
if ( !grid_.gen_prepared ) setGen();
Information( Form( "%d events will be generated", input_params_->generation.maxgen ) );
unsigned int i = 0;
while ( i < input_params_->generation.maxgen ) {
if ( generateOneEvent() ) i++;
}
Information( Form( "%d events generated", i ) );
}
bool
Integrator::generateOneEvent()
{
if ( !grid_.gen_prepared ) setGen();
const unsigned int ndim = function_->dim, max = pow( mbin_, ndim );
std::vector<double> x( ndim, 0. );
//--- correction cycles
if ( ps_bin_ != 0 ) {
bool has_correction = false;
while ( !correctionCycle( x, has_correction ) ) {}
if ( has_correction ) return storeEvent( x );
}
double weight;
double y = -1.;
//--- normal generation cycle
//----- select a Integrator bin and reject if fmax is too small
do {
do {
// ...
ps_bin_ = uniform() * max;
y = uniform() * grid_.f_max_global;
grid_.nm[ps_bin_] += 1;
} while ( y > grid_.f_max[ps_bin_] );
// Select x values in this Integrator bin
int jj = ps_bin_;
for ( unsigned int i=0; i<ndim; i++ ) {
int jjj = jj / mbin_;
grid_.n[i] = jj - jjj * mbin_;
x[i] = ( uniform() + grid_.n[i] ) / mbin_;
jj = jjj;
}
// Get weight for selected x value
weight = F( x );
} while ( y > weight );
if ( weight <= grid_.f_max[ps_bin_] ) ps_bin_ = 0;
// Init correction cycle if weight is higher than fmax or ffmax
else if ( weight <= grid_.f_max_global ) {
grid_.f_max_old = grid_.f_max[ps_bin_];
grid_.f_max[ps_bin_] = weight;
grid_.f_max_diff = weight-grid_.f_max_old;
grid_.correc = ( grid_.nm[ps_bin_] - 1. ) * grid_.f_max_diff / grid_.f_max_global - 1.;
}
else {
grid_.f_max_old = grid_.f_max[ps_bin_];
grid_.f_max[ps_bin_] = weight;
grid_.f_max_diff = weight-grid_.f_max_old;
grid_.f_max_global = weight;
grid_.correc = ( grid_.nm[ps_bin_] - 1. ) * grid_.f_max_diff / grid_.f_max_global * weight / grid_.f_max_global - 1.;
}
Debugging( Form( "Correction applied: %f, phase space bin = %d", grid_.correc, ps_bin_ ) );
// Return with an accepted event
if ( weight > 0. ) return storeEvent( x );
return false;
}
bool
Integrator::correctionCycle( std::vector<double>& x, bool& has_correction )
{
double weight;
const unsigned int ndim = function_->dim;
Debugging( Form( "Correction cycles are started.\n\t"
"j = %f"
"correc = %f"
"corre2 = %f", ps_bin_, grid_.correc2 ) );
if ( grid_.correc >= 1. ) grid_.correc -= 1.;
if ( uniform() < grid_.correc ) {
grid_.correc = -1.;
std::vector<double> xtmp( ndim, 0. );
// Select x values in phase space bin
for ( unsigned int k=0; k<ndim; k++ ) {
xtmp[k] = ( uniform() + grid_.n[k] ) * inv_mbin_;
}
// Compute weight for x value
weight = F( xtmp );
// Parameter for correction of correction
if ( weight > grid_.f_max[ps_bin_] ) {
if ( weight > grid_.f_max2 ) grid_.f_max2 = weight;
grid_.correc2 -= 1.;
grid_.correc += 1.;
}
// Accept event
if ( weight >= grid_.f_max_diff*uniform() + grid_.f_max_old ) { // FIXME!!!!
//Error("Accepting event!!!");
//return storeEvent(x);
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_.f_max[ps_bin_] ) {
grid_.f_max_old = grid_.f_max[ps_bin_];
grid_.f_max[ps_bin_] = grid_.f_max2;
grid_.f_max_diff = grid_.f_max2-grid_.f_max_old;
if ( grid_.f_max2 < grid_.f_max_global ) {
grid_.correc = ( grid_.nm[ps_bin_] - 1. ) * grid_.f_max_diff / grid_.f_max_global - grid_.correc2;
}
else {
grid_.f_max_global = grid_.f_max2;
grid_.correc = ( grid_.nm[ps_bin_] - 1. ) * grid_.f_max_diff / grid_.f_max_global * grid_.f_max2 / grid_.f_max_global - grid_.correc2;
}
grid_.correc2 = 0.;
grid_.f_max2 = 0.;
return false;
}
return true;
}
bool
Integrator::storeEvent( const std::vector<double>& x )
{
input_params_->setStorage( true );
F( x );
input_params_->generation.ngen += 1;
input_params_->setStorage( false );
if ( input_params_->generation.ngen % input_params_->generation.gen_print_every == 0 ) {
Debugging( Form( "Generated events: %d", input_params_->generation.ngen ) );
input_params_->generation.last_event->dump();
}
return true;
}
void
Integrator::setGen()
{
Information( Form( "Preparing the grid for the generation of unweighted events: %d points", input_params_->integrator.npoints ) );
// Variables for debugging
std::ostringstream os;
if ( Logger::get().level >= Logger::Debug ) {
Debugging( Form( "MaxGen = %d", input_params_->generation.maxgen ) );
}
const unsigned int ndim = function_->dim,
max = pow( mbin_, ndim ),
npoin = input_params_->integrator.npoints;
const double inv_npoin = 1./npoin;
if ( ndim > max_dimensions_ ) {
FatalError( Form( "Number of dimensions to integrate exceed the maximum number, %d", max_dimensions_ ) );
}
grid_.nm = std::vector<int>( max, 0 );
grid_.f_max = std::vector<double>( max, 0. );
grid_.n = std::vector<int>( ndim, 0 );
std::vector<double> x( ndim, 0. );
input_params_->generation.ngen = 0;
// ...
double sum = 0., sum2 = 0., sum2p = 0.;
//--- main loop
for ( unsigned int i=0; i<max; i++ ) {
int jj = i;
for ( unsigned int j=0; j<ndim; j++ ) {
int jjj = jj*inv_mbin_;
grid_.n[j] = jj-jjj*mbin_;
jj = jjj;
}
double fsum = 0., fsum2 = 0.;
for ( unsigned int j=0; j<npoin; j++ ) {
for ( unsigned int k=0; k<ndim; k++ ) {
x[k] = ( uniform() + grid_.n[k] ) * inv_mbin_;
}
const double z = F( x );
grid_.f_max[i] = std::max( grid_.f_max[i], z );
fsum += z;
fsum2 += z*z;
}
const double av = fsum*inv_npoin, av2 = fsum2*inv_npoin, sig2 = av2 - av*av;
sum += av;
sum2 += av2;
sum2p += sig2;
grid_.f_max_global = std::max( grid_.f_max_global, grid_.f_max[i] );
if ( Logger::get().level >= Logger::DebugInsideLoop ) {
const double sig = sqrt( sig2 );
const double eff = ( grid_.f_max[i] != 0. ) ? grid_.f_max[i]/av : 1.e4;
os.str(""); for ( unsigned int j=0; j<ndim; j++ ) { os << grid_.n[j]; if ( j != ndim-1 ) os << ", "; }
DebuggingInsideLoop( Form( "In iteration #%d:\n\t"
"av = %f\n\t"
"sig = %f\n\t"
"fmax = %f\n\t"
"eff = %f\n\t"
"n = (%s)",
i, av, sig, grid_.f_max[i], eff, os.str().c_str() ) );
}
} // end of main loop
sum = sum/max;
sum2 = sum2/max;
sum2p = sum2p/max;
if ( Logger::get().level >= Logger::Debug ) {
const double sig = sqrt( sum2-sum*sum ), sigp = sqrt( sum2p );
double eff1 = 0.;
for ( unsigned int i=0; i<max; i++ ) eff1 += ( grid_.f_max[i] / ( max*sum ) );
const double eff2 = grid_.f_max_global/sum;
Debugging( Form( "Average function value = sum = %f\n\t"
"Average function value**2 = sum2 = %f\n\t"
"Overall standard deviation = sig = %f\n\t"
"Average standard deviation = sigp = %f\n\t"
"Maximum function value = ffmax = %f\n\t"
"Average inefficiency = eff1 = %f\n\t"
"Overall inefficiency = eff2 = %f\n\t",
sum, sum2, sig, sigp, grid_.f_max_global, eff1, eff2 ) );
}
grid_.gen_prepared = true;
Information( "Grid prepared! Now launching the production." );
}
std::ostream&
operator<<( std::ostream& os, const Integrator::Type& type )
{
switch ( type ) {
case Integrator::Vegas: return os << "Vegas";
case Integrator::MISER: return os << "MISER";
}
return os;
}
}
diff --git a/CepGen/Core/Parameters.cpp b/CepGen/Core/Parameters.cpp
index a4b2f55..bb2d2a6 100644
--- a/CepGen/Core/Parameters.cpp
+++ b/CepGen/Core/Parameters.cpp
@@ -1,114 +1,122 @@
#include "CepGen/Parameters.h"
#include "CepGen/Core/Exception.h"
+#include "CepGen/Processes/GenericProcess.h"
namespace CepGen
{
Parameters::Parameters() :
store_( false )
{}
Parameters::Parameters( Parameters& param ) :
kinematics( param.kinematics ), integrator( param.integrator ), generation( param.generation ),
taming_functions( param.taming_functions ),
process_( std::move( param.process_ ) ),
store_( param.store_ )
{}
Parameters::Parameters( const Parameters& param ) :
kinematics( param.kinematics ), integrator( param.integrator ), generation( param.generation ),
taming_functions( param.taming_functions ),
store_( param.store_ )
{}
Parameters::~Parameters()
{}
void
Parameters::setThetaRange( float thetamin, float thetamax )
{
kinematics.central_cuts[Cuts::eta_single].in( Particle::thetaToEta( thetamax ), Particle::thetaToEta( thetamin ) );
if ( Logger::get().level >= Logger::Debug ) {
std::ostringstream os; os << kinematics.central_cuts[Cuts::eta_single];
Debugging( Form( "eta in range: %s => theta(min) = %5.2f, theta(max) = %5.2f",
os.str().c_str(), thetamin, thetamax ) );
}
}
+ std::string
+ Parameters::processName() const
+ {
+ if ( process_ ) return process_->name();
+ return "no process";
+ }
+
void
Parameters::dump( std::ostream& out, bool pretty ) const
{
std::ostringstream os;
const int wb = 75, wt = 32;
os.str( "" );
os
<< "Parameters dump" << std::left
<< std::endl << std::endl
<< std::setfill('_') << std::setw( wb ) << "_/¯ RUN INFORMATION ¯\\_" << std::setfill( ' ' ) << std::endl
<< std::right << std::setw( wb ) << std::left << std::endl
<< std::setw( wt ) << "Process to generate";
if ( process_ ) {
os << ( pretty ? boldify( process_->name().c_str() ) : process_->name() ) << std::endl
<< std::setw( wt ) << "" << process_->description();
}
else
os << ( pretty ? boldify( "no process!" ) : "no process!" );
os
<< std::endl
<< std::setw( wt ) << "Events generation? " << ( pretty ? yesno( generation.enabled ) : std::to_string( generation.enabled ) ) << std::endl
<< std::setw( wt ) << "Number of events to generate" << ( pretty ? boldify( generation.maxgen ) : std::to_string( generation.maxgen ) ) << std::endl
<< std::setw( wt ) << "Verbosity level " << Logger::get().level << std::endl
<< std::endl
<< std::setfill( '-' ) << std::setw( wb+6 ) << ( pretty ? boldify( " Integration parameters " ) : "Integration parameters" ) << std::setfill( ' ' ) << std::endl
<< std::endl;
std::ostringstream int_algo; int_algo << integrator.type;
os
<< std::setw( wt ) << "Integration algorithm" << ( pretty ? boldify( int_algo.str().c_str() ) : int_algo.str() ) << std::endl
<< std::setw( wt ) << "Maximum number of iterations" << ( pretty ? boldify( integrator.itvg ) : std::to_string( integrator.itvg ) ) << std::endl
<< std::setw( wt ) << "Number of function calls" << integrator.ncvg << std::endl
<< std::setw( wt ) << "Number of points to try per bin" << integrator.npoints << std::endl
<< std::setw( wt ) << "Random number generator seed" << integrator.seed << std::endl
<< std::endl
<< std::setfill('_') << std::setw( wb ) << "_/¯ EVENTS KINEMATICS ¯\\_" << std::setfill( ' ' ) << std::endl
<< std::endl
<< std::setfill( '-' ) << std::setw( wb+6 ) << ( pretty ? boldify( " Incoming particles " ) : "Incoming particles" ) << std::setfill( ' ' ) << std::endl
<< std::endl;
std::ostringstream proc_mode; proc_mode << kinematics.mode;
std::ostringstream ip1, ip2, op; ip1 << kinematics.inpdg.first; ip2 << kinematics.inpdg.second;
for ( std::vector<Particle::ParticleCode>::const_iterator cp = kinematics.central_system.begin(); cp != kinematics.central_system.end(); ++cp )
op << ( cp != kinematics.central_system.begin() ? ", " : "" ) << *cp;
std::ostringstream q2range; q2range << kinematics.initial_cuts.at( Cuts::q2 );
os
<< std::setw( wt ) << "Subprocess mode" << ( pretty ? boldify( proc_mode.str().c_str() ) : proc_mode.str() ) << std::endl
<< std::setw( wt ) << "Incoming particles" << ( pretty ? boldify( ip1.str().c_str() ) : ip1.str() ) << ", " << ( pretty ? boldify( ip2.str().c_str() ) : ip2.str() ) << std::endl
<< std::setw( wt ) << "Momenta (GeV/c)" << kinematics.inp.first << ", " << kinematics.inp.second << std::endl
<< std::setw( wt ) << "Structure functions" << kinematics.structure_functions << std::endl
<< std::endl
<< std::setfill( '-' ) << std::setw( wb+6 ) << ( pretty ? boldify( " Incoming partons " ) : "Incoming partons" ) << std::setfill( ' ' ) << std::endl
<< std::endl
<< std::setw( wt ) << "Virtuality range" << ( pretty ? boldify( q2range.str().c_str() ) : q2range.str().c_str() ) << " GeV**2" << std::endl
<< std::endl
<< std::setfill( '-' ) << std::setw( wb+6 ) << ( pretty ? boldify( " Outgoing central system " ) : "Outgoing central system" ) << std::setfill( ' ' ) << std::endl
<< std::endl
<< std::setw( wt ) << "Central particles" << ( pretty ? boldify( op.str().c_str() ) : op.str() ) << std::endl;
for ( std::map<Cuts::Central, Kinematics::Limits>::const_iterator lim = kinematics.central_cuts.begin(); lim != kinematics.central_cuts.end(); ++lim ) {
os << std::setw( wt ) << lim->first << lim->second << std::endl;
}
/*std::ostringstream ptrange; ptrange << kinematics.pt_single_central;
std::ostringstream erange; erange << kinematics.e_single_central;
std::ostringstream etarange; etarange << kinematics.eta_single_central;
std::ostringstream mxrange; mxrange << kinematics.mass_remnants;
os
<< std::setw( wt ) << "Lepton(s)' pT range" << ( pretty ? boldify( ptrange.str().c_str() ) : ptrange.str().c_str() ) << " GeV/c" << std::endl
<< std::setw( wt ) << "Lepton(s)' energy range" << ( pretty ? boldify( erange.str().c_str() ) : erange.str().c_str() ) << " GeV" << std::endl
<< std::setw( wt ) << "Pseudorapidity range" << ( pretty ? boldify( etarange.str().c_str() ) : etarange.str().c_str() ) << std::endl
//<< std::setw( wt ) << "Polar angle theta in range [deg]" << "[" << std::setw(3) << mintheta << ", " << std::setw( 3 ) << maxtheta << "]" << std::endl
<< std::endl
<< std::setfill( '-' ) << std::setw( wb+6 ) << ( pretty ? boldify( " Outgoing remnants" ) : "Outgoing remnants" ) << std::endl
<< std::endl << std::setfill( ' ' );
os << std::setw( wt ) << "Mass range" << ( pretty ? boldify( mxrange.str().c_str() ) : mxrange.str().c_str() ) << " GeV/c**2" << std::endl;*/
if ( pretty ) { Information( os.str() ); }
else out << os.str();
}
}
diff --git a/CepGen/Parameters.h b/CepGen/Parameters.h
index 3b93449..aba23c6 100644
--- a/CepGen/Parameters.h
+++ b/CepGen/Parameters.h
@@ -1,102 +1,104 @@
#ifndef CepGen_Parameters_h
#define CepGen_Parameters_h
-#include "CepGen/Processes/GenericProcess.h"
#include "CepGen/Core/TamingFunction.h"
#include "CepGen/Core/Integrator.h"
+#include "CepGen/Physics/Kinematics.h"
+#include "CepGen/Processes/GenericProcess.h"
#include <memory>
namespace CepGen
{
- class Kinematics;
+ class Event;
/// List of parameters used to start and run the simulation job
- class Parameters {
+ class Parameters
+ {
public:
Parameters();
/// Copy constructor (transfers ownership to the process!)
Parameters( Parameters& );
/// Const copy constructor (all but the process)
Parameters( const Parameters& );
~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 console
void dump( std::ostream& os=Logger::get().outputStream, bool pretty=true ) const;
//----- process to compute
/// Process for which the cross-section will be computed and the events will be generated
Process::GenericProcess* process() { return process_.get(); }
/// Name of the process considered
- std::string processName() const { return ( process_ ) ? process_->name() : "no process"; }
+ std::string processName() const;
/// Set the process to study
void setProcess( Process::GenericProcess* proc ) { process_.reset( proc ); }
//----- events kinematics
/// Events kinematics for phase space definition
Kinematics kinematics;
//----- VEGAS
/// Collection of integrator parameters
struct IntegratorParameters
{
IntegratorParameters() : type( Integrator::Vegas ), ncvg( 100000 ), itvg( 10 ), npoints( 100 ), first_run( true ), seed( 0 ) {}
Integrator::Type type;
/// Number of function calls to be computed for each point
unsigned int ncvg; // ??
/// Number of iterations for the integration
unsigned int itvg;
/// Number of points to "shoot" in each integration bin by the algorithm
unsigned int npoints;
/// Is it the first time the integrator is run?
bool first_run;
/// Random number generator seed
unsigned long seed;
};
/// Integrator parameters
IntegratorParameters integrator;
//----- events generation
/// Collection of events generation parameters
struct Generation
{
Generation() : enabled( false ), maxgen( 0 ), symmetrise( false ), ngen( 0 ), gen_print_every( 1 ) {}
/// 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;
/// Pointer to the last event produced in this run
std::shared_ptr<Event> last_event;
/// Do we want the events to be symmetrised with respect to the \f$z\f$-axis ?
bool symmetrise;
/// 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;
};
/// 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_; }
//----- taming functions
/// Functionals to be used to account for rescattering corrections (implemented within the process)
TamingFunctionsCollection taming_functions;
private:
std::unique_ptr<Process::GenericProcess> process_;
bool store_;
};
}
#endif
diff --git a/test/test_physics_processes.cpp b/test/test_physics_processes.cpp
index c813961..ac61062 100644
--- a/test/test_physics_processes.cpp
+++ b/test/test_physics_processes.cpp
@@ -1,140 +1,147 @@
#include "CepGen/Generator.h"
#include "CepGen/Parameters.h"
#include "CepGen/Core/Timer.h"
#include "CepGen/Processes/GamGamLL.h"
#include "CepGen/Processes/PPtoLL.h"
#include "CepGen/StructureFunctions/StructureFunctions.h"
#include "abort.h"
#include <unordered_map>
#include <assert.h>
+using namespace std;
+
int
main( int argc, char* argv[] )
{
- typedef std::vector<std::pair<std::string,std::pair<double,double> > > KinematicsMap;
- typedef std::vector<std::pair<float, KinematicsMap> > ValuesAtCutMap;
+ typedef vector<pair<string,pair<double,double> > > KinematicsMap;
+ typedef vector<pair<float, KinematicsMap> > ValuesAtCutMap;
AbortHandler ctrl_c;
// values defined at pt(single lepton)>15 GeV, |eta(single lepton)|<2.5, mX<1000 GeV
// process -> { pt cut -> { kinematics -> ( sigma, delta(sigma) ) } }
- std::vector<std::pair<std::string,ValuesAtCutMap> > values_map = {
+ vector<pair<string,ValuesAtCutMap> > values_map = {
//--- LPAIR values at sqrt(s) = 13 TeV
{ "lpair", {
{ 3.0, { // pt cut
{ "elastic", { 2.0871703e1, 3.542e-2 } },
{ "singlediss", { 1.5042536e1, 3.256e-2 } },
{ "doublediss", { 1.38835e1, 4.03018e-2 } }
} },
{ 15.0, { // pt cut
{ "elastic", { 4.1994803e-1, 8.328e-4 } },
{ "singlediss", { 4.8504819e-1, 1.171e-3 } },
{ "doublediss", { 6.35650e-1, 1.93968e-3 } }
} },
} },
//--- PPTOLL values
{ "pptoll", {
{ 3.0, { // pt cut
{ "elastic", { 2.0936541e1, 1.4096e-2 } },
{ "singlediss_su", { 1.4844881e1, 2.0723e-2 } }, // SU, qt<50
{ "doublediss_su", { 1.3580637e1, 2.2497e-2 } }, // SU, qt<50
} },
{ 15.0, { // pt cut
{ "elastic", { 4.2515888e-1, 3.0351e-4 } },
{ "singlediss_su", { 4.4903253e-1, 5.8970e-4 } }, // SU, qt<50
{ "doublediss_su", { 5.1923819e-1, 9.6549e-4 } }, // SU, qt<50
/*{ "2_singlediss", { 4.6710287e-1, 6.4842e-4 } }, // SU, qt<500
{ "3_doublediss", { 5.6316353e-1, 1.1829e-3 } }, // SU, qt<500*/
} },
} },
};
const double num_sigma = 3.0;
- if ( argc == 1 || strcmp( argv[1], "debug" ) != 0 ) {
+ if ( argc < 3 || strcmp( argv[2], "debug" ) != 0 ) {
CepGen::Logger::get().level = CepGen::Logger::Nothing;
}
Timer tmr;
CepGen::Generator mg;
+ if ( argc > 1 && strcmp( argv[1], "vegas" ) == 0 ) mg.parameters->integrator.type = CepGen::Integrator::Vegas;
+ if ( argc > 1 && strcmp( argv[1], "miser" ) == 0 ) mg.parameters->integrator.type = CepGen::Integrator::MISER;
+ //{ ostringstream os; os << mg.parameters->integrator.type; Information( Form( "Testing with %s", os.str().c_str() ) ); }
+ { cout << "Testing with " << mg.parameters->integrator.type << endl; }
+
mg.parameters->kinematics.setSqrtS( 13.e3 );
mg.parameters->kinematics.central_cuts[CepGen::Cuts::eta_single].in( -2.5, 2.5 );
mg.parameters->kinematics.remnant_cuts[CepGen::Cuts::mass].max() = 1000.;
//mg.parameters->integrator.ncvg = 50000;
mg.parameters->integrator.itvg = 5;
Information( Form( "Initial configuration time: %.3f ms", tmr.elapsed()*1.e3 ) );
tmr.reset();
unsigned short num_tests = 0, num_tests_passed = 0;
- std::vector<std::string> failed_tests, passed_tests;
+ vector<string> failed_tests, passed_tests;
try {
for ( const auto& values_vs_generator : values_map ) { // loop over all generators
if ( values_vs_generator.first == "lpair" ) mg.parameters->setProcess( new CepGen::Process::GamGamLL );
else if ( values_vs_generator.first == "pptoll" ) {
mg.parameters->setProcess( new CepGen::Process::PPtoLL );
mg.parameters->kinematics.initial_cuts[CepGen::Cuts::qt].max() = 50.0;
}
else { InError( Form( "Unrecognized generator mode: %s", values_vs_generator.first.c_str() ) ); break; }
for ( const auto& values_vs_cut : values_vs_generator.second ) { // loop over the single lepton pT cut
mg.parameters->kinematics.central_cuts[CepGen::Cuts::pt_single].min() = values_vs_cut.first;
for ( const auto& values_vs_kin : values_vs_cut.second ) { // loop over all possible kinematics
- if ( values_vs_kin.first.find( "elastic" ) != std::string::npos ) mg.parameters->kinematics.mode = CepGen::Kinematics::ElasticElastic;
- else if ( values_vs_kin.first.find( "singlediss" ) != std::string::npos ) mg.parameters->kinematics.mode = CepGen::Kinematics::InelasticElastic;
- else if ( values_vs_kin.first.find( "doublediss" ) != std::string::npos ) mg.parameters->kinematics.mode = CepGen::Kinematics::InelasticInelastic;
+ if ( values_vs_kin.first.find( "elastic" ) != string::npos ) mg.parameters->kinematics.mode = CepGen::Kinematics::ElasticElastic;
+ else if ( values_vs_kin.first.find( "singlediss" ) != string::npos ) mg.parameters->kinematics.mode = CepGen::Kinematics::InelasticElastic;
+ else if ( values_vs_kin.first.find( "doublediss" ) != string::npos ) mg.parameters->kinematics.mode = CepGen::Kinematics::InelasticInelastic;
else { InError( Form( "Unrecognized kinematics mode: %s", values_vs_kin.first.c_str() ) ); break; }
- if ( values_vs_kin.first.find( "_su" ) != std::string::npos ) mg.parameters->kinematics.structure_functions = CepGen::StructureFunctions::SzczurekUleshchenko;
+ if ( values_vs_kin.first.find( "_su" ) != string::npos ) mg.parameters->kinematics.structure_functions = CepGen::StructureFunctions::SzczurekUleshchenko;
else mg.parameters->kinematics.structure_functions = CepGen::StructureFunctions::SuriYennie;
Information( Form( "Process: %s/%s\n\tConfiguration time: %.3f ms", values_vs_generator.first.c_str(), values_vs_kin.first.c_str(), tmr.elapsed()*1.e3 ) );
tmr.reset();
mg.clearRun();
const double xsec_ref = values_vs_kin.second.first, err_xsec_ref = values_vs_kin.second.second;
double xsec_cepgen, err_xsec_cepgen;
mg.computeXsection( xsec_cepgen, err_xsec_cepgen );
const double sigma = ( fabs( xsec_ref-xsec_cepgen ) ) / sqrt( err_xsec_cepgen*err_xsec_cepgen + err_xsec_ref*err_xsec_ref );
Information( Form( "Computed cross section:\n\tRef. = %.3e +/- %.3e\n\tCepGen = %.3e +/- %.3e\n\tPull: %.6f", xsec_ref, err_xsec_ref, xsec_cepgen, err_xsec_cepgen, sigma ) );
Information( Form( "Computation time: %.3f ms", tmr.elapsed()*1.e3 ) );
tmr.reset();
- std::string test_res = values_vs_generator.first+":"+
+ string test_res = values_vs_generator.first+":"+
Form( "pt-gt-%.1f", values_vs_cut.first )+":"+
values_vs_kin.first+":"
"ref="+Form( "%.3e", xsec_ref )+":"
"got="+Form( "%.3e", xsec_cepgen )+":"
"pull="+Form( "%.3f", sigma );
if ( fabs( sigma )<num_sigma ) {
passed_tests.emplace_back( test_res );
num_tests_passed++;
}
else {
failed_tests.emplace_back( test_res );
}
num_tests++;
- std::cout << "Test " << num_tests_passed << "/" << num_tests << " passed!" << std::endl;
+ cout << "Test " << num_tests_passed << "/" << num_tests << " passed!" << endl;
}
}
}
} catch ( CepGen::Exception& e ) {}
if ( failed_tests.size() != 0 ) {
- std::ostringstream os_failed, os_passed;
- for ( const auto& fail : failed_tests ) os_failed << " (*) " << fail << std::endl;
- for ( const auto& pass : passed_tests ) os_passed << " (*) " << pass << std::endl;
+ ostringstream os_failed, os_passed;
+ for ( const auto& fail : failed_tests ) os_failed << " (*) " << fail << endl;
+ for ( const auto& pass : passed_tests ) os_passed << " (*) " << pass << endl;
throw CepGen::Exception( __PRETTY_FUNCTION__, Form( "Some tests failed!\n%s\nPassed tests:\n%s", os_failed.str().c_str(), os_passed.str().c_str() ), CepGen::FatalError );
}
else Information( "ALL TESTS PASSED!" );
return 0;
}