diff --git a/CepGen/CMakeLists.txt b/CepGen/CMakeLists.txt index 68a1cee..1ca3189 100644 --- a/CepGen/CMakeLists.txt +++ b/CepGen/CMakeLists.txt @@ -1,127 +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 IO/TextHandler.cpp) file(GLOB hadr_sources Hadronisers/GenericHadroniser.cpp) include_directories(${PROJECT_SOURCE_DIR}) #----- check the external dependencies for SFs set(GRV_PATH ${PROJECT_SOURCE_DIR}/External) file(GLOB grv_sources ${GRV_PATH}/grv_*.f) if(grv_sources) message(STATUS "GRV PDFset found in ${grv_sources}!") add_definitions(-DGRVPDF) list(APPEND sf_sources ${grv_sources}) else() message(STATUS "GRV PDFset not found. Will proceed without it") endif() #----- check the external dependencies for physics utilities if(alphas_sources) list(APPEND phys_sources ${alphas_sources}) endif() set(addons_libraries "") set(strf_libraries ${GSL_LIB} ${GSL_CBLAS_LIB}) #--- linking with LHAPDF if(LHAPDF) message(STATUS "LHAPDF found in ${LHAPDF}") list(APPEND strf_libraries ${LHAPDF}) include_directories(${LHAPDF_INCLUDE}) else() file(GLOB partonic_sf StructureFunctions/Partonic.cpp) list(REMOVE_ITEM sf_sources ${partonic_sf}) endif() #--- linking with Pythia 6 if(PYTHIA6) message(STATUS "Pythia 6 found in ${PYTHIA6}") list(APPEND hadr_sources "Hadronisers/Pythia6Hadroniser.cpp") list(APPEND addons_libraries ${PYTHIA6}) if(PYTHIA6DUMMY) message(STATUS "Pythia 6 addons found in ${PYTHIA6DUMMY}") list(APPEND addons_libraries ${PYTHIA6DUMMY}) endif() elseif(EXISTS $ENV{PYTHIA6_SRC}) file(GLOB pythia6_src $ENV{PYTHIA6_SRC}) message(STATUS "Pythia 6 source found in ${pythia6_src}") set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -Wno-tabs -Wno-maybe-uninitialized -Wno-integer-division -Wno-unused-variable -Wno-unused-dummy-argument") add_library(pythia6 SHARED ${pythia6_src}) list(APPEND hadr_sources "Hadronisers/Pythia6Hadroniser.cpp") list(APPEND addons_libraries pythia6) endif() #--- linking with Pythia 8 if(PYTHIA8) message(STATUS "Pythia 8 found in ${PYTHIA8}") message(STATUS "Pythia 8 will be used for LHEF output") set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-misleading-indentation") include_directories(${PYTHIA8_INCLUDE}) add_definitions(-DPYTHIA8) list(APPEND hadr_sources "Hadronisers/Pythia8Hadroniser.cpp") list(APPEND io_sources "IO/PythiaEventInterface.cpp") list(APPEND io_sources "IO/LHEFPythiaHandler.cpp") list(APPEND addons_libraries ${PYTHIA8} dl) endif() #--- linking with HepMC if(HEPMC_LIB) message(STATUS "HepMC found in ${HEPMC_LIB}") if(HEPMC_LIB MATCHES ".*HepMC3.?.so") message(STATUS "HepMC version 3 found") if(HEPMC_ROOT_LIB) message(STATUS "HepMC ROOT I/O library found") add_definitions(-DHEPMC3_ROOTIO) list(APPEND addons_libraries ${HEPMC_ROOT_LIB}) endif() add_definitions(-DHEPMC3) endif() list(APPEND io_sources "IO/HepMCHandler.cpp") if(NOT PYTHIA8) message(STATUS "HepMC will be used for LHEF output") list(APPEND io_sources "IO/LHEFHepMCHandler.cpp") endif() list(APPEND addons_libraries ${HEPMC_LIB}) include_directories(${HEPMC_INCLUDE}) endif() -#--- linking with Delphes - -if(DELPHES) - message(STATUS "Delphes found in ${DELPHES}") - find_package(ROOT QUIET) - if(ROOT_FOUND) - message(STATUS "ROOT found in ${ROOT_INCLUDE_DIRS}") - list(APPEND addons_libraries ${DELPHES} ${ROOT_LIBRARIES} ${TBB}) +#--- 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/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} ${TBB}) list(APPEND io_sources "IO/DelphesHandler.cpp") include_directories(${DELPHES_INCLUDE}) message(STATUS ${DELPHES_INCLUDE}) - include_directories(${ROOT_INCLUDE_DIRS}) - else() - message(STATUS "ROOT not found! Disabling Delphes module.") 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 d9378ab..ffb2c37 100644 --- a/CepGen/Cards/LpairHandler.cpp +++ b/CepGen/Cards/LpairHandler.cpp @@ -1,249 +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( "OUTP", "Output module", &out_mod_name_ ); registerParameter( "OUTF", "Output file name", &out_file_name_ ); - registerParameter( "KMRG", "KMR grid interpolation path", &kmr_grid_path_ ); //------------------------------------------------------------------------------------------- // 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 fef13af..22503f1 100644 --- a/CepGen/Cards/LpairHandler.h +++ b/CepGen/Cards/LpairHandler.h @@ -1,121 +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_, 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/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/Integrator.cpp b/CepGen/Core/Integrator.cpp index f267ae2..ce7a162 100644 --- a/CepGen/Core/Integrator.cpp +++ b/CepGen/Core/Integrator.cpp @@ -1,475 +1,477 @@ #include "CepGen/Core/Integrator.h" #include "CepGen/Core/GridParameters.h" #include "CepGen/Core/utils.h" #include "CepGen/Core/Exception.h" #include "CepGen/Parameters.h" #include "CepGen/Processes/GenericProcess.h" #include "CepGen/Hadronisers/GenericHadroniser.h" #include "CepGen/IO/GenericExportHandler.h" #include "CepGen/Event/Event.h" #include #include #include #define COORD(s,i,j) ((s)->xi[(i)*(s)->dim + (j)]) namespace cepgen { Integrator::Integrator( unsigned int ndim, double integrand( double*, size_t, void* ), Parameters& params ) : ps_bin_( INVALID_BIN ), input_params_( params ), function_( new gsl_monte_function{ integrand, ndim, (void*)&input_params_ } ), rng_( gsl_rng_alloc( input_params_.integration().rng_engine ) ), grid_( new GridParameters( ndim ) ) { //--- initialise the random number generator unsigned long seed = ( input_params_.integration().rng_seed > 0 ) ? input_params_.integration().rng_seed : time( nullptr ); // seed with time gsl_rng_set( rng_.get(), seed ); //--- a bit of printout for debugging CG_DEBUG( "Integrator:build" ) << "Number of integration dimensions: " << function_->dim << ",\n\t" << "Number of function calls: " << input_params_.integration().ncvg << ",\n\t" << "Random numbers generator: " << gsl_rng_name( rng_.get() ) << "."; switch ( input_params_.integration().type ) { case IntegratorType::Vegas: CG_DEBUG( "Integrator:build" ) << "Vegas parameters:\n\t" << "Number of iterations in Vegas: " << input_params_.integration().vegas.iterations << ",\n\t" << "α-value: " << input_params_.integration().vegas.alpha << ",\n\t" << "Verbosity: " << input_params_.integration().vegas.verbose << ",\n\t" << "Grid interpolation mode: " << (Integrator::VegasMode)input_params_.integration().vegas.mode << "."; break; case IntegratorType::MISER: CG_DEBUG( "Integrator:build" ) << "MISER parameters:\n\t" << "Number of calls: " << input_params_.integration().miser.min_calls << ", " << "per bisection: " << input_params_.integration().miser.min_calls_per_bisection << ",\n\t" << "Estimate fraction: " << input_params_.integration().miser.estimate_frac << ",\n\t" << "α-value: " << input_params_.integration().miser.alpha << ",\n\t" << "Dither: " << input_params_.integration().miser.dither << "."; break; case IntegratorType::plain: break; } } Integrator::~Integrator() {} //----------------------------------------------------------------------------------------------- // integration part //----------------------------------------------------------------------------------------------- void Integrator::integrate( double& result, double& abserr ) { int res = -1; //--- integration bounds std::vector x_low( function_->dim, 0. ), x_up( function_->dim, 1. ); //--- launch integration switch ( input_params_.integration().type ) { case IntegratorType::plain: { std::unique_ptr pln_state( gsl_monte_plain_alloc( function_->dim ), gsl_monte_plain_free ); res = gsl_monte_plain_integrate( function_.get(), &x_low[0], &x_up[0], function_->dim, input_params_.integration().ncvg, rng_.get(), pln_state.get(), &result, &abserr ); } break; case IntegratorType::Vegas: { //----- warmup (prepare the grid) warmupVegas( x_low, x_up, 25000 ); //----- integration unsigned short it_chisq = 0; do { res = gsl_monte_vegas_integrate( function_.get(), &x_low[0], &x_up[0], function_->dim, 0.2 * input_params_.integration().ncvg, rng_.get(), veg_state_.get(), &result, &abserr ); CG_LOG( "Integrator:integrate" ) << "\t>> at call " << ( ++it_chisq ) << ": " << Form( "average = %10.6f " "sigma = %10.6f chi2 = %4.3f.", result, abserr, gsl_monte_vegas_chisq( veg_state_.get() ) ); } while ( fabs( gsl_monte_vegas_chisq( veg_state_.get() )-1. ) > input_params_.integration().vegas_chisq_cut-1. ); CG_DEBUG( "Integrator:integrate" ) << "Vegas grid information:\n\t" << "ran for " << veg_state_->dim << " dimensions, and generated " << veg_state_->bins_max << " bins.\n\t" << "Integration volume: " << veg_state_->vol << "."; grid_->r_boxes = std::pow( veg_state_->bins, function_->dim ); } break; case IntegratorType::MISER: { std::unique_ptr mis_state( gsl_monte_miser_alloc( function_->dim ), gsl_monte_miser_free ); gsl_monte_miser_params_set( mis_state.get(), &input_params_.integration().miser ); res = gsl_monte_miser_integrate( function_.get(), &x_low[0], &x_up[0], function_->dim, input_params_.integration().ncvg, rng_.get(), mis_state.get(), &result, &abserr ); } break; } input_params_.integration().result = result; input_params_.integration().err_result = abserr; if ( input_params_.hadroniser() ) input_params_.hadroniser()->setCrossSection( result, abserr ); + if ( input_params_.outputModule() ) + input_params_.outputModule()->setCrossSection( result, abserr ); if ( res != GSL_SUCCESS ) throw CG_FATAL( "Integrator:integrate" ) << "Error while performing the integration!\n\t" << "GSL error: " << gsl_strerror( res ) << "."; } void Integrator::warmupVegas( std::vector& x_low, std::vector& x_up, unsigned int ncall ) { // start by preparing the grid/state veg_state_.reset( gsl_monte_vegas_alloc( function_->dim ) ); gsl_monte_vegas_params_set( veg_state_.get(), &input_params_.integration().vegas ); // then perform a first integration with the given calls count double result = 0., abserr = 0.; const int res = gsl_monte_vegas_integrate( function_.get(), &x_low[0], &x_up[0], function_->dim, ncall, rng_.get(), veg_state_.get(), &result, &abserr ); // ensure the operation was successful if ( res != GSL_SUCCESS ) throw CG_ERROR( "Integrator:vegas" ) << "Failed to warm-up the Vegas grid.\n\t" << "GSL error: " << gsl_strerror( res ) << "."; CG_INFO( "Integrator:vegas" ) << "Finished the Vegas warm-up."; } //----------------------------------------------------------------------------------------------- // events generation part //----------------------------------------------------------------------------------------------- void Integrator::generateOne( std::function callback ) { if ( !grid_->gen_prepared ) computeGenerationParameters(); std::vector xtmp; //--- correction cycles if ( ps_bin_ != INVALID_BIN ) { bool has_correction = false; while ( !correctionCycle( xtmp, has_correction ) ) {} if ( has_correction ) { storeEvent( xtmp, callback ); return; } } double weight = 0.; //--- normal generation cycle while ( true ) { double y = -1.; //----- select a and reject if fmax is too small while ( true ) { // ... ps_bin_ = uniform() * grid_->size(); y = uniform() * grid_->globalMax(); grid_->num[ps_bin_] += 1; if ( y <= grid_->maxValue( ps_bin_ ) ) break; } // shoot a point x in this bin grid_->shoot( rng_.get(), ps_bin_, xtmp ); // get weight for selected x value weight = eval( xtmp ); if ( weight <= 0. ) continue; if ( weight > y ) break; } if ( weight <= grid_->maxValue( ps_bin_ ) ) ps_bin_ = INVALID_BIN; else { //--- if weight is higher than local or global maximum, // init correction cycle grid_->f_max_old = grid_->maxValue( ps_bin_ ); grid_->f_max_diff = weight-grid_->f_max_old; grid_->setValue( ps_bin_, weight ); grid_->correc = ( grid_->num[ps_bin_]-1. ) * grid_->f_max_diff / grid_->globalMax() - 1.; CG_DEBUG("Integrator::generateOne") << "Correction " << grid_->correc << " will be applied for phase space bin " << ps_bin_ << "."; } // return with an accepted event if ( weight > 0. ) storeEvent( xtmp, callback ); } void Integrator::generate( unsigned long num_events, std::function callback ) { if ( num_events < 1 ) num_events = input_params_.generation().maxgen; if ( input_params_.outputModule() ) input_params_.outputModule()->initialise( input_params_ ); try { while ( input_params_.numGeneratedEvents() < num_events ) generateOne( callback ); } catch ( const Exception& ) { return; } } bool Integrator::correctionCycle( std::vector& x, bool& has_correction ) { CG_DEBUG_LOOP( "Integrator:correction" ) << "Correction cycles are started.\n\t" << "bin = " << ps_bin_ << "\t" << "correc = " << grid_->correc << "\t" << "corre2 = " << grid_->correc2 << "."; if ( grid_->correc >= 1. ) grid_->correc -= 1.; if ( uniform() < grid_->correc ) { grid_->correc = -1.; std::vector xtmp( function_->dim ); // Select x values in phase space bin grid_->shoot( rng_.get(), ps_bin_, xtmp ); const double weight = eval( xtmp ); // Parameter for correction of correction if ( weight > grid_->maxValue( ps_bin_ ) ) { grid_->f_max2 = std::max( grid_->f_max2, weight ); grid_->correc += 1.; grid_->correc2 -= 1.; } // Accept event if ( weight >= grid_->f_max_diff*uniform() + grid_->f_max_old ) { x = xtmp; has_correction = true; return true; } return false; } // Correction if too big weight is found while correction // (All your bases are belong to us...) if ( grid_->f_max2 > grid_->maxValue( ps_bin_ ) ) { grid_->f_max_old = grid_->maxValue( ps_bin_ ); grid_->f_max_diff = grid_->f_max2-grid_->f_max_old; grid_->correc = ( grid_->num[ps_bin_]-1. ) * grid_->f_max_diff / grid_->globalMax(); if ( grid_->f_max2 >= grid_->globalMax() ) grid_->correc *= grid_->f_max2 / grid_->globalMax(); grid_->setValue( ps_bin_, grid_->f_max2 ); grid_->correc -= grid_->correc2; grid_->correc2 = 0.; grid_->f_max2 = 0.; return false; } return true; } bool Integrator::storeEvent( const std::vector& x, std::function callback ) { //--- start by computing the matrix element for that point const double weight = eval( x ); //--- reject if unphysical if ( weight <= 0. ) return false; { if ( input_params_.numGeneratedEvents() % input_params_.generation().gen_print_every == 0 ) { CG_INFO( "Integrator:store" ) << "Generated events: " << input_params_.numGeneratedEvents(); input_params_.process()->last_event->dump(); } const Event& last_event = *input_params_.process()->last_event; if ( callback ) callback( last_event, input_params_.numGeneratedEvents() ); input_params_.addGenerationTime( last_event.time_total ); if ( input_params_.outputModule() ) *input_params_.outputModule() << last_event; } return true; } //----------------------------------------------------------------------------------------------- // initial preparation run before the generation of unweighted events //----------------------------------------------------------------------------------------------- void Integrator::computeGenerationParameters() { input_params_.setStorage( false ); if ( input_params_.generation().treat && input_params_.integration().type != IntegratorType::Vegas ) { CG_INFO( "Integrator:setGen" ) << "Treat switched on without a proper Vegas grid; running a warm-up beforehand."; std::vector x_low( function_->dim, 0. ), x_up( function_->dim, 1. ); try { warmupVegas( x_low, x_up, 25000 ); } catch ( const Exception& ) { throw CG_FATAL( "Integrator::setGen" ) << "Failed to perform a Vegas warm-up.\n\t" << "Try to re-run while disabling integrand treatment..."; } } CG_INFO( "Integrator:setGen" ) << "Preparing the grid (" << input_params_.generation().num_points << " points/bin) " << "for the generation of unweighted events."; const double inv_num_points = 1./input_params_.generation().num_points; std::vector x( function_->dim, 0. ); std::vector n( function_->dim, 0 );; // ... double sum = 0., sum2 = 0., sum2p = 0.; utils::ProgressBar prog_bar( grid_->size(), 5 ); //--- main loop for ( unsigned int i = 0; i < grid_->size(); ++i ) { double fsum = 0., fsum2 = 0.; for ( unsigned int j = 0; j < input_params_.generation().num_points; ++j ) { grid_->shoot( rng_.get(), i, x ); const double weight = eval( x ); grid_->setValue( i, weight ); fsum += weight; fsum2 += weight*weight; } const double av = fsum*inv_num_points, av2 = fsum2*inv_num_points, sig2 = av2-av*av; sum += av; sum2 += av2; sum2p += sig2; // per-bin debugging loop if ( CG_LOG_MATCH( "Integrator:setGen", debugInsideLoop ) ) { const double sig = sqrt( sig2 ); const double eff = ( grid_->maxValue( i ) != 0. ) ? grid_->maxValue( i )/av : 1.e4; CG_DEBUG_LOOP( "Integrator:setGen" ) << "n-vector for bin " << i << ": " << utils::repr( grid_->n( i ) ) << "\n\t" << "av = " << av << "\n\t" << "sig = " << sig << "\n\t" << "fmax = " << grid_->maxValue( i ) << "\n\t" << "eff = " << eff; } prog_bar.update( i+1 ); } // end of main loop const double inv_max = 1./grid_->size(); sum *= inv_max; sum2 *= inv_max; sum2p *= inv_max; const double sig = sqrt( sum2-sum*sum ), sigp = sqrt( sum2p ); double eff1 = 0.; for ( unsigned int i = 0; i < grid_->size(); ++i ) eff1 += sum/grid_->size()*grid_->maxValue( i ); const double eff2 = sum/grid_->globalMax(); CG_DEBUG( "Integrator:setGen" ) << "Average function value = " << sum << "\n\t" << "Average squared function value = " << sum2 << "\n\t" << "Overall standard deviation = " << sig << "\n\t" << "Average standard deviation = " << sigp << "\n\t" << "Maximum function value = " << grid_->globalMax() << "\n\t" << "Average inefficiency = " << eff1 << "\n\t" << "Overall inefficiency = " << eff2; grid_->gen_prepared = true; input_params_.setStorage( true ); CG_INFO( "Integrator:setGen" ) << "Grid prepared! Now launching the production."; } //------------------------------------------------------------------------------------------------ // helper / alias methods //------------------------------------------------------------------------------------------------ unsigned short Integrator::dimensions() const { if ( !function_ ) return 0; return function_->dim; } double Integrator::eval( const std::vector& x ) { if ( !input_params_.generation().treat ) return function_->f( (double*)&x[0], function_->dim, (void*)&input_params_ ); //--- treatment of the integration grid double w = grid_->r_boxes; std::vector x_new( x.size() ); for ( unsigned short j = 0; j < function_->dim; ++j ) { //--- find surrounding coordinates and interpolate const double z = x[j]*veg_state_->bins; const unsigned int id = z; // coordinate of point before const double rel_pos = z-id; // position between coordinates (norm.) const double bin_width = ( id == 0 ) ? COORD( veg_state_, 1, j ) : COORD( veg_state_, id+1, j )-COORD( veg_state_, id, j ); //--- build new coordinate from linear interpolation x_new[j] = COORD( veg_state_, id+1, j )-bin_width*( 1.-rel_pos ); w *= bin_width; } return w*function_->f( (double*)&x_new[0], function_->dim, (void*)&input_params_ ); } double Integrator::uniform() const { return gsl_rng_uniform( rng_.get() ); } //------------------------------------------------------------------------------------------------ std::ostream& operator<<( std::ostream& os, const IntegratorType& type ) { switch ( type ) { case IntegratorType::plain: return os << "plain"; case IntegratorType::Vegas: return os << "Vegas"; case IntegratorType::MISER: return os << "MISER"; } return os; } std::ostream& operator<<( std::ostream& os, const Integrator::VegasMode& mode ) { switch ( mode ) { case Integrator::VegasMode::importance: return os << "importance"; case Integrator::VegasMode::importanceOnly: return os << "importance-only"; case Integrator::VegasMode::stratified: return os << "stratified"; } return os; } } diff --git a/CepGen/IO/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 index 6c214b9..e4031da 100644 --- a/CepGen/IO/TextHandler.cpp +++ b/CepGen/IO/TextHandler.cpp @@ -1,208 +1,330 @@ #include "CepGen/IO/ExportHandler.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/ParametersList.h" +#include "CepGen/Core/utils.h" + #include "CepGen/Event/Event.h" #include "CepGen/Parameters.h" + #include "CepGen/Version.h" +#include + +#include #include #include namespace cepgen { namespace io { /** * \brief Handler for the generic text file output * \author Laurent Forthomme * \date Jul 2019 */ class TextHandler : public GenericExportHandler { public: explicit TextHandler( const ParametersList& ); ~TextHandler(); void initialise( const Parameters& ) override; + void setCrossSection( double xsec, double ) override { xsec_ = xsec; } void operator<<( const Event& ) override; private: + short extractVariableProperties( const std::string& ); + std::string textHistogram( const std::string&, const gsl_histogram* ) const; /// Retrieve a named variable from a particle double variable( const Particle&, const std::string& ) const; /// Retrieve a named variable from the whole event double variable( const Event&, const std::string& ) const; static const std::regex rgx_select_id_, rgx_select_role_; static constexpr double INVALID_OUTPUT = -999.; + static constexpr size_t PLOT_WIDTH = 50; + static constexpr char PLOT_CHAR = '#'; - std::ofstream file_; + std::ofstream file_, hist_file_; const std::vector variables_; - const bool print_banner_, print_variables_; + const bool save_banner_, save_variables_; + const bool show_hists_, save_hists_; + const ParametersList hist_variables_; const std::string separator_; //--- variables definition + std::unordered_map variables_name_; + std::unordered_map variable_stored_; + typedef std::pair IndexedVariable; std::unordered_map > variables_per_id_; std::unordered_map > variables_per_role_; std::vector variables_for_event_; unsigned short num_vars_; std::ostringstream oss_vars_; + double xsec_; + //--- auxiliary helper maps const std::unordered_map role_str_ = { { "ib1", Particle::Role::IncomingBeam1 }, { "ib2", Particle::Role::IncomingBeam2 }, { "ob1", Particle::Role::OutgoingBeam1 }, { "ob2", Particle::Role::OutgoingBeam2 }, { "pa1", Particle::Role::Parton1 }, { "pa2", Particle::Role::Parton2 }, { "cs", Particle::Role::CentralSystem }, { "int", Particle::Role::Intermediate } }; typedef double( Particle::Momentum::*pMethod )(void) const; /// Mapping of string variables to momentum getter methods const std::unordered_map m_mom_str_ = { { "px", &Particle::Momentum::px }, { "py", &Particle::Momentum::py }, { "pz", &Particle::Momentum::pz }, { "pt", &Particle::Momentum::pt }, { "eta", &Particle::Momentum::eta }, { "phi", &Particle::Momentum::phi }, { "m", &Particle::Momentum::mass }, { "e", &Particle::Momentum::energy }, { "p", &Particle::Momentum::p }, { "pt2", &Particle::Momentum::pt2 }, { "th", &Particle::Momentum::theta }, { "y", &Particle::Momentum::rapidity } }; //--- kinematic variables double sqrts_; unsigned long num_evts_; + struct gsl_histogram_deleter + { + void operator()( gsl_histogram* h ) { + gsl_histogram_free( h ); + } + }; + std::unordered_map > hists_; }; const std::regex TextHandler::rgx_select_id_( "(\\w+)\\((\\d+)\\)" ); const std::regex TextHandler::rgx_select_role_( "(\\w+)\\(([a-z]+\\d?)\\)" ); TextHandler::TextHandler( const ParametersList& params ) : GenericExportHandler( "text" ), - file_ ( params.get( "filename", "output.txt" ) ), - variables_ ( params.get >( "variables" ) ), - print_banner_ ( params.get( "saveBanner", true ) ), - print_variables_( params.get( "saveVariables", true ) ), - separator_ ( params.get( "separator", "\t" ) ), - num_vars_( 0 ) + file_ ( params.get( "filename", "output.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 ) ), + hist_variables_( params.get( "histVariables" ) ), + separator_ ( params.get( "separator", "\t" ) ), + num_vars_( 0 ), xsec_( 1. ) { - std::smatch sm; + //--- first extract list of variables to store in output file oss_vars_.clear(); std::string sep; for ( const auto& var : variables_ ) { - if ( std::regex_match( var, sm, rgx_select_id_ ) ) - variables_per_id_[stod( sm[2].str() )].emplace_back( std::make_pair( num_vars_, sm[1].str() ) ); - else if ( std::regex_match( var, sm, rgx_select_role_ ) ) { - const auto& str_role = sm[2].str(); - if ( role_str_.count( str_role ) == 0 ) { - CG_WARNING( "TextHandler" ) - << "Invalid particle role retrieved from configuration: \"" << str_role << "\".\n\t" - << "Skipping the variable \"" << var << "\" in the output module."; - continue; - } - variables_per_role_[role_str_.at( str_role )].emplace_back( std::make_pair( num_vars_, sm[1].str() ) ); + auto id = extractVariableProperties( var ); + if ( id >= 0 ) { + oss_vars_ << sep << var, sep = separator_; + variable_stored_[id] = true; } - else // event-level variables - variables_for_event_.emplace_back( std::make_pair( num_vars_, var ) ); - oss_vars_ << sep << var, sep = separator_; - ++num_vars_; } + //--- then extract list of variables to be plotted in histogram + for ( const auto& var : hist_variables_.keys() ) { + auto id = extractVariableProperties( var ); + if ( id < 0 ) + continue; + const auto& hvar = hist_variables_.get( var ); + const int nbins = hvar.get( "nbins", 10 ); + const double min = hvar.get( "low", 0. ), max = hvar.get( "high", 1. ); + hists_[id].reset( gsl_histogram_alloc( nbins ) ); + gsl_histogram_set_ranges_uniform( hists_[id].get(), min, max ); + CG_INFO( "TextHandler" ) + << "Booking a histogram with " << nbins << " bin" << utils::s( nbins ) + << " between " << min << " and " << max << " for \"" << var << "\"."; + } + if ( save_hists_ && !hists_.empty() ) + hist_file_.open( "lastrun.hists.txt" ); } TextHandler::~TextHandler() { + //--- histograms printout + for ( const auto& var : hist_variables_.keys() ) { + const auto& vn = std::find_if( variables_name_.begin(), variables_name_.end(), + [&var]( auto&& p ) { return p.second == var; } ); + if ( vn == variables_name_.end() ) { + CG_WARNING( "TextHandler" ) + << "Failed to retrieve variable \"" << var << "\" for plotting."; + continue; + } + const auto& hist = hists_.at( vn->first ).get(); + gsl_histogram_scale( hist, xsec_/( num_evts_+1 ) ); + if ( show_hists_ ) + CG_INFO( "TextHandler" ) + << textHistogram( var, hist ); + if ( save_hists_ ) + hist_file_ << "\n" << textHistogram( var, hist ) << "\n"; + } + //--- finalisation of the output file file_.close(); } void TextHandler::initialise( const Parameters& params ) { sqrts_ = params.kinematics.sqrtS(); num_evts_ = 0ul; - if ( print_banner_ ) + if ( save_banner_ ) file_ << banner( params, "#" ) << "\n"; - if ( print_variables_ ) + if ( save_variables_ ) file_ << "# " << oss_vars_.str() << "\n"; + if ( save_hists_ && !hists_.empty() ) + hist_file_ << banner( params, "#" ) << "\n"; } void TextHandler::operator<<( const Event& ev ) { std::vector vars( num_vars_ ); //--- extract and order the variables to be retrieved //--- particle-level variables (indexed by integer id) for ( const auto& id_vars : variables_per_id_ ) { const auto& part = ev[id_vars.first]; //--- loop over the list of variables for this particle for ( const auto& var : id_vars.second ) vars[var.first] = variable( part, var.second ); } //--- particle-level variables (indexed by role) for ( const auto& role_vars : variables_per_role_ ) { const auto& part = ev[role_vars.first][0]; //--- loop over the list of variables for this particle for ( const auto& var : role_vars.second ) vars[var.first] = variable( part, var.second ); } //--- event-level variables for ( const auto& var : variables_for_event_ ) vars[var.first] = variable( ev, var.second ); //--- write down the variables list in the file std::string sep; - for ( const auto& var : vars ) - file_ << sep << var, sep = separator_; + unsigned short i = 0; + for ( const auto& var : vars ) { + if ( variable_stored_.count( i ) > 0 && variable_stored_.at( i ) ) + file_ << sep << var, sep = separator_; + if ( hists_.count( i ) > 0 ) + gsl_histogram_increment( hists_.at( i ).get(), var ); + ++i; + } file_ << "\n"; ++num_evts_; } double TextHandler::variable( const Particle& part, const std::string& var ) const { if ( m_mom_str_.count( var ) ) { auto meth = m_mom_str_.at( var ); return ( part.momentum().*meth )(); } - if ( var == "xi" ) return 1.-part.energy()*2./sqrts_; + if ( var == "xi" ) return 1.-part.momentum().energy()*2./sqrts_; if ( var == "pdg" ) return (double)part.integerPdgId(); if ( var == "charge" ) return part.charge(); if ( var == "status" ) return (double)part.status(); CG_WARNING( "TextHandler" ) << "Failed to retrieve variable \"" << var << "\"."; return INVALID_OUTPUT; } double TextHandler::variable( const Event& ev, const std::string& var ) const { if ( var == "np" ) return (double)ev.size(); if ( var == "nev" ) return (double)num_evts_+1; if ( var == "nob1" || var == "nob2" ) { unsigned short out = 0.; for ( const auto& part : ev[ var == "nob1" ? Particle::Role::OutgoingBeam1 : Particle::Role::OutgoingBeam2 ] ) if ( (int)part.status() > 0 ) out++; return (double)out; } if ( var == "tgen" ) return ev.time_generation; if ( var == "ttot" ) return ev.time_total; CG_WARNING( "TextHandler" ) << "Failed to retrieve the event-level variable \"" << var << "\"."; return INVALID_OUTPUT; } + + short + TextHandler::extractVariableProperties( const std::string& var ) + { + const auto& vn = std::find_if( variables_name_.begin(), variables_name_.end(), + [&var]( auto&& p ) { return p.second == var; } ); + if ( vn != variables_name_.end() ) + return vn->first; + std::smatch sm; + if ( std::regex_match( var, sm, rgx_select_id_ ) ) + variables_per_id_[std::stod( sm[2].str() )].emplace_back( std::make_pair( num_vars_, sm[1].str() ) ); + else if ( std::regex_match( var, sm, rgx_select_role_ ) ) { + const auto& str_role = sm[2].str(); + if ( role_str_.count( str_role ) == 0 ) { + CG_WARNING( "TextHandler" ) + << "Invalid particle role retrieved from configuration: \"" << str_role << "\".\n\t" + << "Skipping the variable \"" << var << "\" in the output module."; + return -1; + } + variables_per_role_[role_str_.at( str_role )].emplace_back( std::make_pair( num_vars_, sm[1].str() ) ); + } + else // event-level variables + variables_for_event_.emplace_back( std::make_pair( num_vars_, var ) ); + variables_name_[num_vars_] = var; + return num_vars_++; + } + + 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( 15, ' ' ); + os + << "plot of \"" << var << "\"\n" + << sep << std::string( PLOT_WIDTH-16-var.size(), ' ' ) + << "d(sig)/d" << var << " (pb/bin)\n" + << sep << Form( "%-5.2f", gsl_histogram_min_val( hist ) ) + << std::string( PLOT_WIDTH-9, ' ' ) + << Form( "%5.2f", 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( "[%6.2f,%6.2f):", min, max ) + << std::string( val, PLOT_CHAR ) << std::string( PLOT_WIDTH-val, ' ' ) + << ": " << Form( "%6.2f", value ); + } + os + << "\n" + << Form( "%15s", var.c_str() ) << ":" << std::string( PLOT_WIDTH, '.' ) << ":\n" // 2nd abscissa axis + << "\t(" + << "bin width=" << ( gsl_histogram_max( hist )-gsl_histogram_min( hist ) )/nbins << ", " + << "mean=" << gsl_histogram_mean( hist ) << ", " + << "st.dev.=" << gsl_histogram_sigma( hist ) + << ")"; + return os.str(); + } } } REGISTER_IO_MODULE( text, TextHandler ) 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/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 0b3531c..d484063 100644 --- a/cmake/UseEnvironment.cmake +++ b/cmake/UseEnvironment.cmake @@ -1,91 +1,93 @@ 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") 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) #--- 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 ${DELPHES_DIR} PATH_SUFFIXES lib) find_path(DELPHES_INCLUDE NAMES modules classes HINTS ${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/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-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; -} -