diff --git a/CepGen/Core/Integrator.cpp b/CepGen/Core/Integrator.cpp index f267ae2..ce7a162 100644 --- a/CepGen/Core/Integrator.cpp +++ b/CepGen/Core/Integrator.cpp @@ -1,475 +1,477 @@ #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/IO/GenericExportHandler.h" #include "CepGen/Event/Event.h" #include #include #include #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 x_low( function_->dim, 0. ), x_up( function_->dim, 1. ); //--- launch integration switch ( input_params_.integration().type ) { case IntegratorType::plain: { std::unique_ptr 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 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 ( input_params_.outputModule() ) + input_params_.outputModule()->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& x_low, std::vector& 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 callback ) { if ( !grid_->gen_prepared ) computeGenerationParameters(); std::vector 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 callback ) { if ( num_events < 1 ) num_events = input_params_.generation().maxgen; if ( input_params_.outputModule() ) input_params_.outputModule()->initialise( input_params_ ); try { while ( input_params_.numGeneratedEvents() < num_events ) generateOne( callback ); } catch ( const Exception& ) { return; } } bool Integrator::correctionCycle( std::vector& 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 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& x, std::function 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 ); if ( input_params_.outputModule() ) *input_params_.outputModule() << last_event; } 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 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 x( function_->dim, 0. ); std::vector n( function_->dim, 0 );; // ... double sum = 0., sum2 = 0., sum2p = 0.; 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 = fsum*inv_num_points, av2 = fsum2*inv_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& 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 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 w*function_->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/IO/TextHandler.cpp b/CepGen/IO/TextHandler.cpp index 1cbc9f6..64eaabe 100644 --- a/CepGen/IO/TextHandler.cpp +++ b/CepGen/IO/TextHandler.cpp @@ -1,245 +1,267 @@ #include "CepGen/IO/ExportHandler.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/ParametersList.h" #include "CepGen/Core/utils.h" #include "CepGen/Event/Event.h" #include "CepGen/Parameters.h" #include "CepGen/Version.h" +#include + +#include #include #include -#include - namespace cepgen { namespace io { /** * \brief Handler for the generic text file output * \author Laurent Forthomme * \date Jul 2019 */ class TextHandler : public GenericExportHandler { public: explicit TextHandler( const ParametersList& ); ~TextHandler(); void initialise( const Parameters& ) override; void setCrossSection( double xsec, double ) override { xsec_ = xsec; } void operator<<( const Event& ) override; private: /// Retrieve a named variable from a particle double variable( const Particle&, const std::string& ) const; /// Retrieve a named variable from the whole event double variable( const Event&, const std::string& ) const; static const std::regex rgx_select_id_, rgx_select_role_; static constexpr double INVALID_OUTPUT = -999.; + static constexpr size_t PLOT_WIDTH = 50; std::ofstream file_; const std::vector variables_; const bool print_banner_, print_variables_; const ParametersList hist_variables_; const std::string separator_; //--- variables definition + std::unordered_map variables_name_; typedef std::pair IndexedVariable; std::unordered_map > variables_per_id_; std::unordered_map > variables_per_role_; std::vector variables_for_event_; unsigned short num_vars_; std::ostringstream oss_vars_; double xsec_; //--- auxiliary helper maps const std::unordered_map role_str_ = { { "ib1", Particle::Role::IncomingBeam1 }, { "ib2", Particle::Role::IncomingBeam2 }, { "ob1", Particle::Role::OutgoingBeam1 }, { "ob2", Particle::Role::OutgoingBeam2 }, { "pa1", Particle::Role::Parton1 }, { "pa2", Particle::Role::Parton2 }, { "cs", Particle::Role::CentralSystem }, { "int", Particle::Role::Intermediate } }; typedef double( Particle::Momentum::*pMethod )(void) const; /// Mapping of string variables to momentum getter methods const std::unordered_map m_mom_str_ = { { "px", &Particle::Momentum::px }, { "py", &Particle::Momentum::py }, { "pz", &Particle::Momentum::pz }, { "pt", &Particle::Momentum::pt }, { "eta", &Particle::Momentum::eta }, { "phi", &Particle::Momentum::phi }, { "m", &Particle::Momentum::mass }, { "e", &Particle::Momentum::energy }, { "p", &Particle::Momentum::p }, { "pt2", &Particle::Momentum::pt2 }, { "th", &Particle::Momentum::theta }, { "y", &Particle::Momentum::rapidity } }; //--- kinematic variables double sqrts_; unsigned long num_evts_; struct gsl_histogram_deleter { void operator()( gsl_histogram* h ) { gsl_histogram_free( h ); } }; std::unordered_map > hists_; }; const std::regex TextHandler::rgx_select_id_( "(\\w+)\\((\\d+)\\)" ); const std::regex TextHandler::rgx_select_role_( "(\\w+)\\(([a-z]+\\d?)\\)" ); TextHandler::TextHandler( const ParametersList& params ) : GenericExportHandler( "text" ), file_ ( params.get( "filename", "output.txt" ) ), variables_ ( params.get >( "variables" ) ), print_banner_ ( params.get( "saveBanner", true ) ), print_variables_( params.get( "saveVariables", true ) ), hist_variables_ ( params.get( "histVariables" ) ), separator_ ( params.get( "separator", "\t" ) ), num_vars_( 0 ), xsec_( 1. ) { std::smatch sm; oss_vars_.clear(); std::string sep; for ( const auto& var : variables_ ) { if ( std::regex_match( var, sm, rgx_select_id_ ) ) variables_per_id_[stod( sm[2].str() )].emplace_back( std::make_pair( num_vars_, sm[1].str() ) ); else if ( std::regex_match( var, sm, rgx_select_role_ ) ) { const auto& str_role = sm[2].str(); if ( role_str_.count( str_role ) == 0 ) { CG_WARNING( "TextHandler" ) << "Invalid particle role retrieved from configuration: \"" << str_role << "\".\n\t" << "Skipping the variable \"" << var << "\" in the output module."; continue; } variables_per_role_[role_str_.at( str_role )].emplace_back( std::make_pair( num_vars_, sm[1].str() ) ); } else // event-level variables variables_for_event_.emplace_back( std::make_pair( num_vars_, var ) ); oss_vars_ << sep << var, sep = separator_; if ( hist_variables_.has( var ) ) { const auto& hvar = hist_variables_.get( var ); const int nbins = hvar.get( "nbins" ); const double min = hvar.get( "low" ), max = hvar.get( "high" ); hists_[num_vars_].reset( gsl_histogram_alloc( nbins ) ); gsl_histogram_set_ranges_uniform( hists_[num_vars_].get(), min, max ); CG_INFO( "TextHandler" ) << "Booking a histogram with " << nbins << " bin" << utils::s( nbins ) << " between " << min << " and " << max << " for \"" << var << "\"."; } + variables_name_[num_vars_] = var; ++num_vars_; } } TextHandler::~TextHandler() { for ( const auto& hs : hists_ ) { const auto& hist = hs.second.get(); gsl_histogram_scale( hist, xsec_/( num_evts_+1 ) ); - gsl_histogram_fprintf( stdout, hist, "%g", "%g" ); + const double inv_max_bin = 1./gsl_histogram_max_val( hist ); + CG_INFO( "TextHandler" ) + << "plot of \"" << variables_name_.at( hs.first ) << "\"\n" + << std::string( 11, ' ' ) + << Form( "%-5.2f", gsl_histogram_min_val( hist ) ) + << std::string( PLOT_WIDTH-12, ' ' ) + << Form( "%5.2f", gsl_histogram_max_val( hist ) ) << " pb\n" + << std::string( 11, ' ' ) + << std::string( PLOT_WIDTH+1, '.' ); + for ( size_t i = 0; i < gsl_histogram_bins( hist ); ++i ) { + double min, max; + gsl_histogram_get_range( hist, i, &min, &max ); + const int val = gsl_histogram_get( hist, i )*PLOT_WIDTH*inv_max_bin; + CG_LOG( "TextHandler" ) << "[" + << std::setw( 4 ) << std::setprecision( 4 ) << min << "," + << std::setw( 4 ) << std::setprecision( 4 ) << max << "):" + << std::string( val, '*' ); + } + //gsl_histogram_fprintf( stdout, hist, "%g", "%g" ); } file_.close(); } void TextHandler::initialise( const Parameters& params ) { sqrts_ = params.kinematics.sqrtS(); num_evts_ = 0ul; if ( print_banner_ ) file_ << banner( params, "#" ) << "\n"; if ( print_variables_ ) file_ << "# " << oss_vars_.str() << "\n"; } void TextHandler::operator<<( const Event& ev ) { std::vector vars( num_vars_ ); //--- extract and order the variables to be retrieved //--- particle-level variables (indexed by integer id) for ( const auto& id_vars : variables_per_id_ ) { const auto& part = ev[id_vars.first]; //--- loop over the list of variables for this particle for ( const auto& var : id_vars.second ) vars[var.first] = variable( part, var.second ); } //--- particle-level variables (indexed by role) for ( const auto& role_vars : variables_per_role_ ) { const auto& part = ev[role_vars.first][0]; //--- loop over the list of variables for this particle for ( const auto& var : role_vars.second ) vars[var.first] = variable( part, var.second ); } //--- event-level variables for ( const auto& var : variables_for_event_ ) vars[var.first] = variable( ev, var.second ); //--- write down the variables list in the file std::string sep; unsigned short i = 0; for ( const auto& var : vars ) { file_ << sep << var, sep = separator_; if ( hists_.count( i ) > 0 ) gsl_histogram_increment( hists_.at( i ).get(), var ); ++i; } file_ << "\n"; ++num_evts_; } double TextHandler::variable( const Particle& part, const std::string& var ) const { if ( m_mom_str_.count( var ) ) { auto meth = m_mom_str_.at( var ); return ( part.momentum().*meth )(); } if ( var == "xi" ) return 1.-part.energy()*2./sqrts_; if ( var == "pdg" ) return (double)part.integerPdgId(); if ( var == "charge" ) return part.charge(); if ( var == "status" ) return (double)part.status(); CG_WARNING( "TextHandler" ) << "Failed to retrieve variable \"" << var << "\"."; return INVALID_OUTPUT; } double TextHandler::variable( const Event& ev, const std::string& var ) const { if ( var == "np" ) return (double)ev.size(); if ( var == "nev" ) return (double)num_evts_+1; if ( var == "nob1" || var == "nob2" ) { unsigned short out = 0.; for ( const auto& part : ev[ var == "nob1" ? Particle::Role::OutgoingBeam1 : Particle::Role::OutgoingBeam2 ] ) if ( (int)part.status() > 0 ) out++; return (double)out; } if ( var == "tgen" ) return ev.time_generation; if ( var == "ttot" ) return ev.time_total; CG_WARNING( "TextHandler" ) << "Failed to retrieve the event-level variable \"" << var << "\"."; return INVALID_OUTPUT; } } } REGISTER_IO_MODULE( text, TextHandler )