diff --git a/CepGen/Cards/Handler.h b/CepGen/Cards/Handler.h index ea1f331..20749c1 100644 --- a/CepGen/Cards/Handler.h +++ b/CepGen/Cards/Handler.h @@ -1,36 +1,36 @@ #ifndef CepGen_Cards_Handler_h #define CepGen_Cards_Handler_h #include "CepGen/Parameters.h" namespace CepGen { class Parameters; /// Location for all steering card parsers/writers - namespace Cards + namespace cards { /// Generic steering card handler class Handler { public: /// Build a configuration from an external steering card Handler() {} ~Handler() {} /// Retrieve a configuration from a parsed steering cart Parameters& parameters() { return params_; } /// Small utility to retrieve the extension of a filename /// (naive approach) static std::string getExtension( const char* filename ) { const std::string file( filename ); return file.substr( file.find_last_of( "." )+1 ); } protected: /// List of parameters parsed from a card handler Parameters params_; }; } } #endif diff --git a/CepGen/Cards/LpairHandler.cpp b/CepGen/Cards/LpairHandler.cpp index fc81d0c..9b886a8 100644 --- a/CepGen/Cards/LpairHandler.cpp +++ b/CepGen/Cards/LpairHandler.cpp @@ -1,217 +1,217 @@ #include "CepGen/Cards/LpairHandler.h" #include "CepGen/Core/ParametersList.h" #include "CepGen/Core/Exception.h" #include "CepGen/Physics/PDG.h" #include "CepGen/StructureFunctions/StructureFunctions.h" #include "CepGen/StructureFunctions/LHAPDF.h" #include "CepGen/Processes/ProcessesHandler.h" #include "CepGen/Hadronisers/Pythia8Hadroniser.h" #include <fstream> namespace CepGen { - namespace Cards + namespace cards { const int LpairHandler::kInvalid = 99999; //----- specialization for LPAIR input cards LpairHandler::LpairHandler( const char* file ) : proc_params_( new ParametersList ), str_fun_( 11 ), 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( ¶ms_ ); //--- parse all fields std::unordered_map<std::string, std::string> 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 ) << ")"; } f.close(); //--- parse the process name auto proc = CepGen::ProcessesHandler::get().build( proc_name_, *proc_params_ ); params_.setProcess( std::move( proc ) ); //--- parse the structure functions code const unsigned long kLHAPDFCodeDec = 10000000, kLHAPDFPartDec = 1000000; if ( str_fun_ / kLHAPDFCodeDec == 1 ) { // SF from parton params_.kinematics.structure_functions = sf::Parameterisation::build( sf::Type::LHAPDF ); auto sf = dynamic_cast<sf::LHAPDF*>( params_.kinematics.structure_functions.get() ); const unsigned long icode = str_fun_ % kLHAPDFCodeDec; sf->params.pdf_code = icode % kLHAPDFPartDec; sf->params.mode = (sf::LHAPDF::Parameters::Mode)( icode / kLHAPDFPartDec ); // 0, 1, 2 } else params_.kinematics.structure_functions = sf::Parameterisation::build( (sf::Type)str_fun_ ); //--- parse the integration algorithm name if ( integr_type_ == "plain" ) params_.integrator.type = Integrator::Type::plain; else if ( integr_type_ == "Vegas" ) params_.integrator.type = Integrator::Type::Vegas; else if ( integr_type_ == "MISER" ) params_.integrator.type = Integrator::Type::MISER; else if ( integr_type_ != "" ) throw CG_FATAL( "LpairHandler" ) << "Unrecognized integrator type: " << integr_type_ << "!"; //--- parse the hadronisation algorithm name if ( hadr_name_ == "pythia8" ) params_.setHadroniser( new hadroniser::Pythia8Hadroniser( params_, ParametersList() ) ); if ( m_params.count( "IEND" ) ) setValue<bool>( "IEND", ( std::stoi( m_params["IEND"] ) > 1 ) ); //--- 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( Parameters* params ) { //------------------------------------------------------------------------------------------- // Process/integration/hadronisation parameters //------------------------------------------------------------------------------------------- registerParameter<std::string>( "PROC", "Process name to simulate", &proc_name_ ); registerParameter<std::string>( "ITYP", "Integration algorithm", &integr_type_ ); registerParameter<std::string>( "HADR", "Hadronisation algorithm", &hadr_name_ ); registerParameter<std::string>( "KMRG", "KMR grid interpolation path", ¶ms_.kinematics.kmr_grid_path ); //------------------------------------------------------------------------------------------- // General parameters //------------------------------------------------------------------------------------------- registerParameter<bool>( "IEND", "Generation type", ¶ms->generation.enabled ); registerParameter<bool>( "NTRT", "Smoothen the integrand", ¶ms->generation.treat ); registerParameter<int>( "DEBG", "Debugging verbosity", (int*)&Logger::get().level ); registerParameter<int>( "NCVG", "Number of function calls", (int*)¶ms->integrator.ncvg ); registerParameter<int>( "ITVG", "Number of integration iterations", (int*)¶ms->integrator.vegas.iterations ); registerParameter<int>( "SEED", "Random generator seed", (int*)¶ms->integrator.rng_seed ); registerParameter<int>( "NTHR", "Number of threads to use for events generation", (int*)¶ms->generation.num_threads ); registerParameter<int>( "MODE", "Subprocess' mode", (int*)¶ms->kinematics.mode ); registerParameter<int>( "NCSG", "Number of points to probe", (int*)¶ms->generation.num_points ); registerParameter<int>( "NGEN", "Number of events to generate", (int*)¶ms->generation.maxgen ); registerParameter<int>( "NPRN", "Number of events before printout", (int*)¶ms->generation.gen_print_every ); //------------------------------------------------------------------------------------------- // Process-specific parameters //------------------------------------------------------------------------------------------- registerParameter<int>( "METH", "Computation method (kT-factorisation)", &proc_params_->operator[]<int>( "method" ) ); registerParameter<int>( "IPOL", "Polarisation states to consider", &proc_params_->operator[]<int>( "polarisationStates" ) ); //------------------------------------------------------------------------------------------- // Process kinematics parameters //------------------------------------------------------------------------------------------- registerParameter<int>( "PMOD", "Outgoing primary particles' mode", &str_fun_ ); registerParameter<int>( "EMOD", "Outgoing primary particles' mode", &str_fun_ ); registerParameter<int>( "PAIR", "Outgoing particles' PDG id", (int*)&proc_params_->operator[]<int>( "pair" ) ); registerParameter<int>( "INA1", "Heavy ion atomic weight (1st incoming beam)", (int*)&hi_1_.first ); registerParameter<int>( "INZ1", "Heavy ion atomic number (1st incoming beam)", (int*)&hi_1_.second ); registerParameter<int>( "INA2", "Heavy ion atomic weight (1st incoming beam)", (int*)&hi_2_.first ); registerParameter<int>( "INZ2", "Heavy ion atomic number (1st incoming beam)", (int*)&hi_2_.second ); registerParameter<double>( "INP1", "Momentum (1st primary particle)", ¶ms->kinematics.incoming_beams.first.pz ); registerParameter<double>( "INP2", "Momentum (2nd primary particle)", ¶ms->kinematics.incoming_beams.second.pz ); registerParameter<double>( "INPP", "Momentum (1st primary particle)", ¶ms->kinematics.incoming_beams.first.pz ); registerParameter<double>( "INPE", "Momentum (2nd primary particle)", ¶ms->kinematics.incoming_beams.second.pz ); registerParameter<double>( "PTCT", "Minimal transverse momentum (single central outgoing particle)", ¶ms->kinematics.cuts.central.pt_single.min() ); registerParameter<double>( "MSCT", "Minimal central system mass", ¶ms->kinematics.cuts.central.mass_sum.min() ); registerParameter<double>( "ECUT", "Minimal energy (single central outgoing particle)", ¶ms->kinematics.cuts.central.energy_single.min() ); registerParameter<double>( "ETMN", "Minimal pseudo-rapidity (central outgoing particles)", ¶ms->kinematics.cuts.central.eta_single.min() ); registerParameter<double>( "ETMX", "Maximal pseudo-rapidity (central outgoing particles)", ¶ms->kinematics.cuts.central.eta_single.max() ); registerParameter<double>( "YMIN", "Minimal rapidity (central outgoing particles)", ¶ms->kinematics.cuts.central.rapidity_single.min() ); registerParameter<double>( "YMAX", "Maximal rapidity (central outgoing particles)", ¶ms->kinematics.cuts.central.rapidity_single.max() ); registerParameter<double>( "Q2MN", "Minimal Q² = -q² (exchanged parton)", ¶ms->kinematics.cuts.initial.q2.min() ); registerParameter<double>( "Q2MX", "Maximal Q² = -q² (exchanged parton)", ¶ms->kinematics.cuts.initial.q2.max() ); registerParameter<double>( "MXMN", "Minimal invariant mass of proton remnants", ¶ms->kinematics.cuts.remnants.mass_single.min() ); registerParameter<double>( "MXMX", "Maximal invariant mass of proton remnants", ¶ms->kinematics.cuts.remnants.mass_single.max() ); } void LpairHandler::store( const char* file ) { std::ofstream f( file, std::fstream::out | std::fstream::trunc ); if ( !f.is_open() ) { 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<double>( key.c_str(), std::stod( value ) ); } catch ( std::invalid_argument& ) {} try { setValue<int>( key.c_str(), std::stoi( value ) ); } catch ( std::invalid_argument& ) {} //setValue<bool>( key.c_str(), std::stoi( value ) ); setValue<std::string>( key.c_str(), value ); } std::string LpairHandler::getParameter( std::string key ) const { double dd = getValue<double>( key.c_str() ); if ( dd != -999. ) return std::to_string( dd ); int ui = getValue<int>( key.c_str() ); if ( ui != 999 ) return std::to_string( ui ); //if ( out = getValue<bool>( key.c_str() ) ); return getValue<std::string>( key.c_str() ); } std::string LpairHandler::getDescription( 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 6425106..bd5fc89 100644 --- a/CepGen/Cards/LpairHandler.h +++ b/CepGen/Cards/LpairHandler.h @@ -1,115 +1,114 @@ #ifndef CepGen_Cards_LpairReader_h #define CepGen_Cards_LpairReader_h #include "CepGen/Cards/Handler.h" #include <unordered_map> using std::string; namespace CepGen { class ParametersList; - namespace Cards + namespace cards { /// 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: template<class T> 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<class T> void registerParameter( const char* key, const char* description, T* def ) {} /// Set a parameter value template<class T> void setValue( const char* key, const T& value ) {} /// Retrieve a parameter value template<class T> 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; static const int kInvalid; std::unordered_map<std::string, Parameter<std::string> > p_strings_; std::unordered_map<std::string, Parameter<double> > p_doubles_; std::unordered_map<std::string, Parameter<int> > p_ints_; std::unordered_map<std::string, Parameter<bool> > p_bools_; void init( Parameters* ); std::shared_ptr<ParametersList> proc_params_; int str_fun_; std::string proc_name_, hadr_name_, integr_type_; std::pair<unsigned short,unsigned short> hi_1_, hi_2_; }; //----- specialised registerers /// Register a string parameter template<> inline void LpairHandler::registerParameter<std::string>( const char* key, const char* description, std::string* def ) { p_strings_.insert( std::make_pair( key, Parameter<std::string>( key, description, def ) ) ); } /// Register a double floating point parameter template<> inline void LpairHandler::registerParameter<double>( const char* key, const char* description, double* def ) { p_doubles_.insert( std::make_pair( key, Parameter<double>( key, description, def ) ) ); } /// Register an integer parameter template<> inline void LpairHandler::registerParameter<int>( const char* key, const char* description, int* def ) { p_ints_.insert( std::make_pair( key, Parameter<int>( key, description, def ) ) ); } /// Register a boolean parameter template<> inline void LpairHandler::registerParameter<bool>( const char* key, const char* description, bool* def ) { p_bools_.insert( std::make_pair( key, Parameter<bool>( key, description, def ) ) ); } //----- specialised setters template<> inline void LpairHandler::setValue<std::string>( 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<double>( 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<int>( 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<bool>( const char* key, const bool& value ) { auto it = p_bools_.find( key ); if ( it != p_bools_.end() ) *it->second.value = value; } //----- specialised getters /// Retrieve a string parameter value template<> inline std::string LpairHandler::getValue( const char* key ) const { const auto& it = p_strings_.find( key ); if ( it != p_strings_.end() ) return *it->second.value; return "null"; } /// Retrieve a floating point parameter value template<> inline double LpairHandler::getValue( const char* key ) const { const auto& it = p_doubles_.find( key ); if ( it != p_doubles_.end() ) return *it->second.value; return -999.; } /// Retrieve an integer parameter value template<> inline int LpairHandler::getValue( const char* key ) const { const auto& it = p_ints_.find( key ); if ( it != p_ints_.end() ) return *it->second.value; return 999; } /// Retrieve a boolean parameter value template<> inline bool LpairHandler::getValue( const char* key ) const { const auto& it = p_bools_.find( key ); if ( it != p_bools_.end() ) return *it->second.value; return true; } } } #endif - diff --git a/CepGen/Cards/PythonHandler.cpp b/CepGen/Cards/PythonHandler.cpp index 83a7f77..a31b208 100644 --- a/CepGen/Cards/PythonHandler.cpp +++ b/CepGen/Cards/PythonHandler.cpp @@ -1,404 +1,404 @@ #include "CepGen/Cards/PythonHandler.h" #include "CepGen/Core/Exception.h" #ifdef PYTHON #include "CepGen/Core/TamingFunction.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/ParametersList.h" #include "CepGen/Processes/ProcessesHandler.h" #include "CepGen/StructureFunctions/StructureFunctions.h" #include "CepGen/StructureFunctions/LHAPDF.h" #include "CepGen/StructureFunctions/MSTWGrid.h" #include "CepGen/StructureFunctions/Schaefer.h" #include "CepGen/Hadronisers/Pythia8Hadroniser.h" #include <algorithm> #if PY_MAJOR_VERSION < 3 # define PYTHON2 #endif namespace CepGen { - namespace Cards + namespace cards { //----- specialization for CepGen input cards PythonHandler::PythonHandler( const char* file ) { setenv( "PYTHONPATH", ".:..:Cards", 1 ); std::string filename = getPythonPath( file ); const size_t fn_len = filename.length()+1; //Py_DebugFlag = 1; //Py_VerboseFlag = 1; #ifdef PYTHON2 char* sfilename = new char[fn_len]; snprintf( sfilename, fn_len, "%s", filename.c_str() ); #else wchar_t* sfilename = new wchar_t[fn_len]; swprintf( sfilename, fn_len, L"%s", filename.c_str() ); #endif if ( sfilename ) Py_SetProgramName( sfilename ); Py_InitializeEx( 1 ); if ( sfilename ) delete [] sfilename; if ( !Py_IsInitialized() ) throw CG_FATAL( "PythonHandler" ) << "Failed to initialise the Python cards parser!"; CG_INFO( "PythonHandler" ) << "Initialised the Python cards parser\n\t" << "Python version: " << Py_GetVersion() << "\n\t" << "Platform: " << Py_GetPlatform() << "."; PyObject* cfg = PyImport_ImportModule( filename.c_str() ); // new if ( !cfg ) throwPythonError( Form( "Failed to parse the configuration card %s", file ) ); PyObject* process = PyObject_GetAttrString( cfg, PROCESS_NAME ); // new if ( !process ) throwPythonError( Form( "Failed to extract a \"%s\" keyword from the configuration card %s", PROCESS_NAME, file ) ); //--- list of process-specific parameters ParametersList proc_params; fillParameter( process, "processParameters", proc_params ); //--- type of process to consider PyObject* pproc_name = getElement( process, MODULE_NAME ); // borrowed if ( !pproc_name ) throwPythonError( Form( "Failed to extract the process name from the configuration card %s", file ) ); const std::string proc_name = get<std::string>( pproc_name ); //--- process mode params_.kinematics.mode = (KinematicsMode)proc_params.get<int>( "mode", (int)KinematicsMode::invalid ); auto proc = CepGen::ProcessesHandler::get().build( proc_name, proc_params ); params_.setProcess( std::move( proc ) ); //--- process kinematics PyObject* pin_kinematics = getElement( process, "inKinematics" ); // borrowed if ( pin_kinematics ) parseIncomingKinematics( pin_kinematics ); PyObject* pout_kinematics = getElement( process, "outKinematics" ); // borrowed if ( pout_kinematics ) parseOutgoingKinematics( pout_kinematics ); //--- taming functions PyObject* ptam = getElement( process, "tamingFunctions" ); // borrowed if ( ptam ) parseTamingFunctions( ptam ); Py_CLEAR( process ); PyObject* plog = PyObject_GetAttrString( cfg, "logger" ); // new if ( plog ) { parseLogging( plog ); Py_CLEAR( plog ); } //--- hadroniser parameters PyObject* phad = PyObject_GetAttrString( cfg, "hadroniser" ); // new if ( phad ) { parseHadroniser( phad ); Py_CLEAR( phad ); } //--- generation parameters PyObject* pint = PyObject_GetAttrString( cfg, "integrator" ); // new if ( pint ) { parseIntegrator( pint ); Py_CLEAR( pint ); } PyObject* pgen = PyObject_GetAttrString( cfg, "generator" ); // new if ( pgen ) { parseGenerator( pgen ); Py_CLEAR( pgen ); } //--- finalisation Py_CLEAR( cfg ); } PythonHandler::~PythonHandler() { if ( Py_IsInitialized() ) Py_Finalize(); } void PythonHandler::parseIncomingKinematics( PyObject* kin ) { //--- retrieve the beams PDG ids std::vector<double> beams_pz; fillParameter( kin, "pz", beams_pz ); if ( beams_pz.size() == 2 ) { params_.kinematics.incoming_beams.first.pz = beams_pz.at( 0 ); params_.kinematics.incoming_beams.second.pz = beams_pz.at( 1 ); } //--- retrieve the beams longitudinal momentum std::vector<int> beams_pdg; fillParameter( kin, "pdgIds", beams_pdg ); if ( beams_pdg.size() == 2 ) { params_.kinematics.incoming_beams.first.pdg = (PDG)beams_pdg.at( 0 ); params_.kinematics.incoming_beams.second.pdg = (PDG)beams_pdg.at( 1 ); } double sqrt_s = -1.; fillParameter( kin, "cmEnergy", sqrt_s ); fillParameter( kin, "kmrGridPath", params_.kinematics.kmr_grid_path ); if ( sqrt_s != -1. ) params_.kinematics.setSqrtS( sqrt_s ); PyObject* psf = getElement( kin, "structureFunctions" ); // borrowed if ( psf ) parseStructureFunctions( psf, params_.kinematics.structure_functions ); std::vector<int> kt_fluxes; fillParameter( kin, "ktFluxes", kt_fluxes ); if ( kt_fluxes.size() > 0 ) params_.kinematics.incoming_beams.first.kt_flux = (KTFlux)kt_fluxes.at( 0 ); if ( kt_fluxes.size() > 1 ) params_.kinematics.incoming_beams.second.kt_flux = (KTFlux)kt_fluxes.at( 1 ); std::vector<int> hi_beam1, hi_beam2; fillParameter( kin, "heavyIonA", hi_beam1 ); if ( hi_beam1.size() == 2 ) params_.kinematics.incoming_beams.first.pdg = HeavyIon{ (unsigned short)hi_beam1[0], (Element)hi_beam1[1] }; fillParameter( kin, "heavyIonB", hi_beam2 ); if ( hi_beam2.size() == 2 ) params_.kinematics.incoming_beams.second.pdg = HeavyIon{ (unsigned short)hi_beam2[0], (Element)hi_beam2[1] }; } void PythonHandler::parseStructureFunctions( PyObject* psf, std::shared_ptr<sf::Parameterisation>& sf_handler ) { int str_fun = 0; fillParameter( psf, "id", str_fun ); sf_handler = sf::Parameterisation::build( (sf::Type)str_fun ); switch( (sf::Type)str_fun ) { case sf::Type::LHAPDF: { auto sf = std::dynamic_pointer_cast<sf::LHAPDF>( params_.kinematics.structure_functions ); fillParameter( psf, "pdfSet", sf->params.pdf_set ); fillParameter( psf, "numFlavours", (unsigned int&)sf->params.num_flavours ); fillParameter( psf, "pdfMember", (unsigned int&)sf->params.pdf_member ); fillParameter( psf, "mode", (unsigned int&)sf->params.mode ); } break; case sf::Type::MSTWgrid: { auto sf = std::dynamic_pointer_cast<mstw::Grid>( params_.kinematics.structure_functions ); fillParameter( psf, "gridPath", sf->params.grid_path ); } break; case sf::Type::Schaefer: { auto sf = std::dynamic_pointer_cast<sf::Schaefer>( params_.kinematics.structure_functions ); fillParameter( psf, "Q2cut", sf->params.q2_cut ); std::vector<double> w2_lims; fillParameter( psf, "W2limits", w2_lims ); if ( w2_lims.size() != 0 ) { if ( w2_lims.size() != 2 ) throwPythonError( Form( "Invalid size for W2limits attribute: %d != 2!", w2_lims.size() ) ); else { sf->params.w2_lo = *std::min_element( w2_lims.begin(), w2_lims.end() ); sf->params.w2_hi = *std::max_element( w2_lims.begin(), w2_lims.end() ); } } PyObject* pcsf = getElement( psf, "continuumSF" ); // borrowed if ( pcsf ) parseStructureFunctions( pcsf, sf->params.continuum_model ); PyObject* ppsf = getElement( psf, "perturbativeSF" ); // borrowed if ( ppsf ) parseStructureFunctions( ppsf, sf->params.perturbative_model ); PyObject* prsf = getElement( psf, "resonancesSF" ); // borrowed if ( prsf ) parseStructureFunctions( prsf, sf->params.resonances_model ); fillParameter( psf, "higherTwist", (bool&)sf->params.higher_twist ); } break; default: break; } } void PythonHandler::parseOutgoingKinematics( PyObject* kin ) { PyObject* pparts = getElement( kin, "minFinalState" ); // borrowed if ( pparts && PyTuple_Check( pparts ) ) for ( unsigned short i = 0; i < PyTuple_Size( pparts ); ++i ) params_.kinematics.minimum_final_state.emplace_back( (PDG)get<int>( PyTuple_GetItem( pparts, i ) ) ); PyObject* pcuts = getElement( kin, "cuts" ); // borrowed if ( pcuts ) parseParticlesCuts( pcuts ); // for LPAIR/collinear matrix elements fillLimits( kin, "q2", params_.kinematics.cuts.initial.q2 ); // for the kT factorised matrix elements fillLimits( kin, "qt", params_.kinematics.cuts.initial.qt ); fillLimits( kin, "phiqt", params_.kinematics.cuts.initial.phi_qt ); fillLimits( kin, "ptdiff", params_.kinematics.cuts.central.pt_diff ); fillLimits( kin, "phiptdiff", params_.kinematics.cuts.central.phi_pt_diff ); fillLimits( kin, "rapiditydiff", params_.kinematics.cuts.central.rapidity_diff ); // generic phase space limits fillLimits( kin, "rapidity", params_.kinematics.cuts.central.rapidity_single ); fillLimits( kin, "eta", params_.kinematics.cuts.central.eta_single ); fillLimits( kin, "pt", params_.kinematics.cuts.central.pt_single ); fillLimits( kin, "ptsum", params_.kinematics.cuts.central.pt_sum ); fillLimits( kin, "invmass", params_.kinematics.cuts.central.mass_sum ); fillLimits( kin, "mx", params_.kinematics.cuts.remnants.mass_single ); } void PythonHandler::parseParticlesCuts( PyObject* cuts ) { if ( !PyDict_Check( cuts ) ) throwPythonError( "Particle cuts object should be a dictionary!" ); PyObject* pkey = nullptr, *pvalue = nullptr; Py_ssize_t pos = 0; while ( PyDict_Next( cuts, &pos, &pkey, &pvalue ) ) { const PDG pdg = (PDG)get<int>( pkey ); fillLimits( pvalue, "pt", params_.kinematics.cuts.central_particles[pdg].pt_single ); fillLimits( pvalue, "energy", params_.kinematics.cuts.central_particles[pdg].energy_single ); fillLimits( pvalue, "eta", params_.kinematics.cuts.central_particles[pdg].eta_single ); fillLimits( pvalue, "rapidity", params_.kinematics.cuts.central_particles[pdg].rapidity_single ); } } void PythonHandler::parseLogging( PyObject* log ) { fillParameter( log, "level", (int&)Logger::get().level ); std::vector<std::string> enabled_modules; fillParameter( log, "enabledModules", enabled_modules ); for ( const auto& mod : enabled_modules ) Logger::get().addExceptionRule( mod ); } void PythonHandler::parseIntegrator( PyObject* integr ) { if ( !PyDict_Check( integr ) ) throwPythonError( "Integrator object should be a dictionary!" ); PyObject* palgo = getElement( integr, MODULE_NAME ); // borrowed if ( !palgo ) throwPythonError( "Failed to retrieve the integration algorithm name!" ); std::string algo = get<std::string>( palgo ); if ( algo == "plain" ) params_.integrator.type = Integrator::Type::plain; else if ( algo == "Vegas" ) { params_.integrator.type = Integrator::Type::Vegas; fillParameter( integr, "alpha", (double&)params_.integrator.vegas.alpha ); fillParameter( integr, "iterations", params_.integrator.vegas.iterations ); fillParameter( integr, "mode", (int&)params_.integrator.vegas.mode ); fillParameter( integr, "verbosity", (int&)params_.integrator.vegas.verbose ); std::string vegas_logging_output = "cerr"; fillParameter( integr, "loggingOutput", vegas_logging_output ); if ( vegas_logging_output == "cerr" ) // redirect all debugging information to the error stream params_.integrator.vegas.ostream = stderr; else if ( vegas_logging_output == "cout" ) // redirect all debugging information to the standard stream params_.integrator.vegas.ostream = stdout; else params_.integrator.vegas.ostream = fopen( vegas_logging_output.c_str(), "w" ); } else if ( algo == "MISER" ) { params_.integrator.type = Integrator::Type::MISER; fillParameter( integr, "estimateFraction", (double&)params_.integrator.miser.estimate_frac ); fillParameter( integr, "minCalls", params_.integrator.miser.min_calls ); fillParameter( integr, "minCallsPerBisection", params_.integrator.miser.min_calls_per_bisection ); fillParameter( integr, "alpha", (double&)params_.integrator.miser.alpha ); fillParameter( integr, "dither", (double&)params_.integrator.miser.dither ); } else throwPythonError( Form( "Invalid integration algorithm: %s", algo.c_str() ) ); fillParameter( integr, "numFunctionCalls", params_.integrator.ncvg ); fillParameter( integr, "seed", (unsigned long&)params_.integrator.rng_seed ); unsigned int rng_engine; fillParameter( integr, "rngEngine", rng_engine ); switch ( rng_engine ) { case 0: default: params_.integrator.rng_engine = (gsl_rng_type*)gsl_rng_mt19937; break; case 1: params_.integrator.rng_engine = (gsl_rng_type*)gsl_rng_taus2; break; case 2: params_.integrator.rng_engine = (gsl_rng_type*)gsl_rng_gfsr4; break; case 3: params_.integrator.rng_engine = (gsl_rng_type*)gsl_rng_ranlxs0; break; } fillParameter( integr, "chiSqCut", params_.integrator.vegas_chisq_cut ); } void PythonHandler::parseGenerator( PyObject* gen ) { if ( !PyDict_Check( gen ) ) throwPythonError( "Generation information object should be a dictionary!" ); params_.generation.enabled = true; fillParameter( gen, "treat", params_.generation.treat ); fillParameter( gen, "numEvents", params_.generation.maxgen ); fillParameter( gen, "printEvery", params_.generation.gen_print_every ); fillParameter( gen, "numThreads", params_.generation.num_threads ); fillParameter( gen, "numPoints", params_.generation.num_points ); } void PythonHandler::parseTamingFunctions( PyObject* tf ) { if ( !PyList_Check( tf ) ) throwPythonError( "Taming functions list should be a list!" ); for ( Py_ssize_t i = 0; i < PyList_Size( tf ); ++i ) { PyObject* pit = PyList_GetItem( tf, i ); // borrowed if ( !pit ) continue; if ( !PyDict_Check( pit ) ) throwPythonError( Form( "Item %d has invalid type %s", i, pit->ob_type->tp_name ) ); PyObject* pvar = getElement( pit, "variable" ), *pexpr = getElement( pit, "expression" ); // borrowed params_.taming_functions->add( get<std::string>( pvar ).c_str(), get<std::string>( pexpr ).c_str() ); } } void PythonHandler::parseHadroniser( PyObject* hadr ) { if ( !PyDict_Check( hadr ) ) throwPythonError( "Hadroniser object should be a dictionary!" ); PyObject* pname = getElement( hadr, MODULE_NAME ); // borrowed if ( !pname ) throwPythonError( "Hadroniser name is required!" ); std::string hadr_name = get<std::string>( pname ); //--- list of module-specific parameters ParametersList mod_params; fillParameter( hadr, "moduleParameters", mod_params ); if ( hadr_name == "pythia8" ) params_.setHadroniser( new hadroniser::Pythia8Hadroniser( params_, mod_params ) ); else throwPythonError( Form( "Unrecognised hadronisation algorithm: \"%s\"!", hadr_name.c_str() ) ); auto h = params_.hadroniser(); { //--- before calling the init() method std::vector<std::string> config; fillParameter( hadr, "preConfiguration", config ); h->readStrings( config ); } h->init(); { //--- after init() has been called std::vector<std::string> config; fillParameter( hadr, "processConfiguration", config ); for ( const auto& block : config ) { std::vector<std::string> config_blk; fillParameter( hadr, block.c_str(), config_blk ); h->readStrings( config_blk ); } } } } } #endif diff --git a/CepGen/Cards/PythonHandler.h b/CepGen/Cards/PythonHandler.h index 3497746..272a230 100644 --- a/CepGen/Cards/PythonHandler.h +++ b/CepGen/Cards/PythonHandler.h @@ -1,71 +1,71 @@ #ifndef CepGen_Cards_PythonHandler_h #define CepGen_Cards_PythonHandler_h #ifdef PYTHON #include <Python.h> #include "Handler.h" namespace CepGen { namespace sf { class Parameterisation; } class ParametersList; - namespace Cards + namespace cards { /// CepGen Python configuration cards reader/writer class PythonHandler : public Handler { public: /// Read a standard configuration card explicit PythonHandler( const char* file ); ~PythonHandler(); static PyObject* getElement( PyObject* obj, const char* key ); static PyObject* encode( const char* str ); private: static constexpr const char* MODULE_NAME = "mod_name"; static constexpr const char* PROCESS_NAME = "process"; static void throwPythonError( const std::string& message ); static std::string getPythonPath( const char* file ); template<typename T> bool is( PyObject* obj ) const; template<typename T> T get( PyObject* obj ) const; void fillLimits( PyObject* obj, const char* key, Limits& lim ); void fillParameter( PyObject* parent, const char* key, bool& out ); void fillParameter( PyObject* parent, const char* key, int& out ); void fillParameter( PyObject* parent, const char* key, unsigned long& out ); void fillParameter( PyObject* parent, const char* key, unsigned int& out ); void fillParameter( PyObject* parent, const char* key, double& out ); void fillParameter( PyObject* parent, const char* key, std::string& out ); void fillParameter( PyObject* parent, const char* key, std::vector<int>& out ); void fillParameter( PyObject* parent, const char* key, std::vector<double>& out ); void fillParameter( PyObject* parent, const char* key, std::vector<std::string>& out ); void fillParameter( PyObject* parent, const char* key, ParametersList& out ); void parseIncomingKinematics( PyObject* ); void parseOutgoingKinematics( PyObject* ); void parseParticlesCuts( PyObject* ); void parseLogging( PyObject* ); void parseIntegrator( PyObject* ); void parseGenerator( PyObject* ); void parseTamingFunctions( PyObject* ); void parseHadroniser( PyObject* ); void parseStructureFunctions( PyObject*, std::shared_ptr<sf::Parameterisation>& sf_handler ); }; template<> bool PythonHandler::is<int>( PyObject* obj ) const; template<> int PythonHandler::get<int>( PyObject* obj ) const; template<> unsigned long PythonHandler::get<unsigned long>( PyObject* obj ) const; template<> bool PythonHandler::is<ParametersList>( PyObject* obj ) const; template<> ParametersList PythonHandler::get<ParametersList>( PyObject* obj ) const; template<> bool PythonHandler::is<double>( PyObject* obj ) const; template<> double PythonHandler::get<double>( PyObject* obj ) const; template<> bool PythonHandler::is<std::string>( PyObject* obj ) const; template<> std::string PythonHandler::get<std::string>( PyObject* obj ) const; } } #endif #endif diff --git a/CepGen/Cards/PythonTypes.cpp b/CepGen/Cards/PythonTypes.cpp index 66361ff..ed2bde2 100644 --- a/CepGen/Cards/PythonTypes.cpp +++ b/CepGen/Cards/PythonTypes.cpp @@ -1,171 +1,170 @@ #include "CepGen/Cards/PythonHandler.h" #include "CepGen/Core/ParametersList.h" #include "CepGen/Core/utils.h" #ifdef PYTHON #if PY_MAJOR_VERSION < 3 # define PYTHON2 #endif namespace CepGen { - namespace Cards + namespace cards { //------------------------------------------------------------------ // typed retrieval helpers //------------------------------------------------------------------ template<> bool PythonHandler::is<int>( PyObject* obj ) const { #ifdef PYTHON2 return ( PyInt_Check( obj ) || PyBool_Check( obj ) ); #else return ( PyLong_Check( obj ) || PyBool_Check( obj ) ); #endif } template<> int PythonHandler::get<int>( PyObject* obj ) const { #ifdef PYTHON2 return PyInt_AsLong( obj ); #else return PyLong_AsLong( obj ); #endif } template<> unsigned long PythonHandler::get<unsigned long>( PyObject* obj ) const { #ifdef PYTHON2 return PyInt_AsUnsignedLongMask( obj ); #else if ( !PyLong_Check( obj ) ) throwPythonError( Form( "Object \"%s\" has invalid type %s", key, obj->ob_type->tp_name ) ); return PyLong_AsUnsignedLong( obj ); #endif } template<> long long PythonHandler::get<long long>( PyObject* obj ) const { return PyLong_AsLongLong( obj ); } template<> bool PythonHandler::is<double>( PyObject* obj ) const { return PyFloat_Check( obj ); } template<> double PythonHandler::get<double>( PyObject* obj ) const { return PyFloat_AsDouble( obj ); } template<> bool PythonHandler::is<std::string>( PyObject* obj ) const { #ifdef PYTHON2 return PyString_Check( obj ); #else return PyUnicode_Check( obj ); #endif } template<> std::string PythonHandler::get<std::string>( PyObject* obj ) const { std::string out; #ifdef PYTHON2 out = PyString_AsString( obj ); // deprecated in python v3+ #else PyObject* pstr = PyUnicode_AsEncodedString( obj, "utf-8", "strict" ); // new if ( !pstr ) throwPythonError( "Failed to decode a Python object!" ); out = PyBytes_AS_STRING( pstr ); Py_CLEAR( pstr ); #endif return out; } template<> bool PythonHandler::is<ParametersList>( PyObject* obj ) const { return PyDict_Check( obj ); } template<> ParametersList PythonHandler::get<ParametersList>( PyObject* obj ) const { ParametersList out; PyObject* pkey = nullptr, *pvalue = nullptr; Py_ssize_t pos = 0; while ( PyDict_Next( obj, &pos, &pkey, &pvalue ) ) { const std::string skey = get<std::string>( pkey ); if ( is<int>( pvalue ) ) out.set<int>( skey, get<int>( pvalue ) ); else if ( is<double>( pvalue ) ) out.set<double>( skey, get<double>( pvalue ) ); else if ( is<std::string>( pvalue ) ) out.set<std::string>( skey, get<std::string>( pvalue ) ); else if ( is<ParametersList>( pvalue ) ) out.set<ParametersList>( skey, get<ParametersList>( pvalue ) ); else if ( PyTuple_Check( pvalue ) || PyList_Check( pvalue ) ) { // vector PyObject* pfirst = PyTuple_GetItem( pvalue, 0 ); PyObject* pit = nullptr; const bool tuple = PyTuple_Check( pvalue ); const Py_ssize_t num_entries = ( tuple ) ? PyTuple_Size( pvalue ) : PyList_Size( pvalue ); if ( is<int>( pfirst ) ) { std::vector<int> vec; for ( Py_ssize_t i = 0; i < num_entries; ++i ) { pit = ( tuple ) ? PyTuple_GetItem( pvalue, i ) : PyList_GetItem( pvalue, i ); if ( pit->ob_type != pfirst->ob_type ) throwPythonError( Form( "Mixed types detected in vector '%s'", skey.c_str() ) ); vec.emplace_back( get<int>( pit ) ); } out.set<std::vector<int> >( skey, vec ); } else if ( is<double>( pfirst ) ) { std::vector<double> vec; for ( Py_ssize_t i = 0; i < num_entries; ++i ) { pit = ( tuple ) ? PyTuple_GetItem( pvalue, i ) : PyList_GetItem( pvalue, i ); if ( pit->ob_type != pfirst->ob_type ) throwPythonError( Form( "Mixed types detected in vector '%s'", skey.c_str() ) ); vec.emplace_back( get<double>( pit ) ); } out.set<std::vector<double> >( skey, vec ); } else if ( is<std::string>( pfirst ) ) { std::vector<std::string> vec; for ( Py_ssize_t i = 0; i < num_entries; ++i ) { pit = ( tuple ) ? PyTuple_GetItem( pvalue, i ) : PyList_GetItem( pvalue, i ); if ( pit->ob_type != pfirst->ob_type ) throwPythonError( Form( "Mixed types detected in vector '%s'", skey.c_str() ) ); vec.emplace_back( get<std::string>( pit ) ); } out.set<std::vector<std::string> >( skey, vec ); } else if ( is<ParametersList>( pfirst ) ) { std::vector<ParametersList> vec; for ( Py_ssize_t i = 0; i < num_entries; ++i ) { pit = ( tuple ) ? PyTuple_GetItem( pvalue, i ) : PyList_GetItem( pvalue, i ); if ( pit->ob_type != pfirst->ob_type ) throwPythonError( Form( "Mixed types detected in vector '%s'", skey.c_str() ) ); vec.emplace_back( get<ParametersList>( pit ) ); } out.set<std::vector<ParametersList> >( skey, vec ); } } } return out; } } } #endif - diff --git a/CepGen/Cards/PythonUtils.cpp b/CepGen/Cards/PythonUtils.cpp index b53a57d..32bca15 100644 --- a/CepGen/Cards/PythonUtils.cpp +++ b/CepGen/Cards/PythonUtils.cpp @@ -1,251 +1,251 @@ #include "CepGen/Cards/PythonHandler.h" #include "CepGen/Core/ParametersList.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/utils.h" #ifdef PYTHON #include <string> #include <algorithm> #include <frameobject.h> #if PY_MAJOR_VERSION < 3 # define PYTHON2 #endif namespace CepGen { - namespace Cards + namespace cards { //------------------------------------------------------------------ // Python API helpers //------------------------------------------------------------------ std::string PythonHandler::getPythonPath( const char* file ) { std::string s_filename = file; s_filename = s_filename.substr( 0, s_filename.find_last_of( "." ) ); // remove the extension std::replace( s_filename.begin(), s_filename.end(), '/', '.' ); // replace all '/' by '.' return s_filename; } void PythonHandler::throwPythonError( const std::string& message ) { PyObject* ptype = nullptr, *pvalue = nullptr, *ptraceback_obj = nullptr; // retrieve error indicator and clear it to handle ourself the error PyErr_Fetch( &ptype, &pvalue, &ptraceback_obj ); PyErr_Clear(); // ensure the objects retrieved are properly normalised and point to compatible objects PyErr_NormalizeException( &ptype, &pvalue, &ptraceback_obj ); std::ostringstream oss; oss << message; if ( ptype != nullptr ) { // we can start the traceback oss << "\n\tError: " #ifdef PYTHON2 << PyString_AsString( PyObject_Str( pvalue ) ); // deprecated in python v3+ #else << PyUnicode_AsUTF8( PyObject_Str( pvalue ) ); #endif PyTracebackObject* ptraceback = (PyTracebackObject*)ptraceback_obj; std::string tabul = "↪ "; if ( ptraceback != nullptr ) { while ( ptraceback->tb_next != nullptr ) { PyFrameObject* pframe = ptraceback->tb_frame; if ( pframe != nullptr ) { int line = PyCode_Addr2Line( pframe->f_code, pframe->f_lasti ); #ifdef PYTHON2 const char* filename = PyString_AsString( pframe->f_code->co_filename ); const char* funcname = PyString_AsString( pframe->f_code->co_name ); #else const char* filename = PyUnicode_AsUTF8( pframe->f_code->co_filename ); const char* funcname = PyUnicode_AsUTF8( pframe->f_code->co_name ); #endif oss << Form( "\n\t%s%s on %s (line %d)", tabul.c_str(), boldify( funcname ).c_str(), filename, line ); } else oss << Form( "\n\t%s issue in line %d", tabul.c_str(), ptraceback->tb_lineno ); tabul = std::string( " " )+tabul; ptraceback = ptraceback->tb_next; } } } Py_Finalize(); throw CG_FATAL( "PythonHandler:error" ) << oss.str(); } PyObject* PythonHandler::encode( const char* str ) { PyObject* obj = PyUnicode_FromString( str ); // new if ( !obj ) throwPythonError( Form( "Failed to encode the following string:\n\t%s", str ) ); return obj; } PyObject* PythonHandler::getElement( PyObject* obj, const char* key ) { PyObject* pout = nullptr, *nink = encode( key ); if ( !nink ) return pout; pout = PyDict_GetItem( obj, nink ); // borrowed Py_CLEAR( nink ); if ( pout ) CG_DEBUG( "PythonHandler:getElement" ) << "retrieved " << pout->ob_type->tp_name << " element \"" << key << "\" " << "from " << obj->ob_type->tp_name << " object\n\t" << "new reference count: " << pout->ob_refcnt; else CG_DEBUG( "PythonHandler:getElement" ) << "did not retrieve a valid element \"" << key << "\""; return pout; } void PythonHandler::fillLimits( PyObject* obj, const char* key, Limits& lim ) { PyObject* pobj = getElement( obj, key ); // borrowed if ( !pobj ) return; if ( !PyTuple_Check( pobj ) ) throw CG_FATAL( "PythonHandler:fillLimits" ) << "Invalid value retrieved for " << key << "."; if ( PyTuple_Size( pobj ) < 1 ) throw CG_FATAL( "PythonHandler:fillLimits" ) << "Invalid number of values unpacked for " << key << "!"; double min = get<double>( PyTuple_GetItem( pobj, 0 ) ); lim.min() = min; if ( PyTuple_Size( pobj ) > 1 ) { double max = get<double>( PyTuple_GetItem( pobj, 1 ) ); if ( max != -1 ) lim.max() = max; } } void PythonHandler::fillParameter( PyObject* parent, const char* key, bool& out ) { PyObject* pobj = getElement( parent, key ); // borrowed if ( !pobj ) return; if ( !PyBool_Check( pobj ) ) throwPythonError( Form( "Object \"%s\" has invalid type %s", key, pobj->ob_type->tp_name ) ); fillParameter( parent, key, (int&)out ); } void PythonHandler::fillParameter( PyObject* parent, const char* key, int& out ) { PyObject* pobj = getElement( parent, key ); // borrowed if ( !pobj ) return; if ( !is<int>( pobj ) ) throwPythonError( Form( "Object \"%s\" has invalid type %s", key, pobj->ob_type->tp_name ) ); out = get<int>( pobj ); } void PythonHandler::fillParameter( PyObject* parent, const char* key, unsigned long& out ) { PyObject* pobj = getElement( parent, key ); // borrowed if ( !pobj ) return; if ( !is<int>( pobj ) ) throwPythonError( Form( "Object \"%s\" has invalid type %s", key, pobj->ob_type->tp_name ) ); out = get<unsigned long>( pobj ); } void PythonHandler::fillParameter( PyObject* parent, const char* key, unsigned int& out ) { PyObject* pobj = getElement( parent, key ); // borrowed if ( !pobj ) return; if ( !is<int>( pobj ) ) throwPythonError( Form( "Object \"%s\" has invalid type %s", key, pobj->ob_type->tp_name ) ); out = get<unsigned long>( pobj ); } void PythonHandler::fillParameter( PyObject* parent, const char* key, double& out ) { PyObject* pobj = getElement( parent, key ); // borrowed if ( !pobj ) return; if ( !is<double>( pobj ) ) throwPythonError( Form( "Object \"%s\" has invalid type %s", key, pobj->ob_type->tp_name ) ); out = get<double>( pobj ); } void PythonHandler::fillParameter( PyObject* parent, const char* key, std::string& out ) { PyObject* pobj = getElement( parent, key ); // borrowed if ( !pobj ) return; if ( !is<std::string>( pobj ) ) throwPythonError( Form( "Object \"%s\" has invalid type %s", key, pobj->ob_type->tp_name ) ); out = get<std::string>( pobj ); } void PythonHandler::fillParameter( PyObject* parent, const char* key, std::vector<double>& out ) { out.clear(); PyObject* pobj = getElement( parent, key ); // borrowed if ( !pobj ) return; if ( !PyTuple_Check( pobj ) ) throwPythonError( Form( "Object \"%s\" has invalid type %s", key, pobj->ob_type->tp_name ) ); for ( Py_ssize_t i = 0; i < PyTuple_Size( pobj ); ++i ) { PyObject* pit = PyTuple_GetItem( pobj, i ); // borrowed if ( is<double>( pit ) ) out.emplace_back( get<double>( pit ) ); } } void PythonHandler::fillParameter( PyObject* parent, const char* key, std::vector<std::string>& out ) { out.clear(); PyObject* pobj = getElement( parent, key ); // borrowed if ( !pobj ) return; if ( !PyTuple_Check( pobj ) ) throwPythonError( Form( "Object \"%s\" has invalid type %s", key, pobj->ob_type->tp_name ) ); for ( Py_ssize_t i = 0; i < PyTuple_Size( pobj ); ++i ) { PyObject* pit = PyTuple_GetItem( pobj, i ); // borrowed out.emplace_back( get<std::string>( pit ) ); } } void PythonHandler::fillParameter( PyObject* parent, const char* key, std::vector<int>& out ) { out.clear(); PyObject* pobj = getElement( parent, key ); // borrowed if ( !pobj ) return; if ( !PyTuple_Check( pobj ) ) throwPythonError( Form( "Object \"%s\" has invalid type", key ) ); for ( Py_ssize_t i = 0; i < PyTuple_Size( pobj ); ++i ) { PyObject* pit = PyTuple_GetItem( pobj, i ); if ( !is<int>( pit ) ) throwPythonError( Form( "Object %d has invalid type", i ) ); out.emplace_back( get<int>( pit ) ); } } void PythonHandler::fillParameter( PyObject* parent, const char* key, ParametersList& out ) { PyObject* pobj = getElement( parent, key ); // borrowed if ( !pobj ) return; if ( !PyDict_Check( pobj ) ) throwPythonError( Form( "Object \"%s\" has invalid type", key ) ); out += get<ParametersList>( pobj ); } } } #endif diff --git a/CepGen/Core/FortranInterface.cpp b/CepGen/Core/FortranInterface.cpp index c5f6f8a..b02d21b 100644 --- a/CepGen/Core/FortranInterface.cpp +++ b/CepGen/Core/FortranInterface.cpp @@ -1,66 +1,66 @@ #include "CepGen/StructureFunctions/StructureFunctions.h" #include "CepGen/Physics/KTFlux.h" #include "CepGen/Physics/HeavyIon.h" #include "CepGen/Physics/ParticleProperties.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; sf::Type sf_mode = (sf::Type)sfmode; CG_DEBUG( "cepgen_structure_functions" ) << sf_mode; static auto& val = sf::Parameterisation::build( sf_mode )->operator()( xbj, q2 ); f2 = val.F2; fl = val.FL; } double cepgen_kt_flux_( int& fmode, double& x, double& kt2, int& sfmode, double& mx ) { using namespace CepGen; static auto sf = sf::Parameterisation::build( (sf::Type)sfmode ); return ktFlux( (KTFlux)fmode, x, kt2, *sf, mx ); } 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 } ); } double cepgen_particle_mass_( int& pdg_id ) { try { - return CepGen::ParticleProperties::mass( (CepGen::PDG)pdg_id ); + return CepGen::part::mass( (CepGen::PDG)pdg_id ); } catch ( const CepGen::Exception& e ) { e.dump(); exit( 0 ); } } double cepgen_particle_charge_( int& pdg_id ) { try { - return CepGen::ParticleProperties::charge( pdg_id ); + return CepGen::part::charge( pdg_id ); } catch ( const CepGen::Exception& e ) { e.dump(); exit( 0 ); } } #ifdef __cplusplus } #endif diff --git a/CepGen/Core/Generator.cpp b/CepGen/Core/Generator.cpp index a402891..d132743 100644 --- a/CepGen/Core/Generator.cpp +++ b/CepGen/Core/Generator.cpp @@ -1,171 +1,170 @@ #include "CepGen/Generator.h" #include "CepGen/Parameters.h" #include "CepGen/Version.h" #include "CepGen/Core/Integrator.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/Timer.h" #include "CepGen/Physics/PDG.h" #include "CepGen/Processes/GenericProcess.h" #include "CepGen/Event/Event.h" #include <fstream> #include <chrono> namespace CepGen { volatile int gSignal; Generator::Generator() : parameters( std::unique_ptr<Parameters>( new Parameters ) ), result_( -1. ), result_error_( -1. ) { CG_DEBUG( "Generator:init" ) << "Generator initialized"; try { printHeader(); } catch ( Exception& e ) { e.dump(); } // Random number initialization std::chrono::system_clock::time_point time = std::chrono::system_clock::now(); srandom( time.time_since_epoch().count() ); } Generator::Generator( Parameters* ip ) : parameters( ip ), result_( -1. ), result_error_( -1. ) {} Generator::~Generator() { if ( parameters->generation.enabled && parameters->process() && parameters->numGeneratedEvents() > 0 ) { CG_INFO( "Generator" ) << "Mean generation time / event: " << parameters->totalGenerationTime()*1.e3/parameters->numGeneratedEvents() << " ms."; } } size_t Generator::numDimensions() const { if ( !parameters->process() ) return 0; parameters->process()->addEventContent(); parameters->process()->setKinematics( parameters->kinematics ); return parameters->process()->numDimensions(); } void Generator::clearRun() { integrator_.reset(); parameters->process()->first_run = true; result_ = result_error_ = -1.; } void Generator::setParameters( Parameters& ip ) { parameters = std::unique_ptr<Parameters>( new Parameters( ip ) ); // copy constructor } void Generator::printHeader() { std::string tmp; std::ostringstream os; os << "version " << version() << std::endl; std::ifstream hf( "README" ); if ( !hf.good() ) throw CG_WARNING( "Generator" ) << "Failed to open README file."; while ( true ) { if ( !hf.good() ) break; getline( hf, tmp ); os << "\n " << tmp; } hf.close(); CG_INFO( "Generator" ) << os.str(); } double Generator::computePoint( double* x ) { - double res = Integrand::eval( x, numDimensions(), (void*)parameters.get() ); + double res = integrand::eval( x, numDimensions(), (void*)parameters.get() ); std::ostringstream os; for ( unsigned int i = 0; i < numDimensions(); ++i ) os << x[i] << " "; CG_DEBUG( "Generator:computePoint" ) << "Result for x[" << numDimensions() << "] = { " << os.str() << "}:\n\t" << res << "."; return res; } void Generator::computeXsection( double& xsec, double& err ) { CG_INFO( "Generator" ) << "Starting the computation of the process cross-section."; integrate(); xsec = result_; err = result_error_; if ( xsec < 1.e-2 ) CG_INFO( "Generator" ) << "Total cross section: " << xsec*1.e3 << " +/- " << err*1.e3 << " fb."; else if ( xsec > 5.e2 ) CG_INFO( "Generator" ) << "Total cross section: " << xsec*1.e-3 << " +/- " << err*1.e-3 << " nb."; else CG_INFO( "Generator" ) << "Total cross section: " << xsec << " +/- " << err << " pb."; } void Generator::integrate() { // first destroy and recreate the integrator instance if ( !integrator_ ) - integrator_ = std::unique_ptr<Integrator>( new Integrator( numDimensions(), Integrand::eval, parameters.get() ) ); + integrator_ = std::unique_ptr<Integrator>( new Integrator( numDimensions(), integrand::eval, parameters.get() ) ); else if ( integrator_->dimensions() != numDimensions() ) - integrator_.reset( new Integrator( numDimensions(), Integrand::eval, parameters.get() ) ); + integrator_.reset( new Integrator( numDimensions(), integrand::eval, parameters.get() ) ); CG_DEBUG( "Generator:newInstance" ) << "New integrator instance created\n\t" << "Considered topology: " << parameters->kinematics.mode << " case\n\t" << "Will proceed with " << numDimensions() << "-dimensional integration."; const int res = integrator_->integrate( result_, result_error_ ); if ( res != 0 ) throw CG_FATAL( "Generator" ) << "Error while computing the cross-section!\n\t" << "GSL error: " << gsl_strerror( res ) << "."; } std::shared_ptr<Event> Generator::generateOneEvent() { integrator_->generateOne(); parameters->addGenerationTime( parameters->process()->last_event->time_total ); return parameters->process()->last_event; } void Generator::generate( std::function<void( const Event&, unsigned long )> callback ) { const Timer tmr; CG_INFO( "Generator" ) << parameters->generation.maxgen << " events will be generated."; integrator_->generate( parameters->generation.maxgen, callback, &tmr ); const double gen_time_s = tmr.elapsed(); CG_INFO( "Generator" ) << parameters->generation.ngen << " events generated " << "in " << gen_time_s << " s " << "(" << gen_time_s/parameters->generation.ngen*1.e3 << " ms/event)."; } } - diff --git a/CepGen/Core/Integrand.cpp b/CepGen/Core/Integrand.cpp index 2e4f9ad..f73e8fb 100644 --- a/CepGen/Core/Integrand.cpp +++ b/CepGen/Core/Integrand.cpp @@ -1,201 +1,201 @@ #include "CepGen/Core/Timer.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/TamingFunction.h" #include "CepGen/Event/Event.h" #include "CepGen/Event/Particle.h" #include "CepGen/Physics/Kinematics.h" #include "CepGen/Physics/PDG.h" #include "CepGen/Processes/GenericProcess.h" #include "CepGen/Hadronisers/GenericHadroniser.h" #include "CepGen/Parameters.h" #include <sstream> #include <fstream> namespace CepGen { - namespace Integrand + namespace integrand { Logger::Level log_level; Timer tmr; double eval( double* x, size_t ndim, void* params ) { log_level = Logger::get().level; std::shared_ptr<Event> ev; Parameters* p = static_cast<Parameters*>( params ); if ( !p ) throw CG_FATAL( "Integrand" ) << "Failed to retrieve the run parameters!"; process::GenericProcess* proc = p->process(); if ( !proc ) throw CG_FATAL( "Integrand" ) << "Failed to retrieve the process!"; //============================================================================================= // start the timer //============================================================================================= tmr.reset(); //============================================================================================= // prepare the event content prior to the process generation //============================================================================================= if ( proc->hasEvent() ) { ev = proc->event(); if ( proc->first_run ) { CG_DEBUG( "Integrand" ) << "Computation launched for " << p->processName() << " process " << "0x" << std::hex << p->process() << std::dec << ".\n\t" << "Process mode considered: " << p->kinematics.mode << "\n\t" << " pz(p1) = " << p->kinematics.incoming_beams.first.pz << "\n\t" << " pz(p2) = " << p->kinematics.incoming_beams.second.pz << "\n\t" << " structure functions: " << p->kinematics.structure_functions; p->clearRunStatistics(); proc->first_run = false; } // passed the first-run preparation proc->clearEvent(); } // event is not empty //============================================================================================= // specify the phase space point to probe //============================================================================================= proc->setPoint( ndim, x ); //============================================================================================= // from this step on, the phase space point is supposed to be set //============================================================================================= p->process()->beforeComputeWeight(); double integrand = p->process()->computeWeight(); //============================================================================================= // invalidate any unphysical behaviour //============================================================================================= if ( integrand <= 0. ) return 0.; //============================================================================================= // speed up the integration process if no event needs to be generated //============================================================================================= if ( !p->storage() && !p->taming_functions && !p->hadroniser() && p->kinematics.cuts.central_particles.size() == 0 ) return integrand; //============================================================================================= // fill in the process' Event object //============================================================================================= p->process()->fillKinematics(); //============================================================================================= // once the kinematics variables have been populated, can apply the collection of taming functions //============================================================================================= if ( p->taming_functions ) { if ( p->taming_functions->has( "m_central" ) || p->taming_functions->has( "pt_central" ) ) { // build the kinematics of the central system Particle::Momentum central_system; for ( const auto& part : ev->getByRole( Particle::CentralSystem ) ) central_system += part.momentum(); // tame the cross-section by the reweighting function if ( p->taming_functions->has( "m_central" ) ) integrand *= p->taming_functions->eval( "m_central", central_system.mass() ); if ( p->taming_functions->has( "pt_central" ) ) integrand *= p->taming_functions->eval( "pt_central", central_system.pt() ); } if ( p->taming_functions->has( "q2" ) ) { integrand *= p->taming_functions->eval( "q2", -ev->getOneByRole( Particle::Parton1 ).momentum().mass() ); integrand *= p->taming_functions->eval( "q2", -ev->getOneByRole( Particle::Parton2 ).momentum().mass() ); } } if ( integrand <= 0. ) return 0.; //============================================================================================= // set the CepGen part of the event generation //============================================================================================= if ( p->storage() ) ev->time_generation = tmr.elapsed(); //============================================================================================= // event hadronisation and resonances decay //============================================================================================= if ( p->hadroniser() ) { double br = -1.; if ( !p->hadroniser()->run( *ev, br, p->storage() ) || br == 0. ) return 0.; integrand *= br; // branching fraction for all decays } //============================================================================================= // apply cuts on final state system (after hadronisation!) // (watch out your cuts, as this might be extremely time-consuming...) //============================================================================================= if ( p->kinematics.cuts.central_particles.size() > 0 ) { for ( const auto& part : ev->getByRole( Particle::CentralSystem ) ) { // retrieve all cuts associated to this final state particle if ( p->kinematics.cuts.central_particles.count( part.pdgId() ) == 0 ) continue; const auto& cuts_pdgid = p->kinematics.cuts.central_particles.at( part.pdgId() ); // apply these cuts on the given particle if ( !cuts_pdgid.pt_single.passes( part.momentum().pt() ) ) return 0.; if ( !cuts_pdgid.energy_single.passes( part.momentum().energy() ) ) return 0.; if ( !cuts_pdgid.eta_single.passes( part.momentum().eta() ) ) return 0.; if ( !cuts_pdgid.rapidity_single.passes( part.momentum().rapidity() ) ) return 0.; } } //============================================================================================= // store the last event in the parameters for its usage by the end user //============================================================================================= if ( p->storage() ) { p->process()->last_event = ev; p->process()->last_event->time_total = tmr.elapsed(); CG_DEBUG( "Integrand" ) << "[process 0x" << std::hex << p->process() << std::dec << "] " << "Individual time (gen+hadr+cuts): " << p->process()->last_event->time_total*1.e3 << " ms"; } //============================================================================================= // a bit of useful debugging //============================================================================================= if ( CG_EXCEPT_MATCH( "Integrand", debugInsideLoop ) ) { std::ostringstream oss; for ( unsigned short i = 0; i < ndim; ++i ) oss << Form( "%10.8f ", x[i] ); CG_DEBUG( "Integrand" ) << "f value for dim-" << ndim << " point ( " << oss.str() << "): " << integrand; } return integrand; } } } diff --git a/CepGen/Event/Particle.cpp b/CepGen/Event/Particle.cpp index 40a8306..69e0fa4 100644 --- a/CepGen/Event/Particle.cpp +++ b/CepGen/Event/Particle.cpp @@ -1,264 +1,264 @@ #include "CepGen/Event/Particle.h" #include "CepGen/Physics/PDG.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/utils.h" #include "CepGen/Physics/Constants.h" namespace CepGen { Particle::Particle() : id_( -1 ), charge_sign_( 1 ), mass_( -1. ), helicity_( 0. ), role_( UnknownRole ), status_( Status::Undefined ), pdg_id_( PDG::invalid ) {} Particle::Particle( Role role, PDG pdgId, Status st ) : id_( -1 ), charge_sign_( 1 ), mass_( -1. ), helicity_( 0. ), role_( role ), status_( st ), pdg_id_( pdgId ) { if ( pdg_id_ != PDG::invalid ) computeMass(); } Particle::Particle( const Particle& part ) : id_( part.id_ ), charge_sign_( part.charge_sign_ ), momentum_( part.momentum_ ), mass_( part.mass_ ), helicity_( part.helicity_ ), role_( part.role_ ), status_( part.status_ ), mothers_( part.mothers_ ), daughters_( part.daughters_ ), pdg_id_( part.pdg_id_ ) {} bool Particle::operator<( const Particle& rhs ) const { return id_ >= 0 && rhs.id_ > 0 && id_ < rhs.id_; } double Particle::thetaToEta( double theta ) { return -log( tan( 0.5*theta*M_PI/180. ) ); } double Particle::etaToTheta( double eta ) { return 2.*atan( exp( -eta ) )*180.*M_1_PI; } bool Particle::valid() { if ( pdg_id_ == PDG::invalid ) return false; if ( momentum_.p() == 0. && mass_ == 0. ) return false; return true; } void Particle::computeMass( bool off_shell ) { if ( !off_shell && pdg_id_ != PDG::invalid ) { // retrieve the mass from the on-shell particle's properties - mass_ = ParticleProperties::mass( pdg_id_ ); + mass_ = part::mass( pdg_id_ ); } else if ( momentum_.energy() >= 0. ) { mass_ = sqrt( energy2() - momentum_.p2() ); } //--- finish by setting the energy accordingly if ( momentum_.energy() < 0. ) { // invalid energy momentum_.setEnergy( sqrt( momentum_.p2() + mass2() ) ); } } void Particle::setMass( double m ) { if ( m >= 0. ) mass_ = m; else computeMass(); } void Particle::addMother( Particle& part ) { mothers_.insert( part.id() ); CG_DEBUG_LOOP( "Particle" ) << "Particle " << id() << " (pdgId=" << part.integerPdgId() << ") " << "is the new mother of " << id_ << " (pdgId=" << (int)pdg_id_ << ")."; part.addDaughter( *this ); } void Particle::addDaughter( Particle& part ) { const auto ret = daughters_.insert( part.id() ); if ( CG_EXCEPT_MATCH( "Particle", debugInsideLoop ) ) { std::ostringstream os; for ( const auto& daugh : daughters_ ) os << Form( "\n\t * id=%d", daugh ); CG_DEBUG_LOOP( "Particle" ) << "Particle " << role_ << " (pdgId=" << (int)pdg_id_ << ") " << "has now " << daughters_.size() << " daughter(s):" << os.str(); } if ( ret.second ) { CG_DEBUG_LOOP( "Particle" ) << "Particle " << part.role() << " (pdgId=" << part.integerPdgId() << ") " << "is a new daughter of " << role_ << " (pdgId=" << (int)pdg_id_ << "%4d)."; if ( part.mothers().find( id_ ) == part.mothers().end() ) part.addMother( *this ); } } void Particle::setMomentum( const Momentum& mom, bool offshell ) { momentum_ = mom; if ( !offshell && mom.mass() > 0. ) mass_ = momentum_.mass(); else computeMass(); } void Particle::setMomentum( double px, double py, double pz, double e ) { momentum_.setP( px, py, pz ); setEnergy( e ); if ( fabs( e-momentum_.energy() ) > 1.e-6 ) // more than 1 eV difference CG_ERROR( Form( "Energy difference: %.5e", e-momentum_.energy() ) ); } double Particle::energy() const { return ( momentum_.energy() < 0. ? std::hypot( mass_, momentum_.p() ) : momentum_.energy() ); } void Particle::setEnergy( double e ) { if ( e < 0. && mass_ >= 0. ) e = std::hypot( mass_, momentum_.p() ); momentum_.setEnergy( e ); } void Particle::setPdgId( short pdg ) { pdg_id_ = (PDG)abs( pdg ); switch ( pdg_id_ ) { case PDG::electron: case PDG::muon: case PDG::tau: charge_sign_ = -pdg/abs( pdg ); break; default: charge_sign_ = pdg/abs( pdg ); break; } } void Particle::setPdgId( const PDG& pdg, short ch ) { pdg_id_ = pdg; switch ( pdg_id_ ) { case PDG::electron: case PDG::muon: case PDG::tau: charge_sign_ = -ch; break; default: charge_sign_ = ch; break; } } int Particle::integerPdgId() const { - const float ch = ParticleProperties::charge( pdg_id_ ); + const float ch = part::charge( pdg_id_ ); if ( ch == 0 ) return static_cast<int>( pdg_id_ ); return static_cast<int>( pdg_id_ ) * charge_sign_ * ( ch/fabs( ch ) ); } void Particle::dump() const { std::ostringstream osm, osd; if ( !primary() ) { osm << ": mother(s): "; unsigned short i = 0; for ( const auto& moth : mothers_ ) { osm << ( i > 0 ? ", " : "" ) << moth; ++i; } } const ParticlesIds daughters_list = daughters(); if ( daughters_list.size() > 0 ) { osd << ": id = "; unsigned short i = 0; for ( const auto& daugh : daughters_list ) { osm << ( i > 0 ? ", " : "" ) << daugh; ++i; } } CG_INFO( "Particle" ) << "Dumping a particle with id=" << id_ << ", role=" << role_ << ", status=" << (int)status_ << "\n\t" << "Particle id: " << integerPdgId() << " (" << pdg_id_ << "), mass = " << mass_ << " GeV\n\t" << "Momentum: " << momentum_ << " GeV\t" << "(|P| = p = " << momentum_.p() << " GeV)\n\t" << " p⟂ = " << momentum_.pt() << " GeV, eta = " << momentum_.eta() << ", phi = " << momentum_.phi() << "\n\t" << "Primary? " << yesno( primary() ) << osm.str() << "\n\t" << numDaughters() << " daughter(s)" << osd.str(); } double Particle::etaToY( double eta_, double m_, double pt_ ) { const double m2 = m_*m_, mt = std::hypot( m_, pt_ ); return asinh( sqrt( ( ( mt*mt-m2 )*cosh( 2.*eta_ )+m2 )/ mt*mt - 1. )*M_SQRT1_2 ); } std::ostream& operator<<( std::ostream& os, const Particle::Role& rl ) { switch ( rl ) { case Particle::UnknownRole: return os << "unknown"; case Particle::IncomingBeam1: return os << "i.beam 1"; case Particle::IncomingBeam2: return os << "i.beam 2"; case Particle::OutgoingBeam1: return os << "o.beam 1"; case Particle::OutgoingBeam2: return os << "o.beam 2"; case Particle::Parton1: return os << "parton 1"; case Particle::Parton2: return os << "parton 2"; case Particle::Parton3: return os << "parton 3"; case Particle::Intermediate: return os << "hard pr."; case Particle::CentralSystem: return os << "central"; } return os; } double CMEnergy( const Particle& p1, const Particle& p2 ) { if ( p1.mass()*p2.mass() < 0. || p1.energy()*p2.energy() < 0. ) return 0.; return sqrt( p1.mass2()+p2.mass2() + 2.*p1.energy()*p2.energy() - 2.*( p1.momentum()*p2.momentum() ) ); } double CMEnergy( const Particle::Momentum& m1, const Particle::Momentum& m2 ) { if ( m1.mass()*m2.mass() < 0. || m1.energy()*m2.energy() < 0. ) return 0.; return sqrt( m1.mass2()+m2.mass2() + 2.*m1.energy()*m2.energy() - 2.*( m1*m2 ) ); } } diff --git a/CepGen/Event/Particle.h b/CepGen/Event/Particle.h index d76ec6b..298206a 100644 --- a/CepGen/Event/Particle.h +++ b/CepGen/Event/Particle.h @@ -1,346 +1,346 @@ #ifndef CepGen_Event_Particle_h #define CepGen_Event_Particle_h #include "CepGen/Physics/ParticleProperties.h" #include <set> #include <map> #include <vector> namespace CepGen { /// A set of integer-type particle identifiers typedef std::set<int> ParticlesIds; /// Kinematic information for one particle class Particle { public: /// Internal status code for a particle enum class Status { PrimordialIncoming = -9, DebugResonance = -5, Resonance = -4, Fragmented = -3, Propagator = -2, Incoming = -1, Undefined = 0, FinalState = 1, Undecayed = 2, Unfragmented = 3 }; /// Role of the particle in the process enum Role { UnknownRole = -1, IncomingBeam1 = 1, IncomingBeam2 = 2, OutgoingBeam1 = 3, OutgoingBeam2 = 5, CentralSystem = 6, Intermediate = 4, Parton1 = 41, Parton2 = 42, Parton3 = 43 }; /** * Container for a particle's 4-momentum, along with useful methods to ease the development of any matrix element level generator * \brief 4-momentum for a particle * \date Dec 2015 * \author Laurent Forthomme <laurent.forthomme@cern.ch> */ class Momentum { public: /// Build a 4-momentum at rest with an invalid energy (no mass information known) Momentum(); /// Build a 4-momentum using its 3-momentum coordinates and its energy Momentum( double x, double y, double z, double t = -1. ); /// Build a 4-momentum using its 3-momentum coordinates and its energy Momentum( double* p ); // --- static definitions /// Build a 3-momentum from its three pseudo-cylindric coordinates static Momentum fromPtEtaPhi( double pt, double eta, double phi, double e = -1. ); /// Build a 4-momentum from its scalar momentum, and its polar and azimuthal angles static Momentum fromPThetaPhi( double p, double theta, double phi, double e = -1. ); /// Build a 4-momentum from its four momentum and energy coordinates static Momentum fromPxPyPzE( double px, double py, double pz, double e ); /// Build a 4-momentum from its transverse momentum, rapidity and mass static Momentum fromPxPyYM( double px, double py, double rap, double m ); // --- vector and scalar operators /// Scalar product of the 3-momentum with another 3-momentum double threeProduct( const Momentum& ) const; /// Scalar product of the 4-momentum with another 4-momentum double fourProduct( const Momentum& ) const; /// Vector product of the 3-momentum with another 3-momentum double crossProduct( const Momentum& ) const; /// Add a 4-momentum through a 4-vector sum Momentum& operator+=( const Momentum& ); /// Subtract a 4-momentum through a 4-vector sum Momentum& operator-=( const Momentum& ); /// Scalar product of the 3-momentum with another 3-momentum double operator*=( const Momentum& ); /// Multiply all 4-momentum coordinates by a scalar Momentum& operator*=( double c ); /// Equality operator bool operator==( const Momentum& ) const; /// Human-readable format for a particle's momentum friend std::ostream& operator<<( std::ostream& os, const Particle::Momentum& mom ); Momentum& betaGammaBoost( double gamma, double betagamma ); /// Forward Lorentz boost Momentum& lorentzBoost( const Particle::Momentum& p ); // --- setters and getters /// Set all the components of the 4-momentum (in GeV) void setP( double px, double py, double pz, double e ); /// Set all the components of the 3-momentum (in GeV) void setP( double px, double py, double pz ); /// Set the energy (in GeV) inline void setEnergy( double e ) { energy_ = e; } /// Compute the energy from the mass inline void setMass( double m ) { setMass2( m*m ); } /// Compute the energy from the mass void setMass2( double m2 ); /// Get one component of the 4-momentum (in GeV) double operator[]( const unsigned int i ) const; /// Get one component of the 4-momentum (in GeV) double& operator[]( const unsigned int i ); /// Momentum along the \f$x\f$-axis (in GeV) inline double px() const { return px_; } /// Momentum along the \f$y\f$-axis (in GeV) inline double py() const { return py_; } /// Longitudinal momentum (in GeV) inline double pz() const { return pz_; } /// Transverse momentum (in GeV) double pt() const; /// Squared transverse momentum (in GeV\f${}^2\f$) double pt2() const; /// 4-vector of double precision floats (in GeV) const std::vector<double> pVector() const; /// 3-momentum norm (in GeV) inline double p() const { return p_; } /// Squared 3-momentum norm (in GeV\f${}^2\f$) inline double p2() const { return p_*p_; } /// Energy (in GeV) inline double energy() const { return energy_; } /// Squared energy (in GeV\f${}^2\f$) inline double energy2() const { return energy_*energy_; } /// Squared mass (in GeV\f${}^2\f$) as computed from its energy and momentum inline double mass2() const { return energy2()-p2(); } /// Mass (in GeV) as computed from its energy and momentum /// \note Returns \f$-\sqrt{|E^2-\mathbf{p}^2|}<0\f$ if \f$\mathbf{p}^2>E^2\f$ double mass() const; /// Polar angle (angle with respect to the longitudinal direction) double theta() const; /// Azimutal angle (angle in the transverse plane) double phi() const; /// Pseudo-rapidity double eta() const; /// Rapidity double rapidity() const; void truncate( double tolerance = 1.e-10 ); /// Rotate the transverse components by an angle phi (and reflect the y coordinate) Momentum& rotatePhi( double phi, double sign ); /// Rotate the particle's momentum by a polar/azimuthal angle Momentum& rotateThetaPhi( double theta_, double phi_ ); /// Apply a \f$ z\rightarrow -z\f$ transformation inline Momentum& mirrorZ() { pz_ = -pz_; return *this; } private: /// Compute the 3-momentum's norm void computeP(); /// Momentum along the \f$x\f$-axis double px_; /// Momentum along the \f$y\f$-axis double py_; /// Momentum along the \f$z\f$-axis double pz_; /// 3-momentum's norm (in GeV/c) double p_; /// Energy (in GeV) double energy_; }; /// Human-readable format for a particle's PDG code friend std::ostream& operator<<( std::ostream& os, const PDG& pc ); /// Human-readable format for a particle's role in the event friend std::ostream& operator<<( std::ostream& os, const Particle::Role& rl ); /// Compute the 4-vector sum of two 4-momenta friend Particle::Momentum operator+( const Particle::Momentum& mom1, const Particle::Momentum& mom2 ); /// Compute the 4-vector difference of two 4-momenta friend Particle::Momentum operator-( const Particle::Momentum& mom1, const Particle::Momentum& mom2 ); /// Compute the inverse per-coordinate 4-vector friend Particle::Momentum operator-( const Particle::Momentum& mom ); /// Scalar product of two 3-momenta friend double operator*( const Particle::Momentum& mom1, const Particle::Momentum& mom2 ); /// Multiply all components of a 4-momentum by a scalar friend Particle::Momentum operator*( const Particle::Momentum& mom, double c ); /// Multiply all components of a 4-momentum by a scalar friend Particle::Momentum operator*( double c, const Particle::Momentum& mom ); //----- static getters /// Convert a polar angle to a pseudo-rapidity static double thetaToEta( double theta ); /// Convert a pseudo-rapidity to a polar angle static double etaToTheta( double eta ); /// Convert a pseudo-rapidity to a rapidity static double etaToY( double eta_, double m_, double pt_ ); Particle(); /// Build using the role of the particle in the process and its PDG id /// \param[in] pdgId PDG identifier /// \param[in] role Role of the particle in the process /// \param[in] st Current status Particle( Role role, PDG pdgId, Status st = Status::Undefined ); /// Copy constructor Particle( const Particle& ); inline ~Particle() {} /// Comparison operator (from unique identifier) bool operator<( const Particle& rhs ) const; /// Comparison operator (from their reference's unique identifier) //bool operator<( Particle *rhs ) const { return ( id < rhs->id ); } // --- general particle properties /// Unique identifier (in a Event object context) int id() const { return id_; } //void setId( int id ) { id_ = id; } /// Set the particle unique identifier in an event void setId( int id ) { id_ = id; } /// Electric charge (given as a float number, for the quarks and bound states) - float charge() const { return charge_sign_ * ParticleProperties::charge( pdg_id_ ); } + float charge() const { return charge_sign_ * part::charge( pdg_id_ ); } /// Set the electric charge sign (+-1 for charged or 0 for neutral particles) void setChargeSign( int sign ) { charge_sign_ = sign; } /// Role in the considered process Role role() const { return role_; } /// Set the particle role in the process void setRole( const Role& role ) { role_ = role; } /** * Codes 1-10 correspond to currently existing partons/particles, and larger codes contain partons/particles which no longer exist, or other kinds of event information * \brief Particle status */ Status status() const { return status_; } /// Set the particle decay/stability status void setStatus( Status status ) { status_ = status; } /// Set the PDG identifier (along with the particle's electric charge) /// \param[in] pdg PDG identifier /// \param[in] ch Electric charge (0, 1, or -1) void setPdgId( const PDG& pdg, short ch = 0 ); /// Set the PDG identifier (along with the particle's electric charge) /// \param[in] pdg_id PDG identifier (incl. electric charge in e) void setPdgId( short pdg_id ); /// Retrieve the objectified PDG identifier inline PDG pdgId() const { return pdg_id_; } /// Retrieve the integer value of the PDG identifier int integerPdgId() const; /// Particle's helicity float helicity() const { return helicity_; } /// Set the helicity of the particle void setHelicity( float heli ) { helicity_ = heli; } /// Particle mass in GeV/c\f${}^2\f$ /// \return Particle's mass inline double mass() const { return mass_; }; /// Compute the particle mass /// \param[in] off_shell Allow the particle to be produced off-shell? /// \note This method ensures that the kinematics is properly set (the mass is set according to the energy and the momentum in priority) void computeMass( bool off_shell = false ); /// Set the particle mass, in GeV/c\f${}^2\f$ /// \param m Mass in GeV/c\f${}^2\f$ /// \note This method ensures that the kinematics is properly set (the mass is set according to the energy and the momentum in priority) void setMass( double m = -1. ); /// Particle squared mass, in GeV\f${}^2\f$/c\f${}^4\f$ inline double mass2() const { return mass_*mass_; }; /// Retrieve the momentum object associated with this particle inline Momentum& momentum() { return momentum_; } /// Retrieve the momentum object associated with this particle inline Momentum momentum() const { return momentum_; } /// Associate a momentum object to this particle void setMomentum( const Momentum& mom, bool offshell = false ); /** * \brief Set the 3- or 4-momentum associated to the particle * \param[in] px Momentum along the \f$x\f$-axis, in GeV/c * \param[in] py Momentum along the \f$y\f$-axis, in GeV/c * \param[in] pz Momentum along the \f$z\f$-axis, in GeV/c * \param[in] e Energy, in GeV */ void setMomentum( double px, double py, double pz, double e = -1. ); /// Set the 4-momentum associated to the particle /// \param[in] p 4-momentum inline void setMomentum( double p[4] ) { setMomentum( p[0], p[1], p[2], p[3] ); } /// Set the particle's energy /// \param[in] e Energy, in GeV void setEnergy( double e = -1. ); /// Get the particle's energy, in GeV double energy() const; /// Get the particle's squared energy, in GeV\f${}^2\f$ inline double energy2() const { return energy()*energy(); }; /// Is this particle a valid particle which can be used for kinematic computations? bool valid(); // --- particle relations /// Is this particle a primary particle? inline bool primary() const { return mothers_.empty(); } /// Set the mother particle /// \param[in] part A Particle object containing all the information on the mother particle void addMother( Particle& part ); /// Get the unique identifier to the mother particle from which this particle arises /// \return An integer representing the unique identifier to the mother of this particle in the event inline ParticlesIds mothers() const { return mothers_; } /** * \brief Add a decay product * \param[in] part The Particle object in which this particle will desintegrate or convert * \return A boolean stating if the particle has been added to the daughters list or if it was already present before */ void addDaughter( Particle& part ); /// Gets the number of daughter particles inline unsigned int numDaughters() const { return daughters_.size(); }; /// Get an identifiers list all daughter particles /// \return An integer vector containing all the daughters' unique identifier in the event inline ParticlesIds daughters() const { return daughters_; } // --- global particle information extraction /// Dump all the information on this particle into the standard output stream void dump() const; private: /// Unique identifier in an event int id_; /// Electric charge (+-1 or 0) short charge_sign_; /// Momentum properties handler Momentum momentum_; /// Mass, in GeV/c\f${}^2\f$ double mass_; /// Helicity float helicity_; /// Role in the process Role role_; /// Decay/stability status Status status_; /// List of mother particles ParticlesIds mothers_; /// List of daughter particles ParticlesIds daughters_; /// PDG id PDG pdg_id_; }; /// Compute the centre of mass energy of two particles (incoming or outgoing states) double CMEnergy( const Particle& p1, const Particle& p2 ); /// Compute the centre of mass energy of two particles (incoming or outgoing states) double CMEnergy( const Particle::Momentum& m1, const Particle::Momentum& m2 ); //bool operator<( const Particle& a, const Particle& b ) { return a.id<b.id; } // --- particle containers /// List of Particle objects typedef std::vector<Particle> Particles; /// List of particles' roles typedef std::vector<Particle::Role> ParticleRoles; /// Map between a particle's role and its associated Particle object typedef std::map<Particle::Role,Particles> ParticlesMap; } #endif diff --git a/CepGen/Generator.h b/CepGen/Generator.h index 3520ac3..b2afd73 100644 --- a/CepGen/Generator.h +++ b/CepGen/Generator.h @@ -1,122 +1,122 @@ #ifndef CepGen_Generator_h #define CepGen_Generator_h #include <sstream> #include <memory> #include <functional> //////////////////////////////////////////////////////////////////////////////// /** * \mainpage Foreword * This Monte Carlo generator was developed as a modern version of the LPAIR code introduced * in the early 1990s by J. Vermaseren *et al*\cite Baranov:1991yq\cite Vermaseren:1982cz. This latter allows to * compute the cross-section and to generate events for the \f$\gamma\gamma\to\ell^{+}\ell^{-}\f$ * process in the scope of high energy physics. * * Soon after the integration of its matrix element, it was extended as a tool to compute and * generate events for any generic 2\f$\rightarrow\f$ 3 central exclusive process. * To do so, the main operation performed here is the integration of the matrix element (given as a * subset of a GenericProcess object) over the full available phase space. * */ //////////////////////////////////////////////////////////////////////////////// /// Common namespace for this Monte Carlo generator namespace CepGen { - namespace Integrand + namespace integrand { /** * Function to be integrated. It returns the value of the weight for one point * of the full phase space (or "event"). This weights includes the matrix element * of the process considered, along with all the kinematic factors, and the cut * restrictions imposed on this phase space. \f$x\f$ is therefore an array of random * numbers defined inside its boundaries (as normalised so that \f$\forall i<\mathrm{ndim}\f$, * \f$0<x_i<1\f$. */ double eval( double*, size_t, void* ); } class Event; class Integrator; class Parameters; //////////////////////////////////////////////////////////////////////////////// /** * This object represents the core of this Monte Carlo generator, with its * capability to generate the events (using the embedded Vegas object) and to * study the phase space in term of the variation of resulting cross section * while scanning the various parameters (point \f${\bf x}\f$ in the * multi-dimensional phase space). * * The phase space is constrained using the Parameters object given as an * argument to the constructor, and the differential cross-sections for each * value of the array \f${\bf x}\f$ are computed in the \a f-function defined * outside (but populated inside) this object. * * This f-function embeds a GenericProcess-inherited object which defines all the * methods to compute this differential cross-section as well as the in- and outgoing * kinematics associated to each particle. * * \author Laurent Forthomme <laurent.forthomme@cern.ch> * \date Feb 2013 * \brief Core of the Monte-Carlo generator * */ class Generator { public: /// Core of the Monte Carlo integrator and events generator Generator(); /// Core of the Monte Carlo integrator and events generator /// \param[in] ip List of input parameters defining the phase space on which to perform the integration Generator( Parameters *ip ); ~Generator(); /// Dump this program's header into the standard output stream void printHeader(); /// Feed the generator with a Parameters object void setParameters( Parameters& ip ); /// Remove all references to a previous generation/run void clearRun(); /** * Compute the cross section for the run parameters defined by this object. * This returns the cross section as well as the absolute error computed along. * \brief Compute the cross-section for the given process * \param[out] xsec The computed cross-section, in pb * \param[out] err The absolute integration error on the computed cross-section, in pb */ void computeXsection( double& xsec, double& err ); /// Integrate the functional over the whole phase space void integrate(); /// Last cross section computed by the generator double crossSection() const { return result_; } /// Last error on the cross section computed by the generator double crossSectionError() const { return result_error_; } //void terminate(); /// Generate one single event given the phase space computed by Vegas in the integration step /// \return A pointer to the Event object generated in this run std::shared_ptr<Event> generateOneEvent(); /// Launch the generation of events void generate( std::function<void( const Event&, unsigned long )> callback = {} ); /// Number of dimensions on which the integration is performed size_t numDimensions() const; /// Compute one single point from the total phase space /// \param[in] x the n-dimensional point to compute /// \return the function value for the given point double computePoint( double* x ); /// Physical Parameters used in the events generation and cross-section computation std::unique_ptr<Parameters> parameters; private: /// Vegas instance which will integrate the function std::unique_ptr<Integrator> integrator_; /// Cross section value computed at the last integration double result_; /// Error on the cross section as computed in the last integration double result_error_; }; } #endif diff --git a/CepGen/Hadronisers/Pythia8Hadroniser.cpp b/CepGen/Hadronisers/Pythia8Hadroniser.cpp index ca6afb6..b188eda 100644 --- a/CepGen/Hadronisers/Pythia8Hadroniser.cpp +++ b/CepGen/Hadronisers/Pythia8Hadroniser.cpp @@ -1,464 +1,464 @@ #include "CepGen/Hadronisers/Pythia8Hadroniser.h" #include "CepGen/Parameters.h" #include "CepGen/Physics/Kinematics.h" #include "CepGen/Physics/Constants.h" #include "CepGen/Physics/PDG.h" #include "CepGen/Core/ParametersList.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/utils.h" #include "CepGen/Event/Event.h" #include "CepGen/Event/Particle.h" #include "CepGen/Version.h" namespace CepGen { namespace hadroniser { #ifdef PYTHIA8 Pythia8::Vec4 momToVec4( const Particle::Momentum& mom ) { return Pythia8::Vec4( mom.px(), mom.py(), mom.pz(), mom.energy() ); } #endif Pythia8Hadroniser::Pythia8Hadroniser( const Parameters& params, const ParametersList& plist ) : GenericHadroniser( "pythia8", plist ), #ifdef PYTHIA8 pythia_( new Pythia8::Pythia ), lhaevt_( new LHAEvent( ¶ms ) ), #endif full_evt_( false ), offset_( 0 ), first_evt_( true ), params_( ¶ms ) { #ifdef PYTHIA8 pythia_->setLHAupPtr( (Pythia8::LHAup*)lhaevt_.get() ); pythia_->settings.parm( "Beams:idA", (short)params.kinematics.incoming_beams.first.pdg ); pythia_->settings.parm( "Beams:idB", (short)params.kinematics.incoming_beams.second.pdg ); // specify we will be using a LHA input pythia_->settings.mode( "Beams:frameType", 5 ); pythia_->settings.parm( "Beams:eCM", params.kinematics.sqrtS() ); #endif for ( const auto& pdgid : params.kinematics.minimum_final_state ) min_ids_.emplace_back( (unsigned short)pdgid ); } Pythia8Hadroniser::~Pythia8Hadroniser() { #ifdef PYTHIA8 pythia_->settings.writeFile( "last_pythia_config.cmd", false ); #endif } void Pythia8Hadroniser::readString( const char* param ) { #ifdef PYTHIA8 if ( !pythia_->readString( param ) ) throw CG_FATAL( "Pythia8Hadroniser" ) << "The Pythia8 core failed to parse the following setting:\n\t" << param; #endif } void Pythia8Hadroniser::init() { #ifdef PYTHIA8 if ( pythia_->settings.flag( "ProcessLevel:all" ) != full_evt_ ) pythia_->settings.flag( "ProcessLevel:all", full_evt_ ); if ( seed_ == -1ll ) pythia_->settings.flag( "Random:setSeed", false ); else { pythia_->settings.flag( "Random:setSeed", true ); pythia_->settings.mode( "Random:seed", seed_ ); } switch ( params_->kinematics.mode ) { case KinematicsMode::ElasticElastic: { pythia_->settings.mode( "BeamRemnants:unresolvedHadron", 3 ); } break; case KinematicsMode::InelasticElastic: { pythia_->settings.mode( "BeamRemnants:unresolvedHadron", 2 ); } break; case KinematicsMode::ElasticInelastic: { pythia_->settings.mode( "BeamRemnants:unresolvedHadron", 1 ); } break; case KinematicsMode::InelasticInelastic: default: { pythia_->settings.mode( "BeamRemnants:unresolvedHadron", 0 ); } break; } if ( !pythia_->init() ) throw CG_FATAL( "Pythia8Hadroniser" ) << "Failed to initialise the Pythia8 core!\n\t" << "See the message above for more details."; #else throw CG_FATAL( "Pythia8Hadroniser" ) << "Pythia8 is not linked to this instance!"; #endif } void Pythia8Hadroniser::setCrossSection( double xsec, double xsec_err ) { #ifdef PYTHIA8 lhaevt_->setCrossSection( 0, xsec, xsec_err ); #endif } bool Pythia8Hadroniser::run( Event& ev, double& weight, bool full ) { //--- initialise the event weight before running any decay algorithm weight = 1.; #ifdef PYTHIA8 if ( !full && !pythia_->settings.flag( "ProcessLevel:resonanceDecays" ) ) return true; //--- switch full <-> partial event if ( full != full_evt_ ) { full_evt_ = full; init(); } //=========================================================================================== // convert our event into a custom LHA format //=========================================================================================== lhaevt_->feedEvent( ev, full, params_->kinematics.mode ); //if ( full ) lhaevt_->listEvent(); //=========================================================================================== // launch the hadronisation / resonances decays, and update the event accordingly //=========================================================================================== ev.num_hadronisation_trials = 0; while ( true ) { if ( ev.num_hadronisation_trials++ > max_trials_ ) return false; //--- run the hadronisation/fragmentation algorithm if ( pythia_->next() ) { //--- hadronisation successful if ( first_evt_ && full ) { offset_ = 0; for ( unsigned short i = 1; i < pythia_->event.size(); ++i ) if ( pythia_->event[i].status() == -12 ) // skip the incoming particles offset_++; first_evt_ = false; } break; } } //=========================================================================================== // update the event content with Pythia's output //=========================================================================================== updateEvent( ev, weight, full ); return true; #else throw CG_FATAL( "Pythia8Hadroniser" ) << "Pythia8 is not linked to this instance!"; #endif } #ifdef PYTHIA8 Particle& Pythia8Hadroniser::addParticle( Event& ev, const Pythia8::Particle& py_part, const Pythia8::Vec4& mom, unsigned short role ) const { Particle& op = ev.addParticle( (Particle::Role)role ); op.setPdgId( static_cast<PDG>( abs( py_part.id() ) ), py_part.charge() ); op.setStatus( py_part.isFinal() ? Particle::Status::FinalState : Particle::Status::Propagator ); op.setMomentum( Particle::Momentum( mom.px(), mom.py(), mom.pz(), mom.e() ) ); op.setMass( mom.mCalc() ); lhaevt_->addCorresp( py_part.index()-offset_, op.id() ); return op; } void Pythia8Hadroniser::updateEvent( Event& ev, double& weight, bool full ) const { for ( unsigned short i = 1+offset_; i < pythia_->event.size(); ++i ) { const Pythia8::Particle& p = pythia_->event[i]; const unsigned short cg_id = lhaevt_->cepgenId( i-offset_ ); if ( cg_id != LHAEvent::invalid_id ) { //----- particle already in the event Particle& cg_part = ev[cg_id]; //--- fragmentation result if ( cg_part.role() == Particle::OutgoingBeam1 || cg_part.role() == Particle::OutgoingBeam2 ) { cg_part.setStatus( Particle::Status::Fragmented ); continue; } //--- particle is not what we expect if ( p.idAbs() != abs( cg_part.integerPdgId() ) ) { CG_INFO( "Pythia8Hadroniser:update" ) << "LHAEVT event content:"; lhaevt_->listEvent(); CG_INFO( "Pythia8Hadroniser:update" ) << "Pythia event content:"; pythia_->event.list(); CG_INFO( "Pythia8Hadroniser:update" ) << "CepGen event content:"; ev.dump(); CG_INFO( "Pythia8Hadroniser:update" ) << "Correspondence:"; lhaevt_->dumpCorresp(); throw CG_FATAL( "Pythia8Hadroniser:update" ) << "Event list corruption detected for (Pythia/CepGen) particle " << i << "/" << cg_id << ":\n\t" << "should be " << abs( p.id() ) << ", " << "got " << cg_part.integerPdgId() << "!"; } //--- resonance decayed; apply branching ratio for this decay if ( p.particleDataEntry().sizeChannels() > 0 ) { weight *= p.particleDataEntry().pickChannel().bRatio(); cg_part.setStatus( Particle::Status::Resonance ); } } else { //----- new particle to be added const unsigned short role = findRole( ev, p ); switch ( (Particle::Role)role ) { default: break; case Particle::OutgoingBeam1: { ev.getByRole( Particle::OutgoingBeam1 )[0].setStatus( Particle::Status::Fragmented ); if ( abs( p.status() ) != 61 ) break; } // no break! case Particle::OutgoingBeam2: { ev.getByRole( Particle::OutgoingBeam2 )[0].setStatus( Particle::Status::Fragmented ); if ( abs( p.status() ) != 61 ) break; } // no break! } // found the role ; now we can add the particle Particle& cg_part = addParticle( ev, p, p.p(), role ); for ( const auto& moth_id : p.motherList() ) { if ( moth_id <= offset_ ) continue; const unsigned short moth_cg_id = lhaevt_->cepgenId( moth_id-offset_ ); if ( moth_cg_id != LHAEvent::invalid_id ) cg_part.addMother( ev[moth_cg_id] ); else cg_part.addMother( addParticle( ev, pythia_->event[moth_id], p.p(), role ) ); if ( !p.isFinal() ) { if ( p.isResonance() || p.daughterList().size() > 0 ) cg_part.setStatus( Particle::Status::Resonance ); else cg_part.setStatus( Particle::Status::Undefined ); } } } } } unsigned short Pythia8Hadroniser::findRole( const Event& ev, const Pythia8::Particle& p ) const { for ( const auto& par_id : p.motherList() ) { if ( par_id == 1 && offset_ > 0 ) return (unsigned short)Particle::OutgoingBeam1; if ( par_id == 2 && offset_ > 0 ) return (unsigned short)Particle::OutgoingBeam2; const unsigned short par_cg_id = lhaevt_->cepgenId( par_id-offset_ ); if ( par_cg_id != LHAEvent::invalid_id ) return (unsigned short)ev.at( par_cg_id ).role(); return findRole( ev, pythia_->event[par_id] ); } return (unsigned short)Particle::UnknownRole; } #endif } //================================================================================================ // Custom LHA event definition //================================================================================================ #ifdef PYTHIA8 - const double LHAEvent::mp_ = ParticleProperties::mass( PDG::proton ); + const double LHAEvent::mp_ = part::mass( PDG::proton ); const double LHAEvent::mp2_ = LHAEvent::mp_*LHAEvent::mp_; LHAEvent::LHAEvent( const Parameters* params ) : LHAup( 3 ), params_( params ) { addProcess( 0, 1., 1., 1.e3 ); if ( params_ ) { setBeamA( (short)params_->kinematics.incoming_beams.first.pdg, params_->kinematics.incoming_beams.first.pz ); setBeamB( (short)params_->kinematics.incoming_beams.second.pdg, params_->kinematics.incoming_beams.second.pz ); } } void LHAEvent::setCrossSection( int id, double xsec, double xsec_err ) { setXSec( id, xsec ); setXErr( id, xsec_err ); //listInit(); } void LHAEvent::feedEvent( const Event& ev, bool full, const KinematicsMode& mode ) { const double scale = ev.getOneByRole( Particle::Intermediate ).mass(); - setProcess( 0, 1., scale, Constants::alphaEM, Constants::alphaQCD ); + setProcess( 0, 1., scale, constants::alphaEM, constants::alphaQCD ); const Particle& part1 = ev.getOneByRole( Particle::Parton1 ), &part2 = ev.getOneByRole( Particle::Parton2 ); const Particle& op1 = ev.getOneByRole( Particle::OutgoingBeam1 ), &op2 = ev.getOneByRole( Particle::OutgoingBeam2 ); const double q2_1 = -part1.momentum().mass2(), q2_2 = -part2.momentum().mass2(); const double x1 = q2_1/( q2_1+op1.mass2()-mp2_ ), x2 = q2_2/( q2_2+op2.mass2()-mp2_ ); unsigned short quark1_id = 0, quark2_id = 0; unsigned short quark1_pdgid = part1.integerPdgId(), quark2_pdgid = part2.integerPdgId(); const Pythia8::Vec4 mom_part1( hadroniser::momToVec4( part1.momentum() ) ), mom_part2( hadroniser::momToVec4( part2.momentum() ) ); if ( !full ) { //============================================================================================= // incoming partons //============================================================================================= addCorresp( sizePart(), part1.id() ); addParticle( quark1_pdgid, -2, quark1_id, 0, 0, 0, mom_part1.px(), mom_part1.py(), mom_part1.pz(), mom_part1.e(), mom_part1.mCalc(), 0., 0. ); addCorresp( sizePart(), part2.id() ); addParticle( quark2_pdgid, -2, quark2_id, 0, 0, 0, mom_part2.px(), mom_part2.py(), mom_part2.pz(), mom_part2.e(), mom_part2.mCalc(), 0., 0. ); } else { // full event content (with collinear partons) const bool inel1 = ( mode == KinematicsMode::InelasticElastic || mode == KinematicsMode::InelasticInelastic ); const bool inel2 = ( mode == KinematicsMode::ElasticInelastic || mode == KinematicsMode::InelasticInelastic ); Pythia8::Vec4 mom_iq1 = mom_part1, mom_iq2 = mom_part2; unsigned short colour_index = 501, quark1_colour = 0, quark2_colour = 0; //FIXME select quark flavours accordingly if ( inel1 ) { quark1_pdgid = 2; quark1_colour = colour_index++; mom_iq1 = hadroniser::momToVec4( x1*ev.getOneByRole( Particle::IncomingBeam1 ).momentum() ); } if ( inel2 ) { quark2_pdgid = 2; quark2_colour = colour_index++; mom_iq2 = hadroniser::momToVec4( x2*ev.getOneByRole( Particle::IncomingBeam2 ).momentum() ); } //--- flavour / x value of hard-process initiators setIdX( part1.integerPdgId(), part2.integerPdgId(), x1, x2 ); //=========================================================================================== // incoming valence quarks //=========================================================================================== quark1_id = sizePart(); addCorresp( quark1_id, op1.id() ); addParticle( quark1_pdgid, -1, 0, 0, quark1_colour, 0, mom_iq1.px(), mom_iq1.py(), mom_iq1.pz(), mom_iq1.e(), mom_iq1.mCalc(), 0., 1. ); quark2_id = sizePart(); addCorresp( quark2_id, op2.id() ); addParticle( quark2_pdgid, -1, 0, 0, quark2_colour, 0, mom_iq2.px(), mom_iq2.py(), mom_iq2.pz(), mom_iq2.e(), mom_iq2.mCalc(), 0., 1. ); //=========================================================================================== // outgoing valence quarks //=========================================================================================== if ( inel1 ) { const Pythia8::Vec4 mom_oq1 = mom_iq1-mom_part1; addParticle( quark1_pdgid, 1, quark1_id, quark2_id, quark1_colour, 0, mom_oq1.px(), mom_oq1.py(), mom_oq1.pz(), mom_oq1.e(), mom_oq1.mCalc(), 0., 1. ); } if ( inel2 ) { const Pythia8::Vec4 mom_oq2 = mom_iq2-mom_part2; addParticle( quark2_pdgid, 1, quark1_id, quark2_id, quark2_colour, 0, mom_oq2.px(), mom_oq2.py(), mom_oq2.pz(), mom_oq2.e(), mom_oq2.mCalc(), 0., 1. ); } } //============================================================================================= // central system //============================================================================================= for ( const auto& p : ev.getByRole( Particle::CentralSystem ) ) { const auto mothers = p.mothers(); unsigned short moth1_id = 1, moth2_id = 2; if ( !full ) { moth1_id = moth2_id = 0; if ( mothers.size() > 0 ) { const unsigned short moth1_cg_id = *mothers.begin(); moth1_id = pythiaId( moth1_cg_id ); if ( moth1_id == invalid_id ) { const Particle& moth = ev.at( moth1_cg_id ); if ( moth.mothers().size() > 0 ) moth1_id = pythiaId( *moth.mothers().begin() ); if ( moth.mothers().size() > 1 ) moth2_id = pythiaId( *moth.mothers().rbegin() ); } if ( mothers.size() > 1 ) { const unsigned short moth2_cg_id = *mothers.rbegin(); moth2_id = pythiaId( moth2_cg_id ); if ( moth2_id == invalid_id ) { const Particle& moth = ev.at( moth2_cg_id ); moth.dump(); moth2_id = pythiaId( *moth.mothers().rbegin() ); } } } } const Pythia8::Vec4 mom_part( p.momentum().px(), p.momentum().py(), p.momentum().pz(), p.momentum().energy() ); addCorresp( sizePart(), p.id() ); addParticle( p.integerPdgId(), 1, moth1_id, moth2_id, 0, 0, mom_part.px(), mom_part.py(), mom_part.pz(), mom_part.e(), mom_part.mCalc(), 0., 0., 0. ); } setPdf( quark1_pdgid, quark2_pdgid, x1, x2, scale, 0., 0., false ); } bool LHAEvent::setInit() { return true; } bool LHAEvent::setEvent( int ) { return true; } void LHAEvent::setProcess( int id, double xsec, double q2_scale, double alpha_qed, double alpha_qcd ) { LHAup::setProcess( id, xsec, q2_scale, alpha_qed, alpha_qcd ); py_cg_corresp_.clear(); } unsigned short LHAEvent::cepgenId( unsigned short py_id ) const { for ( const auto& py_cg : py_cg_corresp_ ) if ( py_cg.first == py_id ) return py_cg.second; return invalid_id; } unsigned short LHAEvent::pythiaId( unsigned short cg_id ) const { for ( const auto& py_cg : py_cg_corresp_ ) if ( py_cg.second == cg_id ) return py_cg.first; return invalid_id; } void LHAEvent::addCorresp( unsigned short py_id, unsigned short cg_id ) { py_cg_corresp_.emplace_back( py_id, cg_id ); } void LHAEvent::dumpCorresp() const { std::ostringstream oss; oss << "List of Pythia <-> CepGen particle ids correspondance"; for ( const auto& py_cg : py_cg_corresp_ ) oss << "\n\t" << py_cg.first << " <-> " << py_cg.second; CG_INFO( "LHAEvent:dump" ) << oss.str(); } #endif } diff --git a/CepGen/IO/ExportHandler.h b/CepGen/IO/ExportHandler.h index dcf77a8..20e9422 100644 --- a/CepGen/IO/ExportHandler.h +++ b/CepGen/IO/ExportHandler.h @@ -1,59 +1,59 @@ #ifndef CepGen_Export_ExportHandler_h #define CepGen_Export_ExportHandler_h #include <iostream> namespace CepGen { class Event; class Parameters; /// Location for all output generators - namespace OutputHandler + namespace output { /** * \brief Output format handler for events export * \author Laurent Forthomme <laurent.forthomme@cern.ch> * \date Sep 2016 */ class ExportHandler { public: /// All types of output available for export enum OutputType { HepMC, ///< HepMC ASCII format LHE ///< LHEF format }; /// Human-readable output name friend std::ostream& operator<<( std::ostream& os, const OutputType& type ) { switch ( type ) { case HepMC: return os << "HepMC ASCII"; case LHE: return os << "LHEF"; } return os; } public: /// \brief Class constructor /// \param[in] type Requested output type explicit ExportHandler( const OutputType& type ) : type_( type ), event_num_( 0. ) {} virtual ~ExportHandler() {} /// Initialise the handler and its inner parameterisation virtual void initialise( const Parameters& ) = 0; /// Set the process cross section and its associated error virtual void setCrossSection( double xsec, double err_xsec ) {} /// Set the event number void setEventNumber( const unsigned int& ev_id ) { event_num_ = ev_id; } /// Writer operator virtual void operator<<( const Event& ) = 0; protected: /// Type of output requested OutputType type_; /// Event index unsigned int event_num_; }; } } #endif diff --git a/CepGen/IO/HepMCHandler.cpp b/CepGen/IO/HepMCHandler.cpp index e8b1948..c9d7990 100644 --- a/CepGen/IO/HepMCHandler.cpp +++ b/CepGen/IO/HepMCHandler.cpp @@ -1,134 +1,133 @@ #include "CepGen/IO/HepMCHandler.h" #include "CepGen/Parameters.h" #include "CepGen/Core/Exception.h" #include "CepGen/Event/Event.h" #include "CepGen/Physics/Constants.h" #ifdef LIBHEPMC # include "HepMC/GenVertex.h" # include "HepMC/GenParticle.h" #endif namespace CepGen { - namespace OutputHandler + namespace output { HepMCHandler::HepMCHandler( const char* filename, const ExportHandler::OutputType& type ) : ExportHandler( type ) #ifdef LIBHEPMC # ifdef HEPMC_VERSION3 , output_( new HepMC::WriterAscii( filename ) ), # else , output_( new HepMC::IO_GenEvent( filename ) ), # endif event_( new HepMC::GenEvent() ) #endif {} void HepMCHandler::operator<<( const Event& evt ) { fillEvent( evt ); #ifdef LIBHEPMC # ifdef HEPMC_VERSION3 output_->write_event( *event_ ); # else output_->write_event( event_.get() ); # endif #endif } void HepMCHandler::setCrossSection( double xsect, double xsect_err ) { #ifdef LIBHEPMC # ifdef HEPMC_VERSION3 xs_->set_cross_section( xsect, xsect_err ); - event_->add_attribute( "AlphaQCD", HepMC::make_shared<HepMC::DoubleAttribute>( Constants::alphaQCD ) ); - event_->add_attribute( "AlphaEM", HepMC::make_shared<HepMC::DoubleAttribute>( Constants::alphaEM ) ); + event_->add_attribute( "AlphaQCD", HepMC::make_shared<HepMC::DoubleAttribute>( constants::alphaQCD ) ); + event_->add_attribute( "AlphaEM", HepMC::make_shared<HepMC::DoubleAttribute>( constants::alphaEM ) ); # else xs_.set_cross_section( xsect, xsect_err ); - event_->set_alphaQCD( Constants::alphaQCD ); - event_->set_alphaQED( Constants::alphaEM ); + event_->set_alphaQCD( constants::alphaQCD ); + event_->set_alphaQED( constants::alphaEM ); # endif #endif } void HepMCHandler::fillEvent( const Event& evt ) { #ifdef LIBHEPMC event_->clear(); // general information event_->set_cross_section( xs_ ); event_->set_event_number( event_num_ ); event_->weights().push_back( 1. ); // unweighted events // filling the particles content const HepMC::FourVector origin( 0., 0., 0., 0. ); Particles part_vec = evt.particles(); int cm_id = 0, idx = 1; # ifdef HEPMC_VERSION3 HepMC::GenVertexPtr v1 = HepMC::make_shared<HepMC::GenVertex>( origin ), v2 = HepMC::make_shared<HepMC::GenVertex>( origin ), vcm = HepMC::make_shared<HepMC::GenVertex>( origin ); # else HepMC::GenVertex* v1 = new HepMC::GenVertex( origin ), *v2 = new HepMC::GenVertex( origin ), *vcm = new HepMC::GenVertex( origin ); # endif for ( unsigned int i = 0; i < part_vec.size(); ++i ) { const Particle part_orig = part_vec.at( i ); HepMC::FourVector pmom( part_orig.momentum().px(), part_orig.momentum().py(), part_orig.momentum().pz(), part_orig.energy() ); # ifdef HEPMC_VERSION3 HepMC::GenParticlePtr part = HepMC::make_shared<HepMC::GenParticle>( pmom, part_orig.integerPdgId(), (int)part_orig.status() ); # else HepMC::GenParticle* part = new HepMC::GenParticle( pmom, part_orig.integerPdgId(), (int)part_orig.status() ); part->suggest_barcode( idx++ ); # endif const ParticlesIds moth = part_orig.mothers(); switch ( part_orig.role() ) { case Particle::IncomingBeam1: { v1->add_particle_in( part ); } break; case Particle::IncomingBeam2: { v2->add_particle_in( part ); } break; case Particle::OutgoingBeam1: { v1->add_particle_out( part ); } break; case Particle::OutgoingBeam2: { v2->add_particle_out( part ); } break; case Particle::Parton1: { v1->add_particle_out( part ); vcm->add_particle_in( part ); } break; case Particle::Parton2: { v2->add_particle_out( part ); vcm->add_particle_in( part ); } break; case Particle::Parton3: { v2->add_particle_out( part ); vcm->add_particle_in( part ); } break; case Particle::Intermediate: { cm_id = i; continue; } break; case Particle::CentralSystem: default: { if ( moth.size() == 0 ) continue; if ( *moth.begin() == cm_id ) vcm->add_particle_out( part ); else throw CG_FATAL( "HepMCHandler:fillEvent" ) << "Other particle requested! Not yet implemented!"; } break; } idx++; } event_->add_vertex( v1 ); event_->add_vertex( v2 ); event_->add_vertex( vcm ); # ifndef HEPMC_VERSION3 event_->set_beam_particles( *v1->particles_in_const_begin(), *v2->particles_in_const_begin() ); event_->set_signal_process_vertex( *v1->vertices_begin() ); event_->set_beam_particles( *v1->particles_in_const_begin(), *v2->particles_in_const_end() ); # endif #endif event_num_++; } } } - diff --git a/CepGen/IO/HepMCHandler.h b/CepGen/IO/HepMCHandler.h index e783477..0d739a6 100644 --- a/CepGen/IO/HepMCHandler.h +++ b/CepGen/IO/HepMCHandler.h @@ -1,67 +1,66 @@ #ifndef CepGen_IO_HepMCHandler_h #define CepGen_IO_HepMCHandler_h #include "CepGen/IO/ExportHandler.h" #ifdef LIBHEPMC # include "HepMC/Version.h" # ifndef HEPMC_VERSION_CODE // HepMC v2 # include "HepMC/IO_GenEvent.h" # include "HepMC/SimpleVector.h" # else // HepMC v3+ # define HEPMC_VERSION3 # include "HepMC/WriterAscii.h" # include "HepMC/FourVector.h" # endif # include "HepMC/GenEvent.h" #endif #include <memory> namespace CepGen { - namespace OutputHandler + namespace output { /** * \brief Handler for the HepMC file output * \author Laurent Forthomme <laurent.forthomme@cern.ch> * \date Sep 2016 */ class HepMCHandler : public ExportHandler { public: /// Class constructor /// \param[in] filename Output file path /// \param[in] type Output type HepMCHandler( const char* filename, const ExportHandler::OutputType& type = ExportHandler::HepMC ); void initialise( const Parameters& params ) override {} /// Writer operator void operator<<( const Event& ) override; void setCrossSection( double, double ) override; private: /// Clear the associated HepMC event content void clearEvent(); /// Populate the associated HepMC event with a Event object void fillEvent( const Event& ); #ifdef LIBHEPMC # ifdef HEPMC_VERSION3 /// Writer object (from HepMC v3+) std::unique_ptr<HepMC::WriterAscii> output_; /// Generator cross section and error HepMC::GenCrossSectionPtr xs_; # else /// Writer object (from HepMC v<3) std::unique_ptr<HepMC::IO_GenEvent> output_; /// Generator cross section and error HepMC::GenCrossSection xs_; # endif /// Associated HepMC event std::shared_ptr<HepMC::GenEvent> event_; #endif }; } } #endif - diff --git a/CepGen/IO/LHEFHandler.cpp b/CepGen/IO/LHEFHandler.cpp index 7429b69..19f44e3 100644 --- a/CepGen/IO/LHEFHandler.cpp +++ b/CepGen/IO/LHEFHandler.cpp @@ -1,252 +1,252 @@ #include "CepGen/IO/LHEFHandler.h" #include "CepGen/StructureFunctions/StructureFunctions.h" #include "CepGen/Event/Event.h" #include "CepGen/Physics/Constants.h" #include "CepGen/Parameters.h" #include "CepGen/Version.h" namespace CepGen { - namespace OutputHandler + namespace output { LHEFHandler::LHEFHandler( const char* filename ) : ExportHandler( ExportHandler::LHE ) #if defined ( HEPMC_LHEF ) , lhe_output_( new LHEF::Writer( filename ) ) #elif defined ( PYTHIA_LHEF ) , pythia_( new Pythia8::Pythia ), lhaevt_( new LHAevent ) #endif { #if defined ( PYTHIA_LHEF ) lhaevt_->openLHEF( filename ); #endif } LHEFHandler::~LHEFHandler() { #if defined ( PYTHIA_LHEF ) if ( lhaevt_ ) lhaevt_->closeLHEF( false ); // we do not want to rewrite the init block #endif } void LHEFHandler::initialise( const Parameters& params ) { std::ostringstream oss_init; oss_init << "<!--\n" << " ***** Sample generated with CepGen v" << version() << " *****\n" << " * process: " << params.processName() << " (" << params.kinematics.mode << ")\n"; if ( params.kinematics.mode != KinematicsMode::ElasticElastic ) { oss_init << " * structure functions: " << params.kinematics.structure_functions->type << "\n"; if ( !params.hadroniserName().empty() ) oss_init << " * hadroniser: " << params.hadroniserName() << "\n"; } oss_init << " *--- incoming state\n"; if ( params.kinematics.cuts.initial.q2.valid() ) oss_init << " * Q2 range (GeV2): " << params.kinematics.cuts.initial.q2.min() << ", " << params.kinematics.cuts.initial.q2.max() << "\n"; if ( params.kinematics.mode != KinematicsMode::ElasticElastic && params.kinematics.cuts.remnants.mass_single.valid() ) oss_init << " * remnants mass range (GeV/c2): " << params.kinematics.cuts.remnants.mass_single.min() << ", " << params.kinematics.cuts.remnants.mass_single.max() << "\n"; oss_init << " *--- central system\n"; if ( params.kinematics.cuts.central.pt_single.valid() ) oss_init << " * single particle pt (GeV/c): " << params.kinematics.cuts.central.pt_single.min() << ", " << params.kinematics.cuts.central.pt_single.max() << "\n"; if ( params.kinematics.cuts.central.energy_single.valid() ) oss_init << " * single particle energy (GeV): " << params.kinematics.cuts.central.energy_single.min() << ", " << params.kinematics.cuts.central.energy_single.max() << "\n"; if ( params.kinematics.cuts.central.eta_single.valid() ) oss_init << " * single particle eta: " << params.kinematics.cuts.central.eta_single.min() << ", " << params.kinematics.cuts.central.eta_single.max() << "\n"; if ( params.kinematics.cuts.central.pt_sum.valid() ) oss_init << " * total pt (GeV/c): " << params.kinematics.cuts.central.mass_sum.min() << ", " << params.kinematics.cuts.central.mass_sum.max() << "\n"; if ( params.kinematics.cuts.central.mass_sum.valid() ) oss_init << " * total invariant mass (GeV/c2): " << params.kinematics.cuts.central.mass_sum.min() << ", " << params.kinematics.cuts.central.mass_sum.max() << "\n"; oss_init << " **************************************************\n" << "-->"; #if defined ( HEPMC_LHEF ) lhe_output_->headerBlock() << oss_init.str(); //params.dump( lhe_output_->initComments(), false ); LHEF::HEPRUP run = lhe_output_->heprup; run.IDBMUP = { (int)params.kinematics.incoming_beams.first.pdg, (int)params.kinematics.incoming_beams.second.pdg }; run.EBMUP = { (double)params.kinematics.incoming_beams.first.pz, (double)params.kinematics.incoming_beams.second.pz }; run.NPRUP = 1; run.resize(); run.XSECUP[0] = params.integrator.result; run.XERRUP[0] = params.integrator.err_result; run.XMAXUP[0] = 1.; run.LPRUP[0] = 1; lhe_output_->heprup = run; lhe_output_->init(); #elif defined ( PYTHIA_LHEF ) oss_init << std::endl; // LHEF is usually not beautifully parsed as a standard XML... lhaevt_->addComments( oss_init.str() ); lhaevt_->initialise( params ); pythia_->settings.mode( "Beams:frameType", 5 ); pythia_->settings.mode( "Next:numberCount", 0 ); // remove some of the Pythia output pythia_->settings.flag( "ProcessLevel:all", false ); // we do not want Pythia to interfere... pythia_->setLHAupPtr( lhaevt_.get() ); pythia_->init(); lhaevt_->initLHEF(); #endif } void LHEFHandler::operator<<( const Event& ev ) { #if defined ( HEPMC_LHEF ) LHEF::HEPEUP out; out.heprup = &lhe_output_->heprup; out.XWGTUP = 1.; out.XPDWUP = std::pair<double,double>( 0., 0. ); out.SCALUP = 0.; - out.AQEDUP = Constants::alphaEM; - out.AQCDUP = Constants::alphaQCD; + out.AQEDUP = constants::alphaEM; + out.AQCDUP = constants::alphaQCD; out.NUP = ev.numParticles(); out.resize(); for ( unsigned short ip = 0; ip < ev.numParticles(); ++ip ) { const Particle part = ev.at( ip ); out.IDUP[ip] = part.integerPdgId(); // PDG id out.ISTUP[ip] = (short)part.status(); // status code out.MOTHUP[ip] = std::pair<int,int>( ( part.mothers().size() > 0 ) ? *part.mothers().begin()+1 : 0, ( part.mothers().size() > 1 ) ? *part.mothers().rbegin()+1 : 0 ); // mothers out.ICOLUP[ip] = std::pair<int,int>( 0, 0 ); out.PUP[ip] = std::vector<double>( { { part.momentum().px(), part.momentum().py(), part.momentum().pz(), part.energy(), part.mass() } } ); // momentum out.VTIMUP[ip] = 0.; // invariant lifetime out.SPINUP[ip] = 0.; } lhe_output_->eventComments() << "haha"; lhe_output_->hepeup = out; lhe_output_->writeEvent(); #elif defined ( PYTHIA_LHEF ) lhaevt_->feedEvent( 0, ev ); pythia_->next(); lhaevt_->eventLHEF(); #endif } void LHEFHandler::setCrossSection( double xsect, double xsect_err ) { #if defined ( PYTHIA_LHEF ) lhaevt_->setCrossSection( 0, xsect, xsect_err ); #endif } //--------------------------------------------------------------------------------------------- // Define LHA event record if one uses Pythia to store the LHE //--------------------------------------------------------------------------------------------- #if defined ( PYTHIA_LHEF ) LHEFHandler::LHAevent::LHAevent() : LHAup( 3 ) {} void LHEFHandler::LHAevent::initialise( const Parameters& params ) { setBeamA( (short)params.kinematics.incoming_beams.first.pdg, params.kinematics.incoming_beams.first.pz ); setBeamB( (short)params.kinematics.incoming_beams.second.pdg, params.kinematics.incoming_beams.second.pz ); addProcess( 0, params.integrator.result, params.integrator.err_result, 100. ); } void LHEFHandler::LHAevent::addComments( const std::string& comments ) { osLHEF << comments; } void LHEFHandler::LHAevent::setCrossSection( unsigned short proc_id, double xsect, double xsect_err ) { setXSec( proc_id, xsect ); setXErr( proc_id, xsect_err ); } void LHEFHandler::LHAevent::feedEvent( unsigned short proc_id, const Event& ev, bool full_event ) { const double scale = ev.getOneByRole( Particle::Intermediate ).mass(); - setProcess( proc_id, 1., scale, Constants::alphaEM, Constants::alphaQCD ); + setProcess( proc_id, 1., scale, constants::alphaEM, constants::alphaQCD ); const Particle& ip1 = ev.getOneByRole( Particle::IncomingBeam1 ), &ip2 = ev.getOneByRole( Particle::IncomingBeam2 ); const Particles& op1 = ev.getByRole( Particle::OutgoingBeam1 ), &op2 = ev.getByRole( Particle::OutgoingBeam2 ); const double q2_1 = -( ip1.momentum()-op1[0].momentum() ).mass2(), q2_2 = -( ip2.momentum()-op2[0].momentum() ).mass2(); const double x1 = q2_1/( q2_1+op1[0].mass2()-ip1.mass2() ), x2 = q2_2/( q2_2+op2[0].mass2()-ip2.mass2() ); setIdX( ip1.integerPdgId(), ip2.integerPdgId(), x1, x2 ); short parton1_pdgid = 0, parton2_pdgid = 0; for ( const auto& part : ev.particles() ) { short pdg_id = part.integerPdgId(), status = 0, moth1 = 0, moth2 = 0; switch ( part.role() ) { case Particle::Parton1: case Particle::Parton2: { if ( part.role() == Particle::Parton1 ) parton1_pdgid = part.integerPdgId(); if ( part.role() == Particle::Parton2 ) parton2_pdgid = part.integerPdgId(); if ( !full_event ) continue; status = -2; // conserving xbj/Q2 } break; case Particle::Intermediate: { if ( !full_event ) continue; status = 2; if ( pdg_id == 0 ) pdg_id = ev.at( *part.mothers().begin() ).integerPdgId(); } break; case Particle::IncomingBeam1: case Particle::IncomingBeam2: { if ( !full_event ) continue; status = -9; } break; case Particle::OutgoingBeam1: case Particle::OutgoingBeam2: case Particle::CentralSystem: { status = (short)part.status(); if ( status != 1 ) continue; } break; default: break; } if ( full_event ) { const auto& mothers = part.mothers(); if ( mothers.size() > 0 ) moth1 = *mothers.begin()+1; if ( mothers.size() > 1 ) moth2 = *mothers.rbegin()+1; } const Particle::Momentum& mom = part.momentum(); addParticle( pdg_id, status, moth1, moth2, 0, 0, mom.px(), mom.py(), mom.pz(), mom.energy(), mom.mass(), 0. ,0., 0. ); } setPdf( parton1_pdgid, parton2_pdgid, x1, x2, scale, 0., 0., true ); } #endif } } diff --git a/CepGen/IO/LHEFHandler.h b/CepGen/IO/LHEFHandler.h index 5e6eb1c..5acf548 100644 --- a/CepGen/IO/LHEFHandler.h +++ b/CepGen/IO/LHEFHandler.h @@ -1,71 +1,70 @@ #ifndef CepGen_Export_LHEFHandler_h #define CepGen_Export_LHEFHandler_h #include "CepGen/IO/ExportHandler.h" #include "CepGen/IO/HepMCHandler.h" #ifndef LIBHEPMC # ifndef PYTHIA8 # pragma message( "HepMC/Pythia8 are not linked to this instance!" ) # endif #else # ifndef HEPMC_VERSION3 # ifdef PYTHIA8 # include "Pythia8/Pythia.h" # define PYTHIA_LHEF 1 # else # pragma message( "HepMC v3 or Pythia8 are required for the LHEF export!" ) # endif # else # include "HepMC/LHEF.h" # define HEPMC_LHEF 1 # endif #endif namespace CepGen { class Event; - namespace OutputHandler + namespace output { /** * \brief Handler for the LHE file output * \author Laurent Forthomme <laurent.forthomme@cern.ch> * \date Sep 2016 */ class LHEFHandler : public ExportHandler { public: /// Class constructor /// \param[in] filename Output file path explicit LHEFHandler( const char* filename ); ~LHEFHandler() override; void initialise( const Parameters& ) override; /// Writer operator void operator<<( const Event& ) override; void setCrossSection( double, double ) override; private: #if defined ( HEPMC_LHEF ) /// Writer object (from HepMC) std::unique_ptr<LHEF::Writer> lhe_output_; LHEF::HEPRUP run_; #elif defined ( PYTHIA_LHEF ) std::unique_ptr<Pythia8::Pythia> pythia_; struct LHAevent : Pythia8::LHAup { explicit LHAevent(); void initialise( const Parameters& ); void setCrossSection( unsigned short, double, double ); void feedEvent( unsigned short proc_id, const Event& ev, bool full_event = false ); bool setInit() override { return true; } bool setEvent( int ) override { return true; } void addComments( const std::string& ); }; std::unique_ptr<LHAevent> lhaevt_; #endif }; } } #endif - diff --git a/CepGen/Physics/Constants.h b/CepGen/Physics/Constants.h index 0c971c6..5182da7 100644 --- a/CepGen/Physics/Constants.h +++ b/CepGen/Physics/Constants.h @@ -1,22 +1,21 @@ #ifndef CepGen_Physics_Constants_h #define CepGen_Physics_Constants_h #include <math.h> namespace CepGen { /// List of physical constants useful that may be used for the matrix element definition - namespace Constants + namespace constants { /// Electromagnetic coupling constant \f$\alpha_{\rm em}=\frac{e^2}{4\pi\epsilon_0\hbar c}\f$ constexpr double alphaEM = 1./137.035; /// Strong coupling constant \f$\alpha_{\rm QCD}\f$ constexpr double alphaQCD = 0.1184; // at the Z pole /// Conversion factor between GeV\f${}^2\f$ and barn constexpr double GeV2toBarn = 0.389351824e9; // 1.e4*(197.3271**2); constexpr double sconstb = 2.1868465e10; // 1.1868465e10; } } #endif - diff --git a/CepGen/Physics/FormFactors.cpp b/CepGen/Physics/FormFactors.cpp index 1a1d42b..34689b7 100644 --- a/CepGen/Physics/FormFactors.cpp +++ b/CepGen/Physics/FormFactors.cpp @@ -1,55 +1,55 @@ #include "CepGen/Physics/FormFactors.h" #include "CepGen/Core/Exception.h" #include "CepGen/Physics/ParticleProperties.h" #include "CepGen/Physics/PDG.h" #include "CepGen/StructureFunctions/SuriYennie.h" namespace CepGen { - const double FormFactors::mp_ = ParticleProperties::mass( PDG::proton ); + const double FormFactors::mp_ = part::mass( PDG::proton ); const double FormFactors::mp2_ = FormFactors::mp_*FormFactors::mp_; FormFactors FormFactors::trivial() { return FormFactors( 1.0, 1.0 ); } FormFactors FormFactors::protonElastic( double q2 ) { const double GE = pow( 1.+q2/0.71, -2. ), GE2 = GE*GE; const double GM = 2.79*GE, GM2 = GM*GM; return FormFactors( ( 4.*mp2_*GE2 + q2*GM2 ) / ( 4.*mp2_ + q2 ), GM2 ); } FormFactors FormFactors::protonInelastic( double q2, double mi2, double mf2, sf::Parameterisation& sf ) { const double xbj = q2 / ( q2 + mf2 - mi2 ); switch ( sf.type ) { case sf::Type::ElasticProton: CG_WARNING( "FormFactors" ) << "Elastic proton form factors requested! Check your process definition!"; return FormFactors::protonElastic( q2 ); case sf::Type::SuriYennie: { sf::SuriYennie suriyennie, sy = (sf::SuriYennie)suriyennie( xbj, q2 ); return FormFactors( sy.F2 * xbj * sqrt( mi2 ) / q2, sy.FM ); //FIXME } break; default: { sf = sf( xbj, q2 ); sf.computeFL( xbj, q2 ); return FormFactors( sf.F2 * xbj / q2, -2.*sf.F1( xbj, q2 ) / q2 ); } break; } } std::ostream& operator<<( std::ostream& os, const FormFactors& ff ) { os << Form( "Form factors: electric: Fe = %.3e ; magnetic: Fm = %.3e", ff.FE, ff.FM ).c_str(); return os; } } diff --git a/CepGen/Physics/KTFlux.cpp b/CepGen/Physics/KTFlux.cpp index c2fe582..fcfef07 100644 --- a/CepGen/Physics/KTFlux.cpp +++ b/CepGen/Physics/KTFlux.cpp @@ -1,103 +1,103 @@ #include "CepGen/Physics/KTFlux.h" #include "CepGen/Physics/FormFactors.h" #include "CepGen/Physics/ParticleProperties.h" #include "CepGen/Physics/HeavyIon.h" #include "CepGen/Physics/GluonGrid.h" #include "CepGen/StructureFunctions/StructureFunctions.h" #include "CepGen/Core/Exception.h" namespace CepGen { const double KTFluxParameters::kMinKTFlux = 1.e-20; - const double KTFluxParameters::kMP = ParticleProperties::mass( PDG::proton ); + const double KTFluxParameters::kMP = part::mass( PDG::proton ); const double KTFluxParameters::kMP2 = KTFluxParameters::kMP*KTFluxParameters::kMP; double ktFlux( const KTFlux& type, double x, double kt2, sf::Parameterisation& sf, double mx ) { double flux = 0.; const double mp2 = KTFluxParameters::kMP2; switch ( type ) { case KTFlux::P_Photon_Elastic: { const double x2 = x*x; const double q2min = x2*mp2/( 1.-x ), q2 = q2min + kt2/( 1.-x ); // electromagnetic form factors const auto& ff = FormFactors::protonElastic( q2 ); const double ela1 = ( 1.-x )*( 1.-q2min/q2 ); //const double ela3 = 1.-( q2-kt2 )/q2; - flux = Constants::alphaEM*M_1_PI*( 1.-x )/q2*( ela1*ff.FE + 0.5*x2*ff.FM ); + flux = constants::alphaEM*M_1_PI*( 1.-x )/q2*( ela1*ff.FE + 0.5*x2*ff.FM ); } break; case KTFlux::P_Photon_Inelastic_Budnev: { const double mx2 = mx*mx, x2 = x*x; const double q2min = ( x*( mx2-mp2 ) + x2*mp2 )/( 1.-x ), q2 = q2min + kt2/( 1.-x ); const double xbj = q2 / ( q2+mx2-mp2 ); // structure functions auto& str_fun = sf( xbj, q2 ); str_fun.computeFL( xbj, q2 ); const double f_D = str_fun.F2/( q2+mx2-mp2 )* ( 1.-x )*( 1.-q2min/q2 ); const double f_C = str_fun.F1( xbj, q2 ) * 2./q2; - flux = Constants::alphaEM*M_1_PI*( 1.-x )/q2*( f_D+0.5*x2*f_C ); + flux = constants::alphaEM*M_1_PI*( 1.-x )/q2*( f_D+0.5*x2*f_C ); } break; case KTFlux::P_Gluon_KMR: { flux = kmr::GluonGrid::get()( log10( x ), log10( kt2 ), 2.*log10( mx ) ); } break; default: throw CG_FATAL( "GenericKTProcess:flux" ) << "Invalid flux type: " << type; } if ( flux < KTFluxParameters::kMinKTFlux ) return 0.; return flux; } double ktFlux( const KTFlux& type, double x, double kt2, const HeavyIon& hi ) { double flux = 0.; switch ( type ) { case KTFlux::HI_Photon_Elastic: { const double r_a = 1.1*std::pow( hi.A, 1./3 ), a0 = 0.7, m_a = hi.A*KTFluxParameters::kMP; const double q2_ela = ( kt2+x*x*m_a*m_a )/( 1.-x ), cons = sqrt( q2_ela )/0.1973; const double tau = cons*r_a, tau1 = cons*a0; // "Realistic nuclear form-factor" as used in STARLIGHT const double ff1 = 3.*( sin( tau )-tau*cos( tau ) )/pow( tau+1.e-10, 3 ); const double ff2 = 1./( 1.+tau1*tau1 ); const double ela1 = pow( kt2/( kt2+x*x*m_a*m_a ), 2 ); const double ela2 = pow( ff1*ff2, 2 )/*, ela3 = 1.-( q2_ela-kt2 )/q2_ela*/; const unsigned int z = (unsigned short)hi.Z; - flux = Constants::alphaEM*M_1_PI*z*z*ela1*ela2/q2_ela; + flux = constants::alphaEM*M_1_PI*z*z*ela1*ela2/q2_ela; } break; default: throw CG_FATAL("GenericKTProcess:flux") << "Invalid flux type: " << type; } if ( flux < KTFluxParameters::kMinKTFlux ) return 0.; return flux; } std::ostream& operator<<( std::ostream& os, const KTFlux& type ) { switch ( type ) { case KTFlux::P_Photon_Elastic: return os << "elastic photon from proton"; case KTFlux::P_Photon_Inelastic: return os << "inelastic photon from proton"; case KTFlux::P_Photon_Inelastic_Budnev: return os << "inelastic photon from proton (Budnev)"; case KTFlux::P_Gluon_KMR: return os << "elastic gluon from proton (KMR)"; case KTFlux::HI_Photon_Elastic: return os << "elastic photon from HI"; case KTFlux::invalid: default: return os << "unrecognized flux (" << (int)type << ")"; } } } diff --git a/CepGen/Physics/ParticleProperties.cpp b/CepGen/Physics/ParticleProperties.cpp index 9560e6a..fd4a1a7 100644 --- a/CepGen/Physics/ParticleProperties.cpp +++ b/CepGen/Physics/ParticleProperties.cpp @@ -1,154 +1,153 @@ #include "CepGen/Physics/ParticleProperties.h" #include "CepGen/Physics/PDG.h" #include "CepGen/Core/Exception.h" namespace CepGen { - namespace ParticleProperties + namespace part { double mass( const PDG& pdg_id ) { switch ( pdg_id ) { case PDG::electron: return 0.510998928e-3; case PDG::muon: return 0.1056583715; case PDG::tau: return 1.77682; case PDG::down: return 0.0048; case PDG::up: return 0.0023; case PDG::strange: return 0.095; case PDG::charm: return 1.29; case PDG::bottom: return 4.18; case PDG::top: return 172.44; case PDG::electronNeutrino: case PDG::muonNeutrino: case PDG::tauNeutrino: return 0.; case PDG::gluon: case PDG::photon: case PDG::pomeron: return 0.; case PDG::Z: return 91.1876; case PDG::W: return 80.385; case PDG::piPlus: return 0.13957018; case PDG::piZero: return 0.1349766; case PDG::KPlus: return 0.49368; case PDG::DPlus: return 1.86962; case PDG::Jpsi: return 3.0969; case PDG::proton: return 0.938272046; case PDG::neutron: return 0.939565346; case PDG::Upsilon1S: return 9.46030; case PDG::Upsilon2S: return 10.02326; case PDG::Upsilon3S: return 10.3552; case PDG::rho770_0: return 0.77526; case PDG::rho1450_0: return 1.465; case PDG::rho1700_0: return 1.720; case PDG::h1380_1: return 1.38619; case PDG::eta: return 0.547862; case PDG::invalid: default: return -1.; } } double charge( const PDG& pdg_id ) { switch ( pdg_id ) { case PDG::proton: case PDG::diffractiveProton: return +1.; case PDG::electron: case PDG::muon: case PDG::tau: return -1.; case PDG::down: case PDG::strange: case PDG::bottom: return -1./3; case PDG::up: case PDG::charm: case PDG::top: return +2./3; case PDG::W: return +1.; case PDG::piPlus: case PDG::KPlus: case PDG::DPlus: return +1.; default: return 0.; } } double charge( int id ) { const short sign = id / abs( id ); return sign * charge( (PDG)abs( id ) ); } unsigned short colours( const PDG& pdg_id ) { switch ( pdg_id ) { case PDG::top: return 3; default: return 1; } } double width( const PDG& pdg_id ) { switch ( pdg_id ) { case PDG::Jpsi: return 92.9e-6; //FIXME case PDG::Z: return 2.4952; case PDG::W: return 2.085; case PDG::Upsilon1S: return 54.02e-6; case PDG::Upsilon2S: return 31.98e-6; case PDG::Upsilon3S: return 20.32e-6; case PDG::rho770_0: return 0.150; // PDG case PDG::rho1450_0: return 0.400; // PDG case PDG::rho1700_0: return 0.250; // PDG default: CG_WARNING( "ParticleProperties:width" ) << "Particle " << pdg_id << " has no registered width."; return -1.; } } } std::ostream& operator<<( std::ostream& os, const PDG& pc ) { switch ( pc ) { case PDG::electron: return os << "e± "; case PDG::electronNeutrino: return os << "ν_e "; case PDG::muon: return os << "µ± "; case PDG::muonNeutrino: return os << "ν_µ "; case PDG::tau: return os << "τ± "; case PDG::tauNeutrino: return os << "ν_τ "; case PDG::gluon: return os << "gluon"; case PDG::photon: return os << "ɣ "; case PDG::Z: return os << "Z"; case PDG::W: return os << "W± "; case PDG::piPlus: return os << "π± "; case PDG::piZero: return os << "π⁰\t"; case PDG::KPlus: return os << "K± "; case PDG::DPlus: return os << "D± "; case PDG::rho770_0: return os << "ρ(770)₀ "; case PDG::rho1450_0: return os << "ρ(1450)₀ "; case PDG::rho1700_0: return os << "ρ(1700)₀ "; case PDG::h1380_1: return os << "h(1380)₁ "; case PDG::eta: return os << "η meson"; case PDG::omega782: return os << "ω(782) "; case PDG::Jpsi: return os << "J/ψ "; case PDG::phi1680: return os << "ɸ(1680) "; case PDG::Upsilon1S: return os << "Υ(1S) "; case PDG::Upsilon2S: return os << "Υ(2S) "; case PDG::Upsilon3S: return os << "Υ(3S) ";; case PDG::proton: return os << "proton"; case PDG::diffractiveProton:return os << "diffr.proton"; case PDG::neutron: return os << "neutron"; case PDG::pomeron: return os << "IP"; case PDG::reggeon: return os << "IR"; case PDG::down: return os << "d"; case PDG::up: return os << "u"; case PDG::strange: return os << "s"; case PDG::charm: return os << "c"; case PDG::bottom: return os << "b"; case PDG::top: return os << "t"; case PDG::invalid: return os << "[...]"; } return os; } } - diff --git a/CepGen/Physics/ParticleProperties.h b/CepGen/Physics/ParticleProperties.h index 44aba8e..5c3ac96 100644 --- a/CepGen/Physics/ParticleProperties.h +++ b/CepGen/Physics/ParticleProperties.h @@ -1,33 +1,33 @@ #ifndef CepGen_Physics_ParticleProperties_h #define CepGen_Physics_ParticleProperties_h namespace CepGen { enum class PDG; - namespace ParticleProperties + /// All useful properties about particles + namespace part { /** \brief Mass of a particle, in GeV/c\f${}^2\f$ * \param pdg_id PDG identifier */ double mass( const PDG& pdg_id ); /** \brief Electric charge of a particle, in \f$e\f$ * \param[in] pdg_id PDG id */ double charge( const PDG& pdg_id ); /** \brief Electric charge of a particle, in \f$e\f$ * \param[in] id integer PDG id */ double charge( int id ); /** \brief Colour factor for a given particle * \param[in] pdg_id PDG id */ unsigned short colours( const PDG& pdg_id ); /** \brief Total decay width of an unstable particle, in GeV * \param[in] pdg_id PDG id */ double width( const PDG& pdg_id ); } } #endif - diff --git a/CepGen/Processes/Fortran/KTStructures.h b/CepGen/Processes/Fortran/KTStructures.h index 2601df2..39f2732 100644 --- a/CepGen/Processes/Fortran/KTStructures.h +++ b/CepGen/Processes/Fortran/KTStructures.h @@ -1,76 +1,76 @@ #ifndef CepGen_Processes_Fortran_KTStructures_h #define CepGen_Processes_Fortran_KTStructures_h namespace CepGen { /// Collection of common blocks for Fortran kT-processes - namespace KTBlock + 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 kT-factorised flux for first incoming parton int iflux2; ///< Type of kT-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 kT-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 { int ipt; ///< Switch for cut on single particle transverse momentum int iene; ///< Switch for cut on single particle energy int ieta; ///< Switch for cut on single particle pseudo-rapidity int iinvm; ///< Switch for cut on central system invariant mass int iptsum; ///< Switch for cut on central system transverse momentum int 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 double pc[10][4]; ///< 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/FortranKTProcess.cpp b/CepGen/Processes/FortranKTProcess.cpp index 9124417..fc14e33 100644 --- a/CepGen/Processes/FortranKTProcess.cpp +++ b/CepGen/Processes/FortranKTProcess.cpp @@ -1,166 +1,166 @@ #include "CepGen/Processes/FortranKTProcess.h" #include "CepGen/Processes/Fortran/KTStructures.h" #include "CepGen/Core/ParametersList.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 params_; - extern CepGen::KTBlock::KTKinematics ktkin_; - extern CepGen::KTBlock::Cuts kincuts_; - extern CepGen::KTBlock::Event evtkin_; + extern CepGen::ktblock::Constants constants_; + extern CepGen::ktblock::Parameters params_; + extern CepGen::ktblock::KTKinematics ktkin_; + extern CepGen::ktblock::Cuts kincuts_; + extern CepGen::ktblock::Event evtkin_; } namespace CepGen { namespace process { FortranKTProcess::FortranKTProcess( const ParametersList& params, const char* name, const char* descr, std::function<void( double& )> func ) : GenericKTProcess( params, name, descr, { { PDG::photon, PDG::photon } }, { PDG::muon, PDG::muon } ), pair_( params.get<int>( "pair", 13 ) ), method_( params.get<int>( "method", 1 ) ), func_( func ) { constants_.m_p = GenericProcess::mp_; - constants_.units = Constants::GeV2toBarn; + constants_.units = constants::GeV2toBarn; constants_.pi = M_PI; - constants_.alpha_em = Constants::alphaEM; + constants_.alpha_em = constants::alphaEM; } void FortranKTProcess::preparePhaseSpace() { mom_ip1_ = event_->getOneByRole( Particle::IncomingBeam1 ).momentum(); mom_ip2_ = event_->getOneByRole( Particle::IncomingBeam2 ).momentum(); registerVariable( y1_, Mapping::linear, cuts_.cuts.central.rapidity_single, { -6., 6. }, "First central particle rapidity" ); registerVariable( y2_, Mapping::linear, cuts_.cuts.central.rapidity_single, { -6., 6. }, "Second central particle rapidity" ); registerVariable( pt_diff_, Mapping::linear, cuts_.cuts.central.pt_diff, { 0., 50. }, "Transverse momentum difference between central particles" ); registerVariable( phi_pt_diff_, Mapping::linear, cuts_.cuts.central.phi_pt_diff, { 0., 2.*M_PI }, "Central particles azimuthal angle difference" ); //=========================================================================================== // feed phase space cuts to the common block //=========================================================================================== cuts_.cuts.central.pt_single.save( (bool&)kincuts_.ipt, kincuts_.pt_min, kincuts_.pt_max ); cuts_.cuts.central.energy_single.save( (bool&)kincuts_.iene, kincuts_.ene_min, kincuts_.ene_max ); cuts_.cuts.central.eta_single.save( (bool&)kincuts_.ieta, kincuts_.eta_min, kincuts_.eta_max ); cuts_.cuts.central.mass_sum.save( (bool&)kincuts_.iinvm, kincuts_.invm_min, kincuts_.invm_max ); cuts_.cuts.central.pt_sum.save( (bool&)kincuts_.iptsum, kincuts_.ptsum_min, kincuts_.ptsum_max ); cuts_.cuts.central.rapidity_diff.save( (bool&)kincuts_.idely, kincuts_.dely_min, kincuts_.dely_max ); //=========================================================================================== // feed run parameters to the common block //=========================================================================================== params_.icontri = (int)cuts_.mode; params_.imethod = method_; params_.sfmod = (int)cuts_.structure_functions->type; params_.pdg_l = pair_; //------------------------------------------------------------------------------------------- // incoming beams information //------------------------------------------------------------------------------------------- params_.inp1 = cuts_.incoming_beams.first.pz; params_.inp2 = cuts_.incoming_beams.second.pz; const HeavyIon in1 = (HeavyIon)cuts_.incoming_beams.first.pdg; if ( in1 ) { params_.a_nuc1 = in1.A; params_.z_nuc1 = (unsigned short)in1.Z; if ( params_.z_nuc1 > 1 ) { event_->getOneByRole( Particle::IncomingBeam1 ).setPdgId( (PDG)in1 ); event_->getOneByRole( Particle::OutgoingBeam1 ).setPdgId( (PDG)in1 ); } } else params_.a_nuc1 = params_.z_nuc1 = 1; const HeavyIon in2 = (HeavyIon)cuts_.incoming_beams.second.pdg; if ( in2 ) { params_.a_nuc2 = in2.A; params_.z_nuc2 = (unsigned short)in2.Z; if ( params_.z_nuc2 > 1 ) { event_->getOneByRole( Particle::IncomingBeam2 ).setPdgId( (PDG)in2 ); event_->getOneByRole( Particle::OutgoingBeam2 ).setPdgId( (PDG)in2 ); } } else params_.a_nuc2 = params_.z_nuc2 = 1; //------------------------------------------------------------------------------------------- // intermediate partons information //------------------------------------------------------------------------------------------- params_.iflux1 = (int)cuts_.incoming_beams.first.kt_flux; params_.iflux2 = (int)cuts_.incoming_beams.second.kt_flux; if ( (KTFlux)params_.iflux1 == KTFlux::P_Gluon_KMR ) event_->getOneByRole( Particle::Parton1 ).setPdgId( PDG::gluon ); if ( (KTFlux)params_.iflux2 == KTFlux::P_Gluon_KMR ) event_->getOneByRole( Particle::Parton2 ).setPdgId( PDG::gluon ); } 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 double weight = 0.; func_( weight ); return weight; } void FortranKTProcess::fillCentralParticlesKinematics() { //=========================================================================================== // outgoing beam remnants //=========================================================================================== PX_ = Particle::Momentum( evtkin_.px ); PY_ = Particle::Momentum( evtkin_.py ); // express these momenta per nucleon PX_ *= 1./params_.a_nuc1; PY_ *= 1./params_.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_->getByRole( Particle::CentralSystem ); for ( int i = 0; i < evtkin_.nout; ++i ) { Particle& p = oc[i]; p.setPdgId( evtkin_.pdg[i] ); p.setStatus( Particle::Status::FinalState ); p.setMomentum( Particle::Momentum( evtkin_.pc[i] ) ); } } } } diff --git a/CepGen/Processes/GamGamLL.cpp b/CepGen/Processes/GamGamLL.cpp index 52c0661..5edfe56 100644 --- a/CepGen/Processes/GamGamLL.cpp +++ b/CepGen/Processes/GamGamLL.cpp @@ -1,1149 +1,1149 @@ #include "CepGen/Processes/GamGamLL.h" #include "CepGen/Event/Event.h" #include "CepGen/Core/Exception.h" #include "CepGen/Physics/Constants.h" #include "CepGen/Physics/FormFactors.h" #include "CepGen/Physics/PDG.h" #include "CepGen/Processes/ProcessesHandler.h" namespace CepGen { namespace process { //--------------------------------------------------------------------------------------------- GamGamLL::Masses::Masses() : MX2_( 0. ), MY2_( 0. ), Ml2_( 0. ), w12_( 0. ), w31_( 0. ), dw31_( 0. ), w52_( 0. ), dw52_( 0. ) {} //--------------------------------------------------------------------------------------------- GamGamLL::GamGamLL( const ParametersList& params ) : GenericProcess( "lpair", "pp → p(*) ( ɣɣ → l⁺l¯ ) p(*)" ), n_opt_( params.get<int>( "nopt", 0 ) ), pair_( params.get<int>( "pair", 0 ) ), ep1_( 0. ), ep2_( 0. ), p_cm_( 0. ), ec4_( 0. ), pc4_( 0. ), mc4_( 0. ), w4_( 0. ), p12_( 0. ), p1k2_( 0. ), p2k1_( 0. ), p13_( 0. ), p14_( 0. ), p25_( 0. ), q1dq_( 0. ), q1dq2_( 0. ), s1_( 0. ), s2_( 0. ), epsi_( 0. ), g5_( 0. ), g6_( 0. ), a5_( 0. ), a6_( 0. ), bb_( 0. ), gram_( 0. ), dd1_( 0. ), dd2_( 0. ), dd3_( 0. ), dd4_( 0. ), dd5_( 0. ), delta_( 0. ), g4_( 0. ), sa1_( 0. ), sa2_( 0. ), sl1_( 0. ), cos_theta4_( 0. ), sin_theta4_( 0. ), al4_( 0. ), be4_( 0. ), de3_( 0. ), de5_( 0. ), pt4_( 0. ), jacobian_( 0. ) {} //--------------------------------------------------------------------------------------------- void GamGamLL::addEventContent() { GenericProcess::setEventContent( { { Particle::IncomingBeam1, PDG::proton }, { Particle::IncomingBeam2, PDG::proton }, { Particle::Parton1, PDG::photon }, { Particle::Parton2, PDG::photon } }, { { Particle::OutgoingBeam1, { PDG::proton } }, { Particle::OutgoingBeam2, { PDG::proton } }, { Particle::CentralSystem, { (PDG)pair_, (PDG)pair_ } } } ); } unsigned int GamGamLL::numDimensions() const { switch ( cuts_.mode ) { case KinematicsMode::ElectronProton: default: throw CG_FATAL( "GamGamLL" ) << "Process mode " << cuts_.mode << " not (yet) supported! " << "Please contact the developers to consider an implementation."; case KinematicsMode::ElasticElastic: return 7; case KinematicsMode::ElasticInelastic: case KinematicsMode::InelasticElastic: return 8; case KinematicsMode::InelasticInelastic: return 9; } } //--------------------------------------------------------------------------------------------- void GamGamLL::setKinematics( const Kinematics& kin ) { GenericProcess::setKinematics( kin ); masses_.Ml2_ = event_->getByRole( Particle::CentralSystem )[0].mass2(); w_limits_ = cuts_.cuts.central.mass_single; if ( !w_limits_.hasMax() ) w_limits_.max() = s_; // The minimal energy for the central system is its outgoing leptons' mass energy (or wmin_ if specified) if ( !w_limits_.hasMin() ) w_limits_.min() = 4.*masses_.Ml2_; // The maximal energy for the central system is its CM energy with the outgoing particles' mass energy substracted (or wmax if specified) w_limits_.max() = std::min( pow( sqs_-MX_-MY_, 2 ), w_limits_.max() ); CG_DEBUG_LOOP( "GamGamLL:setKinematics" ) << "w limits = " << w_limits_ << "\n\t" << "wmax/wmin = " << w_limits_.max()/w_limits_.min(); q2_limits_ = cuts_.cuts.initial.q2; mx_limits_ = cuts_.cuts.remnants.mass_single; } //--------------------------------------------------------------------------------------------- bool GamGamLL::pickin() { CG_DEBUG_LOOP( "GamGamLL" ) << "Optimised mode? " << n_opt_; jacobian_ = 0.; w4_ = mc4_*mc4_; // sig1 = sigma and sig2 = sigma' in [1] const double sig = mc4_+MY_; double sig1 = sig*sig, sig2 = sig1; CG_DEBUG_LOOP( "GamGamLL" ) << "mc4 = " << mc4_ << "\n\t" << "sig1 = " << sig1 << "\n\t" << "sig2 = " << sig2; // Mass difference between the first outgoing particle // and the first incoming particle masses_.w31_ = masses_.MX2_-w1_; // Mass difference between the second outgoing particle // and the second incoming particle masses_.w52_ = masses_.MY2_-w2_; // Mass difference between the two incoming particles masses_.w12_ = w1_-w2_; // Mass difference between the central two-photons system // and the second outgoing particle const double d6 = w4_-masses_.MY2_; CG_DEBUG_LOOP( "GamGamLL" ) << "w1 = " << w1_ << "\n\t" << "w2 = " << w2_ << "\n\t" << "w3 = " << masses_.MX2_ << "\n\t" << "w4 = " << w4_ << "\n\t" << "w5 = " << masses_.MY2_;; CG_DEBUG_LOOP( "GamGamLL" ) << "w31 = " << masses_.w31_ << "\n\t" << "w52 = " << masses_.w52_ << "\n\t" << "w12 = " << masses_.w12_;; const double ss = s_+masses_.w12_; const double rl1 = ss*ss-4.*w1_*s_; // lambda(s, m1**2, m2**2) if ( rl1 <= 0. ) { CG_WARNING( "GamGamLL" ) << "rl1 = " << rl1 << " <= 0"; return false; } sl1_ = sqrt( rl1 ); s2_ = 0.; double ds2 = 0.; if ( n_opt_ == 0 ) { const double smax = s_+masses_.MX2_-2.*MX_*sqs_; map( x(2), Limits( sig1, smax ), s2_, ds2, "s2" ); sig1 = s2_; //FIXME!!!!!!!!!!!!!!!!!!!! } CG_DEBUG_LOOP( "GamGamLL" ) << "s2 = " << s2_; const double sp = s_+masses_.MX2_-sig1, d3 = sig1-w2_; const double rl2 = sp*sp-4.*s_*masses_.MX2_; // lambda(s, m3**2, sigma) if ( rl2 <= 0. ) { CG_DEBUG( "GamGamLL" ) << "rl2 = " << rl2 << " <= 0"; return false; } const double sl2 = sqrt( rl2 ); double t1_max = w1_+masses_.MX2_-( ss*sp+sl1_*sl2 )/( 2.*s_ ), // definition from eq. (A.4) in [1] t1_min = ( masses_.w31_*d3+( d3-masses_.w31_ )*( d3*w1_-masses_.w31_*w2_ )/s_ )/t1_max; // definition from eq. (A.5) in [1] // FIXME dropped in CDF version if ( t1_max > -q2_limits_.min() ) { CG_WARNING( "GamGamLL" ) << "t1max = " << t1_max << " > -q2min = " << ( -q2_limits_.min() ); return false; } if ( t1_min < -q2_limits_.max() && q2_limits_.hasMax() ) { CG_DEBUG( "GamGamLL" ) << "t1min = " << t1_min << " < -q2max = " << -q2_limits_.max(); return false; } if ( t1_max < -q2_limits_.max() && q2_limits_.hasMax() ) t1_max = -q2_limits_.max(); if ( t1_min > -q2_limits_.min() && q2_limits_.hasMin() ) t1_min = -q2_limits_.min(); ///// // t1, the first photon propagator, is defined here t1_ = 0.; double dt1 = 0.; map( x(0), Limits( t1_min, t1_max ), t1_, dt1, "t1" ); // changes wrt mapt1 : dx->-dx dt1 *= -1.; CG_DEBUG_LOOP( "GamGamLL" ) << "Definition of t1 = " << t1_ << " according to\n\t" << "(t1min, t1max) = (" << t1_min << ", " << t1_max << ")"; dd4_ = w4_-t1_; const double d8 = t1_-w2_, t13 = t1_-w1_-masses_.MX2_; sa1_ = -pow( t1_-masses_.w31_, 2 )/4.+w1_*t1_; if ( sa1_ >= 0. ) { CG_WARNING( "GamGamLL" ) << "sa1_ = " << sa1_ << " >= 0"; return false; } const double sl3 = sqrt( -sa1_ ); // one computes splus and (s2x=s2max) double splus, s2max; if ( w1_ != 0. ) { const double inv_w1 = 1./w1_; const double sb = masses_.MX2_ + 0.5 * ( s_*( t1_-masses_.w31_ )+masses_.w12_*t13 )*inv_w1, sd = sl1_*sl3*inv_w1, se =( s_*( t1_*( s_+t13-w2_ )-w2_*masses_.w31_ )+masses_.MX2_*( masses_.w12_*d8+w2_*masses_.MX2_ ) )*inv_w1; if ( fabs( ( sb-sd )/sd ) >= 1. ) { splus = sb-sd; s2max = se/splus; } else { s2max = sb+sd; splus = se/s2max; } } else { // 3 s2max = ( s_*( t1_*( s_+d8-masses_.MX2_ )-w2_*masses_.MX2_ )+w2_*masses_.MX2_*( w2_+masses_.MX2_-t1_ ) )/( ss*t13 ); splus = sig2; } // 4 double s2x = s2max; CG_DEBUG_LOOP( "GamGamLL" ) << "s2x = s2max = " << s2x; if ( n_opt_ < 0 ) { // 5 if ( splus > sig2 ) { sig2 = splus; CG_DEBUG_LOOP( "GamGamLL" ) << "sig2 truncated to splus = " << splus; } if ( n_opt_ < -1 ) map( x(2), Limits( sig2, s2max ), s2_, ds2, "s2" ); else mapla( t1_, w2_, x(2), sig2, s2max, s2_, ds2 ); // n_opt_==-1 s2x = s2_; } else if ( n_opt_ == 0 ) s2x = s2_; // 6 CG_DEBUG_LOOP( "GamGamLL" ) << "s2x = " << s2x; // 7 const double r1 = s2x-d8, r2 = s2x-d6; const double rl4 = ( r1*r1-4.*w2_*s2x )*( r2*r2-4.*masses_.MY2_*s2x ); if ( rl4 <= 0. ) { CG_DEBUG_LOOP( "GamGamLL" ) << "rl4 = " << rl4 << " <= 0"; return false; } const double sl4 = sqrt( rl4 ); // t2max, t2min definitions from eq. (A.12) and (A.13) in [1] const double t2_max = w2_+masses_.MY2_-( r1*r2+sl4 )/s2x * 0.5, t2_min = ( masses_.w52_*dd4_+( dd4_-masses_.w52_ )*( dd4_*w2_-masses_.w52_*t1_ )/s2x )/t2_max; // t2, the second photon propagator, is defined here t2_ = 0.; double dt2 = 0.; map( x(1), Limits( t2_min, t2_max ), t2_, dt2, "t2" ); // changes wrt mapt2 : dx->-dx dt2 *= -1.; // \f$\delta_6=m_4^2-m_5^2\f$ as defined in Vermaseren's paper const double tau = t1_-t2_, r3 = dd4_-t2_, r4 = masses_.w52_-t2_; CG_DEBUG_LOOP( "GamGamLL" ) << "r1 = " << r1 << "\n\t" << "r2 = " << r2 << "\n\t" << "r3 = " << r3 << "\n\t" << "r4 = " << r4; const double b = r3*r4-2.*( t1_+w2_ )*t2_; const double c = t2_*d6*d8+( d6-d8 )*( d6*w2_-d8*masses_.MY2_ ); const double t25 = t2_-w2_-masses_.MY2_; sa2_ = -0.25 * r4*r4 + w2_*t2_; if ( sa2_ >= 0. ) { CG_WARNING( "GamGamLL" ) << "sa2_ = " << sa2_ << " >= 0"; return false; } const double sl6 = 2.*sqrt( -sa2_ ); g4_ = -r3*r3/4.+t1_*t2_; if ( g4_ >= 0. ) { CG_WARNING( "GamGamLL" ) << "g4_ = " << g4_ << " >= 0"; return false; } const double sl7 = 2.*sqrt( -g4_ ), sl5 = sl6*sl7; double s2p, s2min; if ( fabs( ( sl5-b )/sl5 ) >= 1. ) { s2p = 0.5 * ( sl5-b )/t2_; s2min = c/( t2_*s2p ); } else { // 8 s2min = 0.5 * ( -sl5-b )/t2_; s2p = c/( t2_*s2min ); } // 9 if ( n_opt_ > 1 ) map( x( 2 ), Limits( s2min, s2max ), s2_, ds2, "s2" ); else if ( n_opt_ == 1 ) mapla( t1_, w2_, x( 2 ), s2min, s2max, s2_, ds2 ); const double ap = -0.25*pow( s2_+d8, 2 )+s2_*t1_; CG_DEBUG_LOOP( "GamGamLL" ) << "s2 = " << s2_ << ", s2max = " << s2max << ", splus = " << splus; if ( w1_ != 0. ) dd1_ = -0.25 * ( s2_-s2max ) * ( s2_-splus ) * w1_; // 10 else dd1_ = 0.25 * ( s2_-s2max ) * ss * t13; // 11 dd2_ = -0.25 * t2_*( s2_-s2p )*( s2_-s2min ); CG_DEBUG_LOOP( "GamGamLL" ) << "t2 = " << t2_ << "\n\t" << "s2 = " << s2_ << "\n\t" << "s2p = " << s2p << "\n\t" << "s2min = " << s2min << "\n\t" << "dd2 = " << dd2_;; const double yy4 = cos( M_PI*x( 3 ) ); const double dd = dd1_*dd2_; p12_ = 0.5 * ( s_-w1_-w2_ ); const double st = s2_-t1_-w2_; const double delb = ( 2.*w2_*r3+r4*st )*( 4.*p12_*t1_-( t1_-masses_.w31_ )*st )/( 16.*ap ); if ( dd <= 0. ) { CG_DEBUG_LOOP( "GamGamLL" ) << std::scientific << "dd = " << dd << " <= 0\n\t" << "dd1 = " << dd1_ << "\t" << "dd2 = " << dd2_ << std::fixed; return false; } delta_ = delb - yy4*st*sqrt( dd )/ ap * 0.5; s1_ = t2_+w1_+( 2.*p12_*r3-4.*delta_ )/st; if ( ap >= 0. ) { CG_DEBUG_LOOP( "GamGamLL" ) << "ap = " << ap << " >= 0"; return false; } jacobian_ = ds2 * dt1 * dt2 * 0.125 * M_PI*M_PI/( sl1_*sqrt( -ap ) ); CG_DEBUG_LOOP( "GamGamLL" ) << "Jacobian = " << std::scientific << jacobian_ << std::fixed; gram_ = ( 1.-yy4*yy4 )*dd/ap; p13_ = -0.5 * t13; p14_ = 0.5 * ( tau+s1_-masses_.MX2_ ); p25_ = -0.5 * t25; p1k2_ = 0.5 * ( s1_-t2_-w1_ ); p2k1_ = 0.5 * st; if ( w2_ != 0. ) { const double inv_w2 = 1./w2_; const double sbb = 0.5 * ( s_*( t2_-masses_.w52_ )-masses_.w12_*t25 )*inv_w2 + masses_.MY2_, sdd = 0.5 * sl1_*sl6*inv_w2, see = ( s_*( t2_*( s_+t25-w1_ )-w1_*masses_.w52_ )+masses_.MY2_*( w1_*masses_.MY2_-masses_.w12_*( t2_-w1_ ) ) )*inv_w2; double s1m = 0., s1p = 0.; if ( sbb/sdd >= 0. ) { s1p = sbb+sdd; s1m = see/s1p; } else { s1m = sbb-sdd; s1p = see/s1m; } // 12 dd3_ = -0.25 * w2_*( s1p-s1_ )*( s1m-s1_ ); // 13 } else { // 14 const double s1p = ( s_*( t2_*( s_-masses_.MY2_+t2_-w1_ )-w1_*masses_.MY2_ )+w1_*masses_.MY2_*( w1_+masses_.MY2_-t2_ ) )/( t25*( s_-masses_.w12_ ) ); dd3_ = -0.25 * t25*( s_-masses_.w12_ )*( s1p-s1_ ); } // 15 //const double acc3 = (s1p-s1_)/(s1p+s1_); const double ssb = t2_+0.5 * w1_-r3*( masses_.w31_-t1_ )/t1_, ssd = sl3*sl7/t1_, sse = ( t2_-w1_ )*( w4_-masses_.MX2_ )+( t2_-w4_+masses_.w31_ )*( ( t2_-w1_)*masses_.MX2_-(w4_-masses_.MX2_)*w1_)/t1_; double s1pp, s1pm; if ( ssb/ssd >= 0. ) { s1pp = ssb+ssd; s1pm = sse/s1pp; } else { // 16 s1pm = ssb-ssd; s1pp = sse/s1pm; } // 17 dd4_ = -0.25 * t1_*( s1_-s1pp )*( s1_-s1pm ); //const double acc4 = ( s1_-s1pm )/( s1_+s1pm ); dd5_ = dd1_+dd3_+( ( p12_*( t1_-masses_.w31_ )*0.5-w1_*p2k1_ )*( p2k1_*( t2_-masses_.w52_ )-w2_*r3 ) -delta_*( 2.*p12_*p2k1_-w2_*( t1_-masses_.w31_ ) ) ) / p2k1_; return true; } //--------------------------------------------------------------------------------------------- bool GamGamLL::orient() { if ( !pickin() || jacobian_ == 0. ) { CG_DEBUG_LOOP( "GamGamLL" ) << "Pickin failed! Jacobian = " << jacobian_; return false; } const double re = 0.5 / sqs_; ep1_ = re*( s_+masses_.w12_ ); ep2_ = re*( s_-masses_.w12_ ); CG_DEBUG_LOOP( "GamGamLL" ) << std::scientific << " re = " << re << "\n\t" << "w12_ = " << masses_.w12_ << std::fixed; CG_DEBUG_LOOP( "GamGamLL" ) << "Incoming particles' energy = " << ep1_ << ", " << ep2_; p_cm_ = re*sl1_; de3_ = re*(s2_-masses_.MX2_+masses_.w12_); de5_ = re*(s1_-masses_.MY2_-masses_.w12_); // Final state energies const double ep3 = ep1_-de3_, ep5 = ep2_-de5_; ec4_ = de3_+de5_; if ( ec4_ < mc4_ ) { CG_WARNING( "GamGamLL" ) << "ec4_ = " << ec4_ << " < mc4_ = " << mc4_ << "\n\t" << "==> de3 = " << de3_ << ", de5 = " << de5_; return false; } // What if the protons' momenta are not along the z-axis? pc4_ = sqrt( ec4_*ec4_-mc4_*mc4_ ); if ( pc4_ == 0. ) { CG_WARNING( "GamGamLL" ) << "pzc4 is null and should not be..."; return false; } const double pp3 = sqrt( ep3*ep3-masses_.MX2_ ), pt3 = sqrt( dd1_/s_ )/p_cm_, pp5 = sqrt( ep5*ep5-masses_.MY2_ ), pt5 = sqrt( dd3_/s_ )/p_cm_; CG_DEBUG_LOOP( "GamGamLL" ) << "Central system's energy: E4 = " << ec4_ << "\n\t" << " momentum: p4 = " << pc4_ << "\n\t" << " invariant mass: m4 = " << mc4_ << "\n\t" << "Outgoing particles' energy: E3 = " << ep3 << "\n\t" << " E5 = " << ep5; const double sin_theta3 = pt3/pp3, sin_theta5 = pt5/pp5; CG_DEBUG_LOOP( "GamGamLL" ) << std::scientific << "sin_theta3 = " << sin_theta3 << "\n\t" << "sin_theta5 = " << sin_theta5 << std::fixed; if ( sin_theta3 > 1. ) { CG_WARNING( "GamGamLL" ) << "sin_theta3 = " << sin_theta3 << " > 1"; return false; } if ( sin_theta5 > 1. ) { CG_WARNING( "GamGamLL" ) << "sin_theta5 = " << sin_theta5 << " > 1"; return false; } double ct3 = sqrt( 1.-sin_theta3*sin_theta3 ), ct5 = sqrt( 1.-sin_theta5*sin_theta5 ); if ( ep1_*ep3 < p13_ ) ct3 *= -1.; if ( ep2_*ep5 > p25_ ) ct5 *= -1.; CG_DEBUG_LOOP( "GamGamLL" ) << "ct3 = " << ct3 << "\n\t" << "ct5 = " << ct5; if ( dd5_ < 0. ) { CG_WARNING( "GamGamLL" ) << "dd5 = " << dd5_ << " < 0"; return false; } // Centre of mass system kinematics (theta4 and phi4) pt4_ = sqrt( dd5_/s_ )/p_cm_; sin_theta4_ = pt4_/pc4_; if ( sin_theta4_ > 1. ) { CG_WARNING( "GamGamLL" ) << "st4 = " << sin_theta4_ << " > 1"; return false; } cos_theta4_ = sqrt( 1.-sin_theta4_*sin_theta4_ ); if ( ep1_*ec4_ < p14_ ) cos_theta4_ *= -1.; al4_ = 1.-cos_theta4_; be4_ = 1.+cos_theta4_; if ( cos_theta4_ < 0. ) be4_ = sin_theta4_*sin_theta4_/al4_; else al4_ = sin_theta4_*sin_theta4_/be4_; CG_DEBUG_LOOP( "GamGamLL" ) << "ct4 = " << cos_theta4_ << "\n\t" << "al4 = " << al4_ << ", be4 = " << be4_; const double rr = sqrt( -gram_/s_ )/( p_cm_*pt4_ ); const double sin_phi3 = rr / pt3, sin_phi5 = -rr / pt5; if ( fabs( sin_phi3 ) > 1. ) { CG_WARNING( "GamGamLL" ) << "sin(phi_3) = " << sin_phi3 << " while it must be in [-1 ; 1]"; return false; } if ( fabs( sin_phi5 ) > 1. ) { CG_WARNING( "GamGamLL" ) << "sin(phi_5) = " << sin_phi5 << " while it must be in [-1 ; 1]"; return false; } const double cos_phi3 = -sqrt( 1.-sin_phi3*sin_phi3 ), cos_phi5 = -sqrt( 1.-sin_phi5*sin_phi5 ); p3_lab_ = Particle::Momentum( pp3*sin_theta3*cos_phi3, pp3*sin_theta3*sin_phi3, pp3*ct3, ep3 ); p5_lab_ = Particle::Momentum( pp5*sin_theta5*cos_phi5, pp5*sin_theta5*sin_phi5, pp5*ct5, ep5 ); const double a1 = p3_lab_.px()-p5_lab_.px(); CG_DEBUG_LOOP( "GamGamLL" ) << "Kinematic quantities\n\t" << "cos(theta3) = " << ct3 << "\t" << "sin(theta3) = " << sin_theta3 << "\n\t" << "cos( phi3 ) = " << cos_phi3 << "\t" << "sin( phi3 ) = " << sin_phi3 << "\n\t" << "cos(theta4) = " << cos_theta4_ << "\t" << "sin(theta4) = " << sin_theta4_ << "\n\t" << "cos(theta5) = " << ct5 << "\t" << "sin(theta5) = " << sin_theta5 << "\n\t" << "cos( phi5 ) = " << cos_phi5 << "\t" << "sin( phi5 ) = " << sin_phi5 << "\n\t" << "a1 = " << a1; if ( fabs( pt4_+p3_lab_.px()+p5_lab_.px() ) < fabs( fabs( a1 )-pt4_ ) ) { CG_DEBUG_LOOP( "GamGamLL" ) << "|pt4+pt3*cos(phi3)+pt5*cos(phi5)| < | |a1|-pt4 |\n\t" << "pt4 = " << pt4_ << "\t" << "pt5 = " << pt5 << "\n\t" << "cos(phi3) = " << cos_phi3 << "\t" << "cos(phi5) = " << cos_phi5 << "\n\t" << "a1 = " << a1; return true; } if ( a1 < 0. ) p5_lab_[0] = -p5_lab_.px(); else p3_lab_[0] = -p3_lab_.px(); return true; } //--------------------------------------------------------------------------------------------- double GamGamLL::computeOutgoingPrimaryParticlesMasses( double x, double outmass, double lepmass, double& dw ) { - const double mx0 = mp_+ParticleProperties::mass( PDG::piZero ); // 1.07 + const double mx0 = mp_+part::mass( PDG::piZero ); // 1.07 const double wx2min = pow( std::max( mx0, mx_limits_.min() ), 2 ), wx2max = pow( std::min( sqs_-outmass-2.*lepmass, mx_limits_.max() ), 2 ); double mx2 = 0., dmx2 = 0.; map( x, Limits( wx2min, wx2max ), mx2, dmx2, "mx2" ); CG_DEBUG_LOOP( "GamGamLL" ) << "mX^2 in range (" << wx2min << ", " << wx2max << "), x = " << x << "\n\t" << "mX^2 = " << mx2 << ", d(mX^2) = " << dmx2 << "\n\t" << "mX = " << sqrt( mx2 ) << ", d(mX) = " << sqrt( dmx2 ); dw = sqrt( dmx2 ); return sqrt( mx2 ); } //--------------------------------------------------------------------------------------------- void GamGamLL::beforeComputeWeight() { if ( !GenericProcess::is_point_set_ ) return; const Particle& p1 = event_->getOneByRole( Particle::IncomingBeam1 ), &p2 = event_->getOneByRole( Particle::IncomingBeam2 ); ep1_ = p1.energy(); ep2_ = p2.energy(); //std::cout << __PRETTY_FUNCTION__ << ":" << w_limits_ << "|" << q2_limits_ << "|" << mx_limits_ << std::endl; switch ( cuts_.mode ) { case KinematicsMode::ElectronProton: default: CG_ERROR( "GamGamLL" ) << "Case not yet supported!"; break; case KinematicsMode::ElasticElastic: masses_.dw31_ = masses_.dw52_ = 0.; break; case KinematicsMode::InelasticElastic: { const double m = computeOutgoingPrimaryParticlesMasses( x( 7 ), p1.mass(), sqrt( masses_.Ml2_ ), masses_.dw31_ ); event_->getOneByRole( Particle::OutgoingBeam1 ).setMass( m ); - event_->getOneByRole( Particle::OutgoingBeam2 ).setMass( ParticleProperties::mass( p2.pdgId() ) ); + event_->getOneByRole( Particle::OutgoingBeam2 ).setMass( part::mass( p2.pdgId() ) ); } break; case KinematicsMode::ElasticInelastic: { const double m = computeOutgoingPrimaryParticlesMasses( x( 7 ), p2.mass(), sqrt( masses_.Ml2_ ), masses_.dw52_ ); - event_->getOneByRole( Particle::OutgoingBeam1 ).setMass( ParticleProperties::mass( p1.pdgId() ) ); + event_->getOneByRole( Particle::OutgoingBeam1 ).setMass( part::mass( p1.pdgId() ) ); event_->getOneByRole( Particle::OutgoingBeam2 ).setMass( m ); } break; case KinematicsMode::InelasticInelastic: { const double mx = computeOutgoingPrimaryParticlesMasses( x( 7 ), p2.mass(), sqrt( masses_.Ml2_ ), masses_.dw31_ ); event_->getOneByRole( Particle::OutgoingBeam1 ).setMass( mx ); const double my = computeOutgoingPrimaryParticlesMasses( x( 8 ), p1.mass(), sqrt( masses_.Ml2_ ), masses_.dw52_ ); event_->getOneByRole( Particle::OutgoingBeam2 ).setMass( my ); } break; } MX_ = event_->getOneByRole( Particle::OutgoingBeam1 ).mass(); MY_ = event_->getOneByRole( Particle::OutgoingBeam2 ).mass(); masses_.MX2_ = MX_*MX_; masses_.MY2_ = MY_*MY_; } //--------------------------------------------------------------------------------------------- double GamGamLL::computeWeight() { CG_DEBUG_LOOP( "GamGamLL" ) << "sqrt(s) = " << sqs_ << " GeV\n\t" << "m(X1) = " << MX_ << " GeV\t" << "m(X2) = " << MY_ << " GeV"; // compute the two-photon energy for this point w4_ = 0.; double dw4 = 0.; map( x( 4 ), w_limits_, w4_, dw4, "w4" ); mc4_ = sqrt( w4_ ); CG_DEBUG_LOOP( "GamGamLL" ) << "Computed value for w4 = " << w4_ << " → mc4 = " << mc4_; if ( !orient() ) return 0.; if ( jacobian_ == 0. ) { CG_WARNING( "GamGamLL" ) << "dj = " << jacobian_; return 0.; } if ( t1_ > 0. ) { CG_WARNING( "GamGamLL" ) << "t1 = " << t1_ << " > 0"; return 0.; } if ( t2_ > 0. ) { CG_WARNING( "GamGamLL" ) << "t2 = " << t2_ << " > 0"; return 0.; } const double ecm6 = w4_ / ( 2.*mc4_ ), pp6cm = sqrt( ecm6*ecm6-masses_.Ml2_ ); - jacobian_ *= dw4*pp6cm/( mc4_*Constants::sconstb*s_ ); + jacobian_ *= dw4*pp6cm/( mc4_*constants::sconstb*s_ ); // Let the most obscure part of this code begin... const double e1mp1 = w1_ / ( ep1_+p_cm_ ), e3mp3 = masses_.MX2_ / ( p3_lab_.energy()+p3_lab_.p() ); const double al3 = pow( sin( p3_lab_.theta() ), 2 )/( 1.+( p3_lab_.theta() ) ); // 2-photon system kinematics ?! const double eg = ( w4_+t1_-t2_ )/( 2.*mc4_ ); double pg = sqrt( eg*eg-t1_ ); const double pgx = -p3_lab_.px()*cos_theta4_-sin_theta4_*( de3_-e1mp1 + e3mp3 + p3_lab_.p()*al3 ), pgy = -p3_lab_.py(), pgz = mc4_*de3_/( ec4_+pc4_ )-ec4_*de3_*al4_/mc4_-p3_lab_.px()*ec4_*sin_theta4_/mc4_+ec4_*cos_theta4_/mc4_*( p3_lab_.p()*al3+e3mp3-e1mp1 ); CG_DEBUG_LOOP( "GamGamLL" ) << "pg = " << Particle::Momentum( pgx, pgy, pgz ); const double pgp = sqrt( pgx*pgx + pgy*pgy ), // outgoing proton (3)'s transverse momentum pgg = sqrt( pgp*pgp + pgz*pgz ); // outgoing proton (3)'s momentum if ( pgg > pgp*0.9 && pgg > pg ) pg = pgg; //FIXME ??? // angles for the 2-photon system ?! const double cpg = pgx/pgp, spg = pgy/pgp; const double stg = pgp/pg; const int theta_sign = ( pgz>0. ) ? 1 : -1; const double ctg = theta_sign*sqrt( 1.-stg*stg ); double xx6 = x( 5 ); const double amap = 0.5 * ( w4_-t1_-t2_ ), bmap = 0.5 * sqrt( ( pow( w4_-t1_-t2_, 2 )-4.*t1_*t2_ )*( 1.-4.*masses_.Ml2_/w4_ ) ), ymap = ( amap+bmap )/( amap-bmap ), beta = pow( ymap, 2.*xx6-1. ); xx6 = 0.5 * ( 1. + amap/bmap*( beta-1. )/( beta+1. ) ); xx6 = std::max( 0., std::min( xx6, 1. ) ); // xx6 in [0., 1.] CG_DEBUG_LOOP( "GamGamLL" ) << "amap = " << amap << "\n\t" << "bmap = " << bmap << "\n\t" << "ymap = " << ymap << "\n\t" << "beta = " << beta; // 3D rotation of the first outgoing lepton wrt the CM system const double theta6cm = acos( 1.-2.*xx6 ); // match the Jacobian jacobian_ *= ( amap+bmap*cos( theta6cm ) ); jacobian_ *= ( amap-bmap*cos( theta6cm ) ); jacobian_ /= amap; jacobian_ /= bmap; jacobian_ *= log( ymap ); jacobian_ *= 0.5; CG_DEBUG_LOOP( "GamGamLL" ) << "Jacobian = " << jacobian_; CG_DEBUG_LOOP( "GamGamLL" ) << "ctcm6 = " << cos( theta6cm ) << "\n\t" << "stcm6 = " << sin( theta6cm ); const double phi6cm = 2.*M_PI*x( 6 ); // First outgoing lepton's 3-momentum in the centre of mass system Particle::Momentum p6cm = Particle::Momentum::fromPThetaPhi( pp6cm, theta6cm, phi6cm ); CG_DEBUG_LOOP( "GamGamLL" ) << "p3cm6 = " << p6cm; const double h1 = stg*p6cm.pz()+ctg*p6cm.px(); const double pc6z = ctg*p6cm.pz()-stg*p6cm.px(), pc6x = cpg*h1-spg*p6cm.py(); const double qcx = 2.*pc6x, qcz = 2.*pc6z; // qcy == QCY is never defined const double el6 = ( ec4_*ecm6+pc4_*pc6z ) / mc4_, h2 = ( ec4_*pc6z+pc4_*ecm6 ) / mc4_; CG_DEBUG_LOOP( "GamGamLL" ) << "h1 = " << h1 << "\n\th2 = " << h2; // first outgoing lepton's 3-momentum const double p6x = cos_theta4_*pc6x+sin_theta4_*h2, p6y = cpg*p6cm.py()+spg*h1, p6z = cos_theta4_*h2-sin_theta4_*pc6x; // first outgoing lepton's kinematics p6_cm_ = Particle::Momentum( p6x, p6y, p6z, el6 ); CG_DEBUG_LOOP( "GamGamLL" ) << "p6(cm) = " << p6_cm_; const double hq = ec4_*qcz/mc4_; const Particle::Momentum qve( cos_theta4_*qcx+sin_theta4_*hq, 2.*p6y, cos_theta4_*hq-sin_theta4_*qcx, pc4_*qcz/mc4_ // energy ); // Available energy for the second lepton is the 2-photon system's energy with the first lepton's energy removed const double el7 = ec4_-el6; CG_DEBUG_LOOP( "GamGamLL" ) << "Outgoing kinematics\n\t" << " first outgoing lepton: p = " << p6_cm_.p() << ", E = " << p6_cm_.energy() << "\n\t" << "second outgoing lepton: p = " << p7_cm_.p() << ", E = " << p7_cm_.energy();; // Second outgoing lepton's 3-momentum const double p7x = -p6x + pt4_, p7y = -p6y, p7z = -p6z + pc4_*cos_theta4_; // second outgoing lepton's kinematics p7_cm_ = Particle::Momentum( p7x, p7y, p7z, el7 ); //p6_cm_ = Particle::Momentum(pl6*st6*cp6, pl6*st6*sp6, pl6*ct6, el6); //p7_cm_ = Particle::Momentum(pl7*st7*cp7, pl7*st7*sp7, pl7*ct7, el7); q1dq_ = eg*( 2.*ecm6-mc4_ )-2.*pg*p6cm.pz(); q1dq2_ = ( w4_-t1_-t2_ ) * 0.5; const double phi3 = p3_lab_.phi(), cos_phi3 = cos( phi3 ), sin_phi3 = sin( phi3 ), phi5 = p5_lab_.phi(), cos_phi5 = cos( phi5 ), sin_phi5 = sin( phi5 ); bb_ = t1_*t2_+( w4_*pow( sin( theta6cm ), 2 ) + 4.*masses_.Ml2_*pow( cos( theta6cm ), 2 ) )*pg*pg; const double c1 = p3_lab_.pt() * ( qve.px()*sin_phi3 - qve.py()*cos_phi3 ), c2 = p3_lab_.pt() * ( qve.pz()*ep1_ - qve.energy() *p_cm_ ), c3 = ( masses_.w31_*ep1_*ep1_ + 2.*w1_*de3_*ep1_ - w1_*de3_*de3_ + p3_lab_.pt2()*ep1_*ep1_ ) / ( p3_lab_.energy()*p_cm_ + p3_lab_.pz()*ep1_ ); const double b1 = p5_lab_.pt() * ( qve.px()*sin_phi5 - qve.py()*cos_phi5 ), b2 = p5_lab_.pt() * ( qve.pz()*ep2_ + qve.energy() *p_cm_ ), b3 = ( masses_.w52_*ep2_*ep2_ + 2.*w2_*de5_*ep2_ - w2_*de5_*de5_ + p5_lab_.pt2()*ep2_*ep2_ ) / ( ep2_*p5_lab_.pz() - p5_lab_.energy()*p_cm_ ); const double r12 = c2*sin_phi3 + qve.py()*c3, r13 = -c2*cos_phi3 - qve.px()*c3; const double r22 = b2*sin_phi5 + qve.py()*b3, r23 = -b2*cos_phi5 - qve.px()*b3; epsi_ = p12_*c1*b1 + r12*r22 + r13*r23; g5_ = w1_*c1*c1 + r12*r12 + r13*r13; g6_ = w2_*b1*b1 + r22*r22 + r23*r23; const double pt3 = p3_lab_.pt(), pt5 = p5_lab_.pt(); a5_ = -( qve.px()*cos_phi3 + qve.py()*sin_phi3 )*pt3*p1k2_ -( ep1_*qve.energy()-p_cm_*qve.pz() )*( cos_phi3*cos_phi5 + sin_phi3*sin_phi5 )*pt3*pt5 +( de5_*qve.pz()+qve.energy()*( p_cm_+p5_lab_.pz() ) )*c3; a6_ = -( qve.px()*cos_phi5 + qve.py()*sin_phi5 )*pt5*p2k1_ -( ep2_*qve.energy()+p_cm_*qve.pz() )*( cos_phi3*cos_phi5 + sin_phi3*sin_phi5 )*pt3*pt5 +( de3_*qve.pz()-qve.energy()*( p_cm_-p3_lab_.pz() ) )*b3; CG_DEBUG_LOOP( "GamGamLL" ) << "a5 = " << a5_ << "\n\t" << "a6 = " << a6_; //////////////////////////////////////////////////////////////// // END of GAMGAMLL subroutine in the FORTRAN version //////////////////////////////////////////////////////////////// const Particle::Momentum cm = event_->getOneByRole( Particle::IncomingBeam1 ).momentum() + event_->getOneByRole( Particle::IncomingBeam2 ).momentum(); //////////////////////////////////////////////////////////////// // INFO from f.f //////////////////////////////////////////////////////////////// const double gamma = cm.energy() / sqs_, betgam = cm.pz() / sqs_; //--- kinematics computation for both leptons p6_cm_.betaGammaBoost( gamma, betgam ); p7_cm_.betaGammaBoost( gamma, betgam ); //--- cut on mass of final hadronic system (MX/Y) if ( mx_limits_.valid() ) { if ( ( cuts_.mode == KinematicsMode::InelasticElastic || cuts_.mode == KinematicsMode::InelasticInelastic ) && !mx_limits_.passes( MX_ ) ) return 0.; if ( ( cuts_.mode == KinematicsMode::ElasticInelastic || cuts_.mode == KinematicsMode::InelasticInelastic ) && !mx_limits_.passes( MY_ ) ) return 0.; } //--- cut on the proton's Q2 (first photon propagator T1) if ( !cuts_.cuts.initial.q2.passes( -t1_ ) ) return 0.; //--- cuts on outgoing leptons' kinematics if ( !cuts_.cuts.central.mass_sum.passes( ( p6_cm_+p7_cm_ ).mass() ) ) return 0.; //----- cuts on the individual leptons if ( cuts_.cuts.central.pt_single.valid() ) { const Limits& pt_limits = cuts_.cuts.central.pt_single; if ( !pt_limits.passes( p6_cm_.pt() ) || !pt_limits.passes( p7_cm_.pt() ) ) return 0.; } if ( cuts_.cuts.central.energy_single.valid() ) { const Limits& energy_limits = cuts_.cuts.central.energy_single; if ( !energy_limits.passes( p6_cm_.energy() ) || !energy_limits.passes( p7_cm_.energy() ) ) return 0.; } if ( cuts_.cuts.central.eta_single.valid() ) { const Limits& eta_limits = cuts_.cuts.central.eta_single; if ( !eta_limits.passes( p6_cm_.eta() ) || !eta_limits.passes( p7_cm_.eta() ) ) return 0.; } //--- compute the structure functions factors switch ( cuts_.mode ) { // inherited from CDF version case KinematicsMode::ElectronProton: default: jacobian_ *= periPP( 1, 2 ); break; // ep case case KinematicsMode::ElasticElastic: jacobian_ *= periPP( 2, 2 ); break; // elastic case case KinematicsMode::InelasticElastic: jacobian_ *= periPP( 3, 2 )*( masses_.dw31_*masses_.dw31_ ); break; case KinematicsMode::ElasticInelastic: jacobian_ *= periPP( 3, 2 )*( masses_.dw52_*masses_.dw52_ ); break; // single-dissociative case case KinematicsMode::InelasticInelastic: jacobian_ *= periPP( 3, 3 )*( masses_.dw31_*masses_.dw31_ )*( masses_.dw52_*masses_.dw52_ ); break; // double-dissociative case } //--- compute the event weight using the Jacobian - return Constants::GeV2toBarn*jacobian_; + return constants::GeV2toBarn*jacobian_; } //--------------------------------------------------------------------------------------------- void GamGamLL::fillKinematics( bool ) { const Particle::Momentum cm = event_->getOneByRole( Particle::IncomingBeam1 ).momentum() + event_->getOneByRole( Particle::IncomingBeam2 ).momentum(); const double gamma = cm.energy()/sqs_, betgam = cm.pz()/sqs_; Particle::Momentum plab_ip1 = Particle::Momentum( 0., 0., p_cm_, ep1_ ).betaGammaBoost( gamma, betgam ); Particle::Momentum plab_ip2 = Particle::Momentum( 0., 0., -p_cm_, ep2_ ).betaGammaBoost( gamma, betgam ); p3_lab_.betaGammaBoost( gamma, betgam ); p5_lab_.betaGammaBoost( gamma, betgam ); //----- parameterise a random rotation around z-axis const int rany = ( rand() >= 0.5 * RAND_MAX ) ? 1 : -1, ransign = ( rand() >= 0.5 * RAND_MAX ) ? 1 : -1; const double ranphi = ( 0.5 * rand()/RAND_MAX )*2.*M_PI; Particle::Momentum plab_ph1 = ( plab_ip1-p3_lab_ ).rotatePhi( ranphi, rany ); Particle::Momentum plab_ph2 = ( plab_ip2-p5_lab_ ).rotatePhi( ranphi, rany ); p3_lab_.rotatePhi( ranphi, rany ); p5_lab_.rotatePhi( ranphi, rany ); p6_cm_.rotatePhi( ranphi, rany ); p7_cm_.rotatePhi( ranphi, rany ); /*if ( symmetrise_ && rand() >= .5*RAND_MAX ) { p6_cm_.mirrorZ(); p7_cm_.mirrorZ(); }*/ //----- incoming protons event_->getOneByRole( Particle::IncomingBeam1 ).setMomentum( plab_ip1 ); event_->getOneByRole( Particle::IncomingBeam2 ).setMomentum( plab_ip2 ); //----- first outgoing proton Particle& op1 = event_->getOneByRole( Particle::OutgoingBeam1 ); op1.setMomentum( p3_lab_ ); switch ( cuts_.mode ) { case KinematicsMode::ElasticElastic: case KinematicsMode::ElasticInelastic: default: op1.setStatus( Particle::Status::FinalState ); // stable proton break; case KinematicsMode::InelasticElastic: case KinematicsMode::InelasticInelastic: op1.setStatus( Particle::Status::Unfragmented ); // fragmenting remnants op1.setMass( MX_ ); break; } //----- second outgoing proton Particle& op2 = event_->getOneByRole( Particle::OutgoingBeam2 ); op2.setMomentum( p5_lab_ ); switch ( cuts_.mode ) { case KinematicsMode::ElasticElastic: case KinematicsMode::InelasticElastic: default: op2.setStatus( Particle::Status::FinalState ); // stable proton break; case KinematicsMode::ElasticInelastic: case KinematicsMode::InelasticInelastic: op2.setStatus( Particle::Status::Unfragmented ); // fragmenting remnants op2.setMass( MY_ ); break; } //----- first incoming photon Particle& ph1 = event_->getOneByRole( Particle::Parton1 ); ph1.setMomentum( plab_ph1 ); //----- second incoming photon Particle& ph2 = event_->getOneByRole( Particle::Parton2 ); ph2.setMomentum( plab_ph2 ); Particles& central_system = event_->getByRole( Particle::CentralSystem ); //----- first outgoing lepton Particle& ol1 = central_system[0]; ol1.setPdgId( ol1.pdgId(), ransign ); ol1.setMomentum( p6_cm_ ); ol1.setStatus( Particle::Status::FinalState ); //----- second outgoing lepton Particle& ol2 = central_system[1]; ol2.setPdgId( ol2.pdgId(), -ransign ); ol2.setMomentum( p7_cm_ ); ol2.setStatus( Particle::Status::FinalState ); //----- intermediate two-lepton system event_->getOneByRole( Particle::Intermediate ).setMomentum( p6_cm_+p7_cm_ ); } //--------------------------------------------------------------------------------------------- double GamGamLL::periPP( int nup_, int ndown_ ) { CG_DEBUG_LOOP( "GamGamLL" ) << " Nup = " << nup_ << "\n\t" << "Ndown = " << ndown_; FormFactors fp1, fp2; formFactors( -t1_, -t2_, fp1, fp2 ); CG_DEBUG_LOOP( "GamGamLL" ) << "u1 = " << fp1.FM << "\n\t" << "u2 = " << fp1.FE << "\n\t" << "v1 = " << fp2.FM << "\n\t" << "v2 = " << fp2.FE; const double qqq = q1dq_*q1dq_, qdq = 4.*masses_.Ml2_-w4_; const double t11 = 64. *( bb_*( qqq-g4_-qdq*( t1_+t2_+2.*masses_.Ml2_ ) )-2.*( t1_+2.*masses_.Ml2_ )*( t2_+2.*masses_.Ml2_ )*qqq ) * t1_*t2_, // magnetic-magnetic t12 = 128.*( -bb_*( dd2_+g6_ )-2.*( t1_+2.*masses_.Ml2_ )*( sa2_*qqq+a6_*a6_ ) ) * t1_, // electric-magnetic t21 = 128.*( -bb_*( dd4_+g5_ )-2.*( t2_+2.*masses_.Ml2_ )*( sa1_*qqq+a5_*a5_ ) ) * t2_, // magnetic-electric t22 = 512.*( bb_*( delta_*delta_-gram_ )-pow( epsi_-delta_*( qdq+q1dq2_ ), 2 )-sa1_*a6_*a6_-sa2_*a5_*a5_-sa1_*sa2_*qqq ); // electric-electric const double peripp = ( fp1.FM*fp2.FM*t11 + fp1.FE*fp2.FM*t21 + fp1.FM*fp2.FE*t12 + fp1.FE*fp2.FE*t22 ) / pow( 2.*t1_*t2_*bb_, 2 ); CG_DEBUG_LOOP( "GamGamLL" ) << "t11 = " << t11 << "\t" << "t12 = " << t12 << "\n\t" << "t21 = " << t21 << "\t" << "t22 = " << t22 << "\n\t" << "=> PeriPP = " << peripp; return peripp; } void GamGamLL::map( double expo, const Limits& lim, double& out, double& dout, const std::string& var_name_ ) { const double y = lim.max()/lim.min(); out = lim.min()*pow( y, expo ); dout = out*log( y ); CG_DEBUG_LOOP( "GamGamLL:map" ) << "Mapping variable \"" << var_name_ << "\"\n\t" << "limits = " << lim << "\n\t" << "max/min = " << y << "\n\t" << "exponent = " << expo << "\n\t" << "output = " << out << "\n\t" << "d(output) = " << dout; } void GamGamLL::mapla( double y, double z, int u, double xm, double xp, double& x, double& d ) { const double xmb = xm-y-z, xpb = xp-y-z; const double c = -4.*y*z; const double alp = sqrt( xpb*xpb + c ), alm = sqrt( xmb*xmb + c ); const double am = xmb+alm, ap = xpb+alp; const double yy = ap/am, zz = pow( yy, u ); x = y + z + ( am*zz - c / ( am*zz ) ) / 2.; const double ax = sqrt( pow( x-y-z, 2 ) + c ); d = ax*log( yy ); } void GamGamLL::formFactors( double q1, double q2, FormFactors& fp1, FormFactors& fp2 ) const { const double mx2 = MX_*MX_, my2 = MY_*MY_; switch ( cuts_.mode ) { case KinematicsMode::ElectronElectron: default: { fp1 = FormFactors::trivial(); // electron (trivial) form factor fp2 = FormFactors::trivial(); // electron (trivial) form factor } break; case KinematicsMode::ProtonElectron: { fp1 = FormFactors::protonElastic( -t1_ ); // proton elastic form factor fp2 = FormFactors::trivial(); // electron (trivial) form factor } break; case KinematicsMode::ElectronProton: { fp1 = FormFactors::trivial(); // electron (trivial) form factor fp2 = FormFactors::protonElastic( -t2_ ); // proton elastic form factor } break; case KinematicsMode::ElasticElastic: { fp1 = FormFactors::protonElastic( -t1_ ); // proton elastic form factor fp2 = FormFactors::protonElastic( -t2_ ); // proton elastic form factor } break; case KinematicsMode::ElasticInelastic: { fp1 = FormFactors::protonElastic( -t1_ ); fp2 = FormFactors::protonInelastic( -t2_, w2_, my2, *cuts_.structure_functions ); } break; case KinematicsMode::InelasticElastic: { fp1 = FormFactors::protonInelastic( -t1_, w1_, mx2, *cuts_.structure_functions ); fp2 = FormFactors::protonElastic( -t2_ ); } break; case KinematicsMode::InelasticInelastic: { fp1 = FormFactors::protonInelastic( -t1_, w1_, mx2, *cuts_.structure_functions ); fp2 = FormFactors::protonInelastic( -t2_, w2_, my2, *cuts_.structure_functions ); } break; } } // register process and define aliases REGISTER_PROCESS( lpair, GamGamLL ) REGISTER_PROCESS( gamgamll, GamGamLL ) } } diff --git a/CepGen/Processes/GenericProcess.cpp b/CepGen/Processes/GenericProcess.cpp index 7f3f719..a003f82 100644 --- a/CepGen/Processes/GenericProcess.cpp +++ b/CepGen/Processes/GenericProcess.cpp @@ -1,228 +1,228 @@ #include "CepGen/Processes/GenericProcess.h" #include "CepGen/Core/Exception.h" #include "CepGen/Event/Event.h" #include "CepGen/Physics/ParticleProperties.h" #include "CepGen/Physics/Constants.h" #include "CepGen/Physics/FormFactors.h" #include "CepGen/Physics/PDG.h" namespace CepGen { namespace process { - const double GenericProcess::mp_ = ParticleProperties::mass( PDG::proton ); + const double GenericProcess::mp_ = part::mass( PDG::proton ); const double GenericProcess::mp2_ = GenericProcess::mp_*GenericProcess::mp_; GenericProcess::GenericProcess( const std::string& name, const std::string& description, bool has_event ) : name_( name ), description_( description ), first_run( true ), s_( -1. ), sqs_( -1. ), MX_( -1. ), MY_( -1. ), w1_( -1. ), w2_( -1. ), t1_( -1. ), t2_( -1. ), has_event_( has_event ), event_( new Event ), is_point_set_( false ) {} GenericProcess::GenericProcess( const GenericProcess& proc ) : name_( proc.name_ ), description_( proc.description_ ), first_run( proc.first_run ), s_( proc.s_ ), sqs_( proc.sqs_ ), MX_( proc.MX_ ), MY_( proc.MY_ ), w1_( proc.w1_ ), w2_( proc.w2_ ), t1_( -1. ), t2_( -1. ), cuts_( proc.cuts_ ), has_event_( proc.has_event_ ), event_( new Event( *proc.event_.get() ) ), is_point_set_( false ) {} GenericProcess& GenericProcess::operator=( const GenericProcess& proc ) { name_ = proc.name_; description_ = proc.description_; first_run = proc.first_run; s_ = proc.s_; sqs_ = proc.sqs_; MX_ = proc.MX_; MY_ = proc.MY_; w1_ = proc.w1_; w2_ = proc.w2_; cuts_ = proc.cuts_; has_event_ = proc.has_event_; event_.reset( new Event( *proc.event_.get() ) ); is_point_set_ = false; return *this; } void GenericProcess::setPoint( const unsigned int ndim, double* x ) { x_ = std::vector<double>( x, x+ndim ); is_point_set_ = true; if ( CG_EXCEPT_MATCH( "Process:dumpPoint", debugInsideLoop ) ) dumpPoint(); } double GenericProcess::x( unsigned int idx ) const { if ( idx >= x_.size() ) return -1.; return x_[idx]; } void GenericProcess::clearEvent() { event_->restore(); } void GenericProcess::setKinematics( const Kinematics& cuts ) { cuts_ = cuts; prepareKinematics(); } void GenericProcess::prepareKinematics() { if ( !isKinematicsDefined() ) throw CG_FATAL( "GenericProcess" ) << "Kinematics not properly defined for the process."; // at some point introduce non head-on colliding beams? Particle::Momentum p1( 0., 0., cuts_.incoming_beams.first.pz ), p2( 0., 0., -cuts_.incoming_beams.second.pz ); // on-shell beam particles - p1.setMass( ParticleProperties::mass( cuts_.incoming_beams.first.pdg ) ); - p2.setMass( ParticleProperties::mass( cuts_.incoming_beams.second.pdg ) ); + p1.setMass( part::mass( cuts_.incoming_beams.first.pdg ) ); + p2.setMass( part::mass( cuts_.incoming_beams.second.pdg ) ); setIncomingKinematics( p1, p2 ); sqs_ = CMEnergy( p1, p2 ); s_ = sqs_*sqs_; w1_ = p1.mass2(); w2_ = p2.mass2(); CG_DEBUG( "GenericProcess" ) << "Kinematics successfully prepared! sqrt(s) = " << sqs_ << "."; } void GenericProcess::dumpPoint() const { std::ostringstream os; for ( unsigned short i = 0; i < x_.size(); ++i ) { os << Form( " x(%2d) = %8.6f\n\t", i, x_[i] ); } CG_INFO( "GenericProcess" ) << "Number of integration parameters: " << x_.size() << "\n\t" << os.str() << "."; } void GenericProcess::setEventContent( const IncomingState& ini, const OutgoingState& fin ) { if ( !has_event_ ) return; event_->clear(); //----- add the particles in the event //--- incoming state for ( const auto& ip : ini ) { Particle& p = event_->addParticle( ip.first ); - p.setPdgId( ip.second, ParticleProperties::charge( ip.second ) ); - p.setMass( ParticleProperties::mass( ip.second ) ); + p.setPdgId( ip.second, part::charge( ip.second ) ); + p.setMass( part::mass( ip.second ) ); if ( ip.first == Particle::IncomingBeam1 || ip.first == Particle::IncomingBeam2 ) p.setStatus( Particle::Status::PrimordialIncoming ); if ( ip.first == Particle::Parton1 || ip.first == Particle::Parton2 ) p.setStatus( Particle::Status::Incoming ); } //--- central system (if not already there) const auto& central_system = ini.find( Particle::CentralSystem ); if ( central_system == ini.end() ) { Particle& p = event_->addParticle( Particle::Intermediate ); p.setPdgId( PDG::invalid ); p.setStatus( Particle::Status::Propagator ); } //--- outgoing state for ( const auto& opl : fin ) { // pair(role, list of PDGids) for ( const auto& pdg : opl.second ) { Particle& p = event_->addParticle( opl.first ); - p.setPdgId( pdg, ParticleProperties::charge( pdg ) ); - p.setMass( ParticleProperties::mass( pdg ) ); + p.setPdgId( pdg, part::charge( pdg ) ); + p.setMass( part::mass( pdg ) ); } } //----- define the particles parentage const Particles parts = event_->particles(); for ( const auto& p : parts ) { Particle& part = event_->operator[]( p.id() ); switch ( part.role() ) { case Particle::OutgoingBeam1: case Particle::Parton1: part.addMother( event_->getOneByRole( Particle::IncomingBeam1 ) ); break; case Particle::OutgoingBeam2: case Particle::Parton2: part.addMother( event_->getOneByRole( Particle::IncomingBeam2 ) ); break; case Particle::Intermediate: part.addMother( event_->getOneByRole( Particle::Parton1 ) ); part.addMother( event_->getOneByRole( Particle::Parton2 ) ); break; case Particle::CentralSystem: part.addMother( event_->getOneByRole( Particle::Intermediate ) ); break; default: break; } } //----- freeze the event as it is event_->freeze(); last_event = event_; } void GenericProcess::setIncomingKinematics( const Particle::Momentum& p1, const Particle::Momentum& p2 ) { if ( !has_event_ ) return; event_->getOneByRole( Particle::IncomingBeam1 ).setMomentum( p1 ); event_->getOneByRole( Particle::IncomingBeam2 ).setMomentum( p2 ); } bool GenericProcess::isKinematicsDefined() { if ( !has_event_ ) return true; // check the incoming state bool is_incoming_state_set = ( !event_->getByRole( Particle::IncomingBeam1 ).empty() && !event_->getByRole( Particle::IncomingBeam2 ).empty() ); // check the outgoing state bool is_outgoing_state_set = ( !event_->getByRole( Particle::OutgoingBeam1 ).empty() && !event_->getByRole( Particle::OutgoingBeam2 ).empty() && !event_->getByRole( Particle::CentralSystem ).empty() ); // combine both states return is_incoming_state_set && is_outgoing_state_set; } std::ostream& operator<<( std::ostream& os, const GenericProcess& proc ) { return os << proc.name().c_str(); } std::ostream& operator<<( std::ostream& os, const GenericProcess* proc ) { return os << proc->name().c_str(); } } } diff --git a/CepGen/Processes/PPtoFF.cpp b/CepGen/Processes/PPtoFF.cpp index ab21577..671338f 100644 --- a/CepGen/Processes/PPtoFF.cpp +++ b/CepGen/Processes/PPtoFF.cpp @@ -1,378 +1,378 @@ #include "CepGen/Processes/PPtoFF.h" #include "CepGen/Event/Event.h" #include "CepGen/Physics/Constants.h" #include "CepGen/Physics/FormFactors.h" #include "CepGen/Physics/PDG.h" #include "CepGen/Core/Exception.h" #include <iomanip> #include "CepGen/Processes/ProcessesHandler.h" namespace CepGen { namespace process { PPtoFF::PPtoFF( const ParametersList& params ) : GenericKTProcess( params, "pptoff", "ɣɣ → f⁺f¯", { { PDG::photon, PDG::photon } }, { PDG::muon, PDG::muon } ), pair_( params.get<int>( "pair", 0 ) ), method_( params.get<int>( "method", 1 ) ), y1_( 0. ), y2_( 0. ), pt_diff_( 0. ), phi_pt_diff_( 0. ) {} void PPtoFF::preparePhaseSpace() { registerVariable( y1_, Mapping::linear, cuts_.cuts.central.rapidity_single, { -6., 6. }, "First outgoing fermion rapidity" ); registerVariable( y2_, Mapping::linear, cuts_.cuts.central.rapidity_single, { -6., 6. }, "Second outgoing fermion rapidity" ); registerVariable( pt_diff_, Mapping::linear, cuts_.cuts.central.pt_diff, { 0., 50. }, "Fermions transverse momentum difference" ); registerVariable( phi_pt_diff_, Mapping::linear, cuts_.cuts.central.phi_pt_diff, { 0., 2.*M_PI }, "Fermions azimuthal angle difference" ); if ( (PDG)pair_ == PDG::invalid ) throw CG_FATAL( "PPtoFF:prepare" ) << "Invalid fermion pair selected: " << pair_ << "!"; const PDG pdg_f = (PDG)pair_; - mf_ = ParticleProperties::mass( pdg_f ); + mf_ = part::mass( pdg_f ); mf2_ = mf_*mf_; - qf_ = ParticleProperties::charge( pdg_f ); - colf_ = ParticleProperties::colours( pdg_f ); + qf_ = part::charge( pdg_f ); + colf_ = part::colours( pdg_f ); CG_DEBUG( "PPtoFF:prepare" ) << "Produced particles (" << pdg_f << ") " << "with mass = " << mf_ << " GeV, " << "and charge = " << std::setprecision( 2 ) << qf_ << " e"; CG_DEBUG( "PPtoFF:mode" ) << "matrix element computation method: " << method_ << "."; } double PPtoFF::computeKTFactorisedMatrixElement() { //================================================================= // matrix element computation //================================================================= //--- central partons const Particle::Momentum q1t( qt1_*cos( phi_qt1_ ), qt1_*sin( phi_qt1_ ), 0. ); const Particle::Momentum q2t( qt2_*cos( phi_qt2_ ), qt2_*sin( phi_qt2_ ), 0. ); CG_DEBUG_LOOP( "PPtoFF" ) << "q(1/2)t = " << q1t << ", " << q2t << "."; //--- two-parton system const Particle::Momentum ptsum = q1t+q2t; const Particle::Momentum ptdiff( pt_diff_*cos( phi_pt_diff_ ), pt_diff_*sin( phi_pt_diff_ ), 0. ); //--- outgoing fermions const Particle::Momentum p1_cm = 0.5*( ptsum+ptdiff ), p2_cm = 0.5*( ptsum-ptdiff ); //================================================================= // a window in single particle transverse momentum //================================================================= const Limits& pt_limits = cuts_.cuts.central.pt_single; if ( !pt_limits.passes( p1_cm.pt() ) || !pt_limits.passes( p2_cm.pt() ) ) return 0.; //================================================================= // a window in transverse momentum difference //================================================================= if ( !cuts_.cuts.central.pt_diff.passes( fabs( p1_cm.pt()-p2_cm.pt() ) ) ) return 0.; //================================================================= // a window in rapidity distance //================================================================= if ( !cuts_.cuts.central.rapidity_diff.passes( fabs( y1_-y2_ ) ) ) return 0.; //================================================================= // auxiliary quantities //================================================================= // transverse mass for the two fermions const double amt1 = std::hypot( p1_cm.pt(), mf_ ), amt2 = std::hypot( p2_cm.pt(), mf_ ); const double alpha1 = amt1/sqs_*exp( y1_ ), beta1 = amt1/sqs_*exp( -y1_ ), alpha2 = amt2/sqs_*exp( y2_ ), beta2 = amt2/sqs_*exp( -y2_ ); CG_DEBUG_LOOP( "PPtoFF" ) << "Sudakov parameters:\n\t" << " alpha(1/2) = " << alpha1 << ", " << alpha2 << "\n\t" << " beta(1/2) = " << beta1 << ", " << beta2 << "."; const double x1 = alpha1+alpha2, x2 = beta1+beta2; if ( x1 <= 0. || x1 > 1. || x2 <= 0. || x2 > 1. ) return 0.; // sanity check //================================================================= // additional conditions for energy-momentum conservation //================================================================= const double s1_eff = x1*s_-qt1_*qt1_, s2_eff = x2*s_-qt2_*qt2_; const double invm = sqrt( amt1*amt1 + amt2*amt2 + 2.*amt1*amt2*cosh(y1_-y2_) - ptsum.pt2() ); CG_DEBUG_LOOP( "PPtoFF" ) << "s(1/2)eff = " << s1_eff << ", " << s2_eff << " GeV²\n\t" << "central system's invariant mass = " << invm << " GeV."; if ( ( cuts_.mode == KinematicsMode::ElasticInelastic || cuts_.mode == KinematicsMode::InelasticInelastic ) && ( sqrt( s1_eff ) <= ( MY_+invm ) ) ) return 0.; if ( ( cuts_.mode == KinematicsMode::InelasticElastic || cuts_.mode == KinematicsMode::InelasticInelastic ) && ( sqrt( s2_eff ) <= ( MX_+invm ) ) ) return 0.; //================================================================= // four-momenta of the outgoing protons (or remnants) //================================================================= const Particle::Momentum& ak1 = event_->getOneByRole( Particle::IncomingBeam1 ).momentum(), &ak2 = event_->getOneByRole( Particle::IncomingBeam2 ).momentum(); CG_DEBUG_LOOP( "PPtoFF" ) << "incoming particles: p(1/2) = " << ak1 << ", " << ak2 << "."; const double px_plus = ( 1.-x1 )*M_SQRT2*ak1.p(), px_minus = ( MX_*MX_ + q1t.pt2() )*0.5/px_plus; const double py_minus = ( 1.-x2 )*M_SQRT2*ak2.p(), py_plus = ( MY_*MY_ + q2t.pt2() )*0.5/py_minus; CG_DEBUG_LOOP( "PPtoFF" ) << "px± = " << px_plus << ", " << px_minus << "\n\t" << "py± = " << py_plus << ", " << py_minus << "."; PX_ = Particle::Momentum( -q1t.px(), -q1t.py(), ( px_plus-px_minus )*M_SQRT1_2, ( px_plus+px_minus )*M_SQRT1_2 ); PY_ = Particle::Momentum( -q2t.px(), -q2t.py(), ( py_plus-py_minus )*M_SQRT1_2, ( py_plus+py_minus )*M_SQRT1_2 ); CG_DEBUG_LOOP( "PPtoFF" ) << "First remnant: " << PX_ << ", mass = " << PX_.mass() << "\n\t" << "Second remnant: " << PY_ << ", mass = " << PY_.mass() << "."; if ( fabs( PX_.mass()-MX_ ) > 1.e-6 ) throw CG_FATAL( "PPtoFF" ) << "Invalid X system mass: " << PX_.mass() << "/" << MX_ << "."; if ( fabs( PY_.mass()-MY_ ) > 1.e-6 ) throw CG_FATAL( "PPtoFF" ) << "Invalid Y system mass: " << PY_.mass() << "/" << MY_ << "."; //================================================================= // four-momenta of the outgoing l^+ and l^- //================================================================= const Particle::Momentum p1 = p1_cm + alpha1*ak1 + beta1*ak2; const Particle::Momentum p2 = p2_cm + alpha2*ak1 + beta2*ak2; CG_DEBUG_LOOP( "PPtoFF" ) << "unboosted first fermion: " << p1 << ", mass = " << p1.mass() << "\n\t" << " second fermion: " << p2 << ", mass = " << p2.mass() << "."; p_f1_ = Particle::Momentum::fromPxPyYM( p1_cm.px(), p1_cm.py(), y2_, mf_ ); p_f2_ = Particle::Momentum::fromPxPyYM( p2_cm.px(), p2_cm.py(), y1_, mf_ ); CG_DEBUG_LOOP( "PPtoFF" ) << "First fermion: " << p_f1_ << ", mass = " << p_f1_.mass() << "\n\t" << "Second fermion: " << p_f2_ << ", mass = " << p_f2_.mass() << "."; if ( fabs( p_f1_.mass()-mf_ ) > 1.e-4 ) throw CG_FATAL( "PPtoFF" ) << "Invalid fermion 1 mass: " << p_f1_.mass() << "/" << mf_ << "."; if ( fabs( p_f2_.mass()-mf_ ) > 1.e-4 ) throw CG_FATAL( "PPtoFF" ) << "Invalid fermion 2 mass: " << p_f2_.mass() << "/" << mf_ << "."; //================================================================= // matrix elements //================================================================= double amat2 = 0.; //================================================================= // How matrix element is calculated //================================================================= if ( method_ == 0 ) { //=============================================================== // Mendelstam variables //=============================================================== //const double shat = s_*x1*x2; // approximation const double shat = ( q1t+q2t ).mass2(); // exact formula //const double mll = sqrt( shat ); const double that1 = ( q1t-p1 ).mass2(), that2 = ( q2t-p2 ).mass2(); const double uhat1 = ( q1t-p2 ).mass2(), uhat2 = ( q2t-p1 ).mass2(); const double that = 0.5*( that1+that2 ), uhat = 0.5*( uhat1+uhat2 ); amat2 = onShellME( shat, that, uhat ); CG_DEBUG_LOOP( "PPtoFF:onShell" ) << "that(1/2) = " << that1 << " / " << that2 << "\n\t" << "uhat(1/2) = " << uhat1 << " / " << uhat2 << "\n\t" << "squared matrix element: " << amat2 << "."; } else if ( method_ == 1 ) { const double t1abs = ( q1t.pt2() + x1*( MX_*MX_-mp2_ )+x1*x1*mp2_ )/( 1.-x1 ), t2abs = ( q2t.pt2() + x2*( MY_*MY_-mp2_ )+x2*x2*mp2_ )/( 1.-x2 ); const double z1p = alpha1/x1, z1m = alpha2/x1, z2p = beta1 /x2, z2m = beta2 /x2; CG_DEBUG_LOOP( "PPtoFF:offShell" ) << "z(1/2)p = " << z1p << ", " << z2p << "\n\t" << "z(1/2)m = " << z1m << ", " << z2m << "."; amat2 = offShellME( t1abs, t2abs, z1m, z1p, z2m, z2p, q1t, q2t ) * pow( x1*x2*s_, 2 ); } //============================================ // unintegrated photon distributions //============================================ const std::pair<double,double> fluxes = GenericKTProcess::incomingFluxes( x1, q1t.pt2(), x2, q2t.pt2() ); CG_DEBUG_LOOP( "PPtoFF" ) << "Incoming photon fluxes for (x/kt2) = " << "(" << x1 << "/" << q1t.pt2() << "), " << "(" << x2 << "/" << q2t.pt2() << "):\n\t" << fluxes.first << ", " << fluxes.second << "."; //================================================================= // factor 2.*pi from integration over phi_sum // factor 1/4 from jacobian of transformations // factors 1/pi and 1/pi due to integration over // d^2 kappa_1 d^2 kappa_2 instead d kappa_1^2 d kappa_2^2 //================================================================= - const double g_em = 4.*M_PI*Constants::alphaEM*qf_*qf_; + const double g_em = 4.*M_PI*constants::alphaEM*qf_*qf_; const double aintegral = amat2 * colf_ * ( g_em*g_em ) * 1. / pow( 4.*M_PI*( x1*x2*s_ ), 2 ) * fluxes.first*M_1_PI * fluxes.second*M_1_PI * 0.25 - * Constants::GeV2toBarn; + * constants::GeV2toBarn; //================================================================= return aintegral*qt1_*qt2_*pt_diff_; //================================================================= } void PPtoFF::fillCentralParticlesKinematics() { // randomise the charge of the outgoing fermions short sign = ( drand() > 0.5 ) ? +1 : -1; //================================================================= // first outgoing fermion //================================================================= Particle& of1 = event_->getByRole( Particle::CentralSystem )[0]; of1.setPdgId( of1.pdgId(), sign ); of1.setStatus( Particle::Status::FinalState ); of1.setMomentum( p_f1_ ); //================================================================= // second outgoing fermion //================================================================= Particle& of2 = event_->getByRole( Particle::CentralSystem )[1]; of2.setPdgId( of2.pdgId(), -sign ); of2.setStatus( Particle::Status::FinalState ); of2.setMomentum( p_f2_ ); } double PPtoFF::onShellME( double shat, double that, double uhat ) const { CG_DEBUG_LOOP( "PPtoFF:onShell" ) << "shat: " << shat << ", that: " << that << ", uhat: " << uhat << "."; //================================================================= // on-shell formula for M^2 //================================================================= const double ml4 = mf2_*mf2_, ml8 = ml4*ml4; const double term1 = 6. *ml8, term2 = -3. *ml4 *that*that, term3 = -14.*ml4 *that*uhat, term4 = -3. *ml4 *uhat*uhat, term5 = mf2_*that*that*that, term6 = 7.* mf2_*that*that*uhat, term7 = 7.* mf2_*that*uhat*uhat, term8 = mf2_*uhat*uhat*uhat, term9 = -that*that*that*uhat, term10 = -that*uhat*uhat*uhat; return -2.*( term1+term2+term3+term4+term5 +term6+term7+term8+term9+term10 )/( pow( ( mf2_-that )*( mf2_-uhat ), 2) ); } double PPtoFF::offShellME( double t1abs, double t2abs, double z1m, double z1p, double z2m, double z2p, const Particle::Momentum& q1, const Particle::Momentum& q2 ) const { //================================================================= // Wolfgang's formulae //================================================================= const double z1 = z1p*z1m, z2 = z2p*z2m; const double eps12 = mf2_+z1*t1abs, eps22 = mf2_+z2*t2abs; const Particle::Momentum ak1 = ( z1m*p_f1_-z1p*p_f2_ ), ak2 = ( z2m*p_f1_-z2p*p_f2_ ); const Particle::Momentum ph_p1 = ak1+z1p*q2, ph_m1 = ak1-z1m*q2; const Particle::Momentum ph_p2 = ak2+z2p*q1, ph_m2 = ak2-z2m*q1; const Particle::Momentum phi1( ph_p1.px()/( ph_p1.pt2()+eps12 )-ph_m1.px()/( ph_m1.pt2()+eps12 ), ph_p1.py()/( ph_p1.pt2()+eps12 )-ph_m1.py()/( ph_m1.pt2()+eps12 ), 0., 1./( ph_p1.pt2()+eps12 )-1./( ph_m1.pt2()+eps12 ) ); const Particle::Momentum phi2( ph_p2.px()/( ph_p2.pt2()+eps22 )-ph_m2.px()/( ph_m2.pt2()+eps22 ), ph_p2.py()/( ph_p2.pt2()+eps22 )-ph_m2.py()/( ph_m2.pt2()+eps22 ), 0., 1./( ph_p2.pt2()+eps22 )-1./( ph_m2.pt2()+eps22 ) ); const double dot1 = phi1.threeProduct( q1 )/qt1_, cross1 = phi1.crossProduct( q1 )/qt1_; const double dot2 = phi2.threeProduct( q2 )/qt2_, cross2 = phi2.crossProduct( q2 )/qt2_; CG_DEBUG_LOOP( "PPtoFF:offShell" ) << "phi1 = " << phi1 << "\n\t" << "phi2 = " << phi2 << "\n\t" << "(dot): " << dot1 << " / " << dot2 << "\n\t" << "(cross): " << cross1 << " / " << cross2 << "."; //================================================================= // six terms in Wolfgang's formula for // off-shell gamma gamma --> l^+ l^- //================================================================= const unsigned short imat1 = 1, imat2 = 1; const unsigned short itermLL = 1, itermTT = 1, itermLT = 1, itermtt = 1; const double aux2_1 = itermLL * ( mf2_ + 4.*z1*z1*t1abs ) * phi1.energy2() +itermTT * ( ( z1p*z1p + z1m*z1m )*( dot1*dot1 + cross1*cross1 ) ) +itermtt * ( cross1*cross1 - dot1*dot1 ) -itermLT * 4.*z1*( z1p-z1m ) * phi1.energy() * q1.threeProduct( phi1 ); const double aux2_2 = itermLL * ( mf2_ + 4.*z2*z2*t2abs ) * phi2.energy2() +itermTT * ( ( z2p*z2p + z2m*z2m )*( dot2*dot2 + cross2*cross2 ) ) +itermtt * ( cross2*cross2 - dot2*dot2 ) -itermLT * 4.*z2*( z2p-z2m ) * phi2.energy() * q2.threeProduct( phi2 ); //================================================================= // convention of matrix element as in our kt-factorization // for heavy flavours //================================================================= const double amat2_1 = aux2_1*2.*z1*q1.pt2()/( q1.pt2()*q2.pt2() ), amat2_2 = aux2_2*2.*z2*q2.pt2()/( q1.pt2()*q2.pt2() ); //================================================================= // symmetrization //================================================================= CG_DEBUG_LOOP( "PPtoFF:offShell" ) << "aux2(1/2) = " << aux2_1 << " / " << aux2_2 << "\n\t" << "amat2(1/2), amat2 = " << amat2_1 << " / " << amat2_2 << " / " << ( 0.5*( imat1*amat2_1 + imat2*amat2_2 ) ) << "."; return 0.5*( imat1*amat2_1 + imat2*amat2_2 ); } // register process and define aliases REGISTER_PROCESS( pptoll, PPtoFF ) REGISTER_PROCESS( pptoff, PPtoFF ) } } diff --git a/CepGen/Processes/PPtoWW.cpp b/CepGen/Processes/PPtoWW.cpp index 4ccf458..96eddfe 100644 --- a/CepGen/Processes/PPtoWW.cpp +++ b/CepGen/Processes/PPtoWW.cpp @@ -1,393 +1,393 @@ #include "CepGen/Processes/PPtoWW.h" #include "CepGen/Event/Event.h" #include "CepGen/Physics/Constants.h" #include "CepGen/Physics/FormFactors.h" #include "CepGen/Physics/PDG.h" #include "CepGen/Core/Exception.h" #include <assert.h> #include "CepGen/Processes/ProcessesHandler.h" namespace CepGen { namespace process { - const double PPtoWW::mw_ = ParticleProperties::mass( PDG::W ); + const double PPtoWW::mw_ = part::mass( PDG::W ); const double PPtoWW::mw2_ = PPtoWW::mw_*PPtoWW::mw_; PPtoWW::PPtoWW( const ParametersList& params ) : GenericKTProcess( params, "pptoww", "ɣɣ → W⁺W¯", { { PDG::photon, PDG::photon } }, { PDG::W, PDG::W } ), method_( params.get<int>( "method", 1 ) ), pol_state_( (Polarisation)params.get<int>( "polarisationStates", 0 ) ), y1_( 0. ), y2_( 0. ), pt_diff_( 0. ), phi_pt_diff_( 0. ) {} void PPtoWW::preparePhaseSpace() { registerVariable( y1_, Mapping::linear, cuts_.cuts.central.rapidity_single, { -6., 6. }, "First outgoing W rapidity" ); registerVariable( y2_, Mapping::linear, cuts_.cuts.central.rapidity_single, { -6., 6. }, "Second outgoing W rapidity" ); registerVariable( pt_diff_, Mapping::linear, cuts_.cuts.central.pt_diff, { 0., 500. }, "Ws transverse momentum difference" ); registerVariable( phi_pt_diff_, Mapping::linear, cuts_.cuts.central.phi_pt_diff, { 0., 2.*M_PI }, "Ws azimuthal angle difference" ); switch ( pol_state_ ) { case Polarisation::LL: pol_w1_ = pol_w2_ = { 0 }; break; case Polarisation::LT: pol_w1_ = { 0 }; pol_w2_ = { -1, 1 }; break; case Polarisation::TL: pol_w1_ = { -1, 1 }; pol_w2_ = { 0 }; break; case Polarisation::TT: pol_w1_ = pol_w2_ = { -1, 1 }; break; default: case Polarisation::full: pol_w1_ = pol_w2_ = { -1, 0, 1 }; break; } CG_DEBUG( "PPtoWW:mode" ) << "matrix element computation method: " << method_ << "."; } double PPtoWW::computeKTFactorisedMatrixElement() { //================================================================= // matrix element computation //================================================================= //const double stild = s_/2.*(1+sqrt(1.-(4*pow(mp2_, 2))/s_*s_)); // Inner photons const double q1tx = qt1_*cos( phi_qt1_ ), q1ty = qt1_*sin( phi_qt1_ ), q2tx = qt2_*cos( phi_qt2_ ), q2ty = qt2_*sin( phi_qt2_ ); CG_DEBUG_LOOP( "PPtoWW:qt" ) << "q1t(x/y) = " << q1tx << " / " << q1ty << "\n\t" << "q2t(x/y) = " << q2tx << " / " << q2ty << "."; // Two-photon system const double ptsumx = q1tx+q2tx, ptsumy = q1ty+q2ty, ptsum = sqrt( ptsumx*ptsumx+ptsumy*ptsumy ); const double ptdiffx = pt_diff_*cos( phi_pt_diff_ ), ptdiffy = pt_diff_*sin( phi_pt_diff_ ); // Outgoing leptons const double pt1x = ( ptsumx+ptdiffx )*0.5, pt1y = ( ptsumy+ptdiffy )*0.5, pt1 = std::hypot( pt1x, pt1y ), pt2x = ( ptsumx-ptdiffx )*0.5, pt2y = ( ptsumy-ptdiffy )*0.5, pt2 = std::hypot( pt2x, pt2y ); if ( cuts_.cuts.central_particles.count( PDG::W ) > 0 && cuts_.cuts.central_particles.at( PDG::W ).pt_single.valid() ) { const Limits pt_limits = cuts_.cuts.central_particles.at( PDG::W ).pt_single; if ( !pt_limits.passes( pt1 ) || !pt_limits.passes( pt2 ) ) return 0.; } // transverse mass for the two leptons const double amt1 = sqrt( pt1*pt1+mw2_ ), amt2 = sqrt( pt2*pt2+mw2_ ); //================================================================= // a window in two-boson invariant mass //================================================================= const double invm = sqrt( amt1*amt1 + amt2*amt2 + 2.*amt1*amt2*cosh( y1_-y2_ ) - ptsum*ptsum ); if ( !cuts_.cuts.central.mass_sum.passes( invm ) ) return 0.; //================================================================= // a window in transverse momentum difference //================================================================= if ( !cuts_.cuts.central.pt_diff.passes( fabs( pt1-pt2 ) ) ) return 0.; //================================================================= // a window in rapidity distance //================================================================= if ( !cuts_.cuts.central.rapidity_diff.passes( fabs( y1_-y2_ ) ) ) return 0.; //================================================================= // auxiliary quantities //================================================================= const double alpha1 = amt1/sqs_*exp( y1_ ), beta1 = amt1/sqs_*exp( -y1_ ), alpha2 = amt2/sqs_*exp( y2_ ), beta2 = amt2/sqs_*exp( -y2_ ); CG_DEBUG_LOOP( "PPtoWW:sudakov" ) << "Sudakov parameters:\n\t" << " alpha1/2 = " << alpha1 << " / " << alpha2 << "\n\t" << " beta1/2 = " << beta1 << " / " << beta2 << "."; const double q1t2 = q1tx*q1tx+q1ty*q1ty, q2t2 = q2tx*q2tx+q2ty*q2ty; const double x1 = alpha1+alpha2, x2 = beta1+beta2; const double z1p = alpha1/x1, z1m = alpha2/x1, z2p = beta1 /x2, z2m = beta2 /x2; CG_DEBUG_LOOP( "PPtoWW:zeta" ) << "z(1/2)p = " << z1p << " / " << z2p << "\n\t" << "z(1/2)m = " << z1m << " / " << z2m << "."; if ( x1 > 1. || x2 > 1. ) return 0.; // sanity check // FIXME FIXME FIXME const double ak10 = event_->getOneByRole( Particle::IncomingBeam1 ).energy(), ak1z = event_->getOneByRole( Particle::IncomingBeam1 ).momentum().pz(), ak20 = event_->getOneByRole( Particle::IncomingBeam2 ).energy(), ak2z = event_->getOneByRole( Particle::IncomingBeam2 ).momentum().pz(); CG_DEBUG_LOOP( "PPtoWW:incoming" ) << "incoming particles: p1: " << ak1z << " / " << ak10 << "\n\t" << " p2: " << ak2z << " / " << ak20 << "."; //================================================================= // additional conditions for energy-momentum conservation //================================================================= const double s1_eff = x1*s_-qt1_*qt1_, s2_eff = x2*s_-qt2_*qt2_; CG_DEBUG_LOOP( "PPtoWW:central" ) << "s(1/2)_eff = " << s1_eff << " / " << s2_eff << " GeV^2\n\t" << "diboson invariant mass = " << invm << " GeV"; if ( ( cuts_.mode == KinematicsMode::ElasticInelastic || cuts_.mode == KinematicsMode::InelasticInelastic ) && ( sqrt( s1_eff ) <= ( MY_+invm ) ) ) return 0.; if ( ( cuts_.mode == KinematicsMode::InelasticElastic || cuts_.mode == KinematicsMode::InelasticInelastic ) && ( sqrt( s2_eff ) <= ( MX_+invm ) ) ) return 0.; //const double qcaptx = pcaptx, qcapty = pcapty; //================================================================= // four-momenta of the outgoing protons (or remnants) //================================================================= const double px_plus = ( 1.-x1 )*fabs( ak1z )*M_SQRT2, px_minus = ( MX_*MX_ + q1t2 )*0.5/px_plus; const double py_minus = ( 1.-x2 )*fabs( ak2z )*M_SQRT2, // warning! sign of pz?? py_plus = ( MY_*MY_ + q2t2 )*0.5/py_minus; CG_DEBUG_LOOP( "PPtoWW:pxy" ) << "px± = " << px_plus << " / " << px_minus << "\n\t" << "py± = " << py_plus << " / " << py_minus << "."; PX_ = Particle::Momentum( -q1tx, -q1ty, ( px_plus-px_minus )*M_SQRT1_2, ( px_plus+px_minus )*M_SQRT1_2 ); PY_ = Particle::Momentum( -q2tx, -q2ty, ( py_plus-py_minus )*M_SQRT1_2, ( py_plus+py_minus )*M_SQRT1_2 ); CG_DEBUG_LOOP( "PPtoWW:remnants" ) << "First remnant: " << PX_ << ", mass = " << PX_.mass() << "\n\t" << "Second remnant: " << PY_ << ", mass = " << PY_.mass() << "."; /*assert( fabs( PX_.mass()-MX_ ) < 1.e-6 ); assert( fabs( PY_.mass()-MY_ ) < 1.e-6 );*/ //================================================================= // four-momenta squared of the virtual photons //================================================================= const double ww = 0.5 * ( 1.+sqrt( 1.-4.*mp2_/s_ ) ); // FIXME FIXME FIXME ///////////////////// const Particle::Momentum q1( q1tx, q1ty, +0.5 * x1*ww*sqs_*( 1.-q1t2/x1/x1/ww/ww/s_ ), +0.5 * x1*ww*sqs_*( 1.+q1t2/x1/x1/ww/ww/s_ ) ); const Particle::Momentum q2( q2tx, q2ty, -0.5 * x2*ww*sqs_*( 1.-q2t2/x2/x2/ww/ww/s_ ), +0.5 * x2*ww*sqs_*( 1.+q2t2/x2/x2/ww/ww/s_ ) ); ////////////////////////////////////////// CG_DEBUG_LOOP( "PPtoWW:partons" ) << "First photon*: " << q1 << ", mass2 = " << q1.mass2() << "\n\t" << "Second photon*: " << q2 << ", mass2 = " << q2.mass2() << "."; //const double q12 = q1.mass2(), q22 = q2.mass2(); //================================================================= // four-momenta of the outgoing W^+ and W^- //================================================================= p_w1_ = Particle::Momentum( pt1x, pt1y, alpha1*ak1z + beta1*ak2z, alpha1*ak10 + beta1*ak20 ); p_w2_ = Particle::Momentum( pt2x, pt2y, alpha2*ak1z + beta2*ak2z, alpha2*ak10 + beta2*ak20 ); CG_DEBUG_LOOP( "PPtoWW:central" ) << "First W: " << p_w1_ << ", mass = " << p_w1_.mass() << "\n\t" << "Second W: " << p_w2_ << ", mass = " << p_w2_.mass() << "."; //assert( fabs( p_w1_.mass()-event_->getByRole( Particle::CentralSystem )[0].mass() ) < 1.e-6 ); //assert( fabs( p_w2_.mass()-event_->getByRole( Particle::CentralSystem )[1].mass() ) < 1.e-6 ); //================================================================= // Mendelstam variables //================================================================= //const double shat = s_*x1*x2; // ishat = 1 (approximation) const double shat = ( q1+q2 ).mass2(); // ishat = 2 (exact formula) const double that1 = ( q1-p_w1_ ).mass2(), that2 = ( q2-p_w2_ ).mass2(); const double uhat1 = ( q1-p_w2_ ).mass2(), uhat2 = ( q2-p_w1_ ).mass2(); CG_DEBUG_LOOP( "PPtoWW" ) << "that(1/2) = " << that1 << " / " << that2 << "\n\t" << "uhat(1/2) = " << uhat1 << " / " << uhat2 << "."; //const double mll = sqrt( shat ); const double that = 0.5*( that1+that2 ), uhat = 0.5*( uhat1+uhat2 ); //================================================================= // matrix elements //================================================================= double amat2 = 0.; //================================================================= // How matrix element is calculated //================================================================= CG_DEBUG_LOOP( "PPtoWW" ) << "matrix element mode: " << method_ << "."; //================================================================= // matrix element for gamma gamma --> W^+ W^- // (Denner+Dittmaier+Schuster) // (work in collaboration with C. Royon) //================================================================= if ( method_ == 0 ) amat2 = onShellME( shat, that, uhat ); //================================================================= // off-shell Nachtmann formulae //================================================================= else if ( method_ == 1 ) amat2 = offShellME( shat, that, uhat, phi_qt1_+phi_qt2_, phi_qt1_-phi_qt2_ ); if ( amat2 <= 0. ) return 0.; //============================================ // unintegrated photon distributions //============================================ const std::pair<double,double> fluxes = GenericKTProcess::incomingFluxes( x1, q1t2, x2, q2t2 ); CG_DEBUG_LOOP( "PPtoWW:fluxes" ) << "Incoming photon fluxes for (x/kt2) = " << "(" << x1 << "/" << q1t2 << "), " << "(" << x2 << "/" << q2t2 << "):\n\t" << fluxes.first << ", " << fluxes.second << "."; //================================================================= // factor 2.*pi from integration over phi_sum // factor 1/4 from jacobian of transformations // factors 1/pi and 1/pi due to integration over // d^2 kappa_1 d^2 kappa_2 instead d kappa_1^2 d kappa_2^2 //================================================================= const double aintegral = amat2 / ( 16.*M_PI*M_PI*( x1*x2*s_ )*( x1*x2*s_ ) ) * fluxes.first*M_1_PI * fluxes.second*M_1_PI * 0.25 - * Constants::GeV2toBarn; + * constants::GeV2toBarn; /*const double aintegral = amat2 / ( 16.*M_PI*M_PI*x1*x1*x2*x2*s_*s_ ) * fluxes.first*M_1_PI * fluxes.second*M_1_PI - * Constants::GeV2toBarn * 0.25;*/ + * constants::GeV2toBarn * 0.25;*/ //================================================================= return aintegral*qt1_*qt2_*pt_diff_; //================================================================= } void PPtoWW::fillCentralParticlesKinematics() { // randomise the charge of the outgoing leptons short sign = ( drand() > 0.5 ) ? +1 : -1; //================================================================= // first outgoing lepton //================================================================= Particle& ow1 = event_->getByRole( Particle::CentralSystem )[0]; ow1.setPdgId( ow1.pdgId(), sign ); ow1.setStatus( Particle::Status::Undecayed ); ow1.setMomentum( p_w1_ ); //================================================================= // second outgoing lepton //================================================================= Particle& ow2 = event_->getByRole( Particle::CentralSystem )[1]; ow2.setPdgId( ow2.pdgId(), -sign ); ow2.setStatus( Particle::Status::Undecayed ); ow2.setMomentum( p_w2_ ); } double PPtoWW::onShellME( double shat, double that, double uhat ) { const double mw4 = mw2_*mw2_; const double term1 = 2.*shat * ( 2.*shat+3.*mw2_ ) / ( 3.*( mw2_-that )*( mw2_-uhat ) ); const double term2 = 2.*shat*shat * ( shat*shat + 3.*mw4 ) / ( 3.*pow( mw2_-that, 2 )*pow( mw2_-uhat, 2 ) ); const double auxil_gamgam = 1.-term1+term2; const double beta = sqrt( 1.-4.*mw2_/shat ); - return 3.*Constants::alphaEM*Constants::alphaEM*beta / ( 2.*shat ) * auxil_gamgam / ( beta/( 64.*M_PI*M_PI*shat ) ); + return 3.*constants::alphaEM*constants::alphaEM*beta / ( 2.*shat ) * auxil_gamgam / ( beta/( 64.*M_PI*M_PI*shat ) ); } double PPtoWW::offShellME( double shat, double that, double uhat, double phi_sum, double phi_diff ) { - const double e2 = 4.*M_PI*Constants::alphaEM; + const double e2 = 4.*M_PI*constants::alphaEM; double amat2_0 = 0., amat2_1 = 0., amat2_interf = 0.; for ( const auto lam3 : pol_w1_ ) for ( const auto lam4 : pol_w2_ ) { double ampli_pp = amplitudeWW( shat, that, uhat, +1, +1, lam3, lam4 ); double ampli_mm = amplitudeWW( shat, that, uhat, -1, -1, lam3, lam4 ); double ampli_pm = amplitudeWW( shat, that, uhat, +1, -1, lam3, lam4 ); double ampli_mp = amplitudeWW( shat, that, uhat, -1, +1, lam3, lam4 ); amat2_0 += ampli_pp*ampli_pp + ampli_mm*ampli_mm + 2.*cos( 2.*phi_diff )*ampli_pp*ampli_mm; amat2_1 += ampli_pm*ampli_pm + ampli_mp*ampli_mp + 2.*cos( 2.*phi_sum )*ampli_pm*ampli_mp; amat2_interf -= 2.*( cos( phi_sum+phi_diff )*( ampli_pp*ampli_pm+ampli_mm*ampli_mp ) +cos( phi_sum-phi_diff )*( ampli_pp*ampli_mp+ampli_mm*ampli_pm ) ); } return e2*e2*( amat2_0+amat2_1+amat2_interf ); } double PPtoWW::amplitudeWW( double shat, double that, double uhat, short lam1, short lam2, short lam3, short lam4 ) { //--- first compute some kinematic variables const double cos_theta = ( that-uhat ) / shat / sqrt( 1.+1.e-10-4.*mw2_/shat ), cos_theta2 = cos_theta*cos_theta; const double sin_theta2 = 1.-cos_theta2, sin_theta = sqrt( sin_theta2 ); const double beta = sqrt( 1.-4.*mw2_/shat ), beta2 = beta*beta; const double inv_gamma = sqrt( 1.-beta2 ), gamma = 1./inv_gamma, gamma2 = gamma*gamma, inv_gamma2 = inv_gamma*inv_gamma; const double invA = 1./( 1.-beta2*cos_theta2 ); //--- per-helicity amplitude // longitudinal-longitudinal if ( lam3 == 0 && lam4 == 0 ) return invA*inv_gamma2*( ( gamma2+1. )*( 1.-lam1*lam2 )*sin_theta2 - ( 1.+lam1*lam2 ) ); // transverse-longitudinal if ( lam4 == 0 ) return invA*( -M_SQRT2*inv_gamma*( lam1-lam2 )*( 1.+lam1*lam3*cos_theta )*sin_theta ); // longitudinal-transverse if ( lam3 == 0 ) return invA*( -M_SQRT2*inv_gamma*( lam2-lam1 )*( 1.+lam2*lam4*cos_theta )*sin_theta ); // transverse-transverse if ( lam3 != 0 && lam4 != 0 ) return -0.5*invA*( 2.*beta*( lam1+lam2 )*( lam3+lam4 ) -inv_gamma2*( 1.+lam3*lam4 )*( 2.*lam1*lam2+( 1.-lam1*lam2 ) * cos_theta2 ) +( 1.+lam1*lam2*lam3*lam4 )*( 3.+lam1*lam2 ) +2.*( lam1-lam2 )*( lam3-lam4 )*cos_theta +( 1.-lam1*lam2 )*( 1.-lam3*lam4 )*cos_theta2 ); return 0.; } // register process and define aliases REGISTER_PROCESS( pptoww, PPtoWW ) } } diff --git a/CepGen/StructureFunctions/CLAS.cpp b/CepGen/StructureFunctions/CLAS.cpp index c309680..fb448d5 100644 --- a/CepGen/StructureFunctions/CLAS.cpp +++ b/CepGen/StructureFunctions/CLAS.cpp @@ -1,220 +1,220 @@ #include "CepGen/StructureFunctions/CLAS.h" #include "CepGen/Physics/PDG.h" #include "CepGen/Physics/Constants.h" #include "CepGen/Event/Particle.h" #include "CepGen/Core/Exception.h" namespace CepGen { namespace sf { CLAS::Parameters CLAS::Parameters::standard_proton() { Parameters params; params.mode = Parameters::proton; params.mp = mp_; - params.mpi0 = ParticleProperties::mass( PDG::piZero ); + params.mpi0 = part::mass( PDG::piZero ); // SLAC fit parameters params.c_slac = { { 0.25615, 2.1785, 0.89784, -6.7162, 3.7557, 1.6421, 0.37636 } }; // CLAS parameterisation params.x = { { -0.599937, 4.76158, 0.411676 } }; params.b = { { 0.755311, 3.35065, 3.51024, 1.74470 } }; params.alpha = -0.174985; params.beta = 0.00967019; params.mu = -0.0352567; params.mup = 3.51852; Parameters::Resonance r0; r0.amplitude = 1.04; r0.mass = 1.22991; r0.width = 0.106254; r0.angular_momentum = 1; params.resonances.emplace_back( r0 ); Parameters::Resonance r1; r1.amplitude = 0.481327; r1.mass = 1.51015; r1.width = 0.0816620; r1.angular_momentum = 2; params.resonances.emplace_back( r1 ); Parameters::Resonance r2; r2.amplitude = 0.655872; r2.mass = 1.71762; r2.width = 0.125520; r2.angular_momentum = 3; params.resonances.emplace_back( r2 ); Parameters::Resonance r3; r3.amplitude = 0.747338; r3.mass = 1.95381; r3.width = 0.198915; r3.angular_momentum = 2; params.resonances.emplace_back( r3 ); return params; } CLAS::Parameters CLAS::Parameters::standard_neutron() { Parameters params = standard_proton(); params.mode = Parameters::neutron; params.c_slac = { { 0.0640, 0.2250, 4.1060, -7.0790, 3.0550, 1.6421, 0.37636 } }; return params; } CLAS::Parameters CLAS::Parameters::standard_deuteron() { Parameters params = standard_proton(); params.mode = Parameters::deuteron; params.c_slac = { { 0.47709, 2.1602, 3.6274, -10.470, 4.9272, 1.5121, 0.35115 } }; params.x = { { -0.21262, 6.9690, 0.40314 } }; params.b = { { 0.76111, 4.1470, 3.7119, 1.4218 } }; params.alpha = -0.24480; params.beta = 0.014503; params.resonances.clear(); Parameters::Resonance r0; r0.amplitude = 0.74847; r0.mass = 1.2400; r0.width = 0.12115; r0.angular_momentum = 1; params.resonances.emplace_back( r0 ); Parameters::Resonance r1; r1.amplitude = 0.011500; r1.mass = 1.4772; r1.width = 0.0069580; r1.angular_momentum = 2; params.resonances.emplace_back( r1 ); Parameters::Resonance r2; r2.amplitude = 0.12662; r2.mass = 1.5233; r2.width = 0.084095; r2.angular_momentum = 3; params.resonances.emplace_back( r2 ); Parameters::Resonance r3; r3.amplitude = 0.747338; r3.mass = 1.95381; r3.width = 0.198915; r3.angular_momentum = 2; params.resonances.emplace_back( r3 ); return params; } CLAS::CLAS( const Parameters& params ) : Parameterisation( Type::CLAS ), params_( params ) {} CLAS& CLAS::operator()( double xbj, double q2 ) { std::pair<double,double> nv = { xbj, q2 }; if ( nv == old_vals_ ) return *this; old_vals_ = nv; const double mp2 = params_.mp*params_.mp; const double w2 = mp2 + q2*( 1.-xbj )/xbj; const double w_min = params_.mp+params_.mpi0; if ( sqrt( w2 ) < w_min ) { F2 = 0.; return *this; } F2 = f2slac( xbj, q2 ); std::pair<double,double> rb = resbkg( q2, sqrt( w2 ) ); F2 *= ( rb.first+rb.second ); return *this; } double CLAS::f2slac( double xbj, double q2 ) const { if ( xbj >= 1. ) return 0.; const double xsxb = ( q2+params_.c_slac[6] )/( q2+params_.c_slac[5]*xbj ); const double xs = xbj*xsxb; double f2 = 0.; for ( unsigned short i = 0; i < 5; ++i ) f2 += params_.c_slac[i]*pow( 1.-xs, i ); if ( params_.mode == Parameters::deuteron && xbj > 0. ) f2 /= ( 1.-exp( -7.70*( 1./xbj-1.+params_.mp*params_.mp/q2 ) ) ); return f2 * pow( 1.-xs, 3 ) / xsxb; } std::pair<double,double> CLAS::resbkg( double q2, double w ) const { const double mp2 = params_.mp*params_.mp, mpi02 = params_.mpi0*params_.mpi0; const double coef = 6.08974; double wth = params_.mp+params_.mpi0; if ( w < wth ) return std::make_pair( 0., 0. ); if ( w > 4. ) return std::make_pair( 1., 0. ); const double w2 = w*w; double qs = pow( w2+mp2-mpi02, 2 )-4.*mp2*w2; if ( qs <= 0. ) return std::make_pair( 1., 0. ); qs = 0.5 * sqrt( qs )/w; const double omega = 0.5*( w2+q2-mp2 )/params_.mp; const double xn = 0.5*q2/( params_.mp*omega ); const double bkg2 = ( w > params_.b[3] ) ? exp( -params_.b[2]*( w2-params_.b[3]*params_.b[3] ) ) : 1.; double f2bkg = ( params_.b[0] )*( 1.-exp( -params_.b[1]*( w-wth ) ) ) + ( 1.-params_.b[0] )*( 1.-bkg2 ); f2bkg *= ( 1.+( 1.-f2bkg )*( params_.x[0]+params_.x[1]*pow( xn-params_.x[2], 2 ) ) ); double etab = 1., etad = 1.; if ( params_.mode != Parameters::deuteron && q2 <= 2. && w <= 2.5 ) { etab = 1.-2.5*q2*exp( -12.5*q2*q2-50.*( w-1.325 )*( w-1.325 ) ); etad = 1.+2.5*q2*exp( -12.5*q2*q2 ); } f2bkg *= etab; double f2resn = 0.; for ( unsigned short i = 0; i < params_.resonances.size(); ++i ) { const Parameters::Resonance& res = params_.resonances[i]; const double ai = ( i == 0 ) ? etad * ( res.amplitude + q2*std::min( 0., params_.alpha+params_.beta*q2 ) ) : res.amplitude; const double dmi = ( i == 2 ) ? res.mass * ( 1.+params_.mu/( 1.+params_.mup*q2 ) ) : res.mass; double qs0 = pow( dmi*dmi+mp2-mpi02, 2 )-4.*mp2*dmi*dmi; if ( qs0 <= 0. ) break; qs0 = 0.5*sqrt( qs0 )/dmi; int ji = 2*res.angular_momentum; const double dg = 0.5*res.width*pow( qs/qs0, ji+1 )*( 1.+pow( coef*qs0, ji ) )/( 1.+pow( coef*qs, ji ) ); f2resn += ai*dg/( ( w-dmi )*( w-dmi )+dg*dg ); } f2resn *= 0.5*( 1.-params_.b[0] )*bkg2/params_.mp*M_1_PI; return std::make_pair( f2bkg, f2resn ); } } } diff --git a/CepGen/StructureFunctions/ChristyBosted.cpp b/CepGen/StructureFunctions/ChristyBosted.cpp index 9d1fc4e..264b6f4 100644 --- a/CepGen/StructureFunctions/ChristyBosted.cpp +++ b/CepGen/StructureFunctions/ChristyBosted.cpp @@ -1,284 +1,284 @@ #include "CepGen/StructureFunctions/ChristyBosted.h" #include "CepGen/Physics/PDG.h" #include "CepGen/Event/Particle.h" #include "CepGen/Core/Exception.h" namespace CepGen { namespace sf { ChristyBosted::ChristyBosted( const Parameters& params ) : Parameterisation( Type::ChristyBosted ), params_( params ) {} double ChristyBosted::resmod507( char sf, double w2, double q2 ) const { - const double mpi = ParticleProperties::mass( PDG::piZero ), mpi2 = mpi*mpi, - meta = ParticleProperties::mass( PDG::eta ), meta2 = meta*meta; + const double mpi = part::mass( PDG::piZero ), mpi2 = mpi*mpi, + meta = part::mass( PDG::eta ), meta2 = meta*meta; const double w = sqrt( w2 ); const double xb = q2/( q2+w2-mp2_ ); double m0 = 0., q20 = 0.; if ( sf == 'T' ) { // transverse m0 = 0.125; q20 = 0.05; } else if ( sf == 'L' ) { m0 = params_.m0; q20 = 0.125; } else { CG_ERROR( "ChristyBosted" ) << "Invalid direction retrieved! Aborting."; return 0.; } const double norm_q2 = 1./0.330/0.330; const double t = log( log( ( q2+m0 )*norm_q2 )/log( m0*norm_q2 ) ); //--- calculate kinematics needed for threshold relativistic B-W // equivalent photon energies const double k = 0.5 * ( w2 - mp2_ )/mp_; const double kcm = 0.5 * ( w2 - mp2_ )/w; const double epicm = 0.5 * ( w2 + mpi2 - mp2_ )/w, ppicm = sqrt( std::max( 0., epicm* epicm - mpi2 ) ); const double epi2cm = 0.5 * ( w2 + 4.*mpi2 - mp2_ )/w, ppi2cm = sqrt( std::max( 0., epi2cm*epi2cm - 4*mpi2 ) ); const double eetacm = 0.5 * ( w2 + meta2 - mp2_ )/w, petacm = sqrt( std::max( 0., eetacm*eetacm - meta2 ) ); std::array<double,7> width, height, pgam; for ( unsigned short i = 0; i < 7; ++i ) { const Parameters::Resonance& res = params_.resonances[i]; width[i] = res.width; //--- calculate partial widths //----- 1-pion decay mode const double x02 = res.x0*res.x0; const double partial_width_singlepi = pow( ppicm /res.pcmr( mpi2 ), 2.*res.angular_momentum+1. ) * pow( ( res.pcmr( mpi2 )*res.pcmr( mpi2 )+x02 )/( ppicm *ppicm +x02 ), res.angular_momentum ); //----- 2-pion decay mode const double partial_width_doublepi = pow( ppi2cm/res.pcmr( 4.*mpi2 ), 2.*( res.angular_momentum+2. ) ) * pow( ( res.pcmr( 4.*mpi2 )*res.pcmr( 4.*mpi2 )+x02 )/( ppi2cm*ppi2cm+x02 ), res.angular_momentum+2 ) * w / res.mass; //----- eta decay mode (only for S11's) const double partial_width_eta = ( res.br.eta == 0. ) ? 0. : pow( petacm/res.pcmr( meta2 ), 2.*res.angular_momentum+1. ) * pow( ( res.pcmr( meta2 )*res.pcmr( meta2 )+x02 )/( petacm*petacm+x02 ), res.angular_momentum ); // virtual photon width pgam[i] = res.width * pow( kcm/res.kcmr(), 2 ) * ( res.kcmr()*res.kcmr()+x02 )/( kcm*kcm+x02 ); width[i] = ( partial_width_singlepi * res.br.singlepi + partial_width_doublepi * res.br.doublepi + partial_width_eta * res.br.eta ) * res.width; //--- resonance Q^2 dependence calculations if ( sf == 'T' ) height[i] = res.A0_T*( 1.+res.fit_parameters[0]*q2/( 1.+res.fit_parameters[1]*q2 ) )/pow( 1.+q2/0.91, res.fit_parameters[2] ); else if ( sf == 'L' ) height[i] = res.A0_L/( 1.+res.fit_parameters[3]*q2 )*q2*exp( -q2*res.fit_parameters[4] ); height[i] = height[i]*height[i]; } //--- calculate Breit-Wigners for all resonances double sig_res = 0.; for ( unsigned short i = 0; i < 7; ++i ) { const Parameters::Resonance& res = params_.resonances[i]; const double mass2 = res.mass*res.mass, width2 = width[i]*width[i]; const double sigr = height[i]*res.kr()/k*res.kcmr()/kcm/res.width * ( width[i]*pgam[i] / ( pow( w2-mass2, 2 ) + mass2*width2 ) ); sig_res += sigr; } sig_res *= w; //--- non-resonant background calculation const double xpr = 1./( 1.+( w2-pow( mp_+mpi, 2 ) )/( q2+q20 ) ); if ( xpr > 1. ) return 0.; // FIXME double sig_nr = 0.; if ( sf == 'T' ) { // transverse const double wdif = w - ( mp_ + mpi ); if ( wdif >= 0. ) { for ( unsigned short i = 0; i < 2; ++i ) { const double expo = params_.continuum.transverse[i].fit_parameters[1] + params_.continuum.transverse[i].fit_parameters[2]*q2 + params_.continuum.transverse[i].fit_parameters[3]*q2*q2; sig_nr += params_.continuum.transverse[i].sig0 / pow( q2+params_.continuum.transverse[i].fit_parameters[0], expo ) * pow( wdif, i+1.5 ); } } sig_nr *= xpr; } else if ( sf == 'L' ) { // longitudinal for ( unsigned short i = 0; i < 1; ++i ) { const double expo = params_.continuum.longitudinal[i].fit_parameters[0] + params_.continuum.longitudinal[i].fit_parameters[1]; sig_nr += params_.continuum.longitudinal[i].sig0 * pow( 1.-xpr, expo )/( 1.-xb ) * pow( q2/( q2+q20 ), params_.continuum.longitudinal[i].fit_parameters[2] )/( q2+q20 ) * pow( xpr, params_.continuum.longitudinal[i].fit_parameters[3]+params_.continuum.longitudinal[i].fit_parameters[4]*t ); } } return sig_res + sig_nr; } ChristyBosted::Parameters ChristyBosted::Parameters::standard() { Parameters params; params.m0 = 4.2802; params.continuum.transverse = { { Continuum::Direction( 246.06, { { 0.067469, 1.3501, 0.12054, -0.0038495 } } ), Continuum::Direction( -89.360, { { 0.20977, 1.5715, 0.090736, 0.010362 } } ) } }; params.continuum.longitudinal = { { Continuum::Direction( 86.746, { { 0., 4.0294, 3.1285, 0.33403, 4.9623 } } ) } }; { //--- P33(1232) Resonance p33; p33.br = Resonance::BR( 1., 0., 0. ); p33.angular_momentum = 1.; //p33.x0 = 0.15; p33.x0 = 0.14462; p33.mass = 1.2298; p33.width = 0.13573; p33.fit_parameters = { { 4.2291, 1.2598, 2.1242, 19.910, 0.22587 } }; p33.A0_T = 7.7805; p33.A0_L = 29.414; params.resonances.emplace_back( p33 ); } { //--- S11(1535) Resonance s11_1535; s11_1535.br = Resonance::BR( 0.45, 0.1, 0.45 ); s11_1535.angular_momentum = 0.; s11_1535.x0 = 0.215; s11_1535.mass = 1.5304; s11_1535.width = 0.220; s11_1535.fit_parameters = { { 6823.2, 33521., 2.5686, 0., 0. } }; s11_1535.A0_T = 6.3351; s11_1535.A0_L = 0.; params.resonances.emplace_back( s11_1535 ); } { //--- D13(1520) Resonance d13; d13.br = Resonance::BR( 0.65, 0.35, 0. ); d13.angular_momentum = 2.; d13.x0 = 0.215; d13.mass = 1.5057; d13.width = 0.082956; d13.fit_parameters = { { 21.240, 0.055746, 2.4886, 97.046, 0.31042 } }; d13.A0_T = 0.60347; d13.A0_L = 157.92; params.resonances.emplace_back( d13 ); } { //--- F15(1680) Resonance f15; f15.br = Resonance::BR( 0.65, 0.35, 0. ); f15.angular_momentum = 3.; f15.x0 = 0.215; f15.mass = 1.6980; f15.width = 0.095782; f15.fit_parameters = { { -0.28789, 0.18607, 0.063534, 0.038200, 1.2182 } }; f15.A0_T = 2.3305; f15.A0_L = 4.2160; params.resonances.emplace_back( f15 ); } { //--- S11(1650) Resonance s11_1650; s11_1650.br = Resonance::BR( 0.4, 0.5, 0.1 ); s11_1650.angular_momentum = 0.; s11_1650.x0 = 0.215; s11_1650.mass = 1.6650; s11_1650.width = 0.10936; s11_1650.fit_parameters = { { -0.56175, 0.38964, 0.54883, 0.31393, 2.9997 } }; s11_1650.A0_T = 1.9790; s11_1650.A0_L = 13.764; params.resonances.emplace_back( s11_1650 ); } { //--- P11(1440) roper Resonance p11; p11.br = Resonance::BR( 0.65, 0.35, 0. ); p11.angular_momentum = 1.; p11.x0 = 0.215; p11.mass = 1.4333; p11.width = 0.37944; p11.fit_parameters = { { 46.213, 0.19221, 1.9141, 0.053743, 1.3091 } }; p11.A0_T = 0.022506; p11.A0_L = 5.5124; params.resonances.emplace_back( p11 ); } { //--- F37(1950) Resonance f37; f37.br = Resonance::BR( 0.5, 0.5, 0. ); f37.angular_momentum = 3.; f37.x0 = 0.215; f37.mass = 1.9341; f37.width = 0.380; f37.fit_parameters = { { 0., 0., 1., 1.8951, 0.51376 } }; f37.A0_T = 3.4187; f37.A0_L = 1.8951; params.resonances.emplace_back( f37 ); } return params; } double ChristyBosted::Parameters::Resonance::kr() const { return 0.5 * ( mass*mass-mp2_ ) / mp_; } double ChristyBosted::Parameters::Resonance::ecmr( double m2 ) const { if ( mass == 0. ) return 0.; return 0.5 * ( mass*mass+m2-mp2_ ) / mass; } ChristyBosted& ChristyBosted::operator()( double xbj, double q2 ) { std::pair<double,double> nv = { xbj, q2 }; if ( nv == old_vals_ ) return *this; old_vals_ = nv; const double w2 = mp2_ + q2*( 1.-xbj )/xbj; - const double w_min = mp_+ParticleProperties::mass( PDG::piZero ); + const double w_min = mp_+part::mass( PDG::piZero ); if ( sqrt( w2 ) < w_min ) { F2 = 0.; return *this; } //----------------------------- // modification of Christy-Bosted at large q2 as described in the LUXqed paper //----------------------------- const double q21 = 30., q20 = 8.; const double delq2 = q2 - q20; const double qq = q21 - q20; - const double prefac = 1./( 4.*M_PI*M_PI*Constants::alphaEM ) * ( 1.-xbj ); + const double prefac = 1./( 4.*M_PI*M_PI*constants::alphaEM ) * ( 1.-xbj ); //------------------------------ double q2_eff = q2, w2_eff = w2; if ( q2 > q20 ) { q2_eff = q20 + delq2/( 1.+delq2/qq ); w2_eff = mp2_ + q2_eff*( 1.-xbj )/xbj; } const double tau = 4.*xbj*xbj*mp2_/q2_eff; const double sigT = resmod507( 'T', w2_eff, q2_eff ); const double sigL = resmod507( 'L', w2_eff, q2_eff ); - F2 = prefac * q2_eff / ( 1+tau ) * ( sigT+sigL ) / Constants::GeV2toBarn * 1.e6; + F2 = prefac * q2_eff / ( 1+tau ) * ( sigT+sigL ) / constants::GeV2toBarn * 1.e6; if ( q2 > q20 ) F2 *= q21/( q21 + delq2 ); if ( sigT != 0. ) Parameterisation::computeFL( q2_eff, xbj, sigL/sigT ); return *this; } } } diff --git a/CepGen/StructureFunctions/FioreBrasse.cpp b/CepGen/StructureFunctions/FioreBrasse.cpp index a5d3d9d..937a0b5 100644 --- a/CepGen/StructureFunctions/FioreBrasse.cpp +++ b/CepGen/StructureFunctions/FioreBrasse.cpp @@ -1,97 +1,97 @@ #include "CepGen/StructureFunctions/FioreBrasse.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/utils.h" #include "CepGen/Physics/PDG.h" #include "CepGen/Physics/ParticleProperties.h" #include "CepGen/Physics/Constants.h" #include <complex> namespace CepGen { namespace sf { FioreBrasse::Parameters FioreBrasse::Parameters::standard() { Parameters p; p.s0 = 1.14; p.norm = 0.021; p.resonances.emplace_back( Resonance{ -0.8377, 0.95, 0.1473, 1.0, 2.4617, 3./2. } ); // N*(1520) p.resonances.emplace_back( Resonance{ -0.37, 0.95, 0.1471, 0.5399, 2.4617, 5./2. } ); // N*(1680) p.resonances.emplace_back( Resonance{ 0.0038, 0.85, 0.1969, 4.2225, 1.5722, 3./2. } ); // Δ(1236) p.resonances.emplace_back( Resonance{ 0.5645, 0.1126, 1.3086, 19.2694, 4.5259, 1. } ); // exotic return p; } FioreBrasse::Parameters FioreBrasse::Parameters::alternative() { Parameters p; p.s0 = 1.2871; p.norm = 0.0207; p.resonances.emplace_back( Resonance{ -0.8070, 0.9632, 0.1387, 1.0, 2.6066, 3./2. } ); // N*(1520) p.resonances.emplace_back( Resonance{ -0.3640, 0.9531, 0.1239, 0.6086, 2.6066, 5./2. } ); // N*(1680) p.resonances.emplace_back( Resonance{ -0.0065, 0.8355, 0.2320, 4.7279, 1.4828, 3./2. } ); // Δ(1236) p.resonances.emplace_back( Resonance{ 0.5484, 0.1373, 1.3139, 14.7267, 4.6041, 1. } ); // exotic return p; } FioreBrasse::FioreBrasse( const Parameters& params ) : Parameterisation( Type::FioreBrasse ), params( params ) {} FioreBrasse& FioreBrasse::operator()( double xbj, double q2 ) { std::pair<double,double> nv = { xbj, q2 }; if ( nv == old_vals_ ) return *this; old_vals_ = nv; const double akin = 1. + 4.*mp2_ * xbj*xbj/q2; - const double prefactor = q2*( 1.-xbj ) / ( 4.*M_PI*Constants::alphaEM*akin ); + const double prefactor = q2*( 1.-xbj ) / ( 4.*M_PI*constants::alphaEM*akin ); const double s = q2*( 1.-xbj )/xbj + mp2_; double ampli_res = 0., ampli_bg = 0., ampli_tot = 0.; for ( unsigned short i = 0; i < 3; ++i ) { //FIXME 4?? const Parameters::Resonance& res = params.resonances[i]; const double sqrts0 = sqrt( params.s0 ); std::complex<double> alpha; if ( s > params.s0 ) alpha = std::complex<double>( res.alpha0 + res.alpha2*sqrts0 + res.alpha1*s, res.alpha2*sqrt( s-params.s0 ) ); else alpha = std::complex<double>( res.alpha0 + res.alpha1*s + res.alpha2*( sqrts0 - sqrt( params.s0 - s ) ), 0. ); double formfactor = 1./pow( 1. + q2/res.q02, 2 ); double denom = pow( res.spin-std::real( alpha ), 2 ) + pow( std::imag( alpha ), 2 ); double ampli_imag = res.a*formfactor*formfactor*std::imag( alpha )/denom; ampli_res += ampli_imag; } { const Parameters::Resonance& res = params.resonances[3]; double sE = res.alpha2, sqrtsE = sqrt( sE ); std::complex<double> alpha; if ( s > sE ) alpha = std::complex<double>( res.alpha0 + res.alpha1*sqrtsE, res.alpha1*sqrt( s-sE ) ); else alpha = std::complex<double>( res.alpha0 + res.alpha1*( sqrtsE - sqrt( sE-s ) ), 0. ); double formfactor = 1./pow( 1. + q2/res.q02, 2 ); double sp = 1.5*res.spin; double denom = pow( sp-std::real( alpha ), 2 ) + pow( std::imag( alpha ), 2 ); ampli_bg = res.a*formfactor*formfactor*std::imag( alpha )/denom; } ampli_tot = params.norm*( ampli_res+ampli_bg ); CG_DEBUG_LOOP( "FioreBrasse:amplitudes" ) << "Amplitudes:\n\t" << " resonance part: " << ampli_res << ",\n\t" << " background part: " << ampli_bg << ",\n\t" << " total (with norm.): " << ampli_tot << "."; F2 = prefactor*ampli_tot; return *this; } } } diff --git a/CepGen/StructureFunctions/SigmaRatio.cpp b/CepGen/StructureFunctions/SigmaRatio.cpp index 0aa297b..ec4aa98 100644 --- a/CepGen/StructureFunctions/SigmaRatio.cpp +++ b/CepGen/StructureFunctions/SigmaRatio.cpp @@ -1,104 +1,104 @@ #include "CepGen/StructureFunctions/SigmaRatio.h" #include "CepGen/Physics/PDG.h" #include "CepGen/Physics/ParticleProperties.h" #include <math.h> #include <iostream> namespace CepGen { namespace sr { - const double Parameterisation::mp_ = ParticleProperties::mass( PDG::proton ); + const double Parameterisation::mp_ = part::mass( PDG::proton ); const double Parameterisation::mp2_ = Parameterisation::mp_*Parameterisation::mp_; double Parameterisation::theta( double xbj, double q2 ) const { return 1.+12.*( q2/( q2+1. ) )*( 0.125*0.125/( 0.125*0.125+xbj*xbj ) ); } E143::Parameters E143::Parameters::standard() { Parameters out; out.q2_b = 0.34; out.lambda2 = 0.2*0.2; out.a = { { 0.0485, 0.5470, 2.0621, -0.3804, 0.5090, -0.0285 } }; out.b = { { 0.0481, 0.6114, -0.3509, -0.4611, 0.7172, -0.0317 } }; out.c = { { 0.0577, 0.4644, 1.8288, 12.3708, -43.1043, 41.7415 } }; return out; } E143::E143( const Parameters& param ) : params_( param ) {} double E143::operator()( double xbj, double q2, double& err ) const { const double u = q2/params_.q2_b; const double inv_xl = 1./log( q2/params_.lambda2 ); const double pa = ( 1.+params_.a[3]*xbj+params_.a[4]*xbj*xbj )*pow( xbj, params_.a[5] ); const double pb = ( 1.+params_.b[3]*xbj+params_.b[4]*xbj*xbj )*pow( xbj, params_.b[5] ); const double th = theta( xbj, q2 ); const double q2_thr = params_.c[3]*xbj + params_.c[4]*xbj*xbj+params_.c[5]*xbj*xbj*xbj; // here come the three fits const double ra = params_.a[0]*inv_xl*th + params_.a[1]/pow( pow( q2, 4 )+pow( params_.a[2], 4 ), 0.25 )*pa, rb = params_.b[0]*inv_xl*th + ( params_.b[1]/q2+params_.b[2]/( q2*q2+0.3*0.3 ) )*pb, rc = params_.c[0]*inv_xl*th + params_.c[1]*pow( pow( q2-q2_thr, 2 )+pow( params_.c[2], 2 ), -0.5 ); const double r = ( ra+rb+rc ) / 3.; // R is set to be the average of the three fits // numerical safety for low-Q² err = 0.0078-0.013*xbj+( 0.070-0.39*xbj+0.70*xbj*xbj )/( 1.7+q2 ); if ( q2 > params_.q2_b ) return r; return r * 0.5 * ( 3.*u-u*u*u ); } R1990::Parameters R1990::Parameters::standard() { Parameters out; out.lambda2 = 0.2*0.2; out.b = { { 0.0635, 0.5747, -0.3534 } }; return out; } R1990::R1990( const Parameters& param ) : params_( param ) {} double R1990::operator()( double xbj, double q2, double& err ) const { err = 0.; return ( params_.b[0]+theta( xbj, q2 )/log( q2/params_.lambda2 ) + params_.b[1]/q2 + params_.b[2]/( q2*q2+0.09 ) ); } double CLAS::operator()( double xbj, double q2, double& err ) const { // 2 kinematic regions: // - resonances ( w < 2.5 ) // - DIS ( w > 2.5 ) const double w2 = mp2_ + q2*( 1.-xbj )/xbj, w = sqrt( w2 ); const double xth = q2/( q2+2.5*2.5-mp2_ ); // xth = x( W = 2.5 GeV ) const double zeta = log( 25.*q2 ); const double xitmp = ( w < 2.5 ) ? theta( xth, q2 ) : theta( xbj, q2 ); const double tmp = 0.041*xitmp/zeta + 0.592/q2 - 0.331/( 0.09+q2*q2 ); if ( w < 2.5 ) return tmp * pow( ( 1.-xbj )/( 1.-xth ), 3 ); return tmp; } double SibirtsevBlunden::operator()( double xbj, double q2, double& err ) const { err = 0.; return 0.014*q2*( exp( -0.07*q2 )+41.*exp( -0.8*q2 ) ); } } } diff --git a/CepGen/StructureFunctions/StructureFunctions.cpp b/CepGen/StructureFunctions/StructureFunctions.cpp index 272bbfc..5152df0 100644 --- a/CepGen/StructureFunctions/StructureFunctions.cpp +++ b/CepGen/StructureFunctions/StructureFunctions.cpp @@ -1,90 +1,90 @@ #include "CepGen/StructureFunctions/StructureFunctions.h" #include "CepGen/Physics/PDG.h" #include "CepGen/Physics/ParticleProperties.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/utils.h" #include <iostream> namespace CepGen { namespace sf { - const double Parameterisation::mp_ = ParticleProperties::mass( PDG::proton ); + const double Parameterisation::mp_ = part::mass( PDG::proton ); const double Parameterisation::mp2_ = Parameterisation::mp_*Parameterisation::mp_; double Parameterisation::F1( double xbj, double q2 ) const { if ( xbj == 0. || q2 == 0. ) { CG_ERROR( "StructureFunctions:F1" ) << "Invalid range for Q² = " << q2 << " or xBj = " << xbj << "."; return 0.; } const double F1 = 0.5*( ( 1+4.*xbj*xbj*mp2_/q2 )*F2 - FL )/xbj; CG_DEBUG_LOOP( "StructureFunctions:F1" ) << "F1 for Q² = " << q2 << ", xBj = " << xbj << ": " << F1 << "\n\t" << "(F2 = " << F2 << ", FL = " << FL << ")."; return F1; } void Parameterisation::computeFL( double xbj, double q2, const sr::Parameterisation& ratio ) { double r_error = 0.; computeFL( xbj, q2, ratio( xbj, q2, r_error ) ); } void Parameterisation::computeFL( double xbj, double q2, double r ) { const double tau = 4.*xbj*xbj*mp2_/q2; FL = F2 * ( 1.+tau ) * ( r/( 1.+r ) ); } std::string Parameterisation::description() const { std::ostringstream os; os << type; return os.str(); } std::ostream& operator<<( std::ostream& os, const Parameterisation& sf ) { os << sf.description(); if ( sf.old_vals_ != std::pair<double,double>() ) os << " at (" << sf.old_vals_.first << ", " << sf.old_vals_.second << "): " << "F2 = " << sf.F2 << ", FL = " << sf.FL; return os; } } /// Human-readable format of a structure function type std::ostream& operator<<( std::ostream& os, const sf::Type& sf ) { switch ( sf ) { case sf::Type::Invalid: return os << "[INVALID]"; case sf::Type::Electron: return os << "electron"; case sf::Type::ElasticProton: return os << "elastic proton"; case sf::Type::SuriYennie: return os << "Suri-Yennie"; case sf::Type::SzczurekUleshchenko: return os << "Szczurek-Uleshchenko"; case sf::Type::FioreBrasse: return os << "Fiore-Brasse"; case sf::Type::ChristyBosted: return os << "Christy-Bosted"; case sf::Type::CLAS: return os << "CLAS"; case sf::Type::BlockDurandHa: return os << "BDH"; case sf::Type::ALLM91: return os << "ALLM91"; case sf::Type::ALLM97: return os << "ALLM97"; case sf::Type::GD07p: return os << "GD07p"; case sf::Type::GD11p: return os << "GD11p"; case sf::Type::Schaefer: return os << "LUXlike"; case sf::Type::MSTWgrid: return os << "MSTW (grid)"; case sf::Type::LHAPDF: return os << "LHAPDF"; } return os; } } diff --git a/test/TestProcess.h b/test/TestProcess.h index f846601..a4be21a 100644 --- a/test/TestProcess.h +++ b/test/TestProcess.h @@ -1,44 +1,43 @@ #ifndef CepGen_Processes_TestProcess_h #define CepGen_Processes_TestProcess_h #include "CepGen/Processes/GenericProcess.h" #include "CepGen/Core/Functional.h" namespace CepGen { - namespace Process + namespace process { /// Generic process to test the Vegas instance template<size_t N> class TestProcess : public GenericProcess { public: TestProcess() : GenericProcess( "test", ".oO TEST PROCESS Oo.", false ), funct_( "1./(1.-cos(x*_pi)*cos(y*_pi)*cos(z*_pi))", { { "x", "y", "z" } } ) {} TestProcess( const char* formula, std::array<std::string,N> args ) : GenericProcess( "test", Form( ".oO TEST PROCESS (%s) Oo.", formula ), false ), funct_( formula, args ) {} ProcessPtr clone() const override { return ProcessPtr( new TestProcess<N>( *this ) ); } void addEventContent() override {} /// Number of dimensions on which to perform the integration unsigned int numDimensions() const override { return N; } /// Generic formula to compute a weight out of a point in the phase space double computeWeight() override { std::array<double,N> args; std::copy_n( x_.begin(), N, args.begin() ); return funct_.eval( args ); } /// Dummy function to be called on events generation void fillKinematics( bool ) override { return; } private: Functional<N> funct_; }; } } #endif - diff --git a/test/cepgen-ascii.cpp b/test/cepgen-ascii.cpp index 64e4ef1..11d48fa 100644 --- a/test/cepgen-ascii.cpp +++ b/test/cepgen-ascii.cpp @@ -1,86 +1,85 @@ #include "CepGen/Cards/PythonHandler.h" #include "CepGen/Cards/LpairHandler.h" #include "CepGen/IO/HepMCHandler.h" #include "CepGen/IO/LHEFHandler.h" #include "CepGen/Generator.h" #include "CepGen/Core/Exception.h" #include "CepGen/Event/Event.h" #include <iostream> using namespace std; // we use polymorphism here -std::shared_ptr<CepGen::OutputHandler::ExportHandler> writer; +std::shared_ptr<CepGen::output::ExportHandler> writer; void storeEvent( const CepGen::Event& ev, unsigned long ) { if ( !writer ) throw CG_FATAL( "storeEvent" ) << "Failed to retrieve a valid writer!"; *writer << ev; } /** * Main caller for this Monte Carlo generator. Loads the configuration files' * variables if set as an argument to this program, else loads a default * "LHC-like" configuration, then launches the cross-section computation and * the events generation. * \author Laurent Forthomme <laurent.forthomme@cern.ch> */ int main( int argc, char* argv[] ) { if ( argc < 2 ) throw CG_FATAL( "main" ) << "No config file provided!\n\t" << "Usage: " << argv[0] << " config-file [format=lhef,hepmc] [filename=example.dat]"; CepGen::Generator mg; //----------------------------------------------------------------------------------------------- // Steering card readout //----------------------------------------------------------------------------------------------- CG_DEBUG( "main" ) << "Reading config file stored in \"" << argv[1] << "\""; - const string extension = CepGen::Cards::Handler::getExtension( argv[1] ); + const string extension = CepGen::cards::Handler::getExtension( argv[1] ); if ( extension == "card" ) - mg.setParameters( CepGen::Cards::LpairHandler( argv[1] ).parameters() ); + mg.setParameters( CepGen::cards::LpairHandler( argv[1] ).parameters() ); else if ( extension == "py" ) - mg.setParameters( CepGen::Cards::PythonHandler( argv[1] ).parameters() ); + mg.setParameters( CepGen::cards::PythonHandler( argv[1] ).parameters() ); else throw CG_FATAL( "main" ) << "Unrecognized card format: ." << extension; //----------------------------------------------------------------------------------------------- // Output file writer definition //----------------------------------------------------------------------------------------------- const string format = ( argc > 2 ) ? argv[2] : "lhef"; const char* filename = ( argc > 3 ) ? argv[3] : "example.dat"; if ( format == "lhef" ) - writer = std::make_shared<CepGen::OutputHandler::LHEFHandler>( filename ); + writer = std::make_shared<CepGen::output::LHEFHandler>( filename ); else if ( format == "hepmc" ) - writer = std::make_shared<CepGen::OutputHandler::HepMCHandler>( filename ); + writer = std::make_shared<CepGen::output::HepMCHandler>( filename ); else throw CG_FATAL( "main" ) << "Unrecognized output format: " << format; //----------------------------------------------------------------------------------------------- // CepGen run part //----------------------------------------------------------------------------------------------- // We might want to cross-check visually the validity of our run CG_INFO( "main" ) << mg.parameters.get(); // Let there be cross-section... double xsec = 0., err = 0.; mg.computeXsection( xsec, err ); writer->initialise( *mg.parameters ); writer->setCrossSection( xsec, err ); // The events generation starts here! mg.generate( storeEvent ); return 0; } - diff --git a/test/cepgen-root.cxx b/test/cepgen-root.cxx index aac9a87..3fe755d 100644 --- a/test/cepgen-root.cxx +++ b/test/cepgen-root.cxx @@ -1,114 +1,114 @@ #include "CepGen/Cards/PythonHandler.h" #include "CepGen/Cards/LpairHandler.h" #include "CepGen/Generator.h" #include "CepGen/Event/Event.h" #include <iomanip> #include <iostream> #include "TreeInfo.h" #include "abort.h" // ROOT includes #include "TFile.h" #include "TTree.h" #include "TLorentzVector.h" using namespace std; std::unique_ptr<ROOT::CepGenRun> run; std::unique_ptr<ROOT::CepGenEvent> 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; ev->momentum.reserve( event.particles().size() ); for ( const auto& p : event.particles() ) { const CepGen::Particle::Momentum m = p.momentum(); ev->momentum[ev->np].SetPxPyPzE( m.px(), m.py(), m.pz(), p.energy() ); 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 <laurent.forthomme@cern.ch> * @date 27 jan 2014 */ int main( int argc, char* argv[] ) { CepGen::Generator mg; if ( argc < 2 ) throw CG_FATAL( "main" ) << "Usage: " << argv[0] << " input-card [filename=events.root]"; - const std::string extension = CepGen::Cards::Handler::getExtension( argv[1] ); + const std::string extension = CepGen::cards::Handler::getExtension( argv[1] ); if ( extension == "card" ) - mg.setParameters( CepGen::Cards::LpairHandler( argv[1] ).parameters() ); + mg.setParameters( CepGen::cards::LpairHandler( argv[1] ).parameters() ); else if ( extension == "py" ) - mg.setParameters( CepGen::Cards::PythonHandler( argv[1] ).parameters() ); + mg.setParameters( CepGen::cards::PythonHandler( argv[1] ).parameters() ); mg.parameters->generation.enabled = true; CG_INFO( "main" ) << mg.parameters.get(); //----- open the output root file const char* filename = ( argc > 2 ) ? argv[2] : "events.root"; std::unique_ptr<TFile> file( TFile::Open( filename, "recreate" ) ); if ( !file ) throw CG_FATAL( "main" ) << "Failed to create the output file!"; AbortHandler ctrl_c; //----- start by computing the cross section for the list of parameters applied double xsec, err; mg.computeXsection( xsec, err ); //----- then generate the events and the container tree structure std::unique_ptr<TTree> ev_tree( new TTree( "events", "A TTree containing information from the events produced from CepGen" ) ); run.reset( new ROOT::CepGenRun ); run->create(); run->xsect = xsec; run->errxsect = err; run->litigious_events = 0; run->sqrt_s = mg.parameters->kinematics.sqrtS(); run->fill(); ev.reset( new ROOT::CepGenEvent ); ev->create( ev_tree.get() ); // launch the events generation try { mg.generate( fill_event_tree ); } catch ( const CepGen::Exception& ) {} file->Write(); CG_INFO( "main" ) << "Events written on \"" << filename << "\"."; return 0; } diff --git a/test/cepgen.cpp b/test/cepgen.cpp index 9cd8999..a46fee4 100644 --- a/test/cepgen.cpp +++ b/test/cepgen.cpp @@ -1,74 +1,74 @@ //--- steering cards #include "CepGen/Cards/PythonHandler.h" #include "CepGen/Cards/LpairHandler.h" #include "CepGen/Generator.h" #include "CepGen/Core/Exception.h" //--- necessary include to build the default run #include "CepGen/Processes/GamGamLL.h" #include "CepGen/Physics/PDG.h" #include "abort.h" #include <iostream> using namespace std; /** * Main caller for this MC generator. * * loads the configuration files' variables if passed as an argument, * or a default LPAIR-like configuration, * * launches the cross-section computation and the events generation. * \author Laurent Forthomme <laurent.forthomme@cern.ch> */ int main( int argc, char* argv[] ) { //--- first start by defining the generator object CepGen::Generator gen; if ( argc < 2 ) { CG_INFO( "main" ) << "No config file provided. Setting the default parameters."; //--- default run: LPAIR elastic ɣɣ → µ⁺µ¯ at 13 TeV CepGen::ParametersList pgen; pgen.set<int>( "pair", (int)CepGen::PDG::muon ); gen.parameters->setProcess( new CepGen::process::GamGamLL( pgen ) ); gen.parameters->kinematics.mode = CepGen::KinematicsMode::ElasticElastic; gen.parameters->kinematics.cuts.central.pt_single.min() = 15.; gen.parameters->kinematics.cuts.central.eta_single = { -2.5, 2.5 }; gen.parameters->generation.enabled = true; gen.parameters->generation.maxgen = 1e3; } else { CG_INFO( "main" ) << "Reading config file stored in " << argv[1] << "."; - const std::string extension = CepGen::Cards::Handler::getExtension( argv[1] ); + const std::string extension = CepGen::cards::Handler::getExtension( argv[1] ); if ( extension == "card" ) - gen.setParameters( CepGen::Cards::LpairHandler( argv[1] ).parameters() ); + gen.setParameters( CepGen::cards::LpairHandler( argv[1] ).parameters() ); #ifdef PYTHON else if ( extension == "py" ) - gen.setParameters( CepGen::Cards::PythonHandler( argv[1] ).parameters() ); + gen.setParameters( CepGen::cards::PythonHandler( argv[1] ).parameters() ); #endif else throw CG_FATAL( "main" ) << "Unrecognized steering card extension: ." << extension << "!"; } //--- list all parameters CG_INFO( "main" ) << gen.parameters.get(); AbortHandler ctrl_c; try { //--- let there be a cross-section... double xsec = 0., err = 0.; gen.computeXsection( xsec, err ); if ( gen.parameters->generation.enabled ) //--- events generation starts here // (one may use a callback function) gen.generate(); } catch ( const CepGen::RunAbortedException& e ) { CG_INFO( "main" ) << "Run aborted!"; } return 0; } diff --git a/test/test_distributions.cxx b/test/test_distributions.cxx index 3d95514..589666c 100644 --- a/test/test_distributions.cxx +++ b/test/test_distributions.cxx @@ -1,63 +1,63 @@ #include "CepGen/Cards/PythonHandler.h" #include "CepGen/Generator.h" #include "CepGen/Event/Event.h" #include "CepGen/Core/Exception.h" #include "Canvas.h" #include "TH1.h" #include <sstream> using namespace std; void produce_plot( const char* name, TH1* hist ) { CepGen::Canvas c( name, "CepGen Simulation" ); hist->Draw( "hist" ); c.Prettify( hist ); c.SetLogy(); c.Save( "pdf" ); } unique_ptr<TH1D> h_mass, h_ptpair, h_ptsingle, h_etasingle; void process_event( const CepGen::Event& ev, unsigned long event_id ) { const auto central_system = ev.getByRole( CepGen::Particle::CentralSystem ); const auto pl1 = central_system[0].momentum(), pl2 = central_system[1].momentum(); h_mass->Fill( ( pl1+pl2 ).mass() ); h_ptpair->Fill( ( pl1+pl2 ).pt() ); h_ptsingle->Fill( pl1.pt() ); h_etasingle->Fill( pl1.eta() ); } int main( int argc, char* argv[] ) { CepGen::Generator mg; if ( argc < 2 ) throw CG_FATAL( "main" ) << "Usage: " << argv[0] << " [input card]"; - mg.setParameters( CepGen::Cards::PythonHandler( argv[1] ).parameters() ); + mg.setParameters( CepGen::cards::PythonHandler( argv[1] ).parameters() ); h_mass = unique_ptr<TH1D>( new TH1D( "invm", ";Dilepton invariant mass;d#sigma/dM (pb/GeV)", 500, 0., 500. ) ); h_ptpair = unique_ptr<TH1D>( new TH1D( "ptpair", ";Dilepton p_{T};d#sigma/dp_{T} (pb/GeV)", 500, 0., 50. ) ); h_ptsingle = unique_ptr<TH1D>( new TH1D( "pt_single", ";Single lepton p_{T};d#sigma/dp_{T} (pb/GeV)", 100, 0., 100. ) ); h_etasingle = unique_ptr<TH1D>( new TH1D( "eta_single", ";Single lepton #eta;d#sigma/d#eta (pb)\\?.2f", 60, -3., 3. ) ); CG_INFO( "main" ) << "Process name: " << mg.parameters->processName() << "."; //mg.parameters->taming_functions.dump(); mg.generate( process_event ); const double weight = mg.crossSection()/mg.parameters->generation.maxgen; h_mass->Scale( weight, "width" ); h_ptpair->Scale( weight, "width" ); h_ptsingle->Scale( weight, "width" ); h_etasingle->Scale( weight, "width" ); produce_plot( "dilepton_invm", h_mass.get() ); produce_plot( "dilepton_ptpair", h_ptpair.get() ); produce_plot( "singlelepton_pt", h_ptsingle.get() ); produce_plot( "singlelepton_eta", h_etasingle.get() ); return 0; } diff --git a/test/test_event_writer.cpp b/test/test_event_writer.cpp index 90639e7..6846bbf 100644 --- a/test/test_event_writer.cpp +++ b/test/test_event_writer.cpp @@ -1,30 +1,30 @@ #include "CepGen/IO/HepMCHandler.h" #include "CepGen/Physics/PDG.h" #include "CepGen/Event/Event.h" using namespace std; using namespace CepGen; int main() { - OutputHandler::HepMCHandler writer( "example.dat" ); + output::HepMCHandler writer( "example.dat" ); writer.setCrossSection(1., 2.); Event ev; Particle p1( Particle::IncomingBeam1, PDG::proton ); p1.setMomentum( 1., -15., 100. ); p1.setStatus( Particle::Status::Incoming ); ev.addParticle(p1); Particle p2( Particle::IncomingBeam2, PDG::electron ); p2.setMomentum( 10., 5., 3200. ); p2.setStatus( Particle::Status::Incoming ); ev.addParticle(p2); ev.dump(); writer << ev; return 0; } diff --git a/test/test_pythoncfg_parser.cpp b/test/test_pythoncfg_parser.cpp index 8817255..701b4bd 100644 --- a/test/test_pythoncfg_parser.cpp +++ b/test/test_pythoncfg_parser.cpp @@ -1,22 +1,22 @@ #include "CepGen/Cards/PythonHandler.h" #include "CepGen/Core/Exception.h" #include <string> #include <iostream> using namespace std; int main( int argc, char* argv[] ) { if ( argc < 2 ) throw CG_FATAL( "main" ) << "One argument required!"; try { - CepGen::Cards::PythonHandler py( argv[1] ); + CepGen::cards::PythonHandler py( argv[1] ); CG_INFO( "main" ) << py.parameters(); } catch ( CepGen::Exception& e ) { e.dump(); } return 0; }