diff --git a/Cards/legacy/lpair.card b/Cards/legacy/lpair.card index 5c832a2..df55cef 100644 --- a/Cards/legacy/lpair.card +++ b/Cards/legacy/lpair.card @@ -1,22 +1,24 @@ PROC lpair MODE 1 IEND 3 NCVG 100000 ITVG 10 INPP 6500. PMOD 11 INPE 6500. #EMOD 2 PAIR 13 MCUT 2 PTCT 15. ECUT 0. NGEN 100000 ETMN -2.5 ETMX 2.5 MXMX 1000.0 #THMN 0. #THMX 180. XIMN 0.02 XIMX 0.15 +#OUTP lhef +OUTP text diff --git a/Cards/lpair_cfg.py b/Cards/lpair_cfg.py index f6e3da3..3688d02 100644 --- a/Cards/lpair_cfg.py +++ b/Cards/lpair_cfg.py @@ -1,36 +1,43 @@ import Config.Core as cepgen from Config.Integration.vegas_cff import integrator +from Config.PDG_cfi import PDG + +#--- example of a hadronisation algorithm steering #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.), ), + #--- example of a complex taming function definition #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"), #), ) -#--- let the user specify the run conditions +#--- example of an output module parameterisation +#output = cepgen.Module('text', variables = ['nev', 'm(4)', 'tgen'], histVariables={'m(4)': cepgen.Parameters(low=0., high=250.)}) +#output = cepgen.Module('lhef', filename='test.lhe') +#output = cepgen.Module('hepmc', filename='test.hepmc') + +#--- let the user specify the events generation parameters from Config.generator_cff import generator generator = generator.clone( numEvents = 100000, printEvery = 10000, ) diff --git a/Cards/pptoll_cfg.py b/Cards/pptoll_cfg.py index ede3f4f..43dd1c8 100644 --- a/Cards/pptoll_cfg.py +++ b/Cards/pptoll_cfg.py @@ -1,38 +1,38 @@ import Config.Core as cepgen import Config.ktProcess_cfi as kt 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, registerParticle #--- example of an auxiliary particles definition -#registerParticle(1000001, 'sd_l', mass=100., charge=1., fermion=True) +#registerParticle(1000001, 'sd_l', mass=100., charge=1., fermion=True) # right now, only fermionic coupling handled process = kt.process.clone('pptoll', processParameters = cepgen.Parameters( mode = cepgen.ProcessMode.InelasticElastic, pair = PDG.muon, #pair = PDG.up, #pair = PDG.sd_l, # whatever was defined above as "new" particle ), inKinematics = cepgen.Parameters( pz = (6500., 6500.), structureFunctions = cepgen.StructureFunctions.SuriYennie, #structureFunctions = cepgen.StructureFunctions.FioreBrasse, #pdgIds = (PDG.proton, PDG.electron), ), outKinematics = kt.process.outKinematics.clone( pt = (25.,), energy = (0.,), eta = (-2.5, 2.5), mx = (1.07, 1000.), #--- extra cuts on the p1t(l) and p2t(l) plane #ptdiff = (0., 2.5), #--- distance in rapidity between l^+ and l^- #dely = (4., 5.), ), ) #--- events generation from Config.generator_cff import generator generator.numEvents = 10000 diff --git a/CepGen/CMakeLists.txt b/CepGen/CMakeLists.txt index c38f1a4..278cbe7 100644 --- a/CepGen/CMakeLists.txt +++ b/CepGen/CMakeLists.txt @@ -1,110 +1,127 @@ 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) +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") + list(APPEND io_sources "IO/HepMCEventInterface.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() +#--- linking with ROOT/Delphes + +if(ROOT_FOUND) + message(STATUS "ROOT found in ${ROOT_LIBRARY_DIR}") + list(APPEND addons_libraries ${ROOT_LIBRARIES}) + list(APPEND io_sources "IO/ROOTHandler.cpp" "IO/ROOTTreeHandler.cpp") + include_directories(${ROOT_INCLUDE_DIRS}) + link_directories(${ROOT_LIBRARY_DIR}) + if(DELPHES) + message(STATUS "Delphes found in ${DELPHES}") + list(APPEND addons_libraries ${DELPHES}) + list(APPEND io_sources "IO/DelphesHandler.cpp") + include_directories(${DELPHES_INCLUDE} ${DELPHES_INCLUDE}/external) + endif() +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/Cards/LpairHandler.cpp b/CepGen/Cards/LpairHandler.cpp index 4e72012..ffb2c37 100644 --- a/CepGen/Cards/LpairHandler.cpp +++ b/CepGen/Cards/LpairHandler.cpp @@ -1,238 +1,251 @@ #include "CepGen/Cards/LpairHandler.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/ParametersList.h" #include "CepGen/Core/Integrator.h" #include "CepGen/Processes/ProcessesHandler.h" #include "CepGen/Hadronisers/HadronisersHandler.h" +#include "CepGen/IO/ExportHandler.h" #include "CepGen/StructureFunctions/StructureFunctions.h" #include "CepGen/Physics/GluonGrid.h" #include "CepGen/Physics/PDG.h" #include #include namespace cepgen { namespace card { const int LpairHandler::kInvalid = 99999; //----- specialization for LPAIR input cards LpairHandler::LpairHandler( const char* file ) : proc_params_( new ParametersList ), str_fun_( 11 ), sr_type_( 1 ), xi_min_( 0. ), xi_max_( 1. ), hi_1_( { 0, 0 } ), hi_2_( { 0, 0 } ) { std::ifstream f( file, std::fstream::in ); if ( !f.is_open() ) throw CG_FATAL( "LpairHandler" ) << "Failed to parse file \"" << file << "%s\"."; init(); //--- parse all fields std::unordered_map m_params; std::string key, value; std::ostringstream os; while ( f >> key >> value ) { if ( key[0] == '#' ) // FIXME need to ensure there is no extra space before! continue; setParameter( key, value ); m_params.insert( { key, value } ); - if ( getDescription( key ) != "null" ) - os << "\n>> " << key << " = " << std::setw( 15 ) << getParameter( key ) - << " (" << getDescription( key ) << ")"; + if ( description( key ) != "null" ) + os << "\n>> " << key << " = " << std::setw( 15 ) << parameter( key ) + << " (" << description( key ) << ")"; } f.close(); //--- parse the process name params_.setProcess( proc::ProcessesHandler::get().build( proc_name_, *proc_params_ ) ); const Limits lim_xi{ xi_min_, xi_max_ }; if ( lim_xi.valid() ) params_.kinematics.cuts.remnants.energy_single = ( lim_xi+(-1.) )*( -params_.kinematics.incoming_beams.first.pz ); //--- parse the structure functions code auto sf_params = ParametersList() .set( "id", str_fun_ ) .set( "sigmaRatio", ParametersList() .set( "id", sr_type_ ) ); const unsigned long kLHAPDFCodeDec = 10000000, kLHAPDFPartDec = 1000000; if ( str_fun_ / kLHAPDFCodeDec == 1 ) { // SF from parton const unsigned long icode = str_fun_ % kLHAPDFCodeDec; sf_params .set( "id", (int)strfun::Type::Partonic ) .set( "pdfId", icode % kLHAPDFPartDec ) .set( "mode", icode / kLHAPDFPartDec ); // 0, 1, 2 } else if ( str_fun_ == (int)strfun::Type::MSTWgrid ) sf_params .set( "gridPath", mstw_grid_path_ ); params_.kinematics.structure_functions = strfun::StructureFunctionsHandler::get().build( sf_params ); //--- parse the integration algorithm name if ( integr_type_ == "plain" ) params_.integration().type = IntegratorType::plain; else if ( integr_type_ == "Vegas" ) params_.integration().type = IntegratorType::Vegas; else if ( integr_type_ == "MISER" ) params_.integration().type = IntegratorType::MISER; else if ( integr_type_ != "" ) throw CG_FATAL( "LpairHandler" ) << "Unrecognized integrator type: " << integr_type_ << "!"; //--- parse the hadronisation algorithm name if ( !hadr_name_.empty() ) { params_.setHadroniser( cepgen::hadr::HadronisersHandler::get().build( hadr_name_, ParametersList() ) ); params_.hadroniser()->setParameters( params_ ); } + //--- parse the output module name + if ( !out_mod_name_.empty() ) { + ParametersList outm; + if ( !out_file_name_.empty() ) + outm.set( "filename", out_file_name_ ); + params_.setOutputModule( cepgen::io::ExportHandler::get().build( out_mod_name_, outm ) ); + } + if ( m_params.count( "IEND" ) ) setValue( "IEND", ( std::stoi( m_params["IEND"] ) > 1 ) ); if ( m_params.count( "KMRG" ) && !kmr_grid_path_.empty() ) kmr::GluonGrid::get( kmr_grid_path_.c_str() ); //--- check if we are dealing with heavy ions for incoming states HeavyIon hi1{ hi_1_.first, (Element)hi_1_.second }, hi2{ hi_2_.first, (Element)hi_2_.second }; if ( hi1 ) params_.kinematics.incoming_beams.first.pdg = hi1; if ( hi2 ) params_.kinematics.incoming_beams.second.pdg = hi2; CG_INFO( "LpairHandler" ) << "File '" << file << "' succesfully opened!\n\t" << "The following parameters are set:" << os.str(); } void LpairHandler::init() { //------------------------------------------------------------------------------------------- // Process/integration/hadronisation parameters //------------------------------------------------------------------------------------------- registerParameter( "PROC", "Process name to simulate", &proc_name_ ); registerParameter( "ITYP", "Integration algorithm", &integr_type_ ); registerParameter( "HADR", "Hadronisation algorithm", &hadr_name_ ); - registerParameter( "KMRG", "KMR grid interpolation path", &kmr_grid_path_ ); + registerParameter( "OUTP", "Output module", &out_mod_name_ ); + registerParameter( "OUTF", "Output file name", &out_file_name_ ); //------------------------------------------------------------------------------------------- // General parameters //------------------------------------------------------------------------------------------- registerParameter( "IEND", "Generation type", ¶ms_.generation().enabled ); registerParameter( "NTRT", "Smoothen the integrand", ¶ms_.generation().treat ); registerParameter( "DEBG", "Debugging verbosity", (int*)&utils::Logger::get().level ); registerParameter( "NCVG", "Number of function calls", (int*)¶ms_.integration().ncvg ); registerParameter( "ITVG", "Number of integration iterations", (int*)¶ms_.integration().vegas.iterations ); registerParameter( "SEED", "Random generator seed", (int*)¶ms_.integration().rng_seed ); registerParameter( "NTHR", "Number of threads to use for events generation", (int*)¶ms_.generation().num_threads ); registerParameter( "MODE", "Subprocess' mode", (int*)¶ms_.kinematics.mode ); registerParameter( "NCSG", "Number of points to probe", (int*)¶ms_.generation().num_points ); registerParameter( "NGEN", "Number of events to generate", (int*)¶ms_.generation().maxgen ); registerParameter( "NPRN", "Number of events before printout", (int*)¶ms_.generation().gen_print_every ); //------------------------------------------------------------------------------------------- // Process-specific parameters //------------------------------------------------------------------------------------------- registerParameter( "METH", "Computation method (kT-factorisation)", &proc_params_->operator[]( "method" ) ); registerParameter( "IPOL", "Polarisation states to consider", &proc_params_->operator[]( "polarisationStates" ) ); //------------------------------------------------------------------------------------------- // Process kinematics parameters //------------------------------------------------------------------------------------------- + registerParameter( "KMRG", "KMR grid interpolation path", &kmr_grid_path_ ); registerParameter( "MGRD", "MSTW grid interpolation path", &mstw_grid_path_ ); registerParameter( "PMOD", "Outgoing primary particles' mode", &str_fun_ ); registerParameter( "EMOD", "Outgoing primary particles' mode", &str_fun_ ); registerParameter( "RTYP", "R-ratio computation type", &sr_type_ ); registerParameter( "PAIR", "Outgoing particles' PDG id", (int*)&proc_params_->operator[]( "pair" ) ); registerParameter( "INA1", "Heavy ion atomic weight (1st incoming beam)", (int*)&hi_1_.first ); registerParameter( "INZ1", "Heavy ion atomic number (1st incoming beam)", (int*)&hi_1_.second ); registerParameter( "INA2", "Heavy ion atomic weight (1st incoming beam)", (int*)&hi_2_.first ); registerParameter( "INZ2", "Heavy ion atomic number (1st incoming beam)", (int*)&hi_2_.second ); registerParameter( "INP1", "Momentum (1st primary particle)", ¶ms_.kinematics.incoming_beams.first.pz ); registerParameter( "INP2", "Momentum (2nd primary particle)", ¶ms_.kinematics.incoming_beams.second.pz ); registerParameter( "INPP", "Momentum (1st primary particle)", ¶ms_.kinematics.incoming_beams.first.pz ); registerParameter( "INPE", "Momentum (2nd primary particle)", ¶ms_.kinematics.incoming_beams.second.pz ); registerParameter( "PTCT", "Minimal transverse momentum (single central outgoing particle)", ¶ms_.kinematics.cuts.central.pt_single.min() ); registerParameter( "MSCT", "Minimal central system mass", ¶ms_.kinematics.cuts.central.mass_sum.min() ); registerParameter( "ECUT", "Minimal energy (single central outgoing particle)", ¶ms_.kinematics.cuts.central.energy_single.min() ); registerParameter( "ETMN", "Minimal pseudo-rapidity (central outgoing particles)", ¶ms_.kinematics.cuts.central.eta_single.min() ); registerParameter( "ETMX", "Maximal pseudo-rapidity (central outgoing particles)", ¶ms_.kinematics.cuts.central.eta_single.max() ); registerParameter( "YMIN", "Minimal rapidity (central outgoing particles)", ¶ms_.kinematics.cuts.central.rapidity_single.min() ); registerParameter( "YMAX", "Maximal rapidity (central outgoing particles)", ¶ms_.kinematics.cuts.central.rapidity_single.max() ); registerParameter( "Q2MN", "Minimal Q² = -q² (exchanged parton)", ¶ms_.kinematics.cuts.initial.q2.min() ); registerParameter( "Q2MX", "Maximal Q² = -q² (exchanged parton)", ¶ms_.kinematics.cuts.initial.q2.max() ); registerParameter( "MXMN", "Minimal invariant mass of proton remnants", ¶ms_.kinematics.cuts.remnants.mass_single.min() ); registerParameter( "MXMX", "Maximal invariant mass of proton remnants", ¶ms_.kinematics.cuts.remnants.mass_single.max() ); registerParameter( "XIMN", "Minimal fractional momentum loss of outgoing proton (ξ)", &xi_min_ ); registerParameter( "XIMX", "Maximal fractional momentum loss of outgoing proton (ξ)", &xi_max_ ); } void LpairHandler::store( const char* file ) { std::ofstream f( file, std::fstream::out | std::fstream::trunc ); if ( !f.is_open() ) throw CG_ERROR( "LpairHandler" ) << "Failed to open file \"" << file << "%s\" for writing."; for ( const auto& it : p_strings_ ) if ( it.second.value ) f << it.first << " = " << *it.second.value << "\n"; for ( const auto& it : p_ints_ ) if ( it.second.value ) f << it.first << " = " << *it.second.value << "\n"; for ( const auto& it : p_doubles_ ) if ( it.second.value ) f << it.first << " = " << *it.second.value << "\n"; for ( const auto& it : p_bools_ ) if ( it.second.value ) f << it.first << " = " << *it.second.value << "\n"; f.close(); } void LpairHandler::setParameter( const std::string& key, const std::string& value ) { - try { setValue( key.c_str(), std::stod( value ) ); } catch ( const std::invalid_argument& ) {} - try { setValue( key.c_str(), std::stoi( value ) ); } catch ( const std::invalid_argument& ) {} - try { setValue( key.c_str(), value ); } catch ( const std::invalid_argument& ) { - throw CG_FATAL( "LpairHandler:setParameter" ) - << "Failed to add the parameter \"" << key << "\" → \"" << value << "\"!"; + try { setValue( key.c_str(), std::stod( value ) ); } catch ( const std::invalid_argument& ) { + try { setValue( key.c_str(), std::stoi( value ) ); } catch ( const std::invalid_argument& ) { + try { setValue( key.c_str(), value ); } catch ( const std::invalid_argument& ) { + throw CG_FATAL( "LpairHandler:setParameter" ) + << "Failed to add the parameter \"" << key << "\" → \"" << value << "\"!"; + } + } } } std::string - LpairHandler::getParameter( std::string key ) const + LpairHandler::parameter( std::string key ) const { double dd = getValue( key.c_str() ); if ( dd != -999. ) return std::to_string( dd ); int ui = getValue( key.c_str() ); if ( ui != 999 ) return std::to_string( ui ); //if ( out = getValue( key.c_str() ) ); return getValue( key.c_str() ); } std::string - LpairHandler::getDescription( std::string key ) const + LpairHandler::description( std::string key ) const { if ( p_strings_.count( key ) ) return p_strings_.find( key )->second.description; if ( p_ints_.count( key ) ) return p_ints_.find( key )->second.description; if ( p_doubles_.count( key ) ) return p_doubles_.find( key )->second.description; if ( p_bools_.count( key ) ) return p_bools_.find( key )->second.description; return "null"; } } } diff --git a/CepGen/Cards/LpairHandler.h b/CepGen/Cards/LpairHandler.h index 0064dd0..22503f1 100644 --- a/CepGen/Cards/LpairHandler.h +++ b/CepGen/Cards/LpairHandler.h @@ -1,118 +1,121 @@ #ifndef CepGen_Cards_LpairReader_h #define CepGen_Cards_LpairReader_h #include "CepGen/Cards/Handler.h" #include using std::string; namespace cepgen { class ParametersList; namespace card { /// LPAIR-like steering cards parser and writer class LpairHandler : public Handler { public: /// Read a LPAIR steering card explicit LpairHandler( const char* file ); /// Store a configuration into a LPAIR steering card void store( const char* file ); private: /// Single parameter handler /// \tparam T Parameter type template struct Parameter { Parameter( const char* key, const char* descr, T* value ) : key( key ), description( descr ), value( value ) {} std::string key, description; T* value; }; /// Register a parameter to be steered to a configuration variable template void registerParameter( const char* key, const char* description, T* def ) {} /// Set a parameter value template void setValue( const char* key, const T& value ) {} /// Retrieve a parameter value template T getValue( const char* key ) const {} void setParameter( const std::string& key, const std::string& value ); - std::string getParameter( std::string key ) const; - std::string getDescription( std::string key ) const; + std::string parameter( std::string key ) const; + std::string description( std::string key ) const; static const int kInvalid; std::unordered_map > p_strings_; std::unordered_map > p_doubles_; std::unordered_map > p_ints_; std::unordered_map > p_bools_; void init(); std::shared_ptr proc_params_; int str_fun_, sr_type_; double xi_min_, xi_max_; - std::string proc_name_, hadr_name_, integr_type_, kmr_grid_path_, mstw_grid_path_; + std::string proc_name_, hadr_name_, out_mod_name_; + std::string out_file_name_; + std::string integr_type_; + std::string kmr_grid_path_, mstw_grid_path_; std::pair hi_1_, hi_2_; }; //----- specialised registerers /// Register a string parameter template<> inline void LpairHandler::registerParameter( const char* key, const char* description, std::string* def ) { p_strings_.insert( std::make_pair( key, Parameter( key, description, def ) ) ); } /// Register a double floating point parameter template<> inline void LpairHandler::registerParameter( const char* key, const char* description, double* def ) { p_doubles_.insert( std::make_pair( key, Parameter( key, description, def ) ) ); } /// Register an integer parameter template<> inline void LpairHandler::registerParameter( const char* key, const char* description, int* def ) { p_ints_.insert( std::make_pair( key, Parameter( key, description, def ) ) ); } /// Register a boolean parameter template<> inline void LpairHandler::registerParameter( const char* key, const char* description, bool* def ) { p_bools_.insert( std::make_pair( key, Parameter( key, description, def ) ) ); } //----- specialised setters template<> inline void LpairHandler::setValue( const char* key, const std::string& value ) { auto it = p_strings_.find( key ); if ( it != p_strings_.end() ) *it->second.value = value; } template<> inline void LpairHandler::setValue( const char* key, const double& value ) { auto it = p_doubles_.find( key ); if ( it != p_doubles_.end() ) *it->second.value = value; } template<> inline void LpairHandler::setValue( const char* key, const int& value ) { auto it = p_ints_.find( key ); if ( it != p_ints_.end() ) *it->second.value = value; } template<> inline void LpairHandler::setValue( const char* key, const bool& value ) { auto it = p_bools_.find( key ); if ( it != p_bools_.end() ) *it->second.value = value; } //----- specialised getters /// Retrieve a string parameter value template<> inline std::string LpairHandler::getValue( const char* key ) const { const auto& it = p_strings_.find( key ); if ( it != p_strings_.end() ) return *it->second.value; return "null"; } /// Retrieve a floating point parameter value template<> inline double LpairHandler::getValue( const char* key ) const { const auto& it = p_doubles_.find( key ); if ( it != p_doubles_.end() ) return *it->second.value; return -999.; } /// Retrieve an integer parameter value template<> inline int LpairHandler::getValue( const char* key ) const { const auto& it = p_ints_.find( key ); if ( it != p_ints_.end() ) return *it->second.value; return 999; } /// Retrieve a boolean parameter value template<> inline bool LpairHandler::getValue( const char* key ) const { const auto& it = p_bools_.find( key ); if ( it != p_bools_.end() ) return *it->second.value; return true; } } } #endif diff --git a/CepGen/Cards/PythonHandler.cpp b/CepGen/Cards/PythonHandler.cpp index 235e88e..29d9b9c 100644 --- a/CepGen/Cards/PythonHandler.cpp +++ b/CepGen/Cards/PythonHandler.cpp @@ -1,373 +1,393 @@ #include "CepGen/Cards/PythonHandler.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/TamingFunction.h" #include "CepGen/Core/ParametersList.h" #include "CepGen/Core/Integrator.h" #include "CepGen/Processes/ProcessesHandler.h" #include "CepGen/Hadronisers/HadronisersHandler.h" +#include "CepGen/IO/ExportHandler.h" #include "CepGen/StructureFunctions/StructureFunctions.h" #include "CepGen/Physics/GluonGrid.h" #include "CepGen/Physics/PDG.h" #include #if PY_MAJOR_VERSION < 3 # define PYTHON2 #endif namespace cepgen { namespace card { //----- specialization for CepGen input cards PythonHandler::PythonHandler( const char* file ) { setenv( "PYTHONPATH", ".:Cards:test:../Cards", 1 ); setenv( "PYTHONDONTWRITEBYTECODE", "1", 1 ); CG_DEBUG( "PythonHandler" ) << "Python PATH: " << getenv( "PYTHONPATH" ) << "."; std::string filename = pythonPath( file ); const size_t fn_len = filename.length()+1; //Py_DebugFlag = 1; //Py_VerboseFlag = 1; #ifdef PYTHON2 char* sfilename = new char[fn_len]; snprintf( sfilename, fn_len, "%s", filename.c_str() ); #else wchar_t* sfilename = new wchar_t[fn_len]; swprintf( sfilename, fn_len, L"%s", filename.c_str() ); #endif if ( sfilename ) Py_SetProgramName( sfilename ); Py_InitializeEx( 1 ); if ( sfilename ) delete [] sfilename; if ( !Py_IsInitialized() ) throw CG_FATAL( "PythonHandler" ) << "Failed to initialise the Python cards parser!"; CG_DEBUG( "PythonHandler" ) << "Initialised the Python cards parser\n\t" << "Python version: " << Py_GetVersion() << "\n\t" << "Platform: " << Py_GetPlatform() << "."; PyObject* cfg = PyImport_ImportModule( filename.c_str() ); // new if ( !cfg ) throwPythonError( Form( "Failed to parse the configuration card %s", file ) ); //--- additional particles definition PyObject* pextp = PyObject_GetAttrString( cfg, PDGLIST_NAME ); // new if ( pextp ) { parseExtraParticles( pextp ); Py_CLEAR( pextp ); } //--- process definition PyObject* process = PyObject_GetAttrString( cfg, PROCESS_NAME ); // new if ( !process ) throwPythonError( Form( "Failed to extract a \"%s\" keyword from the configuration card %s", PROCESS_NAME, file ) ); //--- list of process-specific parameters ParametersList proc_params; fillParameter( process, "processParameters", proc_params ); //--- type of process to consider PyObject* pproc_name = element( process, MODULE_NAME ); // borrowed if ( !pproc_name ) throwPythonError( Form( "Failed to extract the process name from the configuration card %s", file ) ); const std::string proc_name = get( pproc_name ); //--- process mode params_.kinematics.mode = (KinematicsMode)proc_params.get( "mode", (int)KinematicsMode::invalid ); params_.setProcess( cepgen::proc::ProcessesHandler::get().build( proc_name, proc_params ) ); //--- process kinematics PyObject* pin_kinematics = element( process, "inKinematics" ); // borrowed if ( pin_kinematics ) parseIncomingKinematics( pin_kinematics ); PyObject* pout_kinematics = element( process, "outKinematics" ); // borrowed if ( pout_kinematics ) parseOutgoingKinematics( pout_kinematics ); //--- taming functions PyObject* ptam = element( process, "tamingFunctions" ); // borrowed if ( ptam ) for ( const auto& p : getVector( ptam ) ) params_.taming_functions->add( p.get( "variable" ), p.get( "expression" ) ); Py_CLEAR( process ); PyObject* plog = PyObject_GetAttrString( cfg, LOGGER_NAME ); // new if ( plog ) { parseLogging( plog ); Py_CLEAR( plog ); } //--- hadroniser parameters PyObject* phad = PyObject_GetAttrString( cfg, HADR_NAME ); // new if ( phad ) { parseHadroniser( phad ); Py_CLEAR( phad ); } //--- generation parameters PyObject* pint = PyObject_GetAttrString( cfg, INTEGRATOR_NAME ); // new if ( pint ) { parseIntegrator( pint ); Py_CLEAR( pint ); } PyObject* pgen = PyObject_GetAttrString( cfg, GENERATOR_NAME ); // new if ( pgen ) { parseGenerator( pgen ); Py_CLEAR( pgen ); } + PyObject* pout = PyObject_GetAttrString( cfg, OUTPUT_NAME ); // new + if ( pout ) { + parseOutputModule( pout ); + Py_CLEAR( pout ); + } + //--- finalisation Py_CLEAR( cfg ); } PythonHandler::~PythonHandler() { if ( Py_IsInitialized() ) Py_Finalize(); } void PythonHandler::parseIncomingKinematics( PyObject* kin ) { //--- retrieve the beams PDG ids std::vector beams_pdg; fillParameter( kin, "pdgIds", beams_pdg ); if ( !beams_pdg.empty() ) { if ( beams_pdg.size() != 2 ) throwPythonError( Form( "Invalid list of PDG ids retrieved for incoming beams:\n\t2 PDG ids are expected, %d provided!", beams_pdg.size() ) ); params_.kinematics.incoming_beams. first.pdg = (pdgid_t)beams_pdg.at( 0 ).get( "pdgid" ); params_.kinematics.incoming_beams.second.pdg = (pdgid_t)beams_pdg.at( 1 ).get( "pdgid" ); } //--- incoming beams kinematics std::vector beams_pz; fillParameter( kin, "pz", beams_pz ); if ( !beams_pz.empty() ) { if ( beams_pz.size() != 2 ) throwPythonError( Form( "Invalid list of pz's retrieved for incoming beams:\n\t2 pz's are expected, %d provided!", beams_pz.size() ) ); params_.kinematics.incoming_beams. first.pz = beams_pz.at( 0 ); params_.kinematics.incoming_beams.second.pz = beams_pz.at( 1 ); } double sqrt_s = -1.; fillParameter( kin, "cmEnergy", sqrt_s ); if ( sqrt_s != -1. ) params_.kinematics.setSqrtS( sqrt_s ); //--- structure functions set for incoming beams PyObject* psf = element( kin, "structureFunctions" ); // borrowed if ( psf ) params_.kinematics.structure_functions = strfun::StructureFunctionsHandler::get().build( get( psf ) ); //--- types of parton fluxes for kt-factorisation std::vector kt_fluxes; fillParameter( kin, "ktFluxes", kt_fluxes ); if ( !kt_fluxes.empty() ) { params_.kinematics.incoming_beams.first.kt_flux = (KTFlux)kt_fluxes.at( 0 ); params_.kinematics.incoming_beams.second.kt_flux = ( kt_fluxes.size() > 1 ) ? (KTFlux)kt_fluxes.at( 1 ) : (KTFlux)kt_fluxes.at( 0 ); } //--- specify where to look for the grid path for gluon emission std::string kmr_grid_path; fillParameter( kin, "kmrGridPath", kmr_grid_path ); if ( !kmr_grid_path.empty() ) kmr::GluonGrid::get( kmr_grid_path.c_str() ); //--- parse heavy ions beams std::vector hi_beam1, hi_beam2; fillParameter( kin, "heavyIonA", hi_beam1 ); if ( hi_beam1.size() == 2 ) params_.kinematics.incoming_beams. first.pdg = HeavyIon{ (unsigned short)hi_beam1[0], (Element)hi_beam1[1] }; fillParameter( kin, "heavyIonB", hi_beam2 ); if ( hi_beam2.size() == 2 ) params_.kinematics.incoming_beams.second.pdg = HeavyIon{ (unsigned short)hi_beam2[0], (Element)hi_beam2[1] }; } void PythonHandler::parseOutgoingKinematics( PyObject* kin ) { std::vector parts; fillParameter( kin, "minFinalState", parts ); for ( const auto& pdg : parts ) params_.kinematics.minimum_final_state.emplace_back( (pdgid_t)pdg ); ParametersList part_cuts; fillParameter( kin, "cuts", part_cuts ); for ( const auto& part : part_cuts.keys() ) { const auto pdg = (pdgid_t)stoi( part ); const auto& cuts = part_cuts.get( part ); if ( cuts.has( "pt" ) ) params_.kinematics.cuts.central_particles[pdg].pt_single = cuts.get( "pt" ); if ( cuts.has( "energy" ) ) params_.kinematics.cuts.central_particles[pdg].energy_single = cuts.get( "energy" ); if ( cuts.has( "eta" ) ) params_.kinematics.cuts.central_particles[pdg].eta_single = cuts.get( "eta" ); if ( cuts.has( "rapidity" ) ) params_.kinematics.cuts.central_particles[pdg].rapidity_single = cuts.get( "rapidity" ); } // for LPAIR/collinear matrix elements fillParameter( kin, "q2", params_.kinematics.cuts.initial.q2 ); // for the kT factorised matrix elements fillParameter( kin, "qt", params_.kinematics.cuts.initial.qt ); fillParameter( kin, "phiqt", params_.kinematics.cuts.initial.phi_qt ); fillParameter( kin, "ptdiff", params_.kinematics.cuts.central.pt_diff ); fillParameter( kin, "phiptdiff", params_.kinematics.cuts.central.phi_pt_diff ); fillParameter( kin, "rapiditydiff", params_.kinematics.cuts.central.rapidity_diff ); // generic phase space limits fillParameter( kin, "rapidity", params_.kinematics.cuts.central.rapidity_single ); fillParameter( kin, "eta", params_.kinematics.cuts.central.eta_single ); fillParameter( kin, "pt", params_.kinematics.cuts.central.pt_single ); fillParameter( kin, "ptsum", params_.kinematics.cuts.central.pt_sum ); fillParameter( kin, "invmass", params_.kinematics.cuts.central.mass_sum ); fillParameter( kin, "mx", params_.kinematics.cuts.remnants.mass_single ); fillParameter( kin, "yj", params_.kinematics.cuts.remnants.rapidity_single ); Limits lim_xi; fillParameter( kin, "xi", lim_xi ); if ( lim_xi.valid() ) - params_.kinematics.cuts.remnants.energy_single = ( lim_xi+(-1.) )*( -params_.kinematics.incoming_beams.first.pz ); + //params_.kinematics.cuts.remnants.energy_single = ( lim_xi+(-1.) )*( -params_.kinematics.incoming_beams.first.pz ); + params_.kinematics.cuts.remnants.energy_single = -( lim_xi-1. )*params_.kinematics.incoming_beams.first.pz; } void PythonHandler::parseLogging( PyObject* log ) { int log_level = 0; fillParameter( log, "level", log_level ); utils::Logger::get().level = (utils::Logger::Level)log_level; std::vector enabled_modules; fillParameter( log, "enabledModules", enabled_modules ); for ( const auto& mod : enabled_modules ) utils::Logger::get().addExceptionRule( mod ); } void PythonHandler::parseIntegrator( PyObject* integr ) { if ( !PyDict_Check( integr ) ) throwPythonError( "Integrator object should be a dictionary!" ); PyObject* palgo = element( integr, MODULE_NAME ); // borrowed if ( !palgo ) throwPythonError( "Failed to retrieve the integration algorithm name!" ); std::string algo = get( palgo ); if ( algo == "plain" ) params_.integration().type = IntegratorType::plain; else if ( algo == "Vegas" ) { params_.integration().type = IntegratorType::Vegas; fillParameter( integr, "alpha", (double&)params_.integration().vegas.alpha ); fillParameter( integr, "iterations", params_.integration().vegas.iterations ); fillParameter( integr, "mode", (int&)params_.integration().vegas.mode ); fillParameter( integr, "verbosity", (int&)params_.integration().vegas.verbose ); std::string vegas_logging_output = "cerr"; fillParameter( integr, "loggingOutput", vegas_logging_output ); if ( vegas_logging_output == "cerr" ) // redirect all debugging information to the error stream params_.integration().vegas.ostream = stderr; else if ( vegas_logging_output == "cout" ) // redirect all debugging information to the standard stream params_.integration().vegas.ostream = stdout; else params_.integration().vegas.ostream = fopen( vegas_logging_output.c_str(), "w" ); } else if ( algo == "MISER" ) { params_.integration().type = IntegratorType::MISER; fillParameter( integr, "estimateFraction", (double&)params_.integration().miser.estimate_frac ); fillParameter( integr, "minCalls", params_.integration().miser.min_calls ); fillParameter( integr, "minCallsPerBisection", params_.integration().miser.min_calls_per_bisection ); fillParameter( integr, "alpha", (double&)params_.integration().miser.alpha ); fillParameter( integr, "dither", (double&)params_.integration().miser.dither ); } else throwPythonError( Form( "Invalid integration() algorithm: %s", algo.c_str() ) ); fillParameter( integr, "numFunctionCalls", params_.integration().ncvg ); fillParameter( integr, "seed", (unsigned long&)params_.integration().rng_seed ); unsigned int rng_engine; fillParameter( integr, "rngEngine", rng_engine ); switch ( rng_engine ) { case 0: default: params_.integration().rng_engine = (gsl_rng_type*)gsl_rng_mt19937; break; case 1: params_.integration().rng_engine = (gsl_rng_type*)gsl_rng_taus2; break; case 2: params_.integration().rng_engine = (gsl_rng_type*)gsl_rng_gfsr4; break; case 3: params_.integration().rng_engine = (gsl_rng_type*)gsl_rng_ranlxs0; break; } fillParameter( integr, "chiSqCut", params_.integration().vegas_chisq_cut ); } void PythonHandler::parseGenerator( PyObject* gen ) { if ( !PyDict_Check( gen ) ) throwPythonError( "Generation information object should be a dictionary!" ); params_.generation().enabled = true; fillParameter( gen, "treat", params_.generation().treat ); fillParameter( gen, "numEvents", params_.generation().maxgen ); fillParameter( gen, "printEvery", params_.generation().gen_print_every ); fillParameter( gen, "numThreads", params_.generation().num_threads ); fillParameter( gen, "numPoints", params_.generation().num_points ); } void PythonHandler::parseHadroniser( PyObject* hadr ) { if ( !PyDict_Check( hadr ) ) throwPythonError( "Hadroniser object should be a dictionary!" ); PyObject* pname = element( hadr, MODULE_NAME ); // borrowed if ( !pname ) throwPythonError( "Hadroniser name is required!" ); std::string hadr_name = get( pname ); params_.setHadroniser( cepgen::hadr::HadronisersHandler::get().build( hadr_name, get( hadr ) ) ); auto h = params_.hadroniser(); h->setParameters( params_ ); { //--- before calling the init() method std::vector config; fillParameter( hadr, "preConfiguration", config ); h->readStrings( config ); } h->init(); { //--- after init() has been called std::vector config; fillParameter( hadr, "processConfiguration", config ); for ( const auto& block : config ) { std::vector config_blk; fillParameter( hadr, block.c_str(), config_blk ); h->readStrings( config_blk ); } } } void + PythonHandler::parseOutputModule( PyObject* pout ) + { + if ( !is( pout ) ) + throwPythonError( "Invalid type for output parameters list!" ); + + PyObject* pname = element( pout, MODULE_NAME ); // borrowed + if ( !pname ) + throwPythonError( "Output module name is required!" ); + params_.setOutputModule( io::ExportHandler::get().build( get( pname ), get( pout ) ) ); + } + + void PythonHandler::parseExtraParticles( PyObject* pparts ) { if ( !is( pparts ) ) throwPythonError( "Extra particles definition object should be a parameters list!" ); const auto& parts = get( pparts ); for ( const auto& k : parts.keys() ) { const auto& part = parts.get( k ); if ( part.pdgid == 0 || part.mass < 0. ) continue; CG_DEBUG( "PythonHandler:particles" ) << "Adding a new particle with name \"" << part.name << "\" to the PDG dictionary."; PDG::get().define( part ); } } } } diff --git a/CepGen/Cards/PythonHandler.h b/CepGen/Cards/PythonHandler.h index 977da14..cdc02cd 100644 --- a/CepGen/Cards/PythonHandler.h +++ b/CepGen/Cards/PythonHandler.h @@ -1,82 +1,84 @@ #ifndef CepGen_Cards_PythonHandler_h #define CepGen_Cards_PythonHandler_h #include #include "Handler.h" namespace cepgen { namespace strfun { class Parameterisation; } class Limits; class ParametersList; namespace card { /// CepGen Python configuration cards reader/writer class PythonHandler : public Handler { public: /// Read a standard configuration card explicit PythonHandler( const char* file ); ~PythonHandler(); private: static constexpr const char* MODULE_NAME = "mod_name"; static constexpr const char* PROCESS_NAME = "process"; static constexpr const char* HADR_NAME = "hadroniser"; static constexpr const char* LOGGER_NAME = "logger"; static constexpr const char* INTEGRATOR_NAME = "integrator"; static constexpr const char* GENERATOR_NAME = "generator"; + static constexpr const char* OUTPUT_NAME = "output"; static constexpr const char* PDGLIST_NAME = "PDG"; static void throwPythonError( const std::string& message ); static std::string pythonPath( const char* file ); static PyObject* element( PyObject* obj, const char* key ); static PyObject* encode( const char* str ); template bool is( PyObject* obj ) const; template T get( PyObject* obj ) const; template bool isVector( PyObject* obj ) const; template std::vector getVector( PyObject* obj ) const; void fillParameter( PyObject* parent, const char* key, bool& out ); void fillParameter( PyObject* parent, const char* key, int& out ); void fillParameter( PyObject* parent, const char* key, unsigned long& out ); void fillParameter( PyObject* parent, const char* key, unsigned int& out ); void fillParameter( PyObject* parent, const char* key, double& out ); void fillParameter( PyObject* parent, const char* key, std::string& out ); void fillParameter( PyObject* parent, const char* key, Limits& out ); void fillParameter( PyObject* parent, const char* key, std::vector& out ); void fillParameter( PyObject* parent, const char* key, std::vector& out ); void fillParameter( PyObject* parent, const char* key, std::vector& out ); void fillParameter( PyObject* parent, const char* key, ParametersList& out ); void fillParameter( PyObject* parent, const char* key, std::vector& out ); void parseIncomingKinematics( PyObject* ); void parseOutgoingKinematics( PyObject* ); void parseLogging( PyObject* ); void parseIntegrator( PyObject* ); void parseGenerator( PyObject* ); void parseHadroniser( PyObject* ); - void parseExtraParticles( PyObject* parts ); + void parseOutputModule( PyObject* ); + void parseExtraParticles( PyObject* ); }; template<> bool PythonHandler::is( PyObject* obj ) const; template<> bool PythonHandler::is( PyObject* obj ) const; template<> bool PythonHandler::is( PyObject* obj ) const; template<> int PythonHandler::get( PyObject* obj ) const; template<> unsigned long PythonHandler::get( PyObject* obj ) const; template<> bool PythonHandler::is( PyObject* obj ) const; template<> ParametersList PythonHandler::get( PyObject* obj ) const; template<> bool PythonHandler::is( PyObject* obj ) const; template<> double PythonHandler::get( PyObject* obj ) const; template<> bool PythonHandler::is( PyObject* obj ) const; template<> std::string PythonHandler::get( PyObject* obj ) const; template<> bool PythonHandler::is( PyObject* obj ) const; template<> Limits PythonHandler::get( PyObject* obj ) const; } } #endif diff --git a/CepGen/Core/FortranInterface.cpp b/CepGen/Core/FortranInterface.cpp index 47346ee..26a8673 100644 --- a/CepGen/Core/FortranInterface.cpp +++ b/CepGen/Core/FortranInterface.cpp @@ -1,76 +1,98 @@ #include "CepGen/StructureFunctions/StructureFunctions.h" +#include "CepGen/Processes/FortranKTProcess.h" #include "CepGen/Physics/KTFlux.h" #include "CepGen/Physics/HeavyIon.h" #include "CepGen/Physics/PDG.h" #include "CepGen/Core/ParametersList.h" #include "CepGen/Core/Exception.h" #ifdef __cplusplus extern "C" { #endif /// Expose structure functions calculators to Fortran void cepgen_structure_functions_( int& sfmode, double& xbj, double& q2, double& f2, double& fl ) { using namespace cepgen; static auto sf = strfun::StructureFunctionsHandler::get().build( sfmode ); const auto& val = ( *sf )( xbj, q2 ); f2 = val.F2; fl = val.FL; } /// Compute a \f$k_{\rm T}\f$-dependent flux for single nucleons /// \param[in] fmode Flux mode (see cepgen::KTFlux) /// \param[in] x Fractional momentum loss /// \param[in] kt2 The \f$k_{\rm T}\f$ transverse momentum norm /// \param[in] sfmode Structure functions set for dissociative emission /// \param[in] mx Diffractive state mass for dissociative emission double cepgen_kt_flux_( int& fmode, double& x, double& kt2, int& sfmode, double& mx ) { using namespace cepgen; static auto sf = strfun::StructureFunctionsHandler::get().build( sfmode ); return ktFlux( (KTFlux)fmode, x, kt2, *sf, mx ); } /// Compute a \f$k_{\rm T}\f$-dependent flux for heavy ions /// \param[in] fmode Flux mode (see cepgen::KTFlux) /// \param[in] x Fractional momentum loss /// \param[in] kt2 The \f$k_{\rm T}\f$ transverse momentum norm /// \param[in] a Mass number for the heavy ion /// \param[in] z Atomic number for the heavy ion double cepgen_kt_flux_hi_( int& fmode, double& x, double& kt2, int& a, int& z ) { using namespace cepgen; return ktFlux( (KTFlux)fmode, x, kt2, HeavyIon{ (unsigned short)a, (Element)z } ); } /// Mass of a particle, in GeV/c^2 double cepgen_particle_mass_( int& pdg_id ) { try { return cepgen::PDG::get().mass( (cepgen::pdgid_t)pdg_id ); } catch ( const cepgen::Exception& e ) { e.dump(); exit( 0 ); } } /// Charge of a particle, in e double cepgen_particle_charge_( int& pdg_id ) { try { return cepgen::PDG::get().charge( (cepgen::pdgid_t)pdg_id ); } catch ( const cepgen::Exception& e ) { e.dump(); exit( 0 ); } } + + void + cepgen_list_params_() + { + CG_LOG( "cepgen_list_params" ) << "\t" << cepgen::proc::FortranKTProcess::kProcParameters; + } + + int + cepgen_param_int_( char* pname, int& def ) + { + if ( cepgen::proc::FortranKTProcess::kProcParameters.has( pname ) ) + return cepgen::proc::FortranKTProcess::kProcParameters.get( pname ).pdgid; + return cepgen::proc::FortranKTProcess::kProcParameters.get( pname, def ); + } + + double + cepgen_param_real_( char* pname, double& def ) + { + return cepgen::proc::FortranKTProcess::kProcParameters.get( pname, def ); + } + #ifdef __cplusplus } #endif diff --git a/CepGen/Core/Integrand.cpp b/CepGen/Core/Integrand.cpp index d6e8f12..4109c12 100644 --- a/CepGen/Core/Integrand.cpp +++ b/CepGen/Core/Integrand.cpp @@ -1,201 +1,202 @@ #include "CepGen/Core/Timer.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/TamingFunction.h" #include "CepGen/Event/Event.h" #include "CepGen/Event/Particle.h" #include "CepGen/Physics/Kinematics.h" #include "CepGen/Physics/PDG.h" #include "CepGen/Processes/GenericProcess.h" #include "CepGen/Hadronisers/GenericHadroniser.h" +#include "CepGen/IO/GenericExportHandler.h" #include "CepGen/Parameters.h" #include #include #include namespace cepgen { namespace integrand { utils::Logger::Level log_level; utils::Timer tmr; double eval( double* x, size_t ndim, void* func_params ) { log_level = utils::Logger::get().level; std::shared_ptr ev; Parameters* params = nullptr; if ( !func_params || !( params = static_cast( func_params ) ) ) throw CG_FATAL( "Integrand" ) << "Failed to retrieve the run parameters!"; proc::GenericProcess* proc = params->process(); if ( !proc ) throw CG_FATAL( "Integrand" ) << "Failed to retrieve the process!"; //================================================================ // start the timer //================================================================ tmr.reset(); //================================================================ // prepare the event content prior to the process generation //================================================================ if ( proc->hasEvent() ) // event is not empty ev = proc->event(); params->prepareRun(); //================================================================ // specify the phase space point to probe //================================================================ proc->setPoint( ndim, x ); //================================================================ // from this step on, the phase space point is supposed to be set //================================================================ proc->beforeComputeWeight(); double weight = proc->computeWeight(); //================================================================ // invalidate any unphysical behaviour //================================================================ if ( weight <= 0. ) return 0.; //================================================================ // speed up the integration process if no event is to be generated //================================================================ if ( !ev ) return weight; if ( !params->storage() && !params->taming_functions && !params->hadroniser() && params->kinematics.cuts.central_particles.empty() ) return weight; //================================================================ // fill in the process' Event object //================================================================ proc->fillKinematics(); //================================================================ // once the kinematics variables have been populated, can apply // the collection of taming functions //================================================================ if ( params->taming_functions ) { if ( params->taming_functions->has( "m_central" ) || params->taming_functions->has( "pt_central" ) ) { // build the kinematics of the central system Particle::Momentum central_system; for ( const auto& part : (*ev)[Particle::CentralSystem] ) central_system += part.momentum(); // tame the cross-section by the reweighting function if ( params->taming_functions->has( "m_central" ) ) weight *= params->taming_functions->eval( "m_central", central_system.mass() ); if ( params->taming_functions->has( "pt_central" ) ) weight *= params->taming_functions->eval( "pt_central", central_system.pt() ); } if ( params->taming_functions->has( "q2" ) ) { weight *= params->taming_functions->eval( "q2", -ev->getOneByRole( Particle::Parton1 ).momentum().mass() ); weight *= params->taming_functions->eval( "q2", -ev->getOneByRole( Particle::Parton2 ).momentum().mass() ); } } if ( weight <= 0. ) return 0.; //================================================================ // set the CepGen part of the event generation //================================================================ if ( params->storage() ) ev->time_generation = tmr.elapsed(); //================================================================ // event hadronisation and resonances decay //================================================================ if ( params->hadroniser() ) { double br = -1.; if ( !params->hadroniser()->run( *ev, br, params->storage() ) || br == 0. ) return 0.; weight *= br; // branching fraction for all decays } //================================================================ // apply cuts on final state system (after hadronisation!) // (polish your cuts, as this might be very time-consuming...) //================================================================ if ( !params->kinematics.cuts.central_particles.empty() ) { for ( const auto& part : (*ev)[Particle::CentralSystem] ) { // retrieve all cuts associated to this final state particle if ( params->kinematics.cuts.central_particles.count( part.pdgId() ) == 0 ) continue; const auto& cuts_pdgid = params->kinematics.cuts.central_particles.at( part.pdgId() ); // apply these cuts on the given particle if ( !cuts_pdgid.pt_single.passes( part.momentum().pt() ) ) return 0.; if ( !cuts_pdgid.energy_single.passes( part.momentum().energy() ) ) return 0.; if ( !cuts_pdgid.eta_single.passes( part.momentum().eta() ) ) return 0.; if ( !cuts_pdgid.rapidity_single.passes( part.momentum().rapidity() ) ) return 0.; } } for ( const auto& system : { Particle::OutgoingBeam1, Particle::OutgoingBeam2 } ) for ( const auto& part : (*ev)[system] ) { if ( part.status() != Particle::Status::FinalState ) continue; - if ( !params->kinematics.cuts.remnants.energy_single.passes( fabs( part.momentum().energy() ) ) ) + if ( !params->kinematics.cuts.remnants.energy_single.passes( part.momentum().energy() ) ) return 0.; if ( !params->kinematics.cuts.remnants.rapidity_single.passes( fabs( part.momentum().rapidity() ) ) ) return 0.; } //================================================================ // store the last event in parameters block for a later usage //================================================================ if ( params->storage() ) { proc->last_event = ev; proc->last_event->time_total = tmr.elapsed(); CG_DEBUG( "Integrand" ) << "[process 0x" << std::hex << proc << std::dec << "] " << "Individual time (gen+hadr+cuts): " << proc->last_event->time_total*1.e3 << " ms"; } //================================================================ // a bit of useful debugging //================================================================ if ( CG_LOG_MATCH( "Integrand", debugInsideLoop ) ) { std::ostringstream oss; for ( unsigned short i = 0; i < ndim; ++i ) oss << Form( "%10.8f ", x[i] ); CG_DEBUG( "Integrand" ) << "f value for dim-" << ndim << " point ( " << oss.str() << "): " << weight; } return weight; } } } diff --git a/CepGen/Core/Integrator.cpp b/CepGen/Core/Integrator.cpp index 5071230..ce7a162 100644 --- a/CepGen/Core/Integrator.cpp +++ b/CepGen/Core/Integrator.cpp @@ -1,469 +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& 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 ); + 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/Core/ModuleFactory.h b/CepGen/Core/ModuleFactory.h index 622d57b..5582d78 100644 --- a/CepGen/Core/ModuleFactory.h +++ b/CepGen/Core/ModuleFactory.h @@ -1,99 +1,96 @@ #ifndef CepGen_Core_ModuleFactory_h #define CepGen_Core_ModuleFactory_h #include "CepGen/Core/ParametersList.h" #include #include #include #include #define BUILDERNM( obj ) obj ## Builder #define STRINGIFY( name ) #name namespace cepgen { /// A generic factory to build modules /// \tparam T Base class to build /// \tparam I Indexing variable type template class ModuleFactory { public: /// Retrieve a unique instance of this factory static ModuleFactory& get() { static ModuleFactory instance; return instance; } /// Default destructor ~ModuleFactory() = default; /// Register a named module in the database /// \tparam U Class to register (inherited from T base class) template void registerModule( const I& name, const ParametersList& def_params = ParametersList() ) { - static_assert( std::is_base_of::value, "\n Failed to register an object with improper inheritance into the factory." ); - map_[name] = &create; + static_assert( std::is_base_of::value, "\n\n *** Failed to register an object with improper inheritance into the factory. ***\n" ); + map_[name] = &build; params_map_[name] = def_params; } /// Build one instance of a named module /// \param[in] name Module name to retrieve /// \param[in] params List of parameters to be invoked by the constructor std::unique_ptr build( const I& name, ParametersList params = ParametersList() ) const { if ( name == I() || map_.count( name ) == 0 ) { std::ostringstream oss; - oss << __PRETTY_FUNCTION__ << "\n Failed to retrieve a module with index \"" << name << "\" from factory!"; - throw std::runtime_error( oss.str() ); + oss << __PRETTY_FUNCTION__ << "\n\n *** Failed to build a module with index/name \"" << name << "\" from factory! ***\n"; + throw std::invalid_argument( oss.str() ); } if ( params_map_.count( name ) > 0 ) params += params_map_.at( name ); return map_.at( name )( params ); } /// Build one instance of a named module /// \param[in] params List of parameters to be invoked by the constructor std::unique_ptr build( ParametersList params = ParametersList() ) const { if ( params.has( KEY ) ) { const I& idx = params.get( KEY ); - if ( map_.count( idx ) == 0 ) { - std::ostringstream oss; - oss << __PRETTY_FUNCTION__ << " Failed to retrieve a module with index \"" << idx << "\" from factory!"; - throw std::runtime_error( oss.str() ); - } + if ( map_.count( idx ) == 0 ) + throw std::invalid_argument( std::string( __PRETTY_FUNCTION__ )+"\n\n *** Failed to build a module with index/name \""+std::to_string( idx )+"\" from factory! ***\n" ); if ( params_map_.count( idx ) > 0 ) params += params_map_.at( idx ); return map_.at( idx )( params ); } else - throw std::runtime_error( std::string( __PRETTY_FUNCTION__ )+" Failed to retrieve an indexing key from parameters to build from factory!" ); + throw std::invalid_argument( std::string( __PRETTY_FUNCTION__ )+"\n\n *** Failed to retrieve an indexing key from parameters to build from factory! ***\n" ); } /// List of modules registred in the database std::vector modules() const { std::vector out; for ( const auto& p : map_ ) out.emplace_back( p.first ); return out; } static constexpr const char* KEY = "id"; private: explicit ModuleFactory() = default; /// Construct a module with its parameters set - template static std::unique_ptr create( const ParametersList& params ) { + template static std::unique_ptr build( const ParametersList& params ) { return std::unique_ptr( new U( params ) ); } protected: /// Constructor type for a module typedef std::unique_ptr (*ModCreate)( const ParametersList& ); /// Database of modules handled by this instance std::unordered_map map_; /// Database of default parameters associated to modules std::unordered_map params_map_; public: ModuleFactory( const ModuleFactory& ) = delete; void operator=( const ModuleFactory& ) = delete; }; } #endif diff --git a/CepGen/Core/Parameters.cpp b/CepGen/Core/Parameters.cpp index fc699bf..33ab96a 100644 --- a/CepGen/Core/Parameters.cpp +++ b/CepGen/Core/Parameters.cpp @@ -1,309 +1,332 @@ #include "CepGen/Parameters.h" #include "CepGen/Core/Integrator.h" #include "CepGen/Core/ParametersList.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/TamingFunction.h" #include "CepGen/Physics/PDG.h" #include "CepGen/Processes/GenericProcess.h" #include "CepGen/Hadronisers/GenericHadroniser.h" +#include "CepGen/IO/GenericExportHandler.h" #include "CepGen/StructureFunctions/StructureFunctions.h" #include namespace cepgen { Parameters::Parameters() : general( new ParametersList ), taming_functions( new utils::TamingFunctionsCollection ), store_( false ), total_gen_time_( 0. ), num_gen_events_( 0ul ) {} Parameters::Parameters( Parameters& param ) : general( param.general ), kinematics( param.kinematics ), taming_functions( param.taming_functions ), process_( std::move( param.process_ ) ), hadroniser_( std::move( param.hadroniser_ ) ), + out_module_( std::move( param.out_module_ ) ), store_( false ), total_gen_time_( param.total_gen_time_ ), num_gen_events_( param.num_gen_events_ ), integration_( param.integration_ ), generation_( param.generation_ ) {} Parameters::Parameters( const Parameters& param ) : general( param.general ), kinematics( param.kinematics ), taming_functions( param.taming_functions ), store_( false ), total_gen_time_( param.total_gen_time_ ), num_gen_events_( param.num_gen_events_ ), integration_( param.integration_ ), generation_( param.generation_ ) {} Parameters::~Parameters() // required for unique_ptr initialisation! {} Parameters& Parameters::operator=( Parameters param ) { general = param.general; kinematics = param.kinematics; taming_functions = param.taming_functions; process_ = std::move( param.process_ ); hadroniser_ = std::move( param.hadroniser_ ); + out_module_ = std::move( param.out_module_ ); total_gen_time_ = param.total_gen_time_; num_gen_events_ = param.num_gen_events_; integration_ = param.integration_; generation_ = param.generation_; return *this; } void Parameters::setThetaRange( float thetamin, float thetamax ) { kinematics.cuts.central.eta_single = { Particle::thetaToEta( thetamax ), Particle::thetaToEta( thetamin ) }; CG_DEBUG( "Parameters" ) << "eta in range: " << kinematics.cuts.central.eta_single << " => theta(min) = " << thetamin << ", theta(max) = " << thetamax << "."; } void Parameters::prepareRun() { //--- first-run preparation if ( !process_ || !process_->first_run ) return; CG_DEBUG( "Parameters" ) << "Run started for " << process_->name() << " process " << "0x" << std::hex << process_.get() << std::dec << ".\n\t" << "Process mode considered: " << kinematics.mode << "\n\t" << " first beam: " << kinematics.incoming_beams.first << "\n\t" << " second beam: " << kinematics.incoming_beams.second << "\n\t" << " structure functions: " << kinematics.structure_functions; if ( process_->hasEvent() ) process_->clearEvent(); //--- clear the run statistics total_gen_time_ = 0.; num_gen_events_ = 0ul; process_->first_run = false; } void Parameters::addGenerationTime( double gen_time ) { total_gen_time_ += gen_time; num_gen_events_++; } proc::GenericProcess* Parameters::process() { return process_.get(); } const proc::GenericProcess* Parameters::process() const { return process_.get(); } std::string Parameters::processName() const { if ( !process_ ) return "no process"; return process_->name(); } void Parameters::setProcess( std::unique_ptr proc ) { process_ = std::move( proc ); } void Parameters::setProcess( proc::GenericProcess* proc ) { if ( !proc ) throw CG_FATAL( "Parameters" ) << "Trying to clone an invalid process!"; process_.reset( proc ); } hadr::GenericHadroniser* Parameters::hadroniser() { return hadroniser_.get(); } std::string Parameters::hadroniserName() const { if ( !hadroniser_ ) return ""; return hadroniser_->name(); } void Parameters::setHadroniser( std::unique_ptr hadr ) { hadroniser_ = std::move( hadr ); } void Parameters::setHadroniser( hadr::GenericHadroniser* hadr ) { hadroniser_.reset( hadr ); } + io::GenericExportHandler* + Parameters::outputModule() + { + return out_module_.get(); + } + + void + Parameters::setOutputModule( std::unique_ptr mod ) + { + out_module_ = std::move( mod ); + } + + void + Parameters::setOutputModule( io::GenericExportHandler* mod ) + { + out_module_.reset( mod ); + } + std::ostream& operator<<( std::ostream& os, const Parameters* param ) { const bool pretty = true; const int wb = 90, wt = 33; os << std::left << "\n" << std::setfill('_') << std::setw( wb+3 ) << "_/¯ PROCESS INFORMATION ¯\\_" << std::setfill( ' ' ) << "\n" << std::right << std::setw( wb ) << std::left << std::endl << std::setw( wt ) << "Process to generate" << ( pretty ? boldify( param->processName().c_str() ) : param->processName() ); if ( param->process_ ) { os << ", " << param->process_->description(); for ( const auto& par : param->process()->parameters().keys() ) if ( par != "mode" ) os << "\n" << std::setw( wt ) << "" << par << ": " << param->process_->parameters().getString( par ); std::ostringstream proc_mode; proc_mode << param->kinematics.mode; if ( param->kinematics.mode != KinematicsMode::invalid ) os << "\n" << std::setw( wt ) << "Subprocess mode" << ( pretty ? boldify( proc_mode.str().c_str() ) : proc_mode.str() ) << "\n"; } os << "\n" << std::setfill('_') << std::setw( wb+3 ) << "_/¯ RUN INFORMATION ¯\\_" << std::setfill( ' ' ) << "\n" << std::right << std::setw( wb ) << std::left << std::endl << std::setw( wt ) << "Events generation? " << ( pretty ? yesno( param->generation_.enabled ) : std::to_string( param->generation_.enabled ) ) << "\n" << std::setw( wt ) << "Number of events to generate" << ( pretty ? boldify( param->generation_.maxgen ) : std::to_string( param->generation_.maxgen ) ) << "\n"; + if ( param->out_module_ ) + os << std::setw( wt ) << "Output module" << param->out_module_->name() << "\n"; if ( param->generation_.num_threads > 1 ) os << std::setw( wt ) << "Number of threads" << param->generation_.num_threads << "\n"; os << std::setw( wt ) << "Number of points to try per bin" << param->generation_.num_points << "\n" << std::setw( wt ) << "Integrand treatment" << ( pretty ? yesno( param->generation_.treat ) : std::to_string( param->generation_.treat ) ) << "\n" << std::setw( wt ) << "Verbosity level " << utils::Logger::get().level << "\n"; if ( param->hadroniser_ ) { os << "\n" << std::setfill( '-' ) << std::setw( wb+6 ) << ( pretty ? boldify( " Hadronisation algorithm " ) : "Hadronisation algorithm" ) << std::setfill( ' ' ) << "\n\n" << std::setw( wt ) << "Name" << ( pretty ? boldify( param->hadroniser_->name().c_str() ) : param->hadroniser_->name() ) << "\n" << std::setw( wt ) << "Remnants fragmentation? " << ( pretty ? yesno( param->hadroniser_->fragmentRemnants() ) : std::to_string( param->hadroniser_->fragmentRemnants() ) ) << "\n"; } os << "\n" << std::setfill( '-' ) << std::setw( wb+6 ) << ( pretty ? boldify( " Integration parameters " ) : "Integration parameters" ) << std::setfill( ' ' ) << "\n\n"; std::ostringstream int_algo; int_algo << param->integration_.type; os << std::setw( wt ) << "Integration algorithm" << ( pretty ? boldify( int_algo.str().c_str() ) : int_algo.str() ) << "\n" << std::setw( wt ) << "Number of function calls" << param->integration_.ncvg << "\n" << std::setw( wt ) << "Random number generator seed" << param->integration_.rng_seed << "\n"; if ( param->integration_.rng_engine ) os << std::setw( wt ) << "Random number generator engine" << param->integration_.rng_engine->name << "\n"; os << "\n" << std::setfill('_') << std::setw( wb+3 ) << "_/¯ EVENTS KINEMATICS ¯\\_" << std::setfill( ' ' ) << "\n\n" << std::setw( wt ) << "Incoming particles" << param->kinematics.incoming_beams.first << ",\n" << std::setw( wt ) << "" << param->kinematics.incoming_beams.second << "\n" << std::setw( wt ) << "C.m. energy (GeV)" << param->kinematics.sqrtS() << "\n"; if ( param->kinematics.mode != KinematicsMode::ElasticElastic ) os << std::setw( wt ) << "Structure functions" << *param->kinematics.structure_functions << "\n"; os << "\n" << std::setfill( '-' ) << std::setw( wb+6 ) << ( pretty ? boldify( " Incoming partons " ) : "Incoming partons" ) << std::setfill( ' ' ) << "\n\n"; for ( const auto& lim : param->kinematics.cuts.initial.list() ) // map(particles class, limits) if ( lim.second.valid() ) os << std::setw( wt ) << lim.first << lim.second << "\n"; os << "\n" << std::setfill( '-' ) << std::setw( wb+6 ) << ( pretty ? boldify( " Outgoing central system " ) : "Outgoing central system" ) << std::setfill( ' ' ) << "\n\n"; for ( const auto& lim : param->kinematics.cuts.central.list() ) if ( lim.second.valid() ) os << std::setw( wt ) << lim.first << lim.second << "\n"; if ( param->kinematics.cuts.central_particles.size() > 0 ) { os << std::setw( wt ) << ( pretty ? boldify( ">>> per-particle cuts:" ) : ">>> per-particle cuts:" ) << "\n"; for ( const auto& part_per_lim : param->kinematics.cuts.central_particles ) { os << " * all single " << std::setw( wt-3 ) << PDG::get().name( part_per_lim.first ) << "\n"; for ( const auto& lim : part_per_lim.second.list() ) if ( lim.second.valid() ) os << " - " << std::setw( wt-5 ) << lim.first << lim.second << "\n"; } } os << "\n"; os << std::setfill( '-' ) << std::setw( wb+6 ) << ( pretty ? boldify( " Proton / remnants " ) : "Proton / remnants" ) << std::setfill( ' ' ) << "\n"; for ( const auto& lim : param->kinematics.cuts.remnants.list() ) os << "\n" << std::setw( wt ) << lim.first << lim.second; os << "\n"; return os; } //----------------------------------------------------------------------------------------------- Parameters::Integration::Integration() : type( IntegratorType::Vegas ), ncvg( 50000 ), rng_seed( 0 ), rng_engine( (gsl_rng_type*)gsl_rng_mt19937 ), vegas_chisq_cut( 1.5 ), result( -1. ), err_result( -1. ) { const size_t ndof = 10; // random number of dimensions for VEGAS parameters retrieval { std::shared_ptr tmp_state( gsl_monte_vegas_alloc( ndof ), gsl_monte_vegas_free ); gsl_monte_vegas_params_get( tmp_state.get(), &vegas ); vegas.iterations = 10; } { std::shared_ptr tmp_state( gsl_monte_miser_alloc( ndof ), gsl_monte_miser_free ); gsl_monte_miser_params_get( tmp_state.get(), &miser ); } } Parameters::Integration::Integration( const Integration& rhs ) : type( rhs.type ), ncvg( rhs.ncvg ), rng_seed( rhs.rng_seed ), rng_engine( rhs.rng_engine ), vegas( rhs.vegas ), vegas_chisq_cut( rhs.vegas_chisq_cut ), miser( rhs.miser ), result( -1. ), err_result( -1. ) {} Parameters::Integration::~Integration() { //if ( vegas.ostream && vegas.ostream != stdout && vegas.ostream != stderr ) // fclose( vegas.ostream ); } //----------------------------------------------------------------------------------------------- Parameters::Generation::Generation() : enabled( false ), maxgen( 0 ), symmetrise( false ), treat( true ), gen_print_every( 10000 ), num_threads( 2 ), num_points( 100 ) {} Parameters::Generation::Generation( const Generation& rhs ) : enabled( rhs.enabled ), maxgen( rhs.maxgen ), symmetrise( rhs.symmetrise ), treat( rhs.treat ), gen_print_every( rhs.gen_print_every ), num_threads( rhs.num_threads ), num_points( rhs.num_points ) {} } diff --git a/CepGen/Core/Version.cpp b/CepGen/Core/Version.cpp index 868096f..432ea30 100644 --- a/CepGen/Core/Version.cpp +++ b/CepGen/Core/Version.cpp @@ -1,12 +1,12 @@ #include "CepGen/Version.h" #include "CepGen/Core/utils.h" namespace cepgen { const std::string version() { - return Form( "%02u.%02u.%02u", ( cepgen_version >> 16 ) & 0xff, - ( cepgen_version >> 8 ) & 0xff, - cepgen_version & 0xff ); + return Form( "%u.%u.%u", ( cepgen_version >> 16 ) & 0xff, + ( cepgen_version >> 8 ) & 0xff, + cepgen_version & 0xff ); } } diff --git a/CepGen/Core/utils.h b/CepGen/Core/utils.h index 13f287d..6482326 100644 --- a/CepGen/Core/utils.h +++ b/CepGen/Core/utils.h @@ -1,61 +1,63 @@ #ifndef CepGen_Core_utils_h #define CepGen_Core_utils_h #include #include #include namespace cepgen { /// Format a string using a printf style format descriptor. std::string Form( const std::string fmt, ... ); /// Human-readable boolean printout inline const char* yesno( const bool& test ) { return ( test ) ? "\033[32;1myes\033[0m" : "\033[31;1mno\033[0m"; } //inline const char* boldify( const char* str ) { const std::string out = std::string( "\033[33;1m" ) + std::string( str ) + std::string( "\033[0m" ); return out.c_str(); } /// Boldify a string for TTY-type output streams inline std::string boldify( const std::string& str ) { return Form( "\033[1m%s\033[0m", str.c_str() ); } /// Boldify a string for TTY-type output streams inline std::string boldify( const char* str ) { return boldify( std::string( str ) ); } /// Boldify a double floating point number for TTY-type output streams inline std::string boldify( const double& dbl ) { return boldify( Form("%.2f", dbl ) ); } /// Boldify an integer for TTY-type output streams inline std::string boldify( const int& i ) { return boldify( Form("% d", i ) ); } /// Boldify an unsigned integer for TTY-type output streams inline std::string boldify( const unsigned int& ui ) { return boldify( Form("%d", ui ) ); } /// Boldify an unsigned long integer for TTY-type output streams inline std::string boldify( const unsigned long& ui ) { return boldify( Form("%lu", ui ) ); } /// TTY-type enumeration of colours enum class Colour { gray = 30, red = 31, green = 32, yellow = 33, blue = 34, purple = 35 }; /// Colourise a string for TTY-type output streams inline std::string colourise( const std::string& str, const Colour& col ) { return Form( "\033[%d%s\033[0m", (int)col, str.c_str() ); } /// Replace all occurences of a text by another size_t replace_all( std::string& str, const std::string& from, const std::string& to ); namespace utils { /// Add a trailing "s" when needed - inline const char* s( unsigned short num ) { return ( num > 1 ) ? "s" : ""; } + inline const char* s( size_t num ) { return ( num > 1 ) ? "s" : ""; } + /// Add a trailing "s" when needed + inline std::string s( const std::string& word, size_t num ) { return Form( "%i %s%s", num, word.c_str(), ( num > 1 ) ? "s" : "" ); } /// Helper to print a vector template std::string repr( const std::vector& vec, const std::string& sep = "," ) { return std::accumulate( std::next( vec.begin() ), vec.end(), std::to_string( *vec.begin() ), [&sep]( std::string str, T xv ) { return std::move( str )+sep+std::to_string( xv ); } ); } class ProgressBar { public: ProgressBar( size_t tot, size_t freq = 10 ); void update( size_t iter ) const; private: static constexpr size_t BAR_LENGTH = 50; const std::string bar_pattern_; size_t total_, frequency_; }; } } /// Provide a random number generated along a uniform distribution between 0 and 1 #define drand() (double)rand()/RAND_MAX #endif diff --git a/CepGen/Event/Event.cpp b/CepGen/Event/Event.cpp index 467d8b1..89c48e2 100644 --- a/CepGen/Event/Event.cpp +++ b/CepGen/Event/Event.cpp @@ -1,356 +1,356 @@ #include "CepGen/Event/Event.h" #include "CepGen/Physics/PDG.h" #include "CepGen/Physics/HeavyIon.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/utils.h" #include #include namespace cepgen { Event::Event() : num_hadronisation_trials( 0 ), time_generation( -1. ), time_total( -1. ) {} Event::Event( const Event& rhs ) : num_hadronisation_trials( rhs.num_hadronisation_trials ), time_generation( rhs.time_generation ), time_total( rhs.time_total ), particles_( rhs.particles_ ), evtcontent_( rhs.evtcontent_ ) {} void Event::clear() { particles_.clear(); time_generation = -1.; time_total = -1.; } void Event::freeze() { //--- store a snapshot of the primordial event block if ( particles_.count( Particle::CentralSystem ) > 0 ) evtcontent_.cs = particles_[Particle::CentralSystem].size(); if ( particles_.count( Particle::OutgoingBeam1 ) > 0 ) evtcontent_.op1 = particles_[Particle::OutgoingBeam1].size(); if ( particles_.count( Particle::OutgoingBeam2 ) > 0 ) evtcontent_.op2 = particles_[Particle::OutgoingBeam2].size(); } void Event::restore() { //--- remove all particles after the primordial event block if ( particles_.count( Particle::CentralSystem ) > 0 ) particles_[Particle::CentralSystem].resize( evtcontent_.cs ); if ( particles_.count( Particle::OutgoingBeam1 ) > 0 ) particles_[Particle::OutgoingBeam1].resize( evtcontent_.op1 ); if ( particles_.count( Particle::OutgoingBeam2 ) > 0 ) particles_[Particle::OutgoingBeam2].resize( evtcontent_.op2 ); } double Event::cmEnergy() const { return CMEnergy( getOneByRole( Particle::IncomingBeam1 ), getOneByRole( Particle::IncomingBeam2 ) ); } Particles& Event::operator[]( Particle::Role role ) { //--- retrieve all particles with a given role return particles_[role]; } const Particles& Event::operator[]( Particle::Role role ) const { if ( particles_.count( role ) == 0 ) throw CG_FATAL( "Event" ) << "Failed to retrieve a particle with " << role << " role."; //--- retrieve all particles with a given role return particles_.at( role ); } ParticlesIds Event::ids( Particle::Role role ) const { ParticlesIds out; //--- retrieve all particles ids with a given role if ( particles_.count( role ) == 0 ) return out; for ( const auto& part : particles_.at( role ) ) out.insert( part.id() ); return out; } Particle& Event::getOneByRole( Particle::Role role ) { //--- retrieve the first particle of a given role Particles& parts_by_role = operator[]( role ); if ( parts_by_role.empty() ) throw CG_FATAL( "Event" ) << "No particle retrieved with " << role << " role."; if ( parts_by_role.size() > 1 ) throw CG_FATAL( "Event" ) << "More than one particle with " << role << " role: " << parts_by_role.size() << " particles."; return *parts_by_role.begin(); } const Particle& Event::getOneByRole( Particle::Role role ) const { //--- retrieve the first particle of a given role const Particles& parts_by_role = operator[]( role ); if ( parts_by_role.empty() ) throw CG_FATAL( "Event" ) << "No particle retrieved with " << role << " role."; if ( parts_by_role.size() > 1 ) throw CG_FATAL( "Event" ) << "More than one particle with " << role << " role: " << parts_by_role.size() << " particles"; return *parts_by_role.begin(); } Particle& Event::operator[]( int id ) { for ( auto& role_part : particles_ ) for ( auto& part : role_part.second ) if ( part.id() == id ) return part; throw CG_FATAL( "Event" ) << "Failed to retrieve the particle with id=" << id << "."; } const Particle& Event::operator[]( int id ) const { for ( const auto& role_part : particles_ ) for ( const auto& part : role_part.second ) if ( part.id() == id ) return part; throw CG_FATAL( "Event" ) << "Failed to retrieve the particle with id=" << id << "."; } Particles Event::getByIds( const ParticlesIds& ids ) const { Particles out; for ( const auto& id : ids ) out.emplace_back( operator[]( id ) ); return out; } Particles Event::mothers( const Particle& part ) const { return getByIds( part.mothers() ); } Particles Event::daughters( const Particle& part ) const { return getByIds( part.daughters() ); } ParticleRoles Event::roles() const { ParticleRoles out; for ( const auto& pr : particles_ ) out.emplace_back( pr.first ); return out; } Particle& Event::addParticle( Particle& part, bool replace ) { CG_DEBUG_LOOP( "Event" ) << "Particle with PDGid = " << part.integerPdgId() << " has role " << part.role(); if ( part.role() <= 0 ) throw CG_FATAL( "Event" ) << "Trying to add a particle with role=" << (int)part.role() << "."; //--- retrieve the list of particles with the same role Particles& part_with_same_role = operator[]( part.role() ); //--- specify the id - if ( part_with_same_role.empty() && part.id() < 0 ) part.setId( numParticles() ); // set the id if previously invalid/inexistent + if ( part_with_same_role.empty() && part.id() < 0 ) part.setId( size() ); // set the id if previously invalid/inexistent if ( !part_with_same_role.empty() ) { if ( replace ) part.setId( part_with_same_role[0].id() ); // set the previous id if replacing a particle - else part.setId( numParticles() ); + else part.setId( size() ); } //--- add the particle to the collection if ( replace ) part_with_same_role = Particles( 1, part ); // generate a vector containing only this particle else part_with_same_role.emplace_back( part ); return part_with_same_role.back(); } Particle& Event::addParticle( Particle::Role role, bool replace ) { Particle np( role, PDG::invalid ); return addParticle( np, replace ); } size_t - Event::numParticles() const + Event::size() const { size_t out = 0; for ( const auto& role_part : particles_ ) out += role_part.second.size(); return out; } const Particles Event::particles() const { Particles out; for ( const auto& role_part : particles_ ) out.insert( out.end(), role_part.second.begin(), role_part.second.end() ); std::sort( out.begin(), out.end() ); return out; } const Particles Event::stableParticles() const { Particles out; for ( const auto& role_part : particles_ ) for ( const auto& part : role_part.second ) if ( (short)part.status() > 0 ) out.emplace_back( part ); std::sort( out.begin(), out.end() ); return out; } void Event::checkKinematics() const { // check the kinematics through parentage for ( const auto& part : particles() ) { ParticlesIds daughters = part.daughters(); if ( daughters.empty() ) continue; Particle::Momentum ptot; for ( const auto& daugh : daughters ) { const Particle& d = operator[]( daugh ); const ParticlesIds mothers = d.mothers(); ptot += d.momentum(); if ( mothers.size() < 2 ) continue; for ( const auto& moth : mothers ) if ( moth != part.id() ) ptot -= operator[]( moth ).momentum(); } const double mass_diff = ( ptot-part.momentum() ).mass(); if ( fabs( mass_diff ) > MIN_PRECISION ) { dump(); throw CG_FATAL( "Event" ) << "Error in momentum balance for particle " << part.id() << ": mdiff = " << mass_diff << "."; } } } void Event::dump( bool stable ) const { const Particles parts = ( stable ) ? stableParticles() : particles(); std::ostringstream os; Particle::Momentum p_total; for ( const auto& part : parts ) { const ParticlesIds mothers = part.mothers(); { std::ostringstream oss_pdg; if ( part.pdgId() == PDG::invalid && !mothers.empty() ) { //--- if particles compound std::string delim; for ( unsigned short i = 0; i < mothers.size(); ++i ) try { oss_pdg << delim << PDG::get().name( operator[]( *std::next( mothers.begin(), i ) ).pdgId() ), delim = "/"; } catch ( const Exception& ) { oss_pdg << delim << operator[]( *std::next( mothers.begin(), i ) ).pdgId(), delim = "/"; } os << Form( "\n %2d\t\t %-7s", part.id(), oss_pdg.str().c_str() ); } else { //--- if single particle/HI if ( (HeavyIon)part.pdgId() ) oss_pdg << (HeavyIon)part.pdgId(); else try { oss_pdg << PDG::get().name( part.pdgId() ); } catch ( const Exception& ) { oss_pdg << "?"; } os << Form( "\n %2d\t%-+10d %-7s", part.id(), part.integerPdgId(), oss_pdg.str().c_str() ); } } os << "\t"; if ( part.charge() != (int)part.charge() ) { if ( part.charge()*2 == (int)( part.charge()*2 ) ) os << Form( "%-d/2", (int)( part.charge()*2 ) ); else if ( part.charge()*3 == (int)( part.charge()*3 ) ) os << Form( "%-d/3", (int)( part.charge()*3 ) ); else os << Form( "%-.2f", part.charge() ); } else os << Form( "%-g", part.charge() ); os << "\t"; { std::ostringstream oss; oss << part.role(); os << Form( "%-8s %6d\t", oss.str().c_str(), part.status() ); } if ( !mothers.empty() ) { std::ostringstream oss; unsigned short i = 0; for ( const auto& moth : mothers ) { oss << ( i > 0 ? "+" : "" ) << moth; ++i; } os << Form( "%6s ", oss.str().c_str() ); } else os << " "; const Particle::Momentum mom = part.momentum(); os << Form( "% 9.6e % 9.6e % 9.6e % 9.6e % 12.5f", mom.px(), mom.py(), mom.pz(), part.energy(), part.mass() ); // discard non-primary, decayed particles if ( part.status() >= Particle::Status::Undefined ) { const int sign = ( part.status() == Particle::Status::Undefined ) ? -1 : +1; p_total += sign*mom; } } //--- set a threshold to the computation precision p_total.truncate(); // CG_INFO( "Event" ) << Form( "Dump of event content:\n" " Id\tPDG id\t Name\t\tCharge\tRole\t Status\tMother\tpx py pz E \t M \n" " --\t------\t ----\t\t------\t----\t ------\t------\t----GeV/c--- ----GeV/c--- ----GeV/c--- ----GeV/c---\t --GeV/c²--" "%s\n" " ----------------------------------------------------------------------------------------------------------------------------------\n" "\t\t\t\t\t\t\tBalance% 9.6e % 9.6e % 9.6e % 9.6e", os.str().c_str(), p_total.px(), p_total.py(), p_total.pz(), p_total.energy() ); } //------------------------------------------------------------------------------------------------ Event::NumParticles::NumParticles() : cs( 0 ), op1( 0 ), op2( 0 ) {} Event::NumParticles::NumParticles( const NumParticles& np ) : cs( np.cs ), op1( np.op1 ), op2( np.op2 ) {} } diff --git a/CepGen/Event/Event.h b/CepGen/Event/Event.h index 6ba640c..47bcb2f 100644 --- a/CepGen/Event/Event.h +++ b/CepGen/Event/Event.h @@ -1,124 +1,124 @@ #ifndef CepGen_Event_Event_h #define CepGen_Event_Event_h #include "CepGen/Event/Particle.h" namespace cepgen { /** * Class containing all the information on the in- and outgoing particles' kinematics * \brief Kinematic information on the particles in the event */ class Event { public: /// Build an empty event Event(); /// Copy constructor Event( const Event& ); /// Empty the whole event content void clear(); /// Initialize an "empty" event collection void freeze(); /// Restore the event to its "empty" state void restore(); /// Dump all the known information on every Particle object contained in this Event container in the output stream /// \param[out] os Output stream where to dump the information /// \param[in] stable_ Do we only show the stable particles in this event? void dump( bool stable_ = false ) const; /// Incoming beams centre-of-mass energy, in GeV double cmEnergy() const; //----- particles adders /// \brief Set the information on one particle in the process /// \param[in] part The Particle object to insert or modify in the event /// \param[in] replace Do we replace the particle if already present in the event or do we append another particle with the same role ? Particle& addParticle( Particle& part, bool replace = false ); /// \brief Create a new particle in the event, with no kinematic information but the role it has to play in the process /// \param[in] role The role the particle will play in the process /// \param[in] replace Do we replace the particle if already present in the event or do we append another particle with the same role ? Particle& addParticle( Particle::Role role, bool replace = false ); //----- particles retrievers /// Number of particles in the event - size_t numParticles() const; + size_t size() const; /// Vector of all particles in the event const Particles particles() const; /// Vector of all stable particles in the event const Particles stableParticles() const; /** Get a list of Particle objects corresponding to a certain role in the process kinematics * \param[in] role The role the particles have to play in the process * \return A vector of references to the requested Particle objects */ Particles& operator[]( Particle::Role role ); /// Get a list of constant Particle objects corresponding to a certain role in the process kinematics const Particles& operator[]( Particle::Role role ) const; /// Get a list of particle identifiers in Event corresponding to a certain role in the process kinematics ParticlesIds ids( Particle::Role role ) const; /** \brief Get the first Particle object in the particles list whose role corresponds to the given argument * \param[in] role The role the particle has to play in the event * \return A Particle object corresponding to the first particle with the role */ Particle& getOneByRole( Particle::Role role ); const Particle& getOneByRole( Particle::Role role ) const; /** \brief Get the reference to the Particle object corresponding to a unique identifier in the event * \param[in] id The unique identifier to this particle in the event * \return A reference to the requested Particle object */ Particle& operator[]( int id ); /** \brief Get a const Particle object using its unique identifier * \param[in] id Unique identifier of the particle in the event * \return Constant object to be retrieved */ const Particle& operator[]( int id ) const; /** \brief Get references to the Particle objects corresponding to the unique identifiers in the event * \param[in] ids_ The unique identifiers to the particles to be selected in the event * \return A vector of references to the requested Particle objects */ Particles getByIds( const ParticlesIds& ids_ ) const; /** \brief Get the list of mother particles of any given Particle object in this event * \param[in] part The reference to the Particle object from which we want to extract the mother particles * \return A list of parenting Particle object */ //----- general particles information retriever Particles mothers( const Particle& part ) const; /// Get a vector containing all the daughters from a particle /// \param[in] part The particle for which the daughter particles have to be retrieved /// \return Vector of Particle objects containing all the daughters' kinematic information Particles daughters( const Particle& part ) const; /// Get a list of roles for the given event (really process-dependant for the central system) /// \return Vector of integers corresponding to all the roles the particles can play in the event ParticleRoles roles() const; /// Number of trials before the event was "correctly" hadronised unsigned short num_hadronisation_trials; /// Time needed to generate the event at parton level (in seconds) float time_generation; /// Time needed to generate the hadronised (if needed) event (in seconds) float time_total; private: static constexpr double MIN_PRECISION = 1.e-10; /// Check if the event kinematics is properly defined void checkKinematics() const; /// List of particles in the event, mapped to their role in the process ParticlesMap particles_; /// Typical event indices structure struct NumParticles { NumParticles(); NumParticles( const NumParticles& np ); unsigned short cs; ///< Index of the first central system particle unsigned short op1; ///< Index of the first positive-z outgoing beam state unsigned short op2; ///< Index of the first negative-z outgoing beam state }; /// Event indices structure NumParticles evtcontent_; }; } #endif diff --git a/CepGen/Event/EventBrowser.cpp b/CepGen/Event/EventBrowser.cpp new file mode 100644 index 0000000..365f499 --- /dev/null +++ b/CepGen/Event/EventBrowser.cpp @@ -0,0 +1,96 @@ +#include "CepGen/Event/EventBrowser.h" +#include "CepGen/Event/Event.h" + +#include "CepGen/Core/Exception.h" +#include "CepGen/Core/utils.h" + +namespace cepgen +{ + namespace utils + { + const std::regex EventBrowser::rgx_select_id_( "(\\w+)\\((\\d+)\\)" ); + const std::regex EventBrowser::rgx_select_role_( "(\\w+)\\(([a-z]+\\d?)\\)" ); + + double + EventBrowser::get( const Event& ev, const std::string& var ) const + { + std::smatch sm; + //--- particle-level variables (indexed by integer id) + if ( std::regex_match( var, sm, rgx_select_id_ ) ) { + const auto& var_name = sm[1].str(); + const auto& part = ev[std::stod( sm[2].str() )]; + return variable( ev, part, var_name ); + } + //--- particle-level variables (indexed by role) + else if ( std::regex_match( var, sm, rgx_select_role_ ) ) { + const auto& var_name = sm[1].str(); + 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 INVALID_OUTPUT; + } + const auto& part = ev[role_str_.at( str_role )][0]; + return variable( ev, part, var_name ); + } + //--- event-level variables + else + return variable( ev, var ); + } + + double + EventBrowser::variable( const Event& ev, 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" ) { + const auto& moth = part.mothers(); + if ( moth.empty() ) { + CG_WARNING( "EventBrowser" ) + << "Failed to retrieve parent particle to compute xi " + << "for the following particle:\n" + << part; + return INVALID_OUTPUT; + } + return 1.-part.energy()/ev[*moth.begin()].energy(); + } + if ( var == "pdg" ) return (double)part.integerPdgId(); + if ( var == "charge" ) return part.charge(); + if ( var == "status" ) return (double)part.status(); + CG_WARNING( "EventBrowser" ) + << "Failed to retrieve variable \"" << var << "\"."; + return INVALID_OUTPUT; + } + + double + EventBrowser::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( "EventBrowser" ) + << "Failed to retrieve the event-level variable \"" << var << "\"."; + return INVALID_OUTPUT; + } + } +} + diff --git a/CepGen/Event/EventBrowser.h b/CepGen/Event/EventBrowser.h new file mode 100644 index 0000000..0366bdf --- /dev/null +++ b/CepGen/Event/EventBrowser.h @@ -0,0 +1,61 @@ +#ifndef CepGen_Event_EventBrowser_h +#define CepGen_Event_EventBrowser_h + +#include "CepGen/Event/Particle.h" +#include + +namespace cepgen +{ + class Event; + namespace utils + { + /** + * \brief A user-friendly browser for the Event content + * \author Laurent Forthomme + * \date Jul 2019 + */ + class EventBrowser + { + public: + EventBrowser() = default; + double get( const Event& ev, const std::string& var ) const; + + private: + /// Retrieve a named variable from a particle + double variable( const Event&, 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.; + + //--- 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 } + }; + }; + } +} + +#endif + diff --git a/CepGen/Event/Particle.cpp b/CepGen/Event/Particle.cpp index 4dd99a0..ca67386 100644 --- a/CepGen/Event/Particle.cpp +++ b/CepGen/Event/Particle.cpp @@ -1,286 +1,284 @@ #include "CepGen/Event/Particle.h" #include "CepGen/Physics/PDG.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/utils.h" #include namespace cepgen { Particle::Particle() : id_( -1 ), charge_sign_( 1 ), mass_( -1. ), helicity_( 0. ), role_( UnknownRole ), status_( Status::Undefined ), pdg_id_( PDG::invalid ) {} Particle::Particle( Role role, pdgid_t pdgId, Status st ) : id_( -1 ), charge_sign_( 1 ), mass_( -1. ), helicity_( 0. ), role_( role ), status_( st ), pdg_id_( pdgId ) { try { phys_prop_ = PDG::get()( pdg_id_ ); } catch ( const Exception& ) {} if ( pdg_id_ != PDG::invalid ) computeMass(); } Particle::Particle( const Particle& part ) : id_( part.id_ ), charge_sign_( part.charge_sign_ ), momentum_( part.momentum_ ), mass_( part.mass_ ), helicity_( part.helicity_ ), role_( part.role_ ), status_( part.status_ ), mothers_( part.mothers_ ), daughters_( part.daughters_ ), pdg_id_( part.pdg_id_ ) { try { phys_prop_ = PDG::get()( pdg_id_ ); } catch ( const Exception& ) {} } bool Particle::operator<( const Particle& rhs ) const { return id_ >= 0 && rhs.id_ > 0 && id_ < rhs.id_; } double Particle::thetaToEta( double theta ) { return -log( tan( 0.5*theta*M_PI/180. ) ); } double Particle::etaToTheta( double eta ) { return 2.*atan( exp( -eta ) )*180.*M_1_PI; } bool Particle::valid() { if ( pdg_id_ == PDG::invalid ) return false; if ( momentum_.p() == 0. && mass_ == 0. ) return false; return true; } float Particle::charge() const { return charge_sign_ * phys_prop_.charge/3.; } void Particle::computeMass( bool off_shell ) { if ( !off_shell && pdg_id_ != PDG::invalid ) { // retrieve the mass from the on-shell particle's properties mass_ = phys_prop_.mass; } else if ( momentum_.energy() >= 0. ) { mass_ = sqrt( energy2() - momentum_.p2() ); } //--- finish by setting the energy accordingly if ( momentum_.energy() < 0. ) { // invalid energy momentum_.setEnergy( sqrt( momentum_.p2() + mass2() ) ); } } void Particle::setMass( double m ) { if ( m >= 0. ) mass_ = m; else computeMass(); } void Particle::addMother( Particle& part ) { mothers_.insert( part.id() ); CG_DEBUG_LOOP( "Particle" ) << "Particle " << id() << " (pdgId=" << part.integerPdgId() << ") " << "is the new mother of " << id_ << " (pdgId=" << (int)pdg_id_ << ")."; part.addDaughter( *this ); } void Particle::addDaughter( Particle& part ) { const auto ret = daughters_.insert( part.id() ); if ( CG_LOG_MATCH( "Particle", debugInsideLoop ) ) { std::ostringstream os; for ( const auto& daugh : daughters_ ) os << Form( "\n\t * id=%d", daugh ); CG_DEBUG_LOOP( "Particle" ) << "Particle " << role_ << " (pdgId=" << (int)pdg_id_ << ") " << "has now " << daughters_.size() << " daughter(s):" << os.str(); } if ( ret.second ) { CG_DEBUG_LOOP( "Particle" ) << "Particle " << part.role() << " (pdgId=" << part.integerPdgId() << ") " << "is a new daughter of " << role_ << " (pdgId=" << pdg_id_ << ")."; if ( part.mothers().find( id_ ) == part.mothers().end() ) part.addMother( *this ); } } void Particle::setMomentum( const Momentum& mom, bool offshell ) { momentum_ = mom; if ( !offshell && mom.mass() > 0. ) mass_ = momentum_.mass(); else computeMass(); } void Particle::setMomentum( double px, double py, double pz, double e ) { momentum_.setP( px, py, pz ); setEnergy( e ); if ( fabs( e-momentum_.energy() ) > 1.e-6 ) // more than 1 eV difference CG_WARNING( "Particle" ) << "Energy difference: " << e-momentum_.energy(); } double Particle::energy() const { return ( momentum_.energy() < 0. ? std::hypot( mass_, momentum_.p() ) : momentum_.energy() ); } void Particle::setEnergy( double e ) { if ( e < 0. && mass_ >= 0. ) e = std::hypot( mass_, momentum_.p() ); momentum_.setEnergy( e ); } pdgid_t Particle::pdgId() const { return pdg_id_; } void Particle::setPdgId( short pdg ) { pdg_id_ = abs( pdg ); try { phys_prop_ = PDG::get()( pdg_id_ ); } catch ( const Exception& ) {} switch ( pdg_id_ ) { case PDG::electron: case PDG::muon: case PDG::tau: charge_sign_ = -pdg/abs( pdg ); break; default: charge_sign_ = pdg/abs( pdg ); break; } } void Particle::setPdgId( pdgid_t pdg, short ch ) { pdg_id_ = pdg; switch ( pdg_id_ ) { case PDG::electron: case PDG::muon: case PDG::tau: charge_sign_ = -ch; break; default: charge_sign_ = ch; break; } } int Particle::integerPdgId() const { const float ch = phys_prop_.charge/3.; if ( ch == 0 ) return static_cast( pdg_id_ ); return static_cast( pdg_id_ ) * charge_sign_ * ( ch/fabs( ch ) ); } - void - Particle::dump() const + double + Particle::etaToY( double eta_, double m_, double pt_ ) + { + const double m2 = m_*m_, mt = std::hypot( m_, pt_ ); + return asinh( sqrt( ( ( mt*mt-m2 )*cosh( 2.*eta_ )+m2 )/ mt*mt - 1. )*M_SQRT1_2 ); + } + + std::ostream& + operator<<( std::ostream& os, const Particle& part ) { - std::ostringstream os; os << std::resetiosflags( std::ios::showbase ) - << "Particle[" << id_ << "]{role=" << role_ << ", status=" << (int)status_ << ", " - << "pdg=" << pdg_id_ << ", p4=" << momentum_ << " GeV, m=" << mass_ << " GeV, " - << "p⟂=" << momentum_.pt() << " GeV, eta=" << momentum_.eta() << ", phi=" << momentum_.phi(); - if ( primary() ) + << "Particle[" << part.id_ << "]{role=" << part.role_ << ", status=" << (int)part.status_ << ", " + << "pdg=" << part.pdg_id_ << ", p4=" << part.momentum_ << " GeV, m=" << part.mass_ << " GeV, " + << "p⟂=" << part.momentum_.pt() << " GeV, eta=" << part.momentum_.eta() << ", phi=" << part.momentum_.phi(); + if ( part.primary() ) os << ", primary"; else { - os << ", mother" << utils::s( mothers_.size() ) << "="; + os << ", " << utils::s( "mother", part.mothers_.size() ) << "="; std::string delim; - for ( const auto& moth : mothers_ ) + for ( const auto& moth : part.mothers_ ) os << delim << moth, delim = ","; } - const auto& daughters_list = daughters(); + const auto& daughters_list = part.daughters(); if ( !daughters_list.empty() ) { - os << ", daughter" << utils::s( daughters_list.size() ) << "="; + os << ", " << utils::s( "daughter", daughters_list.size() ) << "="; std::string delim; for ( const auto& daugh : daughters_list ) os << delim << daugh, delim = ","; } - os << "}"; - CG_INFO( "Particle" ) << os.str(); - } - - double - Particle::etaToY( double eta_, double m_, double pt_ ) - { - const double m2 = m_*m_, mt = std::hypot( m_, pt_ ); - return asinh( sqrt( ( ( mt*mt-m2 )*cosh( 2.*eta_ )+m2 )/ mt*mt - 1. )*M_SQRT1_2 ); + return os << "}"; } std::ostream& operator<<( std::ostream& os, const Particle::Role& rl ) { switch ( rl ) { case Particle::UnknownRole: return os << "unknown"; case Particle::IncomingBeam1: return os << "i.beam 1"; case Particle::IncomingBeam2: return os << "i.beam 2"; case Particle::OutgoingBeam1: return os << "o.beam 1"; case Particle::OutgoingBeam2: return os << "o.beam 2"; case Particle::Parton1: return os << "parton 1"; case Particle::Parton2: return os << "parton 2"; case Particle::Intermediate: return os << "hard pr."; case Particle::CentralSystem: return os << "central"; } return os; } double CMEnergy( const Particle& p1, const Particle& p2 ) { if ( p1.mass()*p2.mass() < 0. || p1.energy()*p2.energy() < 0. ) return 0.; return sqrt( p1.mass2()+p2.mass2() + 2.*p1.energy()*p2.energy() - 2.*( p1.momentum()*p2.momentum() ) ); } double CMEnergy( const Particle::Momentum& m1, const Particle::Momentum& m2 ) { if ( m1.mass()*m2.mass() < 0. || m1.energy()*m2.energy() < 0. ) return 0.; return sqrt( m1.mass2()+m2.mass2() + 2.*m1.energy()*m2.energy() - 2.*( m1*m2 ) ); } } diff --git a/CepGen/Event/Particle.h b/CepGen/Event/Particle.h index 70d944d..75c4ecb 100644 --- a/CepGen/Event/Particle.h +++ b/CepGen/Event/Particle.h @@ -1,355 +1,355 @@ #ifndef CepGen_Event_Particle_h #define CepGen_Event_Particle_h #include "CepGen/Core/Hasher.h" #include "CepGen/Physics/Constants.h" #include "CepGen/Physics/ParticleProperties.h" #include #include #include namespace cepgen { /// A set of integer-type particle identifiers typedef std::set ParticlesIds; /// Kinematic information for one particle class Particle { public: /// Internal status code for a particle enum class Status { PrimordialIncoming = -9, ///< Incoming beam particle DebugResonance = -5, ///< Intermediate resonance (for processes developers) Resonance = -4, ///< Already decayed intermediate resonance Fragmented = -3, ///< Already fragmented outgoing beam Propagator = -2, ///< Generic propagator Incoming = -1, ///< Incoming parton Undefined = 0, ///< Undefined particle FinalState = 1, ///< Stable, final state particle Undecayed = 2, ///< Particle to be decayed externally Unfragmented = 3 ///< Particle to be hadronised externally }; /// Role of the particle in the process enum Role { UnknownRole = -1, ///< Undefined role IncomingBeam1 = 1, ///< \f$z>0\f$ incoming beam particle IncomingBeam2 = 2, ///< \f$z<0\f$ incoming beam particle OutgoingBeam1 = 3, ///< \f$z<0\f$ outgoing beam state/particle OutgoingBeam2 = 5, ///< \f$z>0\f$ outgoing beam state/particle CentralSystem = 6, ///< Central particles system Intermediate = 4, ///< Intermediate two-parton system Parton1 = 41, ///< \f$z>0\f$ beam incoming parton Parton2 = 42 ///< \f$z<0\f$ beam incoming parton }; /** * Container for a particle's 4-momentum, along with useful methods to ease the development of any matrix element level generator * \brief 4-momentum for a particle * \date Dec 2015 * \author Laurent Forthomme */ class Momentum { public: /// Build a 4-momentum at rest with an invalid energy (no mass information known) Momentum(); /// Build a 4-momentum using its 3-momentum coordinates and its energy Momentum( double x, double y, double z, double t = -1. ); /// Build a 4-momentum using its 3-momentum coordinates and its energy Momentum( double* p ); // --- static definitions /// Build a 3-momentum from its three pseudo-cylindric coordinates static Momentum fromPtEtaPhi( double pt, double eta, double phi, double e = -1. ); /// Build a 4-momentum from its scalar momentum, and its polar and azimuthal angles static Momentum fromPThetaPhi( double p, double theta, double phi, double e = -1. ); /// Build a 4-momentum from its four momentum and energy coordinates static Momentum fromPxPyPzE( double px, double py, double pz, double e ); /// Build a 4-momentum from its three momentum coordinates and mass static Momentum fromPxPyPzM( double px, double py, double pz, double m ); /// Build a 4-momentum from its transverse momentum, rapidity and mass static Momentum fromPxPyYM( double px, double py, double rap, double m ); // --- vector and scalar operators /// Scalar product of the 3-momentum with another 3-momentum double threeProduct( const Momentum& ) const; /// Scalar product of the 4-momentum with another 4-momentum double fourProduct( const Momentum& ) const; /// Vector product of the 3-momentum with another 3-momentum double crossProduct( const Momentum& ) const; /// Compute the 4-vector sum of two 4-momenta Momentum operator+( const Momentum& ) const; /// Add a 4-momentum through a 4-vector sum Momentum& operator+=( const Momentum& ); /// Unary inverse operator Momentum operator-() const; /// Compute the inverse per-coordinate 4-vector Momentum operator-( const Momentum& ) const; /// Subtract a 4-momentum through a 4-vector sum Momentum& operator-=( const Momentum& ); /// Scalar product of two 3-momenta double operator*( const Momentum& ) const; /// Scalar product of the 3-momentum with another 3-momentum double operator*=( const Momentum& ); /// Multiply all components of a 4-momentum by a scalar Momentum operator*( double c ) const; /// Multiply all 4-momentum coordinates by a scalar Momentum& operator*=( double c ); /// Left-multiply all 4-momentum coordinates by a scalar friend Momentum operator*( double, const Momentum& ); /// Equality operator bool operator==( const Momentum& ) const; /// Human-readable format for a particle's momentum friend std::ostream& operator<<( std::ostream&, const Momentum& ); Momentum& betaGammaBoost( double gamma, double betagamma ); /// Forward Lorentz boost Momentum& lorentzBoost( const Momentum& p ); // --- setters and getters /// Set all the components of the 4-momentum (in GeV) Momentum& setP( double px, double py, double pz, double e ); /// Set all the components of the 3-momentum (in GeV) Momentum& setP( double px, double py, double pz ); /// Set the energy (in GeV) inline Momentum& setEnergy( double e ) { energy_ = e; return *this; } /// Compute the energy from the mass inline Momentum& setMass( double m ) { return setMass2( m*m ); } /// Compute the energy from the mass Momentum& setMass2( double m2 ); /// Get one component of the 4-momentum (in GeV) double operator[]( const unsigned int i ) const; /// Get one component of the 4-momentum (in GeV) double& operator[]( const unsigned int i ); /// Momentum along the \f$x\f$-axis (in GeV) inline double px() const { return px_; } /// Momentum along the \f$y\f$-axis (in GeV) inline double py() const { return py_; } /// Longitudinal momentum (in GeV) inline double pz() const { return pz_; } /// Transverse momentum (in GeV) double pt() const; /// Squared transverse momentum (in GeV\f$^2\f$) double pt2() const; /// 5-vector of double precision floats (in GeV) const std::vector pVector() const; /// 3-momentum norm (in GeV) inline double p() const { return p_; } /// Squared 3-momentum norm (in GeV\f$^2\f$) inline double p2() const { return p_*p_; } /// Energy (in GeV) inline double energy() const { return energy_; } /// Squared energy (in GeV\f$^2\f$) inline double energy2() const { return energy_*energy_; } /// Squared mass (in GeV\f$^2\f$) as computed from its energy and momentum inline double mass2() const { return energy2()-p2(); } /// Mass (in GeV) as computed from its energy and momentum /// \note Returns \f$-\sqrt{|E^2-\mathbf{p}^2|}<0\f$ if \f$\mathbf{p}^2>E^2\f$ double mass() const; /// Polar angle (angle with respect to the longitudinal direction) double theta() const; /// Azimutal angle (angle in the transverse plane) double phi() const; /// Pseudo-rapidity double eta() const; /// Rapidity double rapidity() const; Momentum& truncate( double tolerance = 1.e-10 ); /// Rotate the transverse components by an angle phi (and reflect the y coordinate) Momentum& rotatePhi( double phi, double sign ); /// Rotate the particle's momentum by a polar/azimuthal angle Momentum& rotateThetaPhi( double theta_, double phi_ ); /// Apply a \f$ z\rightarrow -z\f$ transformation inline Momentum& mirrorZ() { pz_ = -pz_; return *this; } private: /// Compute the 3-momentum's norm Momentum& computeP(); /// Momentum along the \f$x\f$-axis double px_; /// Momentum along the \f$y\f$-axis double py_; /// Momentum along the \f$z\f$-axis double pz_; /// 3-momentum's norm (in GeV/c) double p_; /// Energy (in GeV) double energy_; }; /// Human-readable format for a particle's role in the event friend std::ostream& operator<<( std::ostream& os, const Role& rl ); //----- static getters /// Convert a polar angle to a pseudo-rapidity static double thetaToEta( double theta ); /// Convert a pseudo-rapidity to a polar angle static double etaToTheta( double eta ); /// Convert a pseudo-rapidity to a rapidity static double etaToY( double eta_, double m_, double pt_ ); Particle(); /// Build using the role of the particle in the process and its PDG id /// \param[in] role Role of the particle in the process /// \param[in] id PDG identifier /// \param[in] st Current status Particle( Role role, pdgid_t id, Status st = Status::Undefined ); /// Copy constructor Particle( const Particle& ); inline ~Particle() {} /// Comparison operator (from unique identifier) bool operator<( const Particle& rhs ) const; /// Comparison operator (from their reference's unique identifier) //bool operator<( Particle *rhs ) const { return ( id < rhs->id ); } // --- general particle properties /// Unique identifier (in a Event object context) int id() const { return id_; } //void setId( int id ) { id_ = id; } /// Set the particle unique identifier in an event void setId( int id ) { id_ = id; } /// Electric charge (given as a float number, for the quarks and bound states) float charge() const; /// Set the electric charge sign (+-1 for charged or 0 for neutral particles) void setChargeSign( int sign ) { charge_sign_ = sign; } /// Role in the considered process Role role() const { return role_; } /// Set the particle role in the process void setRole( const Role& role ) { role_ = role; } /** * Codes 1-10 correspond to currently existing partons/particles, and larger codes contain partons/particles which no longer exist, or other kinds of event information * \brief Particle status */ Status status() const { return status_; } /// Set the particle decay/stability status void setStatus( Status status ) { status_ = status; } /// Set the PDG identifier (along with the particle's electric charge) /// \param[in] pdg PDG identifier /// \param[in] ch Electric charge (0, 1, or -1) void setPdgId( pdgid_t pdg, short ch = 0 ); /// Set the PDG identifier (along with the particle's electric charge) /// \param[in] pdg_id PDG identifier (incl. electric charge in e) void setPdgId( short pdg_id ); /// Retrieve the objectified PDG identifier pdgid_t pdgId() const; /// Retrieve the integer value of the PDG identifier int integerPdgId() const; /// Particle's helicity float helicity() const { return helicity_; } /// Set the helicity of the particle void setHelicity( float heli ) { helicity_ = heli; } /// Particle mass in GeV/c\f$^2\f$ /// \return Particle's mass inline double mass() const { return mass_; }; /// Compute the particle mass /// \param[in] off_shell Allow the particle to be produced off-shell? /// \note This method ensures that the kinematics is properly set (the mass is set according to the energy and the momentum in priority) void computeMass( bool off_shell = false ); /// Set the particle mass, in GeV/c\f$^2\f$ /// \param m Mass in GeV/c\f$^2\f$ /// \note This method ensures that the kinematics is properly set (the mass is set according to the energy and the momentum in priority) void setMass( double m = -1. ); /// Particle squared mass, in GeV\f$^2\f$/c\f$^4\f$ inline double mass2() const { return mass_*mass_; }; /// Retrieve the momentum object associated with this particle inline Momentum& momentum() { return momentum_; } /// Retrieve the momentum object associated with this particle inline Momentum momentum() const { return momentum_; } /// Associate a momentum object to this particle void setMomentum( const Momentum& mom, bool offshell = false ); /** * \brief Set the 3- or 4-momentum associated to the particle * \param[in] px Momentum along the \f$x\f$-axis, in GeV/c * \param[in] py Momentum along the \f$y\f$-axis, in GeV/c * \param[in] pz Momentum along the \f$z\f$-axis, in GeV/c * \param[in] e Energy, in GeV */ void setMomentum( double px, double py, double pz, double e = -1. ); /// Set the 4-momentum associated to the particle /// \param[in] p 4-momentum inline void setMomentum( double p[4] ) { setMomentum( p[0], p[1], p[2], p[3] ); } /// Set the particle's energy /// \param[in] e Energy, in GeV void setEnergy( double e = -1. ); /// Get the particle's energy, in GeV double energy() const; /// Get the particle's squared energy, in GeV\f$^2\f$ inline double energy2() const { return energy()*energy(); }; /// Is this particle a valid particle which can be used for kinematic computations? bool valid(); // --- particle relations /// Is this particle a primary particle? inline bool primary() const { return mothers_.empty(); } /// Set the mother particle /// \param[in] part A Particle object containing all the information on the mother particle void addMother( Particle& part ); /// Get the unique identifier to the mother particle from which this particle arises /// \return An integer representing the unique identifier to the mother of this particle in the event inline ParticlesIds mothers() const { return mothers_; } /** * \brief Add a decay product * \param[in] part The Particle object in which this particle will desintegrate or convert * \return A boolean stating if the particle has been added to the daughters list or if it was already present before */ void addDaughter( Particle& part ); /// Gets the number of daughter particles inline unsigned int numDaughters() const { return daughters_.size(); }; /// Get an identifiers list all daughter particles /// \return An integer vector containing all the daughters' unique identifier in the event inline ParticlesIds daughters() const { return daughters_; } // --- global particle information extraction /// Dump all the information on this particle into the standard output stream - void dump() const; + friend std::ostream& operator<<( std::ostream&, const Particle& ); - private: + protected: /// Unique identifier in an event int id_; /// Electric charge (+-1 or 0) short charge_sign_; /// Momentum properties handler Momentum momentum_; /// Mass, in GeV/c\f$^2\f$ double mass_; /// Helicity float helicity_; /// Role in the process Role role_; /// Decay/stability status Status status_; /// List of mother particles ParticlesIds mothers_; /// List of daughter particles ParticlesIds daughters_; /// PDG id pdgid_t pdg_id_; /// Collection of standard, bare-level physical properties ParticleProperties phys_prop_; }; /// Compute the centre of mass energy of two particles (incoming or outgoing states) double CMEnergy( const Particle& p1, const Particle& p2 ); /// Compute the centre of mass energy of two particles (incoming or outgoing states) double CMEnergy( const Particle::Momentum& m1, const Particle::Momentum& m2 ); //bool operator<( const Particle& a, const Particle& b ) { return a.id Particles; /// List of particles' roles typedef std::vector ParticleRoles; /// Map between a particle's role and its associated Particle object typedef std::unordered_map > ParticlesMap; } #endif diff --git a/CepGen/Hadronisers/Pythia8Hadroniser.cpp b/CepGen/Hadronisers/Pythia8Hadroniser.cpp index 82043c7..ac17279 100644 --- a/CepGen/Hadronisers/Pythia8Hadroniser.cpp +++ b/CepGen/Hadronisers/Pythia8Hadroniser.cpp @@ -1,348 +1,359 @@ #include "CepGen/Hadronisers/GenericHadroniser.h" #include "CepGen/Hadronisers/HadronisersHandler.h" #include "CepGen/IO/PythiaEventInterface.h" #include "CepGen/Parameters.h" #include "CepGen/Physics/Kinematics.h" #include "CepGen/Physics/Constants.h" #include "CepGen/Physics/PDG.h" #include "CepGen/Core/ParametersList.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/utils.h" #include "CepGen/Event/Event.h" #include "CepGen/Event/Particle.h" #include #include #include #include namespace cepgen { namespace hadr { /** * Full interface to the Pythia8 hadronisation algorithm. It can be used in a single particle decay mode as well as a full event hadronisation using the string model, as in Jetset. * \brief Pythia8 hadronisation algorithm */ class Pythia8Hadroniser : public GenericHadroniser { public: explicit Pythia8Hadroniser( const ParametersList& ); ~Pythia8Hadroniser(); void setParameters( const Parameters& ) override; void readString( const char* param ) override; void init() override; bool run( Event& ev, double& weight, bool full ) override; void setCrossSection( double xsec, double xsec_err ) override; private: static constexpr unsigned short PYTHIA_STATUS_IN_BEAM = 12; static constexpr unsigned short PYTHIA_STATUS_IN_PARTON_KT = 61; std::vector min_ids_; std::unordered_map py_cg_corresp_; unsigned short findRole( const Event& ev, const Pythia8::Particle& p ) const; void updateEvent( Event& ev, double& weight ) const; Particle& addParticle( Event& ev, const Pythia8::Particle&, const Pythia8::Vec4& mom, unsigned short ) const; /// A Pythia8 core to be wrapped std::unique_ptr pythia_; std::unique_ptr cg_evt_; - bool correct_central_; + const bool correct_central_; + const bool debug_lhef_; bool res_decay_; bool enable_hadr_; unsigned short offset_; bool first_evt_; }; Pythia8Hadroniser::Pythia8Hadroniser( const ParametersList& plist ) : GenericHadroniser( plist, "pythia8" ), pythia_( new Pythia8::Pythia ), cg_evt_( new Pythia8::CepGenEvent ), correct_central_( plist.get( "correctCentralSystem", false ) ), + debug_lhef_( plist.get( "debugLHEF", false ) ), res_decay_( true ), enable_hadr_( false ), offset_( 0 ), first_evt_( true ) {} void Pythia8Hadroniser::setParameters( const Parameters& params ) { params_ = ¶ms; cg_evt_->initialise( params ); pythia_->setLHAupPtr( (Pythia8::LHAup*)cg_evt_.get() ); pythia_->settings.parm( "Beams:idA", (short)params_->kinematics.incoming_beams.first.pdg ); pythia_->settings.parm( "Beams:idB", (short)params_->kinematics.incoming_beams.second.pdg ); // specify we will be using a LHA input pythia_->settings.mode( "Beams:frameType", 5 ); pythia_->settings.parm( "Beams:eCM", params_->kinematics.sqrtS() ); min_ids_ = params_->kinematics.minimum_final_state; + if ( debug_lhef_ ) + cg_evt_->openLHEF( "debug.lhe" ); } Pythia8Hadroniser::~Pythia8Hadroniser() { pythia_->settings.writeFile( "last_pythia_config.cmd", false ); + if ( debug_lhef_ ) + cg_evt_->closeLHEF( true ); } void Pythia8Hadroniser::readString( const char* param ) { if ( !pythia_->readString( param ) ) throw CG_FATAL( "Pythia8Hadroniser" ) << "The Pythia8 core failed to parse the following setting:\n\t" << param; } void Pythia8Hadroniser::init() { pythia_->settings.flag( "ProcessLevel:resonanceDecays", res_decay_ ); if ( pythia_->settings.flag( "ProcessLevel:all" ) != enable_hadr_ ) pythia_->settings.flag( "ProcessLevel:all", enable_hadr_ ); if ( seed_ == -1ll ) pythia_->settings.flag( "Random:setSeed", false ); else { pythia_->settings.flag( "Random:setSeed", true ); pythia_->settings.mode( "Random:seed", seed_ ); } #if defined( PYTHIA_VERSION_INTEGER ) && PYTHIA_VERSION_INTEGER >= 8226 switch ( params_->kinematics.mode ) { case KinematicsMode::ElasticElastic: { pythia_->settings.mode( "BeamRemnants:unresolvedHadron", 3 ); pythia_->settings.flag( "PartonLevel:all", false ); } break; case KinematicsMode::InelasticElastic: { pythia_->settings.mode( "BeamRemnants:unresolvedHadron", 2 ); } break; case KinematicsMode::ElasticInelastic: { pythia_->settings.mode( "BeamRemnants:unresolvedHadron", 1 ); } break; case KinematicsMode::InelasticInelastic: default: { pythia_->settings.mode( "BeamRemnants:unresolvedHadron", 0 ); } break; } #else CG_WARNING( "Pythia8Hadroniser" ) << "Beam remnants framework for this version of Pythia " << "(" << Form( "%.3f", PYTHIA_VERSION ) << ")\n\t" << "does not support mixing of unresolved hadron states.\n\t" << "The proton remnants output might hence be wrong.\n\t" << "Please update the Pythia version or disable this part."; #endif if ( correct_central_ && res_decay_ ) CG_WARNING( "Pythia8Hadroniser" ) << "Central system's kinematics correction enabled while resonances are\n\t" << "expected to be decayed. Please check that this is fully intended."; if ( !pythia_->init() ) throw CG_FATAL( "Pythia8Hadroniser" ) << "Failed to initialise the Pythia8 core!\n\t" << "See the message above for more details."; + + if ( debug_lhef_ ) + cg_evt_->initLHEF(); } void Pythia8Hadroniser::setCrossSection( double xsec, double xsec_err ) { cg_evt_->setCrossSection( 0, xsec, xsec_err ); } bool Pythia8Hadroniser::run( Event& ev, double& weight, bool full ) { //--- initialise the event weight before running any decay algorithm weight = 1.; //--- only launch Pythia if: // 1) the full event kinematics (i.e. with remnants) is to be specified, // 2) the remnants are to be fragmented, or // 3) the resonances are to be decayed. if ( full && !remn_fragm_ && !res_decay_ ) return true; if ( !full && !res_decay_ ) return true; //--- switch full <-> partial event if ( full != enable_hadr_ ) { enable_hadr_ = full; init(); } //=========================================================================================== // convert our event into a custom LHA format //=========================================================================================== cg_evt_->feedEvent( ev, full ? Pythia8::CepGenEvent::Type::centralAndBeamRemnants : Pythia8::CepGenEvent::Type::centralAndPartons ); //if ( full ) cg_evt_->listEvent(); + if ( debug_lhef_ && full ) + cg_evt_->eventLHEF(); //=========================================================================================== // launch the hadronisation / resonances decays, and update the event accordingly //=========================================================================================== ev.num_hadronisation_trials = 0; while ( true ) { if ( ev.num_hadronisation_trials++ > max_trials_ ) return false; //--- run the hadronisation/fragmentation algorithm if ( pythia_->next() ) { //--- hadronisation successful if ( first_evt_ && full ) { offset_ = 0; for ( unsigned short i = 1; i < pythia_->event.size(); ++i ) if ( pythia_->event[i].status() == -PYTHIA_STATUS_IN_BEAM ) //--- no incoming particles in further stages offset_++; first_evt_ = false; } break; } } CG_DEBUG( "Pythia8Hadroniser" ) << "Pythia8 hadronisation performed successfully.\n\t" << "Number of trials: " << ev.num_hadronisation_trials << "/" << max_trials_ << ".\n\t" << "Particles multiplicity: " << ev.particles().size() << " → " << pythia_->event.size() << ".\n\t" << " indices offset: " << offset_ << "."; //=========================================================================================== // update the event content with Pythia's output //=========================================================================================== updateEvent( ev, weight ); return true; } Particle& Pythia8Hadroniser::addParticle( Event& ev, const Pythia8::Particle& py_part, const Pythia8::Vec4& mom, unsigned short role ) const { Particle& op = ev.addParticle( (Particle::Role)role ); ParticleProperties prop; const pdgid_t pdg_id = abs( py_part.id() ); try { prop = PDG::get()( pdg_id ); } catch ( const Exception& ) { prop = ParticleProperties{ pdg_id, py_part.name(), py_part.name(), (short)py_part.col(), // colour factor py_part.m0(), py_part.mWidth(), (short)py_part.charge(), // charge py_part.isLepton() }; PDG::get().define( prop ); } op.setPdgId( (short)py_part.id() ); op.setStatus( py_part.isFinal() ? Particle::Status::FinalState : (Particle::Role)role == Particle::Role::CentralSystem ? Particle::Status::Propagator : Particle::Status::Fragmented ); op.setMomentum( Particle::Momentum( mom.px(), mom.py(), mom.pz(), mom.e() ) ); op.setMass( mom.mCalc() ); cg_evt_->addCorresp( py_part.index()-offset_, op.id() ); return op; } void Pythia8Hadroniser::updateEvent( Event& ev, double& weight ) const { std::vector central_parts; for ( unsigned short i = 1+offset_; i < pythia_->event.size(); ++i ) { const Pythia8::Particle& p = pythia_->event[i]; const unsigned short cg_id = cg_evt_->cepgenId( i-offset_ ); if ( cg_id != Pythia8::CepGenEvent::INVALID_ID ) { //----- particle already in the event Particle& cg_part = ev[cg_id]; //--- fragmentation result if ( cg_part.role() == Particle::OutgoingBeam1 || cg_part.role() == Particle::OutgoingBeam2 ) { cg_part.setStatus( Particle::Status::Fragmented ); continue; } //--- resonance decayed; apply branching ratio for this decay if ( cg_part.role() == Particle::CentralSystem && p.status() < 0 ) { if ( res_decay_ ) weight *= p.particleDataEntry().pickChannel().bRatio(); cg_part.setStatus( Particle::Status::Resonance ); central_parts.emplace_back( i ); } //--- particle is not what we expect if ( p.idAbs() != abs( cg_part.integerPdgId() ) ) { CG_INFO( "Pythia8Hadroniser:update" ) << "LHAEVT event content:"; cg_evt_->listEvent(); CG_INFO( "Pythia8Hadroniser:update" ) << "Pythia event content:"; pythia_->event.list(); CG_INFO( "Pythia8Hadroniser:update" ) << "CepGen event content:"; ev.dump(); CG_INFO( "Pythia8Hadroniser:update" ) << "Correspondence:"; cg_evt_->dumpCorresp(); throw CG_FATAL( "Pythia8Hadroniser:update" ) << "Event list corruption detected for (Pythia/CepGen) particle " << i << "/" << cg_id << ":\n\t" << "should be " << abs( p.id() ) << ", " << "got " << cg_part.integerPdgId() << "!"; } } //--- check for messed up particles parentage and discard incoming beam particles /*else if ( p.mother1() > i || p.mother1() <= offset_ ) continue; else if ( p.mother2() > i || p.mother2() <= offset_ ) continue;*/ else { //----- new particle to be added const unsigned short role = findRole( ev, p ); switch ( (Particle::Role)role ) { case Particle::OutgoingBeam1: ev[Particle::OutgoingBeam1][0].setStatus( Particle::Status::Fragmented ); break; case Particle::OutgoingBeam2: ev[Particle::OutgoingBeam2][0].setStatus( Particle::Status::Fragmented ); break; default: break; } // found the role ; now we can add the particle Particle& cg_part = addParticle( ev, p, p.p(), role ); if ( correct_central_ && (Particle::Role)role == Particle::CentralSystem ) { const auto& ip = std::find( central_parts.begin(), central_parts.end(), p.mother1() ); if ( ip != central_parts.end() ) cg_part.setMomentum( ev[cg_evt_->cepgenId( *ip-offset_ )].momentum() ); } for ( const auto& moth_id : p.motherList() ) { if ( moth_id <= offset_ ) continue; const unsigned short moth_cg_id = cg_evt_->cepgenId( moth_id-offset_ ); if ( moth_cg_id != Pythia8::CepGenEvent::INVALID_ID ) cg_part.addMother( ev[moth_cg_id] ); else cg_part.addMother( addParticle( ev, pythia_->event[moth_id], p.p(), role ) ); if ( !p.isFinal() ) { if ( p.isResonance() || !p.daughterList().empty() ) cg_part.setStatus( Particle::Status::Resonance ); else cg_part.setStatus( Particle::Status::Undefined ); } } } } } unsigned short Pythia8Hadroniser::findRole( const Event& ev, const Pythia8::Particle& p ) const { for ( const auto& par_id : p.motherList() ) { if ( par_id == 1 && offset_ > 0 ) return (unsigned short)Particle::OutgoingBeam1; if ( par_id == 2 && offset_ > 0 ) return (unsigned short)Particle::OutgoingBeam2; const unsigned short par_cg_id = cg_evt_->cepgenId( par_id-offset_ ); if ( par_cg_id != Pythia8::CepGenEvent::INVALID_ID ) return (unsigned short)ev[par_cg_id].role(); return findRole( ev, pythia_->event[par_id] ); } return (unsigned short)Particle::UnknownRole; } } } // register hadroniser and define alias REGISTER_HADRONISER( pythia8, Pythia8Hadroniser ) diff --git a/CepGen/IO/DelphesHandler.cpp b/CepGen/IO/DelphesHandler.cpp new file mode 100644 index 0000000..2ea763c --- /dev/null +++ b/CepGen/IO/DelphesHandler.cpp @@ -0,0 +1,136 @@ +#include "CepGen/IO/ExportHandler.h" + +#include "CepGen/Parameters.h" +#include "CepGen/Core/ParametersList.h" +#include "CepGen/Core/Exception.h" + +#include "CepGen/Event/Event.h" + +#include "modules/Delphes.h" +#include "classes/DelphesFactory.h" +#include "classes/DelphesClasses.h" + +#include "ExRootAnalysis/ExRootTreeWriter.h" +#include "ExRootAnalysis/ExRootTreeBranch.h" + +#include "TFile.h" + +namespace cepgen +{ + namespace io + { + /** + * \brief Export handler for Delphes + * \author Laurent Forthomme + * \date Jul 2019 + */ + class DelphesHandler : public GenericExportHandler + { + public: + explicit DelphesHandler( const ParametersList& ); + ~DelphesHandler(); + + void initialise( const Parameters& ) override; + void setCrossSection( double xsec, double /*err_xsec*/ ) override { xsec_ = xsec; } + void operator<<( const Event& ) override; + + private: + std::unique_ptr output_; + const std::string input_card_; + std::unique_ptr delphes_; + std::unique_ptr conf_reader_; + std::unique_ptr tree_writer_; + //--- non-owning + DelphesFactory* factory_; + ExRootTreeBranch* evt_branch_; + TObjArray* out_all_parts_, *out_stab_parts_, *out_partons_; + double xsec_; + }; + + DelphesHandler::DelphesHandler( const ParametersList& params ) : + GenericExportHandler( "delphes" ), + output_( new TFile( params.get( "filename", "output.delphes.root" ).c_str(), "recreate" ) ), + input_card_( params.get( "inputCard", "input.tcl" ) ), + delphes_( new Delphes ), + conf_reader_( new ExRootConfReader ), + tree_writer_( new ExRootTreeWriter( output_.get(), "Delphes") ), + factory_( nullptr ), evt_branch_( nullptr ), + out_all_parts_( nullptr ), out_stab_parts_( nullptr ), out_partons_( nullptr ), + xsec_( -1. ) + { + try { + conf_reader_->ReadFile( input_card_.c_str() ); + } catch ( const std::runtime_error& err ) { + throw CG_FATAL( "DelphesHandler" ) + << "Failed to parse the Delphes configuration card!\n\t" + << err.what(); + } + delphes_->SetTreeWriter( tree_writer_.get() ); + delphes_->SetConfReader( conf_reader_.get() ); + } + + DelphesHandler::~DelphesHandler() + { + delphes_->FinishTask(); + tree_writer_->Write(); + } + + void + DelphesHandler::initialise( const Parameters& params ) + { + factory_ = delphes_->GetFactory(); + if ( !factory_ ) + throw CG_FATAL( "DelphesHandler" ) + << "Failed to retrieve factory object!"; + out_all_parts_ = delphes_->ExportArray( "allParticles" ); + out_stab_parts_ = delphes_->ExportArray( "stableParticles" ); + out_partons_ = delphes_->ExportArray( "partons" ); + evt_branch_ = tree_writer_->NewBranch( "Event", LHEFEvent::Class() ); + delphes_->InitTask(); + } + + void + DelphesHandler::operator<<( const Event& ev ) + { + delphes_->Clear(); + tree_writer_->Clear(); + //--- auxiliary event quantities + auto evt_aux = static_cast( evt_branch_->NewEntry() ); + evt_aux->Number = event_num_++; + evt_aux->ProcessID = 0; + evt_aux->Weight = 1.; // events are unweighted in CepGen + evt_aux->CrossSection = xsec_; + evt_aux->ScalePDF = 0.; // for the time being + evt_aux->AlphaQED = constants::ALPHA_EM; + evt_aux->AlphaQCD = constants::ALPHA_QCD; + evt_aux->ProcTime = ev.time_generation; + //--- particles content + for ( const auto& part : ev.particles() ) { + auto cand = factory_->NewCandidate(); + cand->PID = part.integerPdgId(); + cand->Status = (int)part.status(); + cand->Charge = part.charge(); + //--- kinematics part + cand->Mass = part.mass(); + const auto& mom = part.momentum(); + cand->Momentum.SetPxPyPzE( mom.px(), mom.py(), mom.pz(), mom.energy() ); + // no cand->Position specified (particles produced at origin) + //--- parentage part + cand->M1 = part.primary() ? 0 : *part.mothers().begin(); + cand->M2 = part.mothers().size() < 2 ? 0 : *part.mothers().rbegin(); + cand->D1 = part.daughters().empty() ? -1 : *part.daughters().begin(); + cand->D2 = part.daughters().size() < 2 ? -1 : *part.daughters().rbegin(); + //--- add to the proper collection(s) + out_all_parts_->Add( cand ); + if ( cand->Status == 1 ) + out_stab_parts_->Add( cand ); + else if ( cand->PID <= 5 || cand->PID == 21 || cand->PID == 15 ) + out_partons_->Add( cand ); + } + delphes_->ProcessTask(); + tree_writer_->Fill(); + } + } +} + +REGISTER_IO_MODULE( delphes, DelphesHandler ) diff --git a/CepGen/IO/ExportHandler.h b/CepGen/IO/ExportHandler.h index ed08dcc..6d4352c 100644 --- a/CepGen/IO/ExportHandler.h +++ b/CepGen/IO/ExportHandler.h @@ -1,24 +1,24 @@ #ifndef CepGen_IO_ExportHandler_h #define CepGen_IO_ExportHandler_h #include "CepGen/Core/ModuleFactory.h" #include "CepGen/IO/GenericExportHandler.h" #define REGISTER_IO_MODULE( name, obj ) \ - namespace cepgen { namespace output { \ + namespace cepgen { namespace io { \ struct BUILDERNM( name ) { \ BUILDERNM( name )() { ExportHandler::get().registerModule( STRINGIFY( name ) ); } }; \ static BUILDERNM( name ) g ## name; \ } } namespace cepgen { - namespace output + namespace io { /// A hadroniser modules factory typedef ModuleFactory ExportHandler; } } #endif diff --git a/CepGen/IO/GenericExportHandler.cpp b/CepGen/IO/GenericExportHandler.cpp index 6a16659..5c5cd6b 100644 --- a/CepGen/IO/GenericExportHandler.cpp +++ b/CepGen/IO/GenericExportHandler.cpp @@ -1,69 +1,72 @@ #include "CepGen/IO/GenericExportHandler.h" #include "CepGen/StructureFunctions/StructureFunctions.h" #include "CepGen/Physics/Constants.h" #include "CepGen/Parameters.h" #include "CepGen/Version.h" #include namespace cepgen { - namespace output + namespace io { - GenericExportHandler::GenericExportHandler() : - event_num_( 0. ) + GenericExportHandler::GenericExportHandler( const std::string& name ) : + name_( name ), event_num_( 0. ) + {} + + GenericExportHandler::~GenericExportHandler() {} std::string - GenericExportHandler::banner( const Parameters& params ) + GenericExportHandler::banner( const Parameters& params, const std::string& prep ) { std::ostringstream os; os - << " ***** Sample generated with CepGen v" << version() << " *****\n" - << " * process: " << params.processName() << " (" << params.kinematics.mode << ")\n"; + << prep << " ***** Sample generated with CepGen v" << version() << " *****\n" + << prep << " * process: " << params.processName() << " (" << params.kinematics.mode << ")\n"; if ( params.kinematics.mode != KinematicsMode::ElasticElastic ) { - os << " * structure functions: " << params.kinematics.structure_functions->description() << "\n"; + os << prep << " * structure functions: " << params.kinematics.structure_functions->description() << "\n"; if ( !params.hadroniserName().empty() ) - os << " * hadroniser: " << params.hadroniserName() << "\n"; + os << prep << " * hadroniser: " << params.hadroniserName() << "\n"; } os - << " *--- incoming state\n"; + << prep << " *--- incoming state\n"; if ( params.kinematics.cuts.initial.q2.valid() ) os - << " * Q2 range (GeV2): " + << prep << " * Q2 range (GeV2): " << params.kinematics.cuts.initial.q2 << "\n"; if ( params.kinematics.mode != KinematicsMode::ElasticElastic && params.kinematics.cuts.remnants.mass_single.valid() ) os - << " * remnants mass range (GeV/c2): " + << prep << " * remnants mass range (GeV/c2): " << params.kinematics.cuts.remnants.mass_single << "\n"; - os << " *--- central system\n"; + os << prep << " *--- central system\n"; if ( params.kinematics.cuts.central.pt_single.valid() ) os - << " * single particle pt (GeV/c): " + << prep << " * single particle pt (GeV/c): " << params.kinematics.cuts.central.pt_single << "\n"; if ( params.kinematics.cuts.central.energy_single.valid() ) os - << " * single particle energy (GeV): " + << prep << " * single particle energy (GeV): " << params.kinematics.cuts.central.energy_single << "\n"; if ( params.kinematics.cuts.central.eta_single.valid() ) os - << " * single particle eta: " + << prep << " * single particle eta: " << params.kinematics.cuts.central.eta_single << "\n"; if ( params.kinematics.cuts.central.pt_sum.valid() ) os - << " * total pt (GeV/c): " + << prep << " * total pt (GeV/c): " << params.kinematics.cuts.central.mass_sum << "\n"; if ( params.kinematics.cuts.central.mass_sum.valid() ) os - << " * total invariant mass (GeV/c2): " + << prep << " * total invariant mass (GeV/c2): " << params.kinematics.cuts.central.mass_sum << "\n"; os - << " **************************************************"; + << prep << " **************************************************"; return os.str(); } } } diff --git a/CepGen/IO/GenericExportHandler.h b/CepGen/IO/GenericExportHandler.h index 72261cc..775ce03 100644 --- a/CepGen/IO/GenericExportHandler.h +++ b/CepGen/IO/GenericExportHandler.h @@ -1,41 +1,46 @@ #ifndef CepGen_IO_GenericExportHandler_h #define CepGen_IO_GenericExportHandler_h #include +#include namespace cepgen { class Event; class Parameters; /// Location for all output generators - namespace output + namespace io { /** * \brief Output format handler for events export * \author Laurent Forthomme * \date Sep 2016 */ class GenericExportHandler { public: /// Class constructor - explicit GenericExportHandler(); + explicit GenericExportHandler( const std::string& ); + virtual ~GenericExportHandler(); + const std::string& name() const { return name_; } /// Initialise the handler and its inner parameterisation virtual void initialise( const Parameters& ) = 0; /// Set the process cross section and its associated error virtual void setCrossSection( double /*xsec*/, double /*err_xsec*/ ) {} /// Set the event number void setEventNumber( const unsigned int& ev_id ) { event_num_ = ev_id; } /// Writer operator virtual void operator<<( const Event& ) = 0; protected: - static std::string banner( const Parameters& ); + static std::string banner( const Parameters&, const std::string& prep = "" ); + /// Module unique name + const std::string name_; /// Event index unsigned int event_num_; }; } } #endif diff --git a/CepGen/IO/HepMCEventInterface.cpp b/CepGen/IO/HepMCEventInterface.cpp new file mode 100644 index 0000000..3c21899 --- /dev/null +++ b/CepGen/IO/HepMCEventInterface.cpp @@ -0,0 +1,130 @@ +#include "CepGen/IO/HepMCEventInterface.h" + +#include "CepGen/Parameters.h" +#include "CepGen/Core/Exception.h" + +#include "CepGen/Event/Event.h" +#include "CepGen/Physics/Constants.h" +#include "CepGen/Physics/PDG.h" + +#ifdef HEPMC3 +# include "HepMC3/Version.h" +# include "HepMC3/FourVector.h" +# include "HepMC3/GenEvent.h" +# include "HepMC3/GenRunInfo.h" +# include "HepMC3/GenVertex.h" +# include "HepMC3/GenParticle.h" +#else +# include "HepMC/Version.h" +# if !defined( HEPMC_VERSION_CODE ) // HepMC v2 +# include "HepMC/SimpleVector.h" +# include "HepMC/GenEvent.h" +# include "HepMC/GenVertex.h" +# include "HepMC/GenParticle.h" +# else +# include "HepMC/FourVector.h" +# include "HepMC/GenEvent.h" +# include "HepMC/GenRunInfo.h" +# include "HepMC/GenVertex.h" +# include "HepMC/GenParticle.h" +# define HEPMC3 +# endif +#endif + +namespace HepMC +{ + CepGenEvent::CepGenEvent() : + GenEvent( Units::GEV, Units::MM ) + { +#ifdef HEPMC3 + GenEvent::add_attribute( "AlphaQCD", make_shared( cepgen::constants::ALPHA_QCD ) ); + GenEvent::add_attribute( "AlphaEM", make_shared( cepgen::constants::ALPHA_EM ) ); +#else + GenEvent::set_alphaQCD( cepgen::constants::ALPHA_QCD ); + GenEvent::set_alphaQED( cepgen::constants::ALPHA_EM ); +#endif + } + + void + CepGenEvent::feedEvent( const cepgen::Event& evt ) + { + GenEvent::clear(); + GenEvent::weights().push_back( 1. ); // unweighted events + + // filling the particles content + const FourVector origin( 0., 0., 0., 0. ); + auto& part_vec = evt.particles(); + + int cm_id = 0, idx = 1; + +#ifdef HEPMC3 + GenVertexPtr v1 = make_shared( origin ), v2 = make_shared( origin ); + GenVertexPtr vcm = make_shared( origin ); +#else + GenVertex* v1 = new GenVertex( origin ), *v2 = new GenVertex( origin ); + GenVertex* vcm = new GenVertex( origin ); +#endif + for ( unsigned int i = 0; i < part_vec.size(); ++i ) { + const auto& part_orig = part_vec.at( i ); + const auto& mom_orig = part_orig.momentum(); + FourVector pmom( mom_orig.px(), mom_orig.py(), mom_orig.pz(), part_orig.energy() ); +#ifdef HEPMC3 + auto part = make_shared( pmom, part_orig.integerPdgId(), (int)part_orig.status() ); +#else + auto part = new GenParticle( pmom, part_orig.integerPdgId(), (int)part_orig.status() ); + part->suggest_barcode( idx++ ); +#endif + part->set_generated_mass( cepgen::PDG::get().mass( part_orig.pdgId() ) ); + + switch ( part_orig.role() ) { + case cepgen::Particle::IncomingBeam1: + v1->add_particle_in( part ); + break; + case cepgen::Particle::IncomingBeam2: + v2->add_particle_in( part ); + break; + case cepgen::Particle::OutgoingBeam1: + v1->add_particle_out( part ); + break; + case cepgen::Particle::OutgoingBeam2: + v2->add_particle_out( part ); + break; + case cepgen::Particle::Parton1: + v1->add_particle_out( part ); + vcm->add_particle_in( part ); + break; + case cepgen::Particle::Parton2: + v2->add_particle_out( part ); + vcm->add_particle_in( part ); + break; + case cepgen::Particle::Intermediate: + // skip the two-parton system and propagate the parentage + cm_id = i; + continue; + case cepgen::Particle::CentralSystem: default: { + const auto& moth = part_orig.mothers(); + if ( moth.empty() ) + // skip disconnected lines + continue; + const short m1 = *moth.begin(), m2 = moth.size() > 1 ? *moth.rbegin() : -1; + if ( cm_id == m1 || ( m2 >= 0 && ( m1 < cm_id && cm_id <= m2 ) ) ) // also supports range + // if connected to the central system + vcm->add_particle_out( part ); + else + throw CG_FATAL( "HepMCHandler:fillEvent" ) + << "Other particle requested! Not yet implemented!"; + } break; + } + idx++; + } + GenEvent::add_vertex( v1 ); + GenEvent::add_vertex( v2 ); + GenEvent::add_vertex( vcm ); + +#ifndef HEPMC3 + GenEvent::set_beam_particles( *v1->particles_in_const_begin(), *v2->particles_in_const_begin() ); + GenEvent::set_signal_process_vertex( vcm ); +#endif + } +} + diff --git a/CepGen/IO/HepMCEventInterface.h b/CepGen/IO/HepMCEventInterface.h new file mode 100644 index 0000000..47845f7 --- /dev/null +++ b/CepGen/IO/HepMCEventInterface.h @@ -0,0 +1,29 @@ +#ifndef CepGen_IO_HepMCEventInterface_h +#define CepGen_IO_HepMCEventInterface_h + +#ifdef HEPMC3 +# include "HepMC3/GenEvent.h" +# define HepMC HepMC3 +#else +# include "HepMC/GenEvent.h" +#endif + + +namespace cepgen { class Event; } + +namespace HepMC +{ + /// Interfacing between CepGen and HepMC event definitions + /// \author Laurent Forthomme + /// \date Jul 2019 + class CepGenEvent : public GenEvent + { + public: + explicit CepGenEvent(); + /// Feed a new CepGen event to this conversion object + /// \param[in] ev CepGen event to be fed + void feedEvent( const cepgen::Event& ev ); + }; +} +#endif + diff --git a/CepGen/IO/HepMCHandler.cpp b/CepGen/IO/HepMCHandler.cpp index 2f5dd1c..f6e89a1 100644 --- a/CepGen/IO/HepMCHandler.cpp +++ b/CepGen/IO/HepMCHandler.cpp @@ -1,233 +1,141 @@ #include "CepGen/IO/ExportHandler.h" +#include "CepGen/IO/HepMCEventInterface.h" #include "CepGen/Parameters.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/ParametersList.h" #include "CepGen/Event/Event.h" #include "CepGen/Physics/Constants.h" #ifdef HEPMC3 +# include "HepMC3/Version.h" +# define HEPMC_VERSION HEPMC3_VERSION +# define HEPMC_VERSION_CODE HEPMC3_VERSION_CODE # include "HepMC3/WriterAscii.h" # include "HepMC3/WriterAsciiHepMC2.h" # include "HepMC3/WriterHEPEVT.h" # ifdef HEPMC3_ROOTIO # include "HepMC3/WriterRoot.h" # include "HepMC3/WriterRootTree.h" # endif -# include "HepMC3/FourVector.h" -# include "HepMC3/GenEvent.h" -# include "HepMC3/GenRunInfo.h" -# include "HepMC3/GenVertex.h" -# include "HepMC3/GenParticle.h" -using namespace HepMC3; #else # include "HepMC/Version.h" # if !defined( HEPMC_VERSION_CODE ) // HepMC v2 # include "HepMC/IO_GenEvent.h" -# include "HepMC/SimpleVector.h" -# include "HepMC/GenEvent.h" -# include "HepMC/GenVertex.h" -# include "HepMC/GenParticle.h" +# include "HepMC/IO_AsciiParticles.h" # else # include "HepMC/WriterAscii.h" # include "HepMC/WriterHEPEVT.h" -# include "HepMC/FourVector.h" -# include "HepMC/GenEvent.h" -# include "HepMC/GenRunInfo.h" -# include "HepMC/GenVertex.h" -# include "HepMC/GenParticle.h" # define HEPMC3 # endif -using namespace HepMC; #endif #include +using namespace HepMC; + namespace cepgen { - namespace output + namespace io { /// Handler for the HepMC file output /// \tparam T HepMC writer handler (format-dependent) /// \author Laurent Forthomme /// \date Sep 2016 template class HepMCHandler : public GenericExportHandler { public: /// Class constructor explicit HepMCHandler( const ParametersList& ); - void initialise( const Parameters& /*params*/ ) override; + ~HepMCHandler(); + + void initialise( const Parameters& /*params*/ ) override {} /// Writer operator void operator<<( const Event& ) override; void setCrossSection( double, double ) override; private: - /// Clear the associated HepMC event content - void clearEvent(); - /// Populate the associated HepMC event with a Event object - void fillEvent( const Event& ); - /// Writer object std::unique_ptr output_; /// Generator cross section and error std::shared_ptr xs_; #ifdef HEPMC3 /// Auxiliary information on run std::shared_ptr runinfo_; #endif /// Associated HepMC event - std::shared_ptr event_; + std::shared_ptr event_; }; template HepMCHandler::HepMCHandler( const ParametersList& params ) : - output_( new T( params.get( "filename", "output.hepmc" ) ) ), + GenericExportHandler( "hepmc" ), + output_( new T( params.get( "filename", "output.hepmc" ).c_str() ) ), xs_( new GenCrossSection ), #ifdef HEPMC3 runinfo_( new GenRunInfo ), #endif - event_( new GenEvent( Units::GEV, Units::MM ) ) + event_( new CepGenEvent ) { #ifdef HEPMC3 output_->set_run_info( runinfo_ ); runinfo_->set_weight_names( { "Default" } ); #endif + CG_INFO( "HepMC" ) + << "Interfacing module initialised " + << "for HepMC version " << HEPMC_VERSION << "."; } - template void - HepMCHandler::operator<<( const Event& evt ) - { - fillEvent( evt ); -#ifdef HEPMC3 - output_->write_event( *event_ ); -#else - output_->write_event( event_.get() ); -#endif - } - - template void - HepMCHandler::initialise( const Parameters& ) + template + HepMCHandler::~HepMCHandler() { #ifdef HEPMC3 - event_->add_attribute( "AlphaQCD", make_shared( constants::ALPHA_QCD ) ); - event_->add_attribute( "AlphaEM", make_shared( constants::ALPHA_EM ) ); -#else - event_->set_alphaQCD( constants::ALPHA_QCD ); - event_->set_alphaQED( constants::ALPHA_EM ); + output_->close(); #endif } template void - HepMCHandler::setCrossSection( double xsect, double xsect_err ) - { - xs_->set_cross_section( xsect, xsect_err ); - } - - template void - HepMCHandler::fillEvent( const Event& evt ) + HepMCHandler::operator<<( const Event& evt ) { - event_->clear(); - + event_->feedEvent( evt ); // general information #ifdef HEPMC3 event_->set_cross_section( xs_ ); event_->set_run_info( runinfo_ ); #else event_->set_cross_section( *xs_ ); #endif - - event_->set_event_number( event_num_ ); - event_->weights().push_back( 1. ); // unweighted events - - // filling the particles content - const FourVector origin( 0., 0., 0., 0. ); - Particles part_vec = evt.particles(); - - int cm_id = 0, idx = 1; - + event_->set_event_number( event_num_++ ); #ifdef HEPMC3 - GenVertexPtr v1 = make_shared( origin ), - v2 = make_shared( origin ), - vcm = make_shared( origin ); -#else - GenVertex* v1 = new GenVertex( origin ), - *v2 = new GenVertex( origin ), - *vcm = new GenVertex( origin ); -#endif - for ( unsigned int i = 0; i < part_vec.size(); ++i ) { - const Particle part_orig = part_vec.at( i ); - FourVector pmom( part_orig.momentum().px(), - part_orig.momentum().py(), - part_orig.momentum().pz(), - part_orig.energy() ); -#ifdef HEPMC3 - GenParticlePtr part = make_shared( pmom, part_orig.integerPdgId(), (int)part_orig.status() ); + output_->write_event( *event_ ); #else - GenParticle* part = new GenParticle( pmom, part_orig.integerPdgId(), (int)part_orig.status() ); - part->suggest_barcode( idx++ ); + output_->write_event( event_.get() ); #endif - const ParticlesIds moth = part_orig.mothers(); - - switch ( part_orig.role() ) { - case Particle::IncomingBeam1: { v1->add_particle_in( part ); } break; - case Particle::IncomingBeam2: { v2->add_particle_in( part ); } break; - case Particle::OutgoingBeam1: { v1->add_particle_out( part ); } break; - case Particle::OutgoingBeam2: { v2->add_particle_out( part ); } break; - case Particle::Parton1: { v1->add_particle_out( part ); vcm->add_particle_in( part ); } break; - case Particle::Parton2: { v2->add_particle_out( part ); vcm->add_particle_in( part ); } break; - case Particle::Intermediate: { cm_id = i; continue; } break; - case Particle::CentralSystem: - default: { - if ( moth.empty() ) continue; - if ( *moth.begin() == cm_id ) vcm->add_particle_out( part ); - else - throw CG_FATAL( "HepMCHandler:fillEvent" ) - << "Other particle requested! Not yet implemented!"; - } break; - } - idx++; - } - event_->add_vertex( v1 ); - event_->add_vertex( v2 ); - event_->add_vertex( vcm ); + } -#ifndef HEPMC3 - event_->set_beam_particles( *v1->particles_in_const_begin(), *v2->particles_in_const_begin() ); - event_->set_signal_process_vertex( *v1->vertices_begin() ); - event_->set_beam_particles( *v1->particles_in_const_begin(), *v2->particles_in_const_end() ); -#endif - event_num_++; + template void + HepMCHandler::setCrossSection( double xsect, double xsect_err ) + { + xs_->set_cross_section( xsect, xsect_err ); } -#ifdef HEPMC3 -# if HEPMC_VERSION_CODE >= 3001000 - typedef HepMCHandler HepMC2Handler; -# endif - typedef HepMCHandler HepMC3Handler; - typedef HepMCHandler HEPEVTHandler; -# ifdef HEPMC3_ROOTIO - typedef HepMCHandler RootHandler; - typedef HepMCHandler RootTreeHandler; -# endif -#else - typedef HepMCHandler HepMC2Handler; -#endif } } #ifdef HEPMC3 -REGISTER_IO_MODULE( hepmc3, HepMC3Handler ) -REGISTER_IO_MODULE( hepmc, HepMC3Handler ) -REGISTER_IO_MODULE( hepevt, HEPEVTHandler ) +REGISTER_IO_MODULE( hepmc, HepMCHandler ) +REGISTER_IO_MODULE( hepmc3, HepMCHandler ) +REGISTER_IO_MODULE( hepevt, HepMCHandler ) +# if HEPMC_VERSION_CODE >= 3001000 +REGISTER_IO_MODULE( hepmc2, HepMCHandler ) +# endif # ifdef HEPMC3_ROOTIO -REGISTER_IO_MODULE( hepmc_root, RootHandler ) -REGISTER_IO_MODULE( hepmc_root_tree, RootTreeHandler ) +REGISTER_IO_MODULE( hepmc_root, HepMCHandler ) +REGISTER_IO_MODULE( hepmc_root_tree, HepMCHandler ) # endif -#else -REGISTER_IO_MODULE( hepmc, HepMC2Handler ) -#endif -#if defined( HEPMC3 ) && HEPMC_VERSION_CODE >= 3001000 -REGISTER_IO_MODULE( hepmc2, HepMC2Handler ) +#else // HepMC version 2 and below +REGISTER_IO_MODULE( hepmc, HepMCHandler ) +REGISTER_IO_MODULE( hepmc_ascii, HepMCHandler ) #endif diff --git a/CepGen/IO/LHEFHepMCHandler.cpp b/CepGen/IO/LHEFHepMCHandler.cpp index e18b0a8..4402e48 100644 --- a/CepGen/IO/LHEFHepMCHandler.cpp +++ b/CepGen/IO/LHEFHepMCHandler.cpp @@ -1,105 +1,106 @@ #include "CepGen/IO/ExportHandler.h" #include "CepGen/Parameters.h" #include "CepGen/Core/ParametersList.h" #include "CepGen/Event/Event.h" #include #ifdef HEPMC3 using namespace std; // account for improper scoping in following includes # include "HepMC3/LHEF.h" #else # include "HepMC/Version.h" # ifdef HEPMC_VERSION_CODE // HepMC v3+ # include "HepMC/LHEF.h" # else # define NO_LHEF # endif #endif #ifndef NO_LHEF namespace cepgen { - namespace output + namespace io { /** * \brief Handler for the LHE file output * \author Laurent Forthomme * \date Sep 2016 */ class LHEFHepMCHandler : public GenericExportHandler { public: /// Class constructor explicit LHEFHepMCHandler( const ParametersList& ); void initialise( const Parameters& ) override; /// Writer operator void operator<<( const Event& ) override; void setCrossSection( double, double ) override {} private: /// Writer object (from HepMC) std::unique_ptr lhe_output_; LHEF::HEPRUP run_; }; LHEFHepMCHandler::LHEFHepMCHandler( const ParametersList& params ) : + GenericExportHandler( "lhef" ), lhe_output_( new LHEF::Writer( params.get( "filename", "output.lhe" ) ) ) {} void LHEFHepMCHandler::initialise( const Parameters& params ) { lhe_output_->headerBlock() << ""; //--- first specify information about the run LHEF::HEPRUP run = lhe_output_->heprup; run.IDBMUP = { (int)params.kinematics.incoming_beams.first.pdg, (int)params.kinematics.incoming_beams.second.pdg }; run.EBMUP = { (double)params.kinematics.incoming_beams.first.pz, (double)params.kinematics.incoming_beams.second.pz }; run.NPRUP = 1; run.resize(); run.XSECUP[0] = params.integration().result; run.XERRUP[0] = params.integration().err_result; run.XMAXUP[0] = 1.; run.LPRUP[0] = 1; lhe_output_->heprup = run; //--- ensure everything is properly parsed lhe_output_->init(); } void LHEFHepMCHandler::operator<<( const Event& ev ) { LHEF::HEPEUP out; out.heprup = &lhe_output_->heprup; out.XWGTUP = 1.; out.XPDWUP = std::pair( 0., 0. ); out.SCALUP = 0.; out.AQEDUP = constants::ALPHA_EM; out.AQCDUP = constants::ALPHA_QCD; - out.NUP = ev.numParticles(); + out.NUP = ev.size(); out.resize(); - for ( unsigned short ip = 0; ip < ev.numParticles(); ++ip ) { + for ( unsigned short ip = 0; ip < ev.size(); ++ip ) { const Particle part = ev[ip]; out.IDUP[ip] = part.integerPdgId(); // PDG id out.ISTUP[ip] = (short)part.status(); // status code out.PUP[ip] = part.momentum().pVector(); // momentum out.MOTHUP[ip] = { // mothers part.mothers().size() > 0 ? *part.mothers(). begin()+1 : 0, part.mothers().size() > 1 ? *part.mothers().rbegin()+1 : 0 }; out.ICOLUP[ip] = { 0, 0 }; out.VTIMUP[ip] = 0.; // invariant lifetime out.SPINUP[ip] = 0.; } //lhe_output_->eventComments() << "haha"; lhe_output_->hepeup = out; lhe_output_->writeEvent(); } } } REGISTER_IO_MODULE( lhef, LHEFHepMCHandler ) #endif diff --git a/CepGen/IO/LHEFPythiaHandler.cpp b/CepGen/IO/LHEFPythiaHandler.cpp index cc0c2bb..8e6eadd 100644 --- a/CepGen/IO/LHEFPythiaHandler.cpp +++ b/CepGen/IO/LHEFPythiaHandler.cpp @@ -1,86 +1,87 @@ #include "CepGen/IO/ExportHandler.h" #include "CepGen/IO/PythiaEventInterface.h" #include "CepGen/Core/ParametersList.h" #include "CepGen/Event/Event.h" #include "Pythia8/Pythia.h" #include namespace cepgen { - namespace output + namespace io { /** * \brief Handler for the LHE file output * \author Laurent Forthomme * \date Sep 2016 */ class LHEFPythiaHandler : public GenericExportHandler { public: /// Class constructor explicit LHEFPythiaHandler( const ParametersList& ); ~LHEFPythiaHandler(); void initialise( const Parameters& ) override; /// Writer operator void operator<<( const Event& ) override; void setCrossSection( double, double ) override; private: std::unique_ptr pythia_; std::unique_ptr lhaevt_; }; LHEFPythiaHandler::LHEFPythiaHandler( const ParametersList& params ) : + GenericExportHandler( "lhef" ), pythia_( new Pythia8::Pythia ), lhaevt_( new Pythia8::CepGenEvent ) { lhaevt_->openLHEF( params.get( "filename", "output.lhe" ) ); } LHEFPythiaHandler::~LHEFPythiaHandler() { if ( lhaevt_ ) lhaevt_->closeLHEF( false ); // we do not want to rewrite the init block } void LHEFPythiaHandler::initialise( const Parameters& params ) { std::ostringstream oss_init; oss_init << ""; oss_init << std::endl; // LHEF is usually not as beautifully parsed as a standard XML... // we're physicists, what do you expect? lhaevt_->addComments( oss_init.str() ); lhaevt_->initialise( params ); pythia_->setLHAupPtr( lhaevt_.get() ); pythia_->settings.flag( "ProcessLevel:all", false ); // we do not want Pythia to interfere... pythia_->settings.flag( "PartonLevel:all", false ); // we do not want Pythia to interfere... pythia_->settings.flag( "HadronLevel:all", false ); // we do not want Pythia to interfere... pythia_->settings.mode( "Beams:frameType", 5 ); // LHEF event readout pythia_->settings.mode( "Next:numberCount", 0 ); // remove some of the Pythia output //pythia_->settings.flag( "Check:event", false ); pythia_->init(); lhaevt_->initLHEF(); } void LHEFPythiaHandler::operator<<( const Event& ev ) { lhaevt_->feedEvent( ev, Pythia8::CepGenEvent::Type::centralAndFullBeamRemnants ); pythia_->next(); lhaevt_->eventLHEF(); } void LHEFPythiaHandler::setCrossSection( double xsect, double xsect_err ) { lhaevt_->setCrossSection( 0, xsect, xsect_err ); } } } REGISTER_IO_MODULE( lhef, LHEFPythiaHandler ) diff --git a/CepGen/IO/ROOTHandler.cpp b/CepGen/IO/ROOTHandler.cpp new file mode 100644 index 0000000..fba85b5 --- /dev/null +++ b/CepGen/IO/ROOTHandler.cpp @@ -0,0 +1,82 @@ +#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/Event/EventBrowser.h" +#include "CepGen/Parameters.h" + +#include "CepGen/Version.h" + +#include "TFile.h" +#include "TH1.h" + +namespace cepgen +{ + namespace io + { + /** + * \brief Handler for the generic ROOT file output + * \author Laurent Forthomme + * \date Jul 2019 + */ + class ROOTHandler : public GenericExportHandler + { + public: + explicit ROOTHandler( const ParametersList& ); + ~ROOTHandler(); + + void initialise( const Parameters& ) override {} + void setCrossSection( double xsec, double ) override { xsec_ = xsec; } + void operator<<( const Event& ) override; + + private: + std::unique_ptr file_; + std::vector > hists_; + const ParametersList variables_; + + double xsec_; + const utils::EventBrowser browser_; + }; + + ROOTHandler::ROOTHandler( const ParametersList& params ) : + GenericExportHandler( "text" ), + file_ ( TFile::Open( params.get( "filename", "output.root" ).c_str(), "recreate" ) ), + variables_( params.get( "variables" ) ), + xsec_( 1. ) + { + //--- extract list of variables to be plotted in histograms + for ( const auto& var : variables_.keys() ) { + const auto& hvar = variables_.get( var ); + const int nbins = hvar.get( "nbins", 10 ); + const double min = hvar.get( "low", 0. ), max = hvar.get( "high", 1. ); + const auto title = Form( "%s;%s;d#sigma/d(%s) (pb/bin)", var.c_str(), var.c_str(), var.c_str() ); + hists_.emplace_back( std::make_pair( var, new TH1D( var.c_str(), title.c_str(), nbins, min, max ) ) ); + CG_INFO( "ROOTHandler" ) + << "Booking a histogram with " << nbins << " bin" << utils::s( nbins ) + << " between " << min << " and " << max << " for \"" << var << "\"."; + } + } + + ROOTHandler::~ROOTHandler() + { + //--- finalisation of the output file + for ( const auto& hist : hists_ ) + hist.second->Write( hist.first.c_str() ); + // ROOT and its sumptuous memory management disallows the "delete" here + file_->Close(); + } + + void + ROOTHandler::operator<<( const Event& ev ) + { + //--- increment the corresponding histograms + for ( const auto& h_var : hists_ ) + h_var.second->Fill( browser_.get( ev, h_var.first ), xsec_ ); + } + } +} + +REGISTER_IO_MODULE( root, ROOTHandler ) diff --git a/CepGen/IO/ROOTTreeHandler.cpp b/CepGen/IO/ROOTTreeHandler.cpp new file mode 100644 index 0000000..2b4d51f --- /dev/null +++ b/CepGen/IO/ROOTTreeHandler.cpp @@ -0,0 +1,103 @@ +#include "CepGen/IO/ExportHandler.h" +#include "CepGen/IO/ROOTTreeInfo.h" + +#include "CepGen/Core/Exception.h" +#include "CepGen/Core/ParametersList.h" + +#include "CepGen/Event/Event.h" +#include "CepGen/Parameters.h" + +// ROOT includes +#include "TFile.h" + +#include + +namespace cepgen +{ + namespace io + { + /** + * Handler for the storage of events in a ROOT format + * \author Laurent Forthomme + * \date 27 Jan 2014 + */ + class ROOTTreeHandler : public GenericExportHandler + { + public: + /// Class constructor + explicit ROOTTreeHandler( const ParametersList& ); + ~ROOTTreeHandler(); + + void initialise( const Parameters& ) override; + /// Writer operator + void operator<<( const Event& ) override; + void setCrossSection( double, double ) override; + + private: + std::unique_ptr file_; + std::unique_ptr run_tree_; + std::unique_ptr evt_tree_; + }; + + ROOTTreeHandler::ROOTTreeHandler( const ParametersList& params ) : + GenericExportHandler( "root" ), + file_( TFile::Open( params.get( "filename", "output.root" ).c_str(), "recreate" ) ), + run_tree_( new ROOT::CepGenRun ), evt_tree_( new ROOT::CepGenEvent ) + { + if ( !file_->IsOpen() ) + throw CG_FATAL( "ROOTTreeHandler" ) << "Failed to create the output file!"; + run_tree_->create(); + evt_tree_->create(); + } + + ROOTTreeHandler::~ROOTTreeHandler() + { + run_tree_->fill(); + file_->Write(); + } + + void + ROOTTreeHandler::initialise( const Parameters& params ) + { + run_tree_->litigious_events = 0; + run_tree_->sqrt_s = params.kinematics.sqrtS(); + } + + void + ROOTTreeHandler::operator<<( const Event& ev ) + { + evt_tree_->gen_time = ev.time_generation; + evt_tree_->tot_time = ev.time_total; + evt_tree_->np = 0; + for ( const auto& p : ev.particles() ) { + const cepgen::Particle::Momentum m = p.momentum(); + evt_tree_->rapidity[evt_tree_->np] = m.rapidity(); + evt_tree_->pt[evt_tree_->np] = m.pt(); + evt_tree_->eta[evt_tree_->np] = m.eta(); + evt_tree_->phi[evt_tree_->np] = m.phi(); + evt_tree_->E[evt_tree_->np] = p.energy(); + evt_tree_->m[evt_tree_->np] = p.mass(); + evt_tree_->pdg_id[evt_tree_->np] = p.integerPdgId(); + evt_tree_->parent1[evt_tree_->np] = ( p.mothers().size() > 0 ) ? *p.mothers().begin() : -1; + evt_tree_->parent2[evt_tree_->np] = ( p.mothers().size() > 1 ) ? *p.mothers().rbegin() : -1; + evt_tree_->status[evt_tree_->np] = (int)p.status(); + evt_tree_->stable[evt_tree_->np] = ( (short)p.status() > 0 ); + evt_tree_->charge[evt_tree_->np] = p.charge(); + evt_tree_->role[evt_tree_->np] = p.role(); + + evt_tree_->np++; + } + run_tree_->num_events += 1; + evt_tree_->fill(); + } + + void + ROOTTreeHandler::setCrossSection( double xsect, double xsect_err ) + { + run_tree_->xsect = xsect; + run_tree_->errxsect = xsect_err; + } + } +} + +REGISTER_IO_MODULE( root, ROOTTreeHandler ) diff --git a/test/TreeInfo.h b/CepGen/IO/ROOTTreeInfo.h similarity index 98% rename from test/TreeInfo.h rename to CepGen/IO/ROOTTreeInfo.h index 7d95027..b7082c7 100644 --- a/test/TreeInfo.h +++ b/CepGen/IO/ROOTTreeInfo.h @@ -1,201 +1,200 @@ -#ifndef Test_TreeInfo_h -#define Test_TreeInfo_h +#ifndef CepGen_IO_ROOTTreeInfo_h +#define CepGen_IO_ROOTTreeInfo_h #include "TFile.h" #include "TTree.h" -#include "Math/Vector3D.h" -#include "Math/Vector4D.h" #include #include +#include namespace ROOT { /// All useful information about a generation run class CepGenRun { public: static constexpr const char* TREE_NAME = "run"; ///< Output tree name double sqrt_s; ///< Centre of mass energy for beam particles double xsect; ///< Process cross section, in pb double errxsect; ///< Uncertainty on process cross section, in pb unsigned int num_events; ///< Number of events generated in run unsigned int litigious_events; ///< Number of litigious events in run CepGenRun() { clear(); } /// Reinitialise the run tree void clear() { sqrt_s = -1.; xsect = errxsect = -1.; num_events = litigious_events = 0; } /// Populate the run tree void create() { tree_ = std::make_shared( TREE_NAME, "a tree containing information on the previous run" ); if ( !tree_ ) throw std::runtime_error( "Failed to create the run TTree!" ); tree_->Branch( "xsect", &xsect, "xsect/D" ); tree_->Branch( "errxsect", &errxsect, "errxsect/D" ); tree_->Branch( "num_events", &num_events, "num_events/i" ); tree_->Branch( "litigious_events", &litigious_events, "litigious_events/i" ); tree_->Branch( "sqrt_s", &sqrt_s, "sqrt_s/D" ); } /// Retrieve the ROOT tree TTree* tree() { return tree_.get(); } /// Fill the run tree void fill() { tree_->Fill(); } /// Attach the run tree reader to a given file void attach( const char* filename, const char* run_tree = TREE_NAME ) { attach( TFile::Open( filename ), run_tree ); } /// Attach the run tree reader to a given tree void attach( TFile* file, const char* run_tree = TREE_NAME ) { //--- special constructor to avoid the memory to be cleared at destruction time tree_ = std::shared_ptr( dynamic_cast( file->Get( run_tree ) ), [=]( TTree* ){} ); if ( !tree_ ) throw std::runtime_error( "Failed to attach to the run TTree!" ); tree_->SetBranchAddress( "xsect", &xsect ); tree_->SetBranchAddress( "errxsect", &errxsect ); tree_->SetBranchAddress( "num_events", &num_events ); tree_->SetBranchAddress( "litigious_events", &litigious_events ); tree_->SetBranchAddress( "sqrt_s", &sqrt_s ); if ( tree_->GetEntriesFast() > 1 ) std::cerr << "The run tree has more than one entry." << std::endl; tree_->GetEntry( 0 ); } private: /// ROOT tree used for storage/retrieval of this run information std::shared_ptr tree_; }; /// All useful information about a generated event class CepGenEvent { public: // book a sufficienly large number to allow the large multiplicity // of excited proton fragmentation products static constexpr unsigned short maxpart = 5000; ///< Maximal number of particles in event static constexpr const char* TREE_NAME = "events"; ///< Output tree name float gen_time; ///< Event generation time float tot_time; ///< Total event generation time int nremn_ch[2], nremn_nt[2]; int np; ///< Number of particles in the event double pt[maxpart]; ///< Particles transverse momentum double eta[maxpart]; ///< Particles pseudo-rapidity double phi[maxpart]; ///< Particles azimutal angle double rapidity[maxpart]; ///< Particles rapidity double E[maxpart]; ///< Particles energy, in GeV double m[maxpart]; ///< Particles mass, in GeV/c\f${}^2\f$ double charge[maxpart]; ///< Particles charges, in e int pdg_id[maxpart]; ///< Integer particles PDG id int parent1[maxpart]; ///< First particles mother int parent2[maxpart]; ///< Last particles mother int stable[maxpart]; ///< Whether the particle must decay or not int role[maxpart]; ///< Particles role in the event int status[maxpart]; ///< Integer status code CepGenEvent() { clear(); } /// Reinitialise the event content void clear() { gen_time = tot_time = 0.; for ( unsigned short i = 0; i < 2; ++i ) nremn_ch[i] = nremn_nt[i] = 0; np = 0; for ( unsigned short i = 0; i < maxpart; ++i ) { pt[i] = eta[i] = phi[i] = rapidity[i] = E[i] = m[i] = charge[i] = 0.; pdg_id[i] = parent1[i] = parent2[i] = stable[i] = role[i] = status[i] = 0; } } /// Retrieve the ROOT tree TTree* tree() { return tree_.get(); } /// Fill the tree with a new event void fill() { if ( !tree_ ) throw std::runtime_error( "CepGenEvent: Trying to fill a non-existent tree!" ); tree_->Fill(); clear(); } /// Populate the tree and all associated branches void create() { tree_ = std::make_shared( TREE_NAME, "a tree containing information on events generated in previous run" ); if ( !tree_ ) throw std::runtime_error( "Failed to create the events TTree!" ); tree_->Branch( "npart", &np, "npart/I" ); tree_->Branch( "nremn_charged", nremn_ch, "nremn_charged[2]/I" ); tree_->Branch( "nremn_neutral", nremn_nt, "nremn_neutral[2]/I" ); tree_->Branch( "role", role, "role[npart]/I" ); tree_->Branch( "pt", pt, "pt[npart]/D" ); tree_->Branch( "eta", eta, "eta[npart]/D" ); tree_->Branch( "phi", phi, "phi[npart]/D" ); tree_->Branch( "rapidity", rapidity, "rapidity[npart]/D" ); tree_->Branch( "E", E, "E[npart]/D" ); tree_->Branch( "m", m, "m[npart]/D" ); tree_->Branch( "charge", charge, "charge[npart]/D" ); tree_->Branch( "pdg_id", pdg_id, "pdg_id[npart]/I" ); tree_->Branch( "parent1", parent1, "parent1[npart]/I" ); tree_->Branch( "parent2", parent2, "parent2[npart]/I" ); tree_->Branch( "stable", stable, "stable[npart]/I" ); tree_->Branch( "status", status, "status[npart]/I" ); tree_->Branch( "generation_time", &gen_time, "generation_time/F" ); tree_->Branch( "total_time", &tot_time, "total_time/F" ); } /// Attach the event tree reader to a given file void attach( const char* filename, const char* events_tree = TREE_NAME ) { file_.reset( TFile::Open( filename ) ); attach( file_.get(), events_tree ); } /// Attach the event tree reader to a given ROOT file void attach( TFile* f, const char* events_tree = TREE_NAME ) { //--- special constructor to avoid the memory to be cleared at destruction time tree_ = std::shared_ptr( dynamic_cast( f->Get( events_tree ) ), [=]( TTree* ){} ); attach(); } /// Attach the event tree reader to a given tree void attach() { if ( !tree_ ) throw std::runtime_error( "Failed to attach to the events TTree!" ); tree_->SetBranchAddress( "npart", &np ); tree_->SetBranchAddress( "nremn_charged", nremn_ch ); tree_->SetBranchAddress( "nremn_neutral", nremn_ch ); tree_->SetBranchAddress( "role", role ); tree_->SetBranchAddress( "pt", pt ); tree_->SetBranchAddress( "eta", eta ); tree_->SetBranchAddress( "phi", phi ); tree_->SetBranchAddress( "rapidity", rapidity ); tree_->SetBranchAddress( "E", E ); tree_->SetBranchAddress( "m", m ); tree_->SetBranchAddress( "charge", charge ); tree_->SetBranchAddress( "pdg_id", pdg_id ); tree_->SetBranchAddress( "parent1", parent1 ); tree_->SetBranchAddress( "parent2", parent2 ); tree_->SetBranchAddress( "stable", stable ); tree_->SetBranchAddress( "status", status ); tree_->SetBranchAddress( "generation_time", &gen_time ); tree_->SetBranchAddress( "total_time", &tot_time ); } private: /// Tree for which the event is booked std::shared_ptr tree_; std::unique_ptr file_; }; } constexpr const char* ROOT::CepGenRun::TREE_NAME; constexpr const char* ROOT::CepGenEvent::TREE_NAME; #endif diff --git a/CepGen/IO/TextHandler.cpp b/CepGen/IO/TextHandler.cpp new file mode 100644 index 0000000..6858607 --- /dev/null +++ b/CepGen/IO/TextHandler.cpp @@ -0,0 +1,193 @@ +#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/Event/EventBrowser.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: + std::string textHistogram( const std::string&, const gsl_histogram* ) const; + + static constexpr size_t PLOT_WIDTH = 50; + static constexpr char PLOT_CHAR = '#'; + + std::ofstream file_, hist_file_; + std::string hist_filename_; + //--- variables definition + const std::vector variables_; + const bool save_banner_, save_variables_; + const bool show_hists_, save_hists_; + const std::string separator_; + + const utils::EventBrowser browser_; + + std::ostringstream oss_vars_; + + double xsec_; + + //--- kinematic variables + double sqrts_; + unsigned long num_evts_; + struct gsl_histogram_deleter + { + void operator()( gsl_histogram* h ) { + gsl_histogram_free( h ); + } + }; + std::vector > > hists_; + }; + + TextHandler::TextHandler( const ParametersList& params ) : + GenericExportHandler( "text" ), + file_ ( params.get( "filename", "output.txt" ) ), + hist_filename_ ( params.get( "histFilename", "output.hists.txt" ) ), + variables_ ( params.get >( "variables" ) ), + save_banner_ ( params.get( "saveBanner", true ) ), + save_variables_( params.get( "saveVariables", true ) ), + show_hists_ ( params.get( "showHistograms", true ) ), + save_hists_ ( params.get( "saveHistograms", false ) ), + separator_ ( params.get( "separator", "\t" ) ), + xsec_( 1. ) + { + //--- first extract list of variables to store in output file + oss_vars_.clear(); + std::string sep; + for ( const auto& var : variables_ ) + oss_vars_ << sep << var, sep = separator_; + //--- then extract list of variables to be plotted in histogram + const auto& hist_vars = params.get( "histVariables" ); + for ( const auto& var : hist_vars.keys() ) { + const auto& hvar = hist_vars.get( var ); + const int nbins = hvar.get( "nbins", 10 ); + const double min = hvar.get( "low", 0. ), max = hvar.get( "high", 1. ); + auto hist = gsl_histogram_alloc( nbins ); + gsl_histogram_set_ranges_uniform( hist, min, max ); + hists_.emplace_back( std::make_pair( var, hist ) ); + CG_INFO( "TextHandler" ) + << "Booking a histogram with " << nbins << " bin" << utils::s( nbins ) + << " between " << min << " and " << max << " for \"" << var << "\"."; + } + if ( save_hists_ && !hists_.empty() ) + hist_file_.open( hist_filename_ ); + } + + TextHandler::~TextHandler() + { + //--- finalisation of the output file + file_.close(); + //--- histograms printout + if ( !show_hists_ && !save_hists_ ) + return; + for ( const auto& h_var : hists_ ) { + const auto& hist = h_var.second.get(); + gsl_histogram_scale( hist, xsec_/( num_evts_+1 ) ); + const auto& h_out = textHistogram( h_var.first, hist ); + if ( show_hists_ ) + CG_INFO( "TextHandler" ) << h_out; + if ( save_hists_ ) + hist_file_ << "\n" << h_out << "\n"; + } + if ( save_hists_ ) + CG_INFO( "TextHandler" ) + << "Saved " << utils::s( "histogram", hists_.size() ) + << " into \"" << hist_filename_ << "\"."; + } + + void + TextHandler::initialise( const Parameters& params ) + { + sqrts_ = params.kinematics.sqrtS(); + num_evts_ = 0ul; + if ( save_banner_ ) + file_ << banner( params, "#" ) << "\n"; + if ( save_variables_ ) + file_ << "# " << oss_vars_.str() << "\n"; + if ( save_hists_ && !hists_.empty() ) + hist_file_ << banner( params, "#" ) << "\n"; + } + + void + TextHandler::operator<<( const Event& ev ) + { + //--- write down the variables list in the file + std::string sep; + for ( const auto& var : variables_ ) + file_ << sep << browser_.get( ev, var ), sep = separator_; + file_ << "\n"; + //--- increment the corresponding histograms + for ( const auto& h_var : hists_ ) + gsl_histogram_increment( h_var.second.get(), browser_.get( ev, h_var.first ) ); + ++num_evts_; + } + + std::string + TextHandler::textHistogram( const std::string& var, const gsl_histogram* hist ) const + { + std::ostringstream os; + const size_t nbins = gsl_histogram_bins( hist ); + const double max_bin = gsl_histogram_max_val( hist ); + const double inv_max_bin = max_bin > 0. ? 1./max_bin : 0.; + const std::string sep( 17, ' ' ); + os + << "plot of \"" << var << "\"\n" + << sep << std::string( PLOT_WIDTH-15-var.size(), ' ' ) + << "d(sig)/d" << var << " (pb/bin)\n" + << sep << Form( "%-5.2f", gsl_histogram_min_val( hist ) ) + << std::string( PLOT_WIDTH-11, ' ' ) + << Form( "%5.2e", gsl_histogram_max_val( hist ) ) << "\n" + << sep << std::string( PLOT_WIDTH+2, '.' ); // abscissa axis + for ( size_t i = 0; i < nbins; ++i ) { + double min, max; + gsl_histogram_get_range( hist, i, &min, &max ); + const double value = gsl_histogram_get( hist, i ); + const int val = value*PLOT_WIDTH*inv_max_bin; + os + << "\n" << Form( "[%7.2f,%7.2f):", min, max ) + << std::string( val, PLOT_CHAR ) << std::string( PLOT_WIDTH-val, ' ' ) + << ": " << Form( "%6.2e", value ); + } + const double bin_width = ( gsl_histogram_max( hist )-gsl_histogram_min( hist ) )/nbins; + os + << "\n" + << Form( "%17s", var.c_str() ) << ":" << std::string( PLOT_WIDTH, '.' ) << ":\n" // 2nd abscissa axis + << "\t(" + << "bin width=" << bin_width << " unit" << utils::s( (int)bin_width ) << ", " + << "mean=" << gsl_histogram_mean( hist ) << ", " + << "st.dev.=" << gsl_histogram_sigma( hist ) + << ")"; + return os.str(); + } + } +} + +REGISTER_IO_MODULE( text, TextHandler ) diff --git a/CepGen/Parameters.h b/CepGen/Parameters.h index 7411797..75d6dd7 100644 --- a/CepGen/Parameters.h +++ b/CepGen/Parameters.h @@ -1,147 +1,157 @@ #ifndef CepGen_Parameters_h #define CepGen_Parameters_h #include "CepGen/Physics/Kinematics.h" #include #include #include namespace cepgen { class Event; class ParametersList; namespace proc { class GenericProcess; } namespace hadr { class GenericHadroniser; } + namespace io { class GenericExportHandler; } namespace utils { class TamingFunctionsCollection; } enum class IntegratorType; /// List of parameters used to start and run the simulation job class Parameters { public: Parameters(); /// Copy constructor (transfers ownership to the process/hadroniser!) Parameters( Parameters& ); /// Const copy constructor (all but the process and the hadroniser) Parameters( const Parameters& ); ~Parameters(); // required for unique_ptr initialisation! /// Assignment operator Parameters& operator=( 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 terminal friend std::ostream& operator<<( std::ostream&, const Parameters* ); std::shared_ptr general; //----- process to compute /// Process for which the cross-section will be computed and the events will be generated proc::GenericProcess* process(); /// Process for which the cross-section will be computed and the events will be generated const proc::GenericProcess* process() const; /// Name of the process considered std::string processName() const; /// Set the process to study void setProcess( std::unique_ptr proc ); /// Set the process to study void setProcess( proc::GenericProcess* proc ); //----- events kinematics /// Events kinematics for phase space definition Kinematics kinematics; //----- VEGAS /// Collection of integrator parameters struct Integration { Integration(); Integration( const Integration& ); ~Integration(); IntegratorType type; unsigned int ncvg; ///< Number of function calls to be computed for each point long rng_seed; ///< Random number generator seed gsl_rng_type* rng_engine; ///< Random number generator engine gsl_monte_vegas_params vegas; double vegas_chisq_cut; gsl_monte_miser_params miser; double result, err_result; }; Integration& integration() { return integration_; } const Integration& integration() const { return integration_; } //----- events generation /// Collection of events generation parameters struct Generation { Generation(); Generation( const Generation& ); bool enabled; ///< Are we generating events ? (true) or are we only computing the cross-section ? (false) unsigned long maxgen; ///< Maximal number of events to generate in this run bool symmetrise; ///< Do we want the events to be symmetrised with respect to the \f$z\f$-axis ? bool treat; ///< Is the integrand to be smoothed for events generation? unsigned int gen_print_every; ///< Frequency at which the events are displayed to the end-user unsigned int num_threads; ///< Number of threads to perform the integration unsigned int num_points; ///< Number of points to "shoot" in each integration bin by the algorithm }; Generation& generation() { return generation_; } const Generation& generation() const { return 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_; } + /// Set a new output module definition + void setOutputModule( std::unique_ptr mod ); + /// Set the pointer to a output module + void setOutputModule( io::GenericExportHandler* mod ); + /// Output module definition + io::GenericExportHandler* outputModule(); + //----- hadronisation algorithm /// Hadronisation algorithm to use for the proton(s) fragmentation hadr::GenericHadroniser* hadroniser(); /// Name of the hadroniser (if applicable) std::string hadroniserName() const; /// Set the hadronisation algorithm void setHadroniser( std::unique_ptr hadr ); /// Set the hadronisation algorithm void setHadroniser( hadr::GenericHadroniser* hadr ); //----- taming functions /// Functionals to be used to account for rescattering corrections (implemented within the process) std::shared_ptr taming_functions; //----- run operations /// Reset the total generation time and the number of events generated for this run, prepare kinematics void prepareRun(); /// Add a new timing into the total generation time /// \param[in] gen_time Time to add (in seconds) void addGenerationTime( double gen_time ); /// Return the total generation time for this run (in seconds) inline double totalGenerationTime() const { return total_gen_time_; } /// Total number of events already generated in this run inline unsigned int numGeneratedEvents() const { return num_gen_events_; } private: std::unique_ptr process_; std::unique_ptr hadroniser_; + /// Storage object + std::unique_ptr out_module_; bool store_; /// Total generation time (in seconds) double total_gen_time_; /// Number of events already generated unsigned long num_gen_events_; /// Integrator parameters Integration integration_; /// Events generation parameters Generation generation_; }; } #endif diff --git a/CepGen/Physics/Limits.cpp b/CepGen/Physics/Limits.cpp index 5ea5805..af80887 100644 --- a/CepGen/Physics/Limits.cpp +++ b/CepGen/Physics/Limits.cpp @@ -1,132 +1,154 @@ #include "CepGen/Physics/Limits.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/utils.h" namespace cepgen { Limits::Limits( double min, double max ) : std::pair( min, max ) {} Limits::Limits( const Limits& rhs ) : std::pair( rhs.first, rhs.second ) {} + Limits + Limits::operator-() const + { + return Limits( -second, -first ); + } + Limits& Limits::operator+=( double c ) { first += c; second += c; return *this; } Limits& + Limits::operator-=( double c ) + { + first -= c; + second -= c; + return *this; + } + + Limits& Limits::operator*=( double c ) { first *= c; second *= c; if ( c < 0. ) { const double tmp = first; // will be optimised by compiler anyway... first = second; second = tmp; } return *this; } void Limits::in( double low, double up ) { first = low; second = up; } double Limits::range() const { if ( !hasMin() || !hasMax() ) return 0.; return second-first; } bool Limits::hasMin() const { return first != INVALID; } bool Limits::hasMax() const { return second != INVALID; } bool Limits::passes( double val ) const { if ( hasMin() && val < min() ) return false; if ( hasMax() && val > max() ) return false; return true; } bool Limits::valid() const { return hasMin() || hasMax(); } void Limits::save( bool& on, double& lmin, double& lmax ) const { on = false; lmin = lmax = 0.; if ( !valid() ) return; on = true; if ( hasMin() ) lmin = min(); if ( hasMax() ) lmax = max(); if ( lmin == lmax ) on = false; } double Limits::x( double v ) const { if ( v < 0. || v > 1. ) throw CG_ERROR( "Limits:shoot" ) << "x must be comprised between 0 and 1; x value = " << v << "."; if ( !valid() ) return INVALID; return first + ( second-first ) * v; } std::ostream& operator<<( std::ostream& os, const Limits& lim ) { if ( !lim.hasMin() && !lim.hasMax() ) return os << "no cuts"; if ( !lim.hasMin() ) return os << Form( "below %g", lim.max() ); if ( !lim.hasMax() ) return os << Form( "above %g", lim.min() ); return os << Form( "%g to %g", lim.min(), lim.max() ); } + Limits operator+( Limits lim, double c ) { lim += c; return lim; } Limits + operator-( Limits lim, double c ) + { + lim -= c; + return lim; + } + + Limits operator*( Limits lim, double c ) { lim *= c; return lim; } } diff --git a/CepGen/Physics/Limits.h b/CepGen/Physics/Limits.h index 7c982b7..1ec246a 100644 --- a/CepGen/Physics/Limits.h +++ b/CepGen/Physics/Limits.h @@ -1,56 +1,59 @@ #ifndef CepGen_Physics_Limits_h #define CepGen_Physics_Limits_h #include #include namespace cepgen { /// Validity interval for a variable class Limits : private std::pair { public: /// Define lower and upper limits on a quantity Limits( double min = INVALID, double max = INVALID ); Limits( const Limits& ); + Limits operator-() const; Limits& operator+=( double c ); + Limits& operator-=( double c ); Limits& operator*=( double c ); friend Limits operator+( Limits lim, double c ); + friend Limits operator-( Limits lim, double c ); friend Limits operator*( Limits lim, double c ); /// Lower limit to apply on the variable double min() const { return first; } /// Lower limit to apply on the variable double& min() { return first; } /// Upper limit to apply on the variable double max() const { return second; } /// Upper limit to apply on the variable double& max() { return second; } /// Export the limits into external variables void save( bool& on, double& lmin, double& lmax ) const; /// Find the [0,1] value scaled between minimum and maximum double x( double v ) const; /// Specify the lower and upper limits on the variable void in( double low, double up ); /// Full variable range allowed double range() const; /// Have a lower limit? bool hasMin() const; /// Have an upper limit? bool hasMax() const; /// Check if the value is inside limits' boundaries bool passes( double val ) const; /// Is there a lower and upper limit? bool valid() const; /// Human-readable expression of the limits friend std::ostream& operator<<( std::ostream&, const Limits& ); private: static constexpr double INVALID = -999.999; }; } #endif diff --git a/CepGen/Processes/CMakeLists.txt b/CepGen/Processes/CMakeLists.txt index 7489f0c..7ed064e 100644 --- a/CepGen/Processes/CMakeLists.txt +++ b/CepGen/Processes/CMakeLists.txt @@ -1,14 +1,15 @@ file(GLOB sources *.cpp Fortran/*.f) include_directories(${PROJECT_SOURCE_DIR}) include_directories(Fortran) set(F77_PROC_PATH ${PROJECT_SOURCE_DIR}/External/Processes) file(GLOB oth_proc_sources ${F77_PROC_PATH}/*.f ${F77_PROC_PATH}/CepGenWrapper.cpp) if(oth_proc_sources) + #message(STATUS "Fortran processes found: ${oth_proc_sources}") list(APPEND sources ${oth_proc_sources}) endif() add_library(CepGenProcesses SHARED ${sources}) install(TARGETS CepGenProcesses DESTINATION lib) diff --git a/CepGen/Processes/Fortran/KTBlocks.inc b/CepGen/Processes/Fortran/KTBlocks.inc index da5c605..dd8424b 100644 --- a/CepGen/Processes/Fortran/KTBlocks.inc +++ b/CepGen/Processes/Fortran/KTBlocks.inc @@ -1,61 +1,67 @@ !> \file KTBblocks.inc c ================================================================ c F77 process-to-CepGen helper c > this header file defines a set of common blocks useful for c > the interfacing of any kt-factorised matrix element computation c > to a mother CepGen instance c ================================================================ c ================================================================= !> collection of fundamental physics constants common/constants/am_p,units,pi,alpha_em double precision am_p,units,pi,alpha_em c ================================================================= c information on the full process c inp1 = proton energy in lab frame c inp2 = nucleus energy **per nucleon** in LAB frame c Collision is along z-axis c ================================================================= - common/genparams/icontri,iflux1,iflux2,imethod,sfmod,pdg_l, - & a_nuc1,z_nuc1,a_nuc2,z_nuc2, - & inp1,inp2 - integer icontri,iflux1,iflux2,imethod,sfmod,pdg_l, + common/genparams/icontri,iflux1,iflux2,sfmod, + & a_nuc1,z_nuc1,a_nuc2,z_nuc2,inp1,inp2 + integer icontri,iflux1,iflux2,sfmod, & a_nuc1,z_nuc1,a_nuc2,z_nuc2 double precision inp1,inp2 c ================================================================= c kt-factorisation kinematics c ================================================================= common/ktkin/q1t,q2t,phiq1t,phiq2t,y1,y2,ptdiff,phiptdiff, & am_x,am_y double precision q1t,q2t,phiq1t,phiq2t,y1,y2,ptdiff,phiptdiff, & am_x,am_y c ================================================================= c phase space cuts c ================================================================= common/kincuts/ipt,iene,ieta,iinvm,iptsum,idely, & pt_min,pt_max,ene_min,ene_max,eta_min,eta_max, & invm_min,invm_max,ptsum_min,ptsum_max, & dely_min,dely_max logical ipt,iene,ieta,iinvm,iptsum,idely double precision pt_min,pt_max,ene_min,ene_max,eta_min,eta_max, & invm_min,invm_max,ptsum_min,ptsum_max, & dely_min,dely_max c ================================================================= c generated event kinematics c ================================================================= common/evtkin/nout,ipdg,idum,pc,px,py integer nout,idum,ipdg(10) double precision pc(10,4),px(4),py(4) c ================================================================= c helpers for the evaluation of unintegrated parton fluxes c ================================================================= external CepGen_kT_flux,CepGen_kT_flux_HI double precision CepGen_kT_flux,CepGen_kT_flux_HI external CepGen_particle_charge,CepGen_particle_mass double precision CepGen_particle_charge,CepGen_particle_mass +c ================================================================= +c helpers for the retrieval of input parameters from the process +c ================================================================= + external CepGen_param_int,CepGen_param_real + integer CepGen_param_int + double precision CepGen_param_real + diff --git a/CepGen/Processes/Fortran/KTStructures.h b/CepGen/Processes/Fortran/KTStructures.h index cad251b..a9d7b3f 100644 --- a/CepGen/Processes/Fortran/KTStructures.h +++ b/CepGen/Processes/Fortran/KTStructures.h @@ -1,77 +1,75 @@ #ifndef CepGen_Processes_Fortran_KTStructures_h #define CepGen_Processes_Fortran_KTStructures_h namespace cepgen { /// Collection of common blocks for Fortran \f$k_{\rm T}\f$-processes namespace ktblock { /// General physics constants struct Constants { double m_p; ///< Proton mass double units; ///< Conversion factor GeV\f$^2\to\f$ barn double pi; ///< \f$\pi\f$ double alpha_em; ///< Electromagnetic coupling constant }; /// Generic run parameters struct Parameters { int icontri; ///< Kinematics mode int iflux1; ///< Type of \f$k_{\rm T}\f$-factorised flux for first incoming parton int iflux2; ///< Type of \f$k_{\rm T}\f$-factorised flux for second incoming parton - int imethod; ///< Computation method for matrix element int sfmod; ///< Structure functions modelling - int pdg_l; ///< Central system PDG id int a_nuc1; ///< First beam mass number int z_nuc1; ///< First beam atomic number int a_nuc2; ///< Second beam mass number int z_nuc2; ///< Second beam atomic number double inp1; ///< First beam momentum, in GeV/c double inp2; ///< Second beam momentum, in GeV/c }; /// Kinematics properties of the \f$k_{\rm T}\f$-factorised process struct KTKinematics { double q1t; ///< Transverse momentum of the first incoming parton double q2t; ///< Transverse momentum of the second incoming parton double phiq1t; ///< Azimutal angle of the first incoming parton double phiq2t; ///< Azimutal angle of the second incoming parton double y1; ///< First incoming parton rapidity double y2; ///< Second incoming parton rapidity double ptdiff; ///< Central system pT balance double phiptdiff; ///< Central system azimutal angle difference double m_x; ///< Invariant mass for the first diffractive state double m_y; ///< Invariant mass for the second diffractive state }; /// Phase space cuts for event kinematics struct Cuts { bool ipt; ///< Switch for cut on single particle transverse momentum bool iene; ///< Switch for cut on single particle energy bool ieta; ///< Switch for cut on single particle pseudo-rapidity bool iinvm; ///< Switch for cut on central system invariant mass bool iptsum; ///< Switch for cut on central system transverse momentum bool idely; ///< Switch for cut on rapididty difference double pt_min; ///< Minimal single particle transverse momentum double pt_max; ///< Maximal single particle transverse momentum double ene_min; ///< Minimal single particle energy double ene_max; ///< Maximal single particle energy double eta_min; ///< Minimal single particle pseudo-rapidity double eta_max; ///< Maximal single particle pseudo-rapidity double invm_min; ///< Minimal central system invariant mass double invm_max; ///< Maximal central system invariant mass double ptsum_min; ///< Minimal central system transverse momentum double ptsum_max; ///< Maximal central system transverse momentum double dely_min; ///< Minimal rapidity difference for central system double dely_max; ///< Maximal rapidity difference for central system }; /// Single event kinematics struct Event { int nout; ///< Number of particles in central system int pdg[10]; ///< PDG ids of all particles in central system int idum; ///< Padding double pc[4][10]; ///< 4-momenta of all particles in central system double px[4]; ///< 4-momentum of first outgoing proton state double py[4]; ///< 4-momentum of second outgoing proton state }; } } #endif diff --git a/CepGen/Processes/Fortran/cepgen_print.f b/CepGen/Processes/Fortran/cepgen_print.f index af7a984..b99fe70 100644 --- a/CepGen/Processes/Fortran/cepgen_print.f +++ b/CepGen/Processes/Fortran/cepgen_print.f @@ -1,42 +1,43 @@ !> \file cepgen_print.f subroutine cepgen_print !> Print useful run information in standard stream implicit none include 'KTBlocks.inc' logical params_shown data params_shown/.false./ save params_shown if(params_shown) return print *,'========================================================' print *,'Parameter value(s)' print *,'--------------------------------------------------------' print 101,'Process mode:',icontri - print 101,'Computation method:',imethod print 101,'Structure functions:',sfmod - print 101,'Central system PDG:',pdg_l print 103,'Beams momenta:',inp1,inp2 print 102,'Fluxes modes:',iflux1,iflux2 print 104,'Beams (A,Z):',a_nuc1,z_nuc1,a_nuc2,z_nuc2 print *,'========================================================' print *,'Cut enabled minimum maximum' print *,'--------------------------------------------------------' print 100,'pt(single)',ipt,pt_min,pt_max print 100,'energy(single)',iene,ene_min,ene_max print 100,'eta(single)',ieta,eta_min,eta_max print 100,'m(sum)',iinvm,invm_min,invm_max print 100,'pt(sum)',iptsum,ptsum_min,ptsum_max print 100,'delta(y)',idely,dely_min,dely_max print *,'========================================================' + print *,'Process-specific parameters' + call cepgen_list_params + print *,'========================================================' params_shown=.true. 100 format(A26,' ',L2,f12.4,f12.4) 101 format(A33,I12) 102 format(A33,I12,I12) 103 format(A33,f12.2,f12.2) 104 format(A33,' (',I3,',',I3,'), (',I3,',',I3,')') end diff --git a/CepGen/Processes/FortranKTProcess.cpp b/CepGen/Processes/FortranKTProcess.cpp index aae8742..eb44546 100644 --- a/CepGen/Processes/FortranKTProcess.cpp +++ b/CepGen/Processes/FortranKTProcess.cpp @@ -1,181 +1,180 @@ #include "CepGen/Processes/FortranKTProcess.h" #include "CepGen/Processes/Fortran/KTStructures.h" #include "CepGen/Core/ParametersList.h" #include "CepGen/Core/Exception.h" #include "CepGen/StructureFunctions/StructureFunctions.h" #include "CepGen/Event/Event.h" #include "CepGen/Physics/KTFlux.h" #include "CepGen/Physics/Constants.h" #include "CepGen/Physics/PDG.h" extern "C" { extern cepgen::ktblock::Constants constants_; extern cepgen::ktblock::Parameters genparams_; extern cepgen::ktblock::KTKinematics ktkin_; extern cepgen::ktblock::Cuts kincuts_; extern cepgen::ktblock::Event evtkin_; } namespace cepgen { namespace proc { + ParametersList + FortranKTProcess::kProcParameters; + FortranKTProcess::FortranKTProcess( const ParametersList& params, const char* name, const char* descr, std::function func ) : GenericKTProcess( params, name, descr, { { PDG::photon, PDG::photon } }, { PDG::muon, PDG::muon } ), - pair_( params.get( "pair" ).pdgid ), - method_( params.get( "method", 1 ) ), func_( func ) { constants_.m_p = GenericProcess::mp_; constants_.units = constants::GEVM2_TO_PB; constants_.pi = M_PI; constants_.alpha_em = constants::ALPHA_EM; } void FortranKTProcess::preparePhaseSpace() { mom_ip1_ = event_->getOneByRole( Particle::IncomingBeam1 ).momentum(); mom_ip2_ = event_->getOneByRole( Particle::IncomingBeam2 ).momentum(); registerVariable( y1_, Mapping::linear, kin_.cuts.central.rapidity_single, { -6., 6. }, "First central particle rapidity" ); registerVariable( y2_, Mapping::linear, kin_.cuts.central.rapidity_single, { -6., 6. }, "Second central particle rapidity" ); registerVariable( pt_diff_, Mapping::linear, kin_.cuts.central.pt_diff, { 0., 50. }, "Transverse momentum difference between central particles" ); registerVariable( phi_pt_diff_, Mapping::linear, kin_.cuts.central.phi_pt_diff, { 0., 2.*M_PI }, "Central particles azimuthal angle difference" ); //=========================================================================================== // feed phase space cuts to the common block //=========================================================================================== kin_.cuts.central.pt_single.save( kincuts_.ipt, kincuts_.pt_min, kincuts_.pt_max ); kin_.cuts.central.energy_single.save( kincuts_.iene, kincuts_.ene_min, kincuts_.ene_max ); kin_.cuts.central.eta_single.save( kincuts_.ieta, kincuts_.eta_min, kincuts_.eta_max ); kin_.cuts.central.mass_sum.save( kincuts_.iinvm, kincuts_.invm_min, kincuts_.invm_max ); kin_.cuts.central.pt_sum.save( kincuts_.iptsum, kincuts_.ptsum_min, kincuts_.ptsum_max ); kin_.cuts.central.rapidity_diff.save( kincuts_.idely, kincuts_.dely_min, kincuts_.dely_max ); //=========================================================================================== // feed run parameters to the common block //=========================================================================================== genparams_.icontri = (int)kin_.mode; - genparams_.imethod = method_; genparams_.sfmod = (int)kin_.structure_functions->type; - genparams_.pdg_l = pair_; //------------------------------------------------------------------------------------------- // incoming beams information //------------------------------------------------------------------------------------------- genparams_.inp1 = kin_.incoming_beams.first.pz; genparams_.inp2 = kin_.incoming_beams.second.pz; const HeavyIon in1 = (HeavyIon)kin_.incoming_beams.first.pdg; if ( in1 ) { genparams_.a_nuc1 = in1.A; genparams_.z_nuc1 = (unsigned short)in1.Z; if ( genparams_.z_nuc1 > 1 ) { event_->getOneByRole( Particle::IncomingBeam1 ).setPdgId( (pdgid_t)in1 ); event_->getOneByRole( Particle::OutgoingBeam1 ).setPdgId( (pdgid_t)in1 ); } } else genparams_.a_nuc1 = genparams_.z_nuc1 = 1; const HeavyIon in2 = (HeavyIon)kin_.incoming_beams.second.pdg; if ( in2 ) { genparams_.a_nuc2 = in2.A; genparams_.z_nuc2 = (unsigned short)in2.Z; if ( genparams_.z_nuc2 > 1 ) { event_->getOneByRole( Particle::IncomingBeam2 ).setPdgId( (pdgid_t)in2 ); event_->getOneByRole( Particle::OutgoingBeam2 ).setPdgId( (pdgid_t)in2 ); } } else genparams_.a_nuc2 = genparams_.z_nuc2 = 1; //------------------------------------------------------------------------------------------- // intermediate partons information //------------------------------------------------------------------------------------------- genparams_.iflux1 = (int)kin_.incoming_beams.first.kt_flux; genparams_.iflux2 = (int)kin_.incoming_beams.second.kt_flux; switch ( (KTFlux)genparams_.iflux1 ) { case KTFlux::P_Gluon_KMR: event_->getOneByRole( Particle::Parton1 ).setPdgId( PDG::gluon ); break; case KTFlux::P_Photon_Elastic: case KTFlux::P_Photon_Inelastic: case KTFlux::P_Photon_Inelastic_Budnev: case KTFlux::HI_Photon_Elastic: event_->getOneByRole( Particle::Parton2 ).setPdgId( PDG::photon ); break; case KTFlux::invalid: throw CG_FATAL( "FortranKTProcess" ) << "Invalid flux for 2nd incoming parton: " << genparams_.iflux2 << "!"; } switch ( (KTFlux)genparams_.iflux2 ) { case KTFlux::P_Gluon_KMR: event_->getOneByRole( Particle::Parton2 ).setPdgId( PDG::gluon ); break; case KTFlux::P_Photon_Elastic: case KTFlux::P_Photon_Inelastic: case KTFlux::P_Photon_Inelastic_Budnev: case KTFlux::HI_Photon_Elastic: event_->getOneByRole( Particle::Parton2 ).setPdgId( PDG::photon ); break; case KTFlux::invalid: throw CG_FATAL( "FortranKTProcess" ) << "Invalid flux for 2nd incoming parton: " << genparams_.iflux2 << "!"; } } double FortranKTProcess::computeKTFactorisedMatrixElement() { ktkin_.q1t = qt1_; ktkin_.q2t = qt2_; ktkin_.phiq1t = phi_qt1_; ktkin_.phiq2t = phi_qt2_; ktkin_.y1 = y1_; ktkin_.y2 = y2_; ktkin_.ptdiff = pt_diff_; ktkin_.phiptdiff = phi_pt_diff_; ktkin_.m_x = MX_; ktkin_.m_y = MY_; //--- compute the event weight return func_(); } void FortranKTProcess::fillCentralParticlesKinematics() { //=========================================================================================== // outgoing beam remnants //=========================================================================================== PX_ = Particle::Momentum( evtkin_.px ); PY_ = Particle::Momentum( evtkin_.py ); // express these momenta per nucleon PX_ *= 1./genparams_.a_nuc1; PY_ *= 1./genparams_.a_nuc2; //=========================================================================================== // intermediate partons //=========================================================================================== const Particle::Momentum mom_par1 = mom_ip1_-PX_, mom_par2 = mom_ip2_-PY_; event_->getOneByRole( Particle::Parton1 ).setMomentum( mom_par1 ); event_->getOneByRole( Particle::Parton2 ).setMomentum( mom_par2 ); event_->getOneByRole( Particle::Intermediate ).setMomentum( mom_par1+mom_par2 ); //=========================================================================================== // central system //=========================================================================================== Particles& oc = (*event_)[Particle::CentralSystem]; for ( int i = 0; i < evtkin_.nout; ++i ) { Particle& p = oc[i]; p.setPdgId( (pdgid_t)evtkin_.pdg[i] ); p.setStatus( Particle::Status::FinalState ); p.setMomentum( Particle::Momentum( evtkin_.pc[i] ) ); } } } } diff --git a/CepGen/Processes/FortranKTProcess.h b/CepGen/Processes/FortranKTProcess.h index ea6aa42..d99c764 100644 --- a/CepGen/Processes/FortranKTProcess.h +++ b/CepGen/Processes/FortranKTProcess.h @@ -1,37 +1,37 @@ #ifndef CepGen_Processes_FortranKTProcess_h #define CepGen_Processes_FortranKTProcess_h #include "CepGen/Processes/GenericKTProcess.h" #include namespace cepgen { namespace proc { /// Compute the matrix element for a generic \f$k_{\rm T}\f$-factorised process defined in a Fortran weighting function class FortranKTProcess : public GenericKTProcess { public: FortranKTProcess( const ParametersList& params, const char* name, const char* descr, std::function func ); ProcessPtr clone( const ParametersList& /*params*/ ) const override { return ProcessPtr( new FortranKTProcess( *this ) ); } + static ParametersList kProcParameters; + private: void preparePhaseSpace() override; double computeKTFactorisedMatrixElement() override; void fillCentralParticlesKinematics() override; - int pair_; ///< Outgoing particles type - int method_; ///< Computation method for the process std::function func_; ///< Function to be called for weight computation double y1_; ///< First outgoing particle rapidity double y2_; ///< Second outgoing particle rapidity double pt_diff_; ///< Transverse momentum balance between outgoing particles double phi_pt_diff_; ///< Azimutal angle difference between outgoing particles Particle::Momentum mom_ip1_; ///< First incoming beam momentum Particle::Momentum mom_ip2_; ///< Second incoming beam momentum }; } } #endif diff --git a/CepGen/Processes/ProcessesHandler.h b/CepGen/Processes/ProcessesHandler.h index 39603a7..1b01160 100644 --- a/CepGen/Processes/ProcessesHandler.h +++ b/CepGen/Processes/ProcessesHandler.h @@ -1,40 +1,42 @@ #ifndef CepGen_Processes_ProcessesHandler_h #define CepGen_Processes_ProcessesHandler_h #include "CepGen/Core/ModuleFactory.h" #include "CepGen/Core/ParametersList.h" #include "CepGen/Processes/GenericProcess.h" #include "CepGen/Processes/FortranKTProcess.h" /** \file */ /// Add a generic process definition to the list of handled processes #define REGISTER_PROCESS( name, obj ) \ namespace cepgen { namespace proc { \ struct BUILDERNM( name ) { \ BUILDERNM( name )() { ProcessesHandler::get().registerModule( STRINGIFY( name ) ); } }; \ static BUILDERNM( name ) g ## name; \ } } /// Declare a Fortran process function name #define DECLARE_FORTRAN_FUNCTION( method ) \ extern "C" { extern double method ## _(); } #define PROCESS_F77_NAME( name ) F77_ ## name /// Add the Fortran process definition to the list of handled processes #define REGISTER_FORTRAN_PROCESS( name, method, description ) \ struct PROCESS_F77_NAME( name ) : public cepgen::proc::FortranKTProcess { \ PROCESS_F77_NAME( name )( const cepgen::ParametersList& params = cepgen::ParametersList() ) : \ - cepgen::proc::FortranKTProcess( params, STRINGIFY( name ), description, method ## _ ) {} }; \ + cepgen::proc::FortranKTProcess( params, STRINGIFY( name ), description, method ## _ ) { \ + cepgen::proc::FortranKTProcess::kProcParameters = params; \ + } }; \ REGISTER_PROCESS( name, PROCESS_F77_NAME( name ) ) namespace cepgen { namespace proc { /// A processes factory typedef ModuleFactory ProcessesHandler; } } #endif diff --git a/CepGen/README b/CepGen/README index 0f5eb38..b6ff04f 100644 --- a/CepGen/README +++ b/CepGen/README @@ -1,26 +1,26 @@ ______ ______ / ____/__ ____ / ____/__ ____ / / / _ \/ __ \/ / __/ _ \/ __ \ / /___/ __/ /_/ / /_/ / __/ / / / \____/\___/ .___/\____/\___/_/ /_/ /_/ Copyright (c) 1976 J. Vermaseren (LPAIR matrix element) (*) - 2014-2018 L. Forthomme + 2014-2019 L. Forthomme Please use https://arxiv.org/abs/1808.06059 for citations. This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see . (*) with contributions from: the ZEUS offline group, M. GVH, D. Bocian, O. Duenger, N. Schul, and others. diff --git a/CepGen/Version.h b/CepGen/Version.h index 0e42102..b9cc0b5 100644 --- a/CepGen/Version.h +++ b/CepGen/Version.h @@ -1,18 +1,18 @@ #ifndef CepGen_Version_h #define CepGen_Version_h #include namespace cepgen { /// CepGen version /// \note Format: 0xMMmmff, with /// - MM = major version /// - mm = minor version /// - ff = feature(s) release - const unsigned int cepgen_version = 0x000906; + const unsigned int cepgen_version = 0x000907; /// Human-readable version number const std::string version(); } #endif diff --git a/External/Processes/nucl_to_ff.f b/External/Processes/nucl_to_ff.f index f5b65ce..0e69d88 100644 --- a/External/Processes/nucl_to_ff.f +++ b/External/Processes/nucl_to_ff.f @@ -1,495 +1,496 @@ function nucl_to_ff() implicit none double precision nucl_to_ff c ================================================================= c CepGen common blocks for kinematics definition c ================================================================= include 'KTBlocks.inc' data iflux1,iflux2,sfmod,pdg_l/10,100,11,13/ data a_nuc1,z_nuc1,a_nuc2,z_nuc2/1,1,208,82/ c ================================================================= c local variables c ================================================================= double precision s,s12 double precision alpha1,alpha2,amt1,amt2 double precision pt1,pt2,eta1,eta2,dely double precision pt1x,pt1y,pt2x,pt2y double precision ak10,ak1x,ak1y,ak1z,ak20,ak2x,ak2y,ak2z double precision beta1,beta2,x1,x2 double precision z1p,z1m,z2p,z2m double precision q1tx,q1ty,q1z,q10,q1t2 double precision q2tx,q2ty,q2z,q20,q2t2 double precision ptsum double precision that1,that2,that,uhat1,uhat2,uhat double precision term1,term2,term3,term4,term5,term6,term7 double precision term8,term9,term10,amat2 double precision ak1_x,ak1_y,ak2_x,ak2_y double precision eps12,eps22 double precision aux2_1,aux2_2,f1,f2 double precision Phi10,Phi102,Phi11_x,Phi11_y,Phi112 double precision Phi20,Phi202,Phi21_x,Phi21_y,Phi212 double precision Phi11_dot_e,Phi11_cross_e double precision Phi21_dot_e,Phi21_cross_e double precision aintegral - integer imat1,imat2 + integer imethod,pdg_l,imat1,imat2 double precision px_plus,px_minus,py_plus,py_minus double precision r1,r2 double precision ptdiffx,ptsumx,ptdiffy,ptsumy double precision invm,invm2,s1_eff,s2_eff double precision t1abs,t2abs double precision am_l,q_l double precision inp_A,inp_B,am_A,am_B double precision p1_plus, p2_minus double precision amat2_1,amat2_2 integer iterm11,iterm22,iterm12,itermtt double precision coupling c ================================================================= c quarks production c ================================================================= #ifdef ALPHA_S double precision t_max,amu2,alphas #endif logical first_init data first_init/.true./ - save first_init,am_l,q_l - - call CepGen_print + save first_init,imethod,pdg_l,am_l,q_l if(first_init) then + call CepGen_print + imethod = CepGen_param_int('method', 1) + pdg_l = CepGen_param_int('pair', 13) am_l = CepGen_particle_mass(pdg_l) ! central particles mass q_l = CepGen_particle_charge(pdg_l) ! central particles charge if(iflux1.ge.20.and.iflux1.lt.40) then if(icontri.eq.3.or.icontri.eq.4) then print *,'Invalid process mode for collinear gluon emission!' stop endif #ifdef ALPHA_S print *,'Initialisation of the alpha(S) evolution algorithm..' call initAlphaS(0,1.d0,1.d0,0.5d0, & CepGen_particle_mass(4), ! charm & CepGen_particle_mass(5), ! bottom & CepGen_particle_mass(6)) ! top #endif endif first_init = .false. endif c ================================================================= c FIXME c ================================================================= nucl_to_ff = 0.d0 amat2 = 0.d0 eps12 = 0.d0 eps22 = 0.d0 q10 = 0.d0 q1z = 0.d0 q20 = 0.d0 q2z = 0.d0 c ================================================================= c go to energy of Nucleus = A*inp in order that generator puts out c proper momenta in LAB frame c ================================================================= inp_A = inp1*a_nuc1 am_A = am_p*a_nuc1 inp_B = inp2*a_nuc2 am_B = am_p*a_nuc2 c ================================================================= c four-momenta for incoming beams in LAB !!!! c ================================================================= r1 = dsqrt(1.d0+am_A**2/inp_A**2) r2 = dsqrt(1.d0+am_B**2/inp_B**2) ak10 = inp_A*r1 ak1x = 0.d0 ak1y = 0.d0 ak1z = inp_A ak20 = inp_B*r2 ak2x = 0.d0 ak2y = 0.d0 ak2z = -inp_B s = 4.*inp_A*inp_B*(1.d0 + r1*r2)/2.d0+(am_A**2+am_B**2) s12 = dsqrt(s) p1_plus = (ak10+ak1z)/dsqrt(2.d0) p2_minus = (ak20-ak2z)/dsqrt(2.d0) c terms in the matrix element c iterm11 = 1 ! LL iterm22 = 1 ! TT iterm12 = 1 ! LT itermtt = 1 ! TT' c ================================================================= c two terms in the Wolfgang's formula for c off-shell gamma gamma --> l^+ l^- c ================================================================= imat1 = 1 imat2 = 1 c ================================================================= c Outgoing proton final state's mass c ================================================================= if((icontri.eq.1).or.(icontri.eq.2)) am_x = am_A if((icontri.eq.1).or.(icontri.eq.3)) am_y = am_B q1tx = q1t*cos(phiq1t) q1ty = q1t*sin(phiq1t) q2tx = q2t*cos(phiq2t) q2ty = q2t*sin(phiq2t) ptsumx = q1tx+q2tx ptsumy = q1ty+q2ty ptsum = sqrt(ptsumx**2+ptsumy**2) c ================================================================= c a window in final state transverse momentum c ================================================================= if(iptsum) then if(ptsum.lt.ptsum_min.or.ptsum.gt.ptsum_max) return endif c ================================================================= c compute the individual central particles momentum c ================================================================= ptdiffx = ptdiff*cos(phiptdiff) ptdiffy = ptdiff*sin(phiptdiff) pt1x = 0.5*(ptsumx+ptdiffx) pt1y = 0.5*(ptsumy+ptdiffy) pt2x = 0.5*(ptsumx-ptdiffx) pt2y = 0.5*(ptsumy-ptdiffy) pt1 = sqrt(pt1x**2+pt1y**2) pt2 = sqrt(pt2x**2+pt2y**2) if(ipt) then if(pt1.lt.pt_min.or.pt2.lt.pt_min) return if(pt_max.gt.0d0.and.(pt2.gt.pt_max.or.pt2.gt.pt_max)) return endif amt1 = dsqrt(pt1**2+am_l**2) amt2 = dsqrt(pt2**2+am_l**2) invm2 = amt1**2 + amt2**2 + 2.d0*amt1*amt2*dcosh(y1-y2) & -ptsum**2 invm = dsqrt(invm2) c ================================================================= c a window in final state invariant mass c ================================================================= if(iinvm) then if(invm.lt.invm_min) return if(invm_max.gt.0d0.and.invm.gt.invm_max) return endif c ================================================================= c a window in rapidity distance c ================================================================= dely = dabs(y1-y2) if(idely) then if(dely.lt.dely_min.or.dely.gt.dely_max) return endif c ================================================================= c auxiliary quantities c ================================================================= alpha1 = amt1/(dsqrt(2.d0)*p1_plus)*dexp( y1) alpha2 = amt2/(dsqrt(2.d0)*p1_plus)*dexp( y2) beta1 = amt1/(dsqrt(2.d0)*p2_minus)*dexp(-y1) beta2 = amt2/(dsqrt(2.d0)*p2_minus)*dexp(-y2) q1t2 = q1tx**2 + q1ty**2 q2t2 = q2tx**2 + q2ty**2 x1 = alpha1 + alpha2 x2 = beta1 + beta2 z1p = alpha1/x1 z1m = alpha2/x1 z2p = beta1/x2 z2m = beta2/x2 c ----------------------------------------------------------------- if(x1.gt.1.0.or.x2.gt.1.0) return c ----------------------------------------------------------------- s1_eff = x1*s - q1t**2 s2_eff = x2*s - q2t**2 c ----------------------------------------------------------------- c Additional conditions for energy-momentum conservation c ----------------------------------------------------------------- if(((icontri.eq.2).or.(icontri.eq.4)) 1 .and.(dsqrt(s1_eff).le.(am_y+invm))) return if(((icontri.eq.3).or.(icontri.eq.4)) 1 .and.(dsqrt(s2_eff).le.(am_x+invm))) return c ================================================================= c >>> TO THE OUTPUT COMMON BLOCK c ================================================================= c ================================================================= c four-momenta of the outgoing protons (or remnants) c ================================================================= px_plus = (1.d0-x1) * p1_plus px_minus = ((a_nuc1*am_x)**2 + q1tx**2 + q1ty**2)/2.d0/px_plus px(1) = - q1tx px(2) = - q1ty px(3) = (px_plus - px_minus)/dsqrt(2.d0) px(4) = (px_plus + px_minus)/dsqrt(2.d0) py_minus = (1.d0-x2) * p2_minus py_plus = ((a_nuc2*am_y)**2 + q2tx**2 + q2ty**2)/2.d0/py_minus py(1) = - q2tx py(2) = - q2ty py(3) = (py_plus - py_minus)/dsqrt(2.d0) py(4) = (py_plus + py_minus)/dsqrt(2.d0) q1t = dsqrt(q1t2) q2t = dsqrt(q2t2) c ================================================================= c four-momenta of the outgoing central particles c ================================================================= nout = 2 ipdg(1) = pdg_l pc(1,1) = pt1x pc(1,2) = pt1y pc(1,3) = alpha1*ak1z + beta1*ak2z pc(1,4) = alpha1*ak10 + beta1*ak20 ipdg(2) = -pdg_l pc(2,1) = pt2x pc(2,2) = pt2y pc(2,3) = alpha2*ak1z + beta2*ak2z pc(2,4) = alpha2*ak10 + beta2*ak20 eta1 = 0.5d0*dlog((dsqrt(amt1**2*(dcosh(y1))**2 - am_l**2) + 2 amt1*dsinh(y1))/(dsqrt(amt1**2*(dcosh(y1))**2 - am_l**2) 3 - amt1*dsinh(y1))) eta2 = 0.5d0*dlog((dsqrt(amt2**2*(dcosh(y2))**2 - am_l**2) + 2 amt2*dsinh(y2))/(dsqrt(amt2**2*(dcosh(y2))**2 - am_l**2) 3 - amt2*dsinh(y2))) if(ieta) then if(eta1.lt.eta_min.or.eta1.gt.eta_max) return if(eta2.lt.eta_min.or.eta2.gt.eta_max) return endif c ================================================================= c matrix element squared c averaged over initial spin polarizations c and summed over final spin polarizations c (--> see Wolfgang's notes c ================================================================= c ================================================================= c Mendelstam variables c ================================================================= that1 = (q10-pc(1,4))**2 & -(q1tx-pc(1,1))**2-(q1ty-pc(2,1))**2-(q1z-pc(3,1))**2 uhat1 = (q10-pc(2,4))**2 & -(q1tx-pc(1,2))**2-(q1ty-pc(2,2))**2-(q1z-pc(3,2))**2 that2 = (q20-pc(2,4))**2 & -(q2tx-pc(1,2))**2-(q2ty-pc(2,2))**2-(q2z-pc(3,2))**2 uhat2 = (q20-pc(1,4))**2 & -(q2tx-pc(1,1))**2-(q2ty-pc(2,1))**2-(q2z-pc(3,1))**2 that = (that1+that2)/2.d0 uhat = (uhat1+uhat2)/2.d0 c ================================================================= c matrix elements c ================================================================= c How matrix element is calculated c imethod = 0: on-shell formula c imethod = 1: off-shell formula c ================================================================= if(imethod.eq.0) then c ================================================================= c on-shell formula for M^2 c ================================================================= term1 = 6.d0*am_l**8 term2 = -3.d0*am_l**4*that**2 term3 = -14.d0*am_l**4*that*uhat term4 = -3.d0*am_l**4*uhat**2 term5 = am_l**2*that**3 term6 = 7.d0*am_l**2*that**2*uhat term7 = 7.d0*am_l**2*that*uhat**2 term8 = am_l**2*uhat**3 term9 = -that**3*uhat term10 = -that*uhat**3 amat2 = -2.d0*( term1+term2+term3+term4+term5 2 +term6+term7+term8+term9+term10 ) 3 / ( (am_l**2-that)**2 * (am_l**2-uhat)**2 ) elseif(imethod.eq.1)then c ================================================================= c Wolfgang's formulae c ================================================================= ak1_x = z1m*pt1x-z1p*pt2x ak1_y = z1m*pt1y-z1p*pt2y ak2_x = z2m*pt1x-z2p*pt2x ak2_y = z2m*pt1y-z2p*pt2y !FIXME FIXME FIXME am_p or am_A/B??? t1abs = (q1t2+x1*(am_x**2-am_p**2)+x1**2*am_p**2)/(1.d0-x1) t2abs = (q2t2+x2*(am_y**2-am_p**2)+x2**2*am_p**2)/(1.d0-x2) eps12 = am_l**2 + z1p*z1m*t1abs eps22 = am_l**2 + z2p*z2m*t2abs Phi10 = 1.d0/((ak1_x+z1p*q2tx)**2+(ak1_y+z1p*q2ty)**2+eps12) 2 - 1.d0/((ak1_x-z1m*q2tx)**2+(ak1_y-z1m*q2ty)**2+eps12) Phi11_x = (ak1_x+z1p*q2tx)/ 2 ((ak1_x+z1p*q2tx)**2+(ak1_y+z1p*q2ty)**2+eps12) 3 - (ak1_x-z1m*q2tx)/ 4 ((ak1_x-z1m*q2tx)**2+(ak1_y-z1m*q2ty)**2+eps12) Phi11_y = (ak1_y+z1p*q2ty)/ 2 ((ak1_x+z1p*q2tx)**2+(ak1_y+z1p*q2ty)**2+eps12) 3 - (ak1_y-z1m*q2ty)/ 4 ((ak1_x-z1m*q2tx)**2+(ak1_y-z1m*q2ty)**2+eps12) Phi102 = Phi10*Phi10 Phi112 = Phi11_x**2+Phi11_y**2 Phi20 = 1.d0/((ak2_x+z2p*q1tx)**2+(ak2_y+z2p*q1ty)**2+eps22) 2 - 1.d0/((ak2_x-z2m*q1tx)**2+(ak2_y-z2m*q1ty)**2+eps22) Phi21_x = (ak2_x+z2p*q1tx)/ 2 ((ak2_x+z2p*q1tx)**2+(ak2_y+z2p*q1ty)**2+eps22) 3 - (ak2_x-z2m*q1tx)/ 4 ((ak2_x-z2m*q1tx)**2+(ak2_y-z2m*q1ty)**2+eps22) Phi21_y = (ak2_y+z2p*q1ty)/ 2 ((ak2_x+z2p*q1tx)**2+(ak2_y+z2p*q1ty)**2+eps22) 3 - (ak2_y-z2m*q1ty)/ 4 ((ak2_x-z2m*q1tx)**2+(ak2_y-z2m*q1ty)**2+eps22) Phi202 = Phi20*Phi20 Phi212 = Phi21_x**2+Phi21_y**2 Phi11_dot_e = (Phi11_x*q1tx + Phi11_y*q1ty)/dsqrt(q1t2) Phi11_cross_e = (Phi11_x*q1ty - Phi11_y*q1tx)/dsqrt(q1t2) Phi21_dot_e = (Phi21_x*q2tx +Phi21_y*q2ty)/dsqrt(q2t2) Phi21_cross_e = (Phi21_x*q2ty -Phi21_y*q2tx)/dsqrt(q2t2) aux2_1 = iterm11*(am_l**2+4.d0*z1p**2*z1m**2*t1abs)*Phi102 1 +iterm22*( (z1p**2 + z1m**2)*(Phi11_dot_e**2 + 2 Phi11_cross_e**2) ) 3 + itermtt*( Phi11_cross_e**2 - Phi11_dot_e**2) 4 - iterm12*4.d0*z1p*z1m*(z1p-z1m)*Phi10 5 *(q1tx*Phi11_x+q1ty*Phi11_y) aux2_2 = iterm11*(am_l**2+4.d0*z2p**2*z2m**2*t2abs)*Phi202 1 +iterm22*( (z2p**2 + z2m**2)*(Phi21_dot_e**2 + 2 Phi21_cross_e**2) ) 3 + itermtt*( Phi21_cross_e**2 - Phi21_dot_e**2) 4 - iterm12*4.d0*z2p*z2m*(z2p-z2m)*Phi20 5 *(q2tx*Phi21_x+q2ty*Phi21_y) c ================================================================= c convention of matrix element as in our kt-factorization c for heavy flavours c ================================================================= amat2_1 = (x1*x2*s)**2 * aux2_1 * 2.*z1p*z1m*q1t2 / (q1t2*q2t2) amat2_2 = (x1*x2*s)**2 * aux2_2 * 2.*z2p*z2m*q2t2 / (q1t2*q2t2) c ================================================================= c symmetrization c ================================================================= amat2 = (imat1*amat2_1 + imat2*amat2_2)/2.d0 endif coupling = 1.d0 c ================================================================= c first parton coupling c ================================================================= if(iflux1.ge.20.and.iflux1.lt.40) then ! at least one gluon exchanged #ifdef ALPHA_S t_max = max(amt1**2,amt2**2) amu2 = max(eps12,t_max) am_x = dsqrt(amu2) coupling = coupling * 4.d0*pi*alphaS(am_x)/2.d0 ! colour flow #else print *,'alphaS not linked to this instance!' stop #endif else ! photon exchange coupling = coupling * 4.d0*pi*alpha_em*q_l**2 endif c ================================================================= c second parton coupling c ================================================================= coupling = coupling * 4.d0*pi*alpha_em*q_l**2 ! photon exchange coupling = coupling * 3.d0 c ============================================ c unintegrated parton distributions c ============================================ if(a_nuc1.le.1) then f1 = CepGen_kT_flux(iflux1,x1,q1t2,sfmod,am_x) else f1 = CepGen_kT_flux_HI(iflux1,x1,q1t2,a_nuc1,z_nuc1) endif if(a_nuc2.le.1) then f2 = CepGen_kT_flux(iflux2,x2,q2t2,sfmod,am_y) else f2 = CepGen_kT_flux_HI(iflux2,x2,q2t2,a_nuc2,z_nuc2) endif c ================================================================= c factor 2.*pi below from integration over phi_sum c factor 1/4 below from jacobian of transformations c factors 1/pi and 1/pi due to integration c over d^2 kappa_1 d^2 kappa_2 instead d kappa_1^2 d kappa_2^2 c ================================================================= aintegral = (2.d0*pi)*1.d0/(16.d0*pi**2*(x1*x2*s)**2) * amat2 & * coupling & * f1/pi * f2/pi * (1.d0/4.d0) * units & * 0.5d0 * 4.0d0 / (4.d0*pi) c ***************************************************************** c ================================================================= nucl_to_ff = aintegral*q1t*q2t*ptdiff c ================================================================= c ***************************************************************** c print *,nucl_to_ff,aintegral,coupling return end diff --git a/cmake/UseEnvironment.cmake b/cmake/UseEnvironment.cmake index 6114b63..f1f94cb 100644 --- a/cmake/UseEnvironment.cmake +++ b/cmake/UseEnvironment.cmake @@ -1,84 +1,94 @@ if(NOT CMAKE_VERSION VERSION_LESS 3.1) set(CMAKE_CXX_STANDARD 14) else() set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++14") endif() set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -Wall -cpp") #--- check if we are at CERN if($ENV{HOSTNAME} MATCHES "^lxplus[0-9]+.cern.ch") set(IS_LXPLUS "yes") endif() #--- ensure a proper version of the compiler is found if(CMAKE_CXX_COMPILER_ID MATCHES "Clang" OR CMAKE_CXX_COMPILER_VERSION VERSION_GREATER 6.1) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wno-long-long -pedantic-errors -g") else() message(STATUS "clang or gcc above 6.1 is required") if(IS_LXPLUS) set(LXPLUS_SRC_ENV "source ${CMAKE_SOURCE_DIR}/source-lxplus.sh") message(STATUS "Compiling on LXPLUS. Did you properly source the environment variables? E.g.\n\n\t${LXPLUS_SRC_ENV}\n") endif() message(FATAL_ERROR "Please clean up this build environment, i.e.\n\trm -rf CMake*\nand try again...") endif() #--- set the default paths for external dependencies if(IS_LXPLUS) set(BASE_DIR "/cvmfs/sft.cern.ch/lcg") list(APPEND CMAKE_PREFIX_PATH "${BASE_DIR}/external/CMake/2.8.9/Linux-i386/share/cmake-2.8/Modules") set(GSL_DIR "${BASE_DIR}/releases/GSL/2.5-32fc5/x86_64-centos7-gcc62-opt") set(HEPMC_DIR "${BASE_DIR}/releases/HepMC/2.06.09-0a23a/x86_64-centos7-gcc62-opt") set(LHAPDF_DIR "${BASE_DIR}/releases/MCGenerators/lhapdf/6.2.2-8a3e6/x86_64-centos7-gcc62-opt") set(PYTHIA6_DIR "${BASE_DIR}/releases/MCGenerators/pythia6/429.2-c4089/x86_64-centos7-gcc62-opt") set(PYTHIA8_DIR "${BASE_DIR}/releases/MCGenerators/pythia8/240p1-ecd34/x86_64-centos7-gcc62-opt") + #--- disabled for the time being (no compatible version found on CVMFS) + #set(DELPHES_DIR "${BASE_DIR}/releases/delphes/3.4.0-03b2c/x86_64-centos7-gcc62-opt") + #set(TBB_DIR "${BASE_DIR}/releases/tbb/2019_U7-ba3eb/x86_64-centos7-gcc62-opt") set(PYTHON_DIR "${BASE_DIR}/releases/Python/2.7.15-075d4/x86_64-centos7-gcc62-opt") set(PYTHON_LIBRARY "${PYTHON_DIR}/lib/libpython2.7.so") set(PYTHON_EXECUTABLE "${PYTHON_DIR}/bin/python") set(PYTHON_INCLUDE_DIR "${PYTHON_DIR}/include/python2.7") message(STATUS "Compiling on LXPLUS. Do not forget to source the environment variables!") message(STATUS "e.g. `${LXPLUS_SRC_ENV}`") endif() #--- searching for GSL find_library(GSL_LIB gsl HINTS ${GSL_DIR} PATH_SUFFIXES lib REQUIRED) find_library(GSL_CBLAS_LIB gslcblas HINTS ${GSL_DIR} PATH_SUFFIXES lib) find_path(GSL_INCLUDE gsl HINTS ${GSL_DIR} PATH_SUFFIXES include) include_directories(${GSL_INCLUDE}) #--- searching for LHAPDF find_library(LHAPDF LHAPDF HINTS ${LHAPDF_DIR} PATH_SUFFIXES lib) find_path(LHAPDF_INCLUDE LHAPDF HINTS ${LHAPDF_DIR} PATH_SUFFIXES include) #--- searching for HepMC -find_library(HEPMC_LIB NAMES HepMC3 HepMC HINTS ${HEPMC_DIR} PATH_SUFFIXES lib) -find_library(HEPMC_ROOT_LIB NAMES HepMC3rootIO PATH_SUFFIXES root) -find_path(HEPMC_INCLUDE NAMES HepMC3 HepMC HINTS ${HEPMC_DIR} PATH_SUFFIXES include) +find_library(HEPMC_LIB NAMES HepMC3 HepMC HINTS $ENV{HEPMC_DIR} ${HEPMC_DIR} PATH_SUFFIXES lib64 lib) +find_library(HEPMC_ROOT_LIB NAMES HepMC3rootIO HepMCrootIO PATH_SUFFIXES root) +find_path(HEPMC_INCLUDE NAMES HepMC3 HepMC HINTS $ENV{HEPMC_DIR} ${HEPMC_DIR} PATH_SUFFIXES include) #--- searching for Pythia 6 set(PYTHIA6_DIRS $ENV{PYTHIA6_DIR} ${PYTHIA6_DIR} /usr /usr/local /opt/pythia6) find_library(PYTHIA6 pythia6 HINTS ${PYTHIA6_DIRS} PATH_SUFFIXES lib) find_library(PYTHIA6DUMMY pythia6_dummy HINTS ${PYTHIA6_DIRS} PATH_SUFFIXES lib) #--- searching for Pythia 8 set(PYTHIA8_DIRS $ENV{PYTHIA8_DIR} ${PYTHIA8_DIR} /usr /usr/local /opt/pythia8) find_library(PYTHIA8 pythia8 HINTS ${PYTHIA8_DIRS} PATH_SUFFIXES lib) find_path(PYTHIA8_INCLUDE Pythia8 HINTS ${PYTHIA8_DIRS} PATH_SUFFIXES include include/Pythia8 include/pythia8) +#--- searching for Delphes +find_library(DELPHES Delphes HINTS $ENV{DELPHES_DIR} ${DELPHES_DIR} PATH_SUFFIXES lib) +find_path(DELPHES_INCLUDE NAMES modules classes HINTS $ENV{DELPHES_DIR} ${DELPHES_DIR} PATH_SUFFIXES include) +#--- searching for tbb +find_library(TBB tbb HINTS ${TBB_DIR} PATH_SUFFIXES lib) +#--- searching for ROOT +find_package(ROOT QUIET) message(STATUS "GSL found in ${GSL_LIB}") list(APPEND CEPGEN_EXTERNAL_CORE_REQS ${GSL_LIB} ${GSL_CBLAS_LIB}) find_package(PythonLibs 2.7) find_library(MUPARSER muparser) if(MUPARSER) message(STATUS "muParser found in ${MUPARSER}") list(APPEND CEPGEN_EXTERNAL_CORE_REQS ${MUPARSER}) add_definitions(-DMUPARSER) else() find_path(EXPRTK exprtk.hpp PATH_SUFFIXES include) if(EXPRTK) message(STATUS "exprtk found in ${EXPRTK}") add_definitions(-DEXPRTK) include_directories(${EXPRTK}) endif() endif() #--- semi-external dependencies set(ALPHAS_PATH ${PROJECT_SOURCE_DIR}/External) file(GLOB alphas_sources ${ALPHAS_PATH}/alphaS.f) if(alphas_sources) message(STATUS "alphaS evolution found in ${alphas_sources}") add_definitions(-DALPHA_S) endif() diff --git a/release-notes.md b/release-notes.md index 81fd31b..e15b8b8 100644 --- a/release-notes.md +++ b/release-notes.md @@ -1,45 +1,54 @@ # Release notes +## v0.9.7 (25 Jul 2019) +* Fortran processes can now be fed a generic set of parameters, thanks to additional getter functions +* Output handlers may now be constructed directly from steering cards, thus enhancing overall modularity. +* Added a helper for the retrieval of events properties through human-readable getters +* New output handlers: text (raw text output, and ASCII histograms), HepMC ASCII output (for HepMC v<3), + ROOT histogram collections and ntuple files writers +* In addition, added an interface to Delphes for the simulation of detectors effects +* Refactored HepMC event builder in preparation for future developments + ## v0.9.6 (11 Jul 2019) * Added support of Pythia6 hadronisation/fragmentation algorithm for legacy tests * Structure functions parameterisation objects polished * New output modes handled for HepMC interfacing module ## v0.9.5 (25 Jun 2019) * Increased flexibility in particles (and their associated properties) definitions * Small corrections in the LPAIR process definition of the output kinematics * Improved Pythia8 interfacing, better handling of the LHEF output format when remnants dissociation is triggered * Better exceptions handling ## v0.9.4 (9 May 2019) * pptoff polarisation terms may be steered for off-shell ME * Functional is now handling exprtk in addition to μParser * Refactored exceptions and tests * Improved accessors for integrator/grid definition * Better memory management for parameters/processes/hadronisers definition ## v0.9.1-3 (29-30 Nov 2018) * Fixes in includes chain * Fix in LPAIR cards parser, simplified cards handler interface * External Fortran processes now defined through a function instead of a subroutine * Fix in Pythia8 interface parton colours definition and resonances decay * Structure functions interface modified ## v0.9 (28-29 Nov 2018) * Fixed memory leaks in Python configuration files handler * First HI process! (and first Fortran process interfaced) * Polarisation states may be specified for the off-shell γγ → W+W- ME * New structure functions definitions: ALLM (+subsequent), MSTW grid, CLAS, LUXlike, LHAPDF * Introduced R-ratio definitions for FL calculation from F2 * New cuts definition * Better memory management ## v0.8 (24 Aug 2017) * Major refactoring of the code * Dropped the hadronisers and other unfinished processes definition for clarity * Several bug fixes and memory management improvements through the usage of smart pointers / C++11 migration * Taming functions introduced ## v0.7 (3 May 2016) * Introduced a generic templating class for the kT factorisation method * Added a placeholder for the γγ → W+W- process computation diff --git a/source-lxplus.sh b/source-lxplus.sh index 5da8659..8967225 100755 --- a/source-lxplus.sh +++ b/source-lxplus.sh @@ -1,8 +1,9 @@ #!/bin/sh source /cvmfs/sft.cern.ch/lcg/releases/gcc/6.2.0/x86_64-centos7/setup.sh -source /cvmfs/sft.cern.ch/lcg/releases/ROOT/6.16.00-4d56b/x86_64-centos7-gcc62-opt/bin/thisroot.sh +source /cvmfs/sft.cern.ch/lcg/releases/ROOT/6.10.06-a3425/x86_64-centos7-gcc62-opt/bin/thisroot.sh #--- for Delphes +#source /cvmfs/sft.cern.ch/lcg/releases/ROOT/6.16.00-4d56b/x86_64-centos7-gcc62-opt/bin/thisroot.sh export PYTHONHOME=/cvmfs/sft.cern.ch/lcg/releases/Python/2.7.15-075d4/x86_64-centos7-gcc62-opt export PYTHIA8_DIR=/cvmfs/sft.cern.ch/lcg/releases/MCGenerators/pythia8/240p1-ecd34/x86_64-centos7-gcc62-opt export PYTHIA8DATA=${PYTHIA8_DIR}/share/Pythia8/xmldoc echo "Environment prepared for LXPLUS" diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 867ae34..16b6510 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,53 +1,52 @@ #----- list all test executables file(GLOB executables_noroot RELATIVE ${PROJECT_SOURCE_DIR}/test *.cpp) file(GLOB executables_root RELATIVE ${PROJECT_SOURCE_DIR}/test *.cxx) file(GLOB executables_fortran RELATIVE ${PROJECT_SOURCE_DIR}/test *.f) #----- include the utilitaries set(tests "") #----- build all tests and link them to the core library foreach(exec_src ${executables_noroot}) string(REPLACE ".cpp" "" exec_bin ${exec_src}) add_executable(${exec_bin} ${exec_src}) set_target_properties(${exec_bin} PROPERTIES EXCLUDE_FROM_ALL true) target_link_libraries(${exec_bin} ${CEPGEN_LIBRARIES}) list(APPEND tests ${exec_bin}) endforeach() foreach(exec_src ${executables_fortran}) string(REPLACE ".f" "" exec_bin ${exec_src}) add_executable(${exec_bin} ${exec_src}) set_target_properties(${exec_bin} PROPERTIES EXCLUDE_FROM_ALL true) target_link_libraries(${exec_bin} ${CEPGEN_LIBRARIES}) list(APPEND tests ${exec_bin}) endforeach() #----- specify the tests requiring ROOT -find_package(ROOT QUIET) if(ROOT_FOUND) include_directories(${ROOT_INCLUDE_DIRS}) link_directories(${ROOT_LIBRARY_DIR}) foreach(exec_src ${executables_root}) string(REPLACE ".cxx" "" exec_bin ${exec_src}) add_executable(${exec_bin} ${exec_src}) set_target_properties(${exec_bin} PROPERTIES EXCLUDE_FROM_ALL true) target_link_libraries(${exec_bin} ${CEPGEN_LIBRARIES} ${ROOT_LIBRARIES}) list(APPEND tests ${exec_bin}) endforeach() else() message(STATUS "ROOT not found! skipping these tests!") endif() message(STATUS "The following tests/executables are available:\n${tests}") file(GLOB_RECURSE test_input_cards RELATIVE ${PROJECT_SOURCE_DIR}/test test_processes/* __init__.py) foreach(_files ${test_input_cards}) configure_file(${_files} ${_files} COPYONLY) endforeach() configure_file(${PROJECT_SOURCE_DIR}/test/test_processes.cfg test_processes.cfg COPYONLY) diff --git a/test/cepgen-event.cpp b/test/cepgen-event.cpp index 7f61022..3485918 100644 --- a/test/cepgen-event.cpp +++ b/test/cepgen-event.cpp @@ -1,82 +1,82 @@ #include "CepGen/Cards/Handler.h" #include "CepGen/IO/ExportHandler.h" #include "CepGen/Generator.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/ParametersList.h" #include "CepGen/Event/Event.h" #include "AbortHandler.h" #include using namespace std; // we use polymorphism here -std::unique_ptr writer; +std::unique_ptr writer; void storeEvent( const cepgen::Event& ev, unsigned long ) { if ( !writer ) throw CG_FATAL( "storeEvent" ) << "Failed to retrieve a valid writer!"; *writer << ev; } /** * Main caller for this Monte Carlo generator. Loads the configuration files' * variables if set as an argument to this program, else loads a default * "LHC-like" configuration, then launches the cross-section computation and * the events generation. * \author Laurent Forthomme */ int main( int argc, char* argv[] ) { ostringstream os; string delim; - for ( const auto& mod : cepgen::output::ExportHandler::get().modules() ) + for ( const auto& mod : cepgen::io::ExportHandler::get().modules() ) os << delim << mod, delim = ","; if ( argc < 2 ) throw CG_FATAL( "main" ) << "No config file provided!\n\t" << "Usage: " << argv[0] << " config-file [format=" << os.str() << "] [filename=example.dat]"; cepgen::Generator mg; //----------------------------------------------------------------------------------------------- // Steering card readout //----------------------------------------------------------------------------------------------- CG_DEBUG( "main" ) << "Reading config file stored in \"" << argv[1] << "\""; mg.setParameters( cepgen::card::Handler::parse( argv[1] ) ); //----------------------------------------------------------------------------------------------- // Output file writer definition //----------------------------------------------------------------------------------------------- const string format = ( argc > 2 ) ? argv[2] : "lhef"; const char* filename = ( argc > 3 ) ? argv[3] : "example.dat"; - writer = cepgen::output::ExportHandler::get().build( format, cepgen::ParametersList() + writer = cepgen::io::ExportHandler::get().build( format, cepgen::ParametersList() .set( "filename", filename ) ); //----------------------------------------------------------------------------------------------- // CepGen run part //----------------------------------------------------------------------------------------------- // We might want to cross-check visually the validity of our run CG_INFO( "main" ) << mg.parametersPtr(); cepgen::utils::AbortHandler ctrl_c; //--- let there be cross-section... double xsec = 0., err = 0.; mg.computeXsection( xsec, err ); writer->initialise( mg.parameters() ); writer->setCrossSection( xsec, err ); // The events generation starts here! mg.generate( storeEvent ); return 0; } diff --git a/test/cepgen-root.cxx b/test/cepgen-root.cxx deleted file mode 100644 index d40ced3..0000000 --- a/test/cepgen-root.cxx +++ /dev/null @@ -1,116 +0,0 @@ -#include "CepGen/Cards/Handler.h" -#include "CepGen/Generator.h" -#include "CepGen/Event/Event.h" - -#include -#include - -#include "TreeInfo.h" -#include "AbortHandler.h" - -// ROOT includes -#include "TFile.h" -#include "TTree.h" -#include "TLorentzVector.h" - -using namespace std; - -std::unique_ptr run; -std::unique_ptr ev; - -void fill_event_tree( const cepgen::Event& event, unsigned long ev_id ) -{ - //if ( ev_id % 10 == 0 ) - // cout << ">> event " << ev_id << " generated" << endl; - - if ( !ev || !run ) - return; - - ev->gen_time = event.time_generation; - ev->tot_time = event.time_total; - ev->np = 0; - //cout << event.particles().size() << endl; - for ( const auto& p : event.particles() ) { - const cepgen::Particle::Momentum m = p.momentum(); - ev->rapidity[ev->np] = m.rapidity(); - ev->pt[ev->np] = m.pt(); - ev->eta[ev->np] = m.eta(); - ev->phi[ev->np] = m.phi(); - ev->E[ev->np] = p.energy(); - ev->m[ev->np] = p.mass(); - ev->pdg_id[ev->np] = p.integerPdgId(); - ev->parent1[ev->np] = ( p.mothers().size() > 0 ) ? *p.mothers().begin() : -1; - ev->parent2[ev->np] = ( p.mothers().size() > 1 ) ? *p.mothers().rbegin() : -1; - ev->status[ev->np] = (int)p.status(); - ev->stable[ev->np] = ( (short)p.status() > 0 ); - ev->charge[ev->np] = p.charge(); - ev->role[ev->np] = p.role(); - - ev->np++; - } - run->num_events += 1; - ev->fill(); -} - -/** - * Generation of events and storage in a ROOT format - * \author Laurent Forthomme - * \date 27 jan 2014 - * \defgroup Executables List of executables - * \addtogroup Executables - */ -int main( int argc, char* argv[] ) -{ - cepgen::utils::AbortHandler ctrl_c; - - cepgen::Generator mg; - - if ( argc < 2 ) - throw CG_FATAL( "main" ) << "Usage: " << argv[0] << " input-card [filename=events.root]"; - - mg.setParameters( cepgen::card::Handler::parse( argv[1] ) ); - CG_INFO( "main" ) << mg.parametersPtr(); - - //----- open the output root file - - const char* filename = ( argc > 2 ) ? argv[2] : "events.root"; - TFile file( filename, "recreate" ); - if ( !file.IsOpen() ) - throw CG_FATAL( "main" ) << "Failed to create the output file!"; - - //----- then generate the events and the container tree structure - - run.reset( new ROOT::CepGenRun ); - run->create(); - ev.reset( new ROOT::CepGenEvent ); - ev->create(); - - try { - //----- start by computing the cross section for the list of parameters applied - double xsec, err; - mg.computeXsection( xsec, err ); - - //----- populate the run tree - - run->xsect = xsec; - run->errxsect = err; - run->litigious_events = 0; - run->sqrt_s = mg.parameters().kinematics.sqrtS(); - - //----- launch the events generation - mg.generate( fill_event_tree ); - run->fill(); - } catch ( const cepgen::utils::RunAbortedException& e ) { - e.dump(); - } catch ( const cepgen::Exception& e ) { - e.dump(); - } - - file.Write(); - CG_INFO( "main" ) - << run->num_events << " event" << cepgen::utils::s( run->num_events ) - << " written in \"" << filename << "\"."; - - return 0; -} - diff --git a/test/cepgen.cpp b/test/cepgen.cpp index becdc19..31b562d 100644 --- a/test/cepgen.cpp +++ b/test/cepgen.cpp @@ -1,49 +1,52 @@ //--- steering cards #include "CepGen/Cards/Handler.h" #include "CepGen/Generator.h" #include "CepGen/Core/Exception.h" #include "AbortHandler.h" using namespace std; /** Example executable for CepGen * - loads the steering card variables into the environment, * - launches the cross-section computation and the events generation (if requested). * \author Laurent Forthomme * \addtogroup Executables */ int main( int argc, char* argv[] ) { if ( argc < 2 ) throw CG_FATAL( "main" ) << "No config file provided!\n\t" << "Usage: " << argv[0] << " config-file"; //--- first start by defining the generator object cepgen::Generator gen; gen.setParameters( cepgen::card::Handler::parse( argv[1] ) ); //--- list all parameters CG_LOG( "main" ) << gen.parametersPtr(); cepgen::utils::AbortHandler ctrl_c; try { //--- let there be a cross-section... double xsec = 0., err = 0.; gen.computeXsection( xsec, err ); if ( gen.parameters().generation().enabled ) //--- events generation starts here // (one may use a callback function) gen.generate(); } catch ( const cepgen::utils::RunAbortedException& e ) { CG_DEBUG( "main" ) << "Run aborted!"; } catch ( const cepgen::Exception& e ) { e.dump(); + } catch ( ... ) { + CG_FATAL( "main" ) << "Other exception caught!"; + throw; } return 0; } diff --git a/test/test_event_writer.cpp b/test/test_event_writer.cpp index 1f2e7e4..bc02f37 100644 --- a/test/test_event_writer.cpp +++ b/test/test_event_writer.cpp @@ -1,30 +1,30 @@ #include "CepGen/IO/ExportHandler.h" #include "CepGen/Physics/PDG.h" #include "CepGen/Event/Event.h" using namespace std; using namespace cepgen; int main() { - auto writer = output::ExportHandler::get().build( "hepmc" ); + auto writer = io::ExportHandler::get().build( "hepmc" ); writer->setCrossSection( 1., 2. ); Event ev; Particle p1( Particle::IncomingBeam1, PDG::proton ); p1.setMomentum( 1., -15., 100. ); p1.setStatus( Particle::Status::Incoming ); ev.addParticle(p1); Particle p2( Particle::IncomingBeam2, PDG::electron ); p2.setMomentum( 10., 5., 3200. ); p2.setStatus( Particle::Status::Incoming ); ev.addParticle(p2); ev.dump(); *writer << ev; return 0; }