diff --git a/CepGen/IO/TextHandler.cpp b/CepGen/IO/TextHandler.cpp index f8dbea9..eb7df45 100644 --- a/CepGen/IO/TextHandler.cpp +++ b/CepGen/IO/TextHandler.cpp @@ -1,302 +1,302 @@ #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 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: short extractVariableProperties( const std::string& ); /// 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_; std::unordered_map variable_stored_; 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. ) { //--- first extract list of variables to store in output file oss_vars_.clear(); std::string sep; for ( const auto& var : variables_ ) { auto id = extractVariableProperties( var ); if ( id >= 0 ) { oss_vars_ << sep << var, sep = separator_; variable_stored_[id] = true; } } //--- then extract list of variables to be plotted in histogram for ( const auto& var : hist_variables_.keys() ) { auto id = extractVariableProperties( var ); if ( id < 0 ) continue; const auto& hvar = hist_variables_.get( var ); const int nbins = hvar.get( "nbins", 10 ); const double min = hvar.get( "low", 0. ), max = hvar.get( "high", 1. ); hists_[id].reset( gsl_histogram_alloc( nbins ) ); gsl_histogram_set_ranges_uniform( hists_[id].get(), min, max ); CG_INFO( "TextHandler" ) << "Booking a histogram with " << nbins << " bin" << utils::s( nbins ) << " between " << min << " and " << max << " for \"" << var << "\"."; } } TextHandler::~TextHandler() { //--- histograms printout for ( const auto& var : hist_variables_.keys() ) { const auto& vn = std::find_if( variables_name_.begin(), variables_name_.end(), [&var]( auto&& p ) { return p.second == var; } ); if ( vn == variables_name_.end() ) { CG_WARNING( "TextHandler" ) << "Failed to retrieve variable \"" << var << "\" for plotting."; continue; } const auto& hist = hists_.at( vn->first ).get(); const size_t nbins = gsl_histogram_bins( hist ); gsl_histogram_scale( hist, xsec_/( num_evts_+1 ) ); const double max_bin = gsl_histogram_max_val( hist ); const double inv_max_bin = max_bin > 0. ? 1./max_bin : 0.; CG_INFO( "TextHandler" ) << "plot of \"" << var << "\" (" << "bin width=" << ( gsl_histogram_max( hist )-gsl_histogram_min( hist ) )/nbins << ", " << "mean=" << gsl_histogram_mean( hist ) << ", " - << "st.err=" << gsl_histogram_sigma( hist ) << ")\n" + << "st.dev=" << gsl_histogram_sigma( hist ) << ")\n" << std::string( 15, ' ' ) << 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( 15, ' ' ) << std::string( PLOT_WIDTH+1, '.' ); for ( size_t i = 0; i < nbins; ++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" ) << Form( "[%6.2f,%6.2f):", min, max ) << std::string( val, '*' ); } } //--- finalisation of the output file 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 ) { if ( variable_stored_.count( i ) > 0 && variable_stored_.at( i ) ) 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; } short TextHandler::extractVariableProperties( const std::string& var ) { const auto& vn = std::find_if( variables_name_.begin(), variables_name_.end(), [&var]( auto&& p ) { return p.second == var; } ); if ( vn != variables_name_.end() ) return vn->first; std::smatch sm; if ( std::regex_match( var, sm, rgx_select_id_ ) ) variables_per_id_[std::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."; return -1; } 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 ) ); variables_name_[num_vars_] = var; return num_vars_++; } } } REGISTER_IO_MODULE( text, TextHandler )