diff --git a/Cards/lpair_cfg.py b/Cards/lpair_cfg.py index 423041d..92304ef 100644 --- a/Cards/lpair_cfg.py +++ b/Cards/lpair_cfg.py @@ -1,38 +1,38 @@ import Config.Core as cepgen from Config.Integration.vegas_cff import integrator #from Config.Hadronisation.pythia6_cff import pythia6 as hadroniser #from Config.Hadronisation.pythia8_cff import pythia8 as hadroniser from Config.PDG_cfi import PDG process = cepgen.Module('lpair', processParameters = cepgen.Parameters( mode = cepgen.ProcessMode.InelasticElastic, pair = PDG.muon, ), inKinematics = cepgen.Parameters( pz = (6500., 6500.), structureFunctions = cepgen.StructureFunctions.SuriYennie, #structureFunctions = cepgen.StructureFunctions.FioreBrasse, #structureFunctions = cepgen.StructureFunctions.LUXlike, ), outKinematics = cepgen.Parameters( pt = (25.,), energy = (0.,), eta = (-2.5, 2.5), mx = (1.07, 1000.), ), #tamingFunctions = cepgen.Parameters( # # example of a complex taming function # cepgen.Parameters(variable = "m_ll", expression = "(m_ll>80.) ? exp(-(m_ll-80)/10) : 1.0"), #), ) -output = cepgen.Module('text', variables = ['m(4)', 'pt(4)']) +output = cepgen.Module('text', variables = ['m(4)', 'pt(4)', 'xi(ob1)']) #--- let the user specify the run conditions from Config.generator_cff import generator generator = generator.clone( numEvents = 100000, printEvery = 10000, ) diff --git a/CepGen/CMakeLists.txt b/CepGen/CMakeLists.txt index b098853..a946120 100644 --- a/CepGen/CMakeLists.txt +++ b/CepGen/CMakeLists.txt @@ -1,110 +1,110 @@ file(GLOB core_sources Core/*.cpp *.cpp) file(GLOB phys_sources Physics/*.cpp) file(GLOB sf_sources StructureFunctions/*.cpp) -file(GLOB io_sources IO/GenericExportHandler.cpp IO/GenericTextHandler.cpp) +file(GLOB io_sources IO/GenericExportHandler.cpp IO/TextHandler.cpp) file(GLOB hadr_sources Hadronisers/GenericHadroniser.cpp) include_directories(${PROJECT_SOURCE_DIR}) #----- check the external dependencies for SFs set(GRV_PATH ${PROJECT_SOURCE_DIR}/External) file(GLOB grv_sources ${GRV_PATH}/grv_*.f) if(grv_sources) message(STATUS "GRV PDFset found in ${grv_sources}!") add_definitions(-DGRVPDF) list(APPEND sf_sources ${grv_sources}) else() message(STATUS "GRV PDFset not found. Will proceed without it") endif() #----- check the external dependencies for physics utilities if(alphas_sources) list(APPEND phys_sources ${alphas_sources}) endif() set(addons_libraries "") set(strf_libraries ${GSL_LIB} ${GSL_CBLAS_LIB}) #--- linking with LHAPDF if(LHAPDF) message(STATUS "LHAPDF found in ${LHAPDF}") list(APPEND strf_libraries ${LHAPDF}) include_directories(${LHAPDF_INCLUDE}) else() file(GLOB partonic_sf StructureFunctions/Partonic.cpp) list(REMOVE_ITEM sf_sources ${partonic_sf}) endif() #--- linking with Pythia 6 if(PYTHIA6) message(STATUS "Pythia 6 found in ${PYTHIA6}") list(APPEND hadr_sources "Hadronisers/Pythia6Hadroniser.cpp") list(APPEND addons_libraries ${PYTHIA6}) if(PYTHIA6DUMMY) message(STATUS "Pythia 6 addons found in ${PYTHIA6DUMMY}") list(APPEND addons_libraries ${PYTHIA6DUMMY}) endif() elseif(EXISTS $ENV{PYTHIA6_SRC}) file(GLOB pythia6_src $ENV{PYTHIA6_SRC}) message(STATUS "Pythia 6 source found in ${pythia6_src}") set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -Wno-tabs -Wno-maybe-uninitialized -Wno-integer-division -Wno-unused-variable -Wno-unused-dummy-argument") add_library(pythia6 SHARED ${pythia6_src}) list(APPEND hadr_sources "Hadronisers/Pythia6Hadroniser.cpp") list(APPEND addons_libraries pythia6) endif() #--- linking with Pythia 8 if(PYTHIA8) message(STATUS "Pythia 8 found in ${PYTHIA8}") message(STATUS "Pythia 8 will be used for LHEF output") set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-misleading-indentation") include_directories(${PYTHIA8_INCLUDE}) add_definitions(-DPYTHIA8) list(APPEND hadr_sources "Hadronisers/Pythia8Hadroniser.cpp") list(APPEND io_sources "IO/PythiaEventInterface.cpp") list(APPEND io_sources "IO/LHEFPythiaHandler.cpp") list(APPEND addons_libraries ${PYTHIA8} dl) endif() #--- linking with HepMC if(HEPMC_LIB) message(STATUS "HepMC found in ${HEPMC_LIB}") if(HEPMC_LIB MATCHES ".*HepMC3.?.so") message(STATUS "HepMC version 3 found") if(HEPMC_ROOT_LIB) message(STATUS "HepMC ROOT I/O library found") add_definitions(-DHEPMC3_ROOTIO) list(APPEND addons_libraries ${HEPMC_ROOT_LIB}) endif() add_definitions(-DHEPMC3) endif() list(APPEND io_sources "IO/HepMCHandler.cpp") if(NOT PYTHIA8) message(STATUS "HepMC will be used for LHEF output") list(APPEND io_sources "IO/LHEFHepMCHandler.cpp") endif() list(APPEND addons_libraries ${HEPMC_LIB}) include_directories(${HEPMC_INCLUDE}) endif() #----- build the objects add_library(CepGenCore SHARED ${core_sources} ${phys_sources} ${sf_sources}) target_link_libraries(CepGenCore ${CEPGEN_EXTERNAL_CORE_REQS}) target_link_libraries(CepGenCore ${strf_libraries}) add_library(CepGenAddOns SHARED ${io_sources} ${hadr_sources}) target_link_libraries(CepGenAddOns ${addons_libraries}) target_link_libraries(CepGenAddOns CepGenEvent) #----- installation rules install(TARGETS CepGenCore DESTINATION lib) install(TARGETS CepGenAddOns DESTINATION lib) diff --git a/CepGen/Core/Integrator.cpp b/CepGen/Core/Integrator.cpp index 5071230..e9110b1 100644 --- a/CepGen/Core/Integrator.cpp +++ b/CepGen/Core/Integrator.cpp @@ -1,469 +1,473 @@ #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 ( 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& e ) { return; } + } 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 ); } 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/GenericTextHandler.cpp b/CepGen/IO/TextHandler.cpp similarity index 78% rename from CepGen/IO/GenericTextHandler.cpp rename to CepGen/IO/TextHandler.cpp index 748bebc..3329d92 100644 --- a/CepGen/IO/GenericTextHandler.cpp +++ b/CepGen/IO/TextHandler.cpp @@ -1,137 +1,151 @@ #include "CepGen/IO/ExportHandler.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/ParametersList.h" #include "CepGen/Event/Event.h" +#include "CepGen/Parameters.h" #include #include namespace cepgen { namespace io { /** * \brief Handler for the generic text file output * \author Laurent Forthomme * \date Jul 2019 */ - class GenericTextHandler : public GenericExportHandler + class TextHandler : public GenericExportHandler { public: - explicit GenericTextHandler( const ParametersList& ); - ~GenericTextHandler(); + explicit TextHandler( const ParametersList& ); + ~TextHandler(); void initialise( const Parameters& ) override; void operator<<( const Event& ) override; private: /// Retrieve a named variable from a particle double variable( const Particle&, const std::string& ) const; static const std::regex rgx_select_id_, rgx_select_role_; + static constexpr double INVALID_OUTPUT = -999.; + std::ofstream file_; std::vector variables_; std::unordered_map > > variables_per_id_; std::unordered_map > > variables_per_role_; unsigned short num_vars_; + + //--- kinematic variables + double sqrts_; }; - const std::regex GenericTextHandler::rgx_select_id_( "(\\w)\\((\\d+)\\)" ); - const std::regex GenericTextHandler::rgx_select_role_( "(\\w)\\(([a-z]+\\d?)\\)" ); + const std::regex TextHandler::rgx_select_id_( "(\\w+)\\((\\d+)\\)" ); + const std::regex TextHandler::rgx_select_role_( "(\\w+)\\(([a-z]+\\d?)\\)" ); - GenericTextHandler::GenericTextHandler( const ParametersList& params ) : + TextHandler::TextHandler( const ParametersList& params ) : GenericExportHandler( "text" ), file_( params.get( "filename", "output.txt" ) ), variables_( params.get >( "variables" ) ), num_vars_( 0 ) { const auto vars_tmp = variables_; variables_.clear(); std::smatch sm; std::string sep; file_ << "# "; for ( const auto& var : vars_tmp ) { 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(); auto role = Particle::Role::UnknownRole; if ( str_role == "ib1" ) role = Particle::Role::IncomingBeam1; else if ( str_role == "ib2" ) role = Particle::Role::IncomingBeam2; else if ( str_role == "ob1" ) role = Particle::Role::OutgoingBeam1; else if ( str_role == "ob2" ) role = Particle::Role::OutgoingBeam2; else if ( str_role == "cs" ) role = Particle::Role::CentralSystem; else if ( str_role == "int" ) role = Particle::Role::Intermediate; else if ( str_role == "pa1" ) role = Particle::Role::Parton1; else if ( str_role == "pa2" ) role = Particle::Role::Parton2; else { - CG_WARNING( "GenericTextHandler" ) + 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].emplace_back( std::make_pair( num_vars_, sm[1].str() ) ); } else { - CG_WARNING( "GenericTextHandler" ) + CG_WARNING( "TextHandler" ) << "Generic variables retrieval not yet supported.\n\t" << "Skipping the variable \"" << var << "\" in the output module."; variables_.emplace_back( var ); } file_ << sep << var, sep = "\t"; ++num_vars_; } file_ << "\n"; } - GenericTextHandler::~GenericTextHandler() + TextHandler::~TextHandler() { file_.close(); } void - GenericTextHandler::initialise( const Parameters& params ) + TextHandler::initialise( const Parameters& params ) { + sqrts_ = params.kinematics.sqrtS(); } void - GenericTextHandler::operator<<( const Event& ev ) + TextHandler::operator<<( const Event& ev ) { std::vector vars( num_vars_ ); //--- extract and order the variables to be retrieved for ( const auto& id_vars : variables_per_id_ ) { + //--- first get the particle const auto& part = ev[id_vars.first]; + //--- then loop on the variables for ( const auto& var : id_vars.second ) vars[var.first] = variable( part, var.second ); } for ( const auto& role_vars : variables_per_role_ ) { + //--- first get the particle const auto& part = ev[role_vars.first][0]; + //--- then loop on the variables for ( const auto& var : role_vars.second ) vars[var.first] = variable( part, var.second ); } //--- write down the variables list in the file std::string separator; for ( const auto& var : vars ) file_ << separator << var, separator = "\t"; file_ << "\n"; } double - GenericTextHandler::variable( const Particle& part, const std::string& var ) const + TextHandler::variable( const Particle& part, const std::string& var ) const { if ( var == "px" ) return part.momentum().px(); else if ( var == "py" ) return part.momentum().py(); else if ( var == "pz" ) return part.momentum().pz(); else if ( var == "pt" ) return part.momentum().pt(); else if ( var == "m" ) return part.mass(); else if ( var == "e" ) return part.energy(); + else if ( var == "xi" ) return 1.-part.energy()*2./sqrts_; else if ( var == "eta" ) return part.momentum().eta(); else if ( var == "phi" ) return part.momentum().phi(); else if ( var == "status" ) return (double)part.status(); - return -999.; + CG_WARNING( "TextHandler" ) + << "Failed to retrieve variable \"" << var << "\"."; + return INVALID_OUTPUT; } } } -REGISTER_IO_MODULE( text, GenericTextHandler ) +REGISTER_IO_MODULE( text, TextHandler )