diff --git a/CMakeLists.txt b/CMakeLists.txt index 72a1cca80f47dce56bb7d86f5434a7bf602e2add..910ecb60adb71547f6a3a5fc70681b3909fc1dad 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -34,6 +34,9 @@ set(CMAKE_INCLUDE_CURRENT_DIR ON) # FIND LIBRARIES =========================================================================== +# find libraries: gsl +find_package(GSL) + # find libraries: Qt4 find_package(Qt4 REQUIRED QtCore QtSql QtXmlPatterns) @@ -57,7 +60,7 @@ find_package(ElementaryUtils REQUIRED) find_package(NumA++ REQUIRED) # directories containing headers -include_directories(${QT_INCLUDE_DIRS} ${SFML_INCLUDE_DIR} ${CLN_INCLUDE_DIR} ${ELEMENTARY_UTILS_INCLUDE_DIR} ${NUMA_INCLUDE_DIR}) +include_directories($(GSL_INCLUDE_DIRS) ${QT_INCLUDE_DIRS} ${SFML_INCLUDE_DIR} ${CLN_INCLUDE_DIR} ${ELEMENTARY_UTILS_INCLUDE_DIR} ${NUMA_INCLUDE_DIR}) include_directories(include) # FINALIZE ================================================================================== @@ -89,6 +92,8 @@ add_library( target_link_libraries( PARTONS + + ${GSL_LIBRARIES} ${QT_LIBRARIES} diff --git a/include/partons/ModuleObjectFactory.h b/include/partons/ModuleObjectFactory.h index 1a652ea7406837a1806f3f7f657935917b2cfce7..293fde01648111855c7c2642675e939dbb3866d5 100644 --- a/include/partons/ModuleObjectFactory.h +++ b/include/partons/ModuleObjectFactory.h @@ -33,6 +33,7 @@ class RadonInverseModule; class RunningAlphaStrongModule; class ScalesModule; class XiConverterModule; +class MWFModule; /** * @class ModuleObjectFactory @@ -236,6 +237,19 @@ public: */ ScalesModule* newScalesModule(const std::string &className); + /** + * Specialization of ModuleObjectFactory::newModuleObject into a MWFModule. + * @param classId Unique identifier of last child class. + * @return ScalesModule pointer. + */ + MWFModule* newMWFModule(unsigned int classId); + /** + * Specialization of ModuleObjectFactory::newModuleObject into a MWFModule. + * @param className Name of last child class. + * @return ScaleModule pointer. + */ + MWFModule* newMWFModule(const std::string &className); + /** * Specialization of ModuleObjectFactory::newModuleObject into a XiConverterModule. * @param classId Unique identifier of last child class. diff --git a/include/partons/beans/MesonType.h b/include/partons/beans/MesonType.h new file mode 100644 index 0000000000000000000000000000000000000000..1aa63b1c273906fdcc9bcc65fd6a3213811439f3 --- /dev/null +++ b/include/partons/beans/MesonType.h @@ -0,0 +1,99 @@ +#ifndef MESON_TYPE_H +#define MESON_TYPE_H + +/** + * + * @file MesonType.h + * @author Pawel Sznajder (NCBJ) + * @date 23 April 2019 + * @version 1.0 + */ + +#include + +namespace PARTONS { + +/** + * @class MesonType + * + * @brief Definition of enumeration values for meson types. + * + * This class defines a set of enumeration values that are used to distinguish between meson types. In addition, a declared object of this class is always associated to one meson type (see MesonType::m_type), so member functions can act on it. + */ +class MesonType { + +public: + + /** + * Definition of enumerate values corresponding to quark flavors. + */ + enum Type { + UNDEFINED = 0, //!< Undefined type. + RHO0 = 1, //!< \f$\rho^{0}\f$ + PHI = 2 //!< \f$\phi\f$ + }; + + /** + * Default constructor. + */ + MesonType(); + + /** + * Copy constructor. + * @param other Object to be copied. + */ + MesonType(const MesonType &other); + + /** + * Assignment constructor. + * @param type Type to be assigned. + */ + MesonType(Type type); + + /** + * Destructor. + */ + virtual ~MesonType(); + + /** + * Automatic cast to enum. + */ + operator Type() const; + + /** + * Get string representation of type being assigned to a declared object of this class. + * @return String representation of assigned type, like "UP" for MesonType::UP. + */ + std::string toString() const; + + /** + * Get short name representation of type being assigned to a declared object of this class. + * @return Short string representation of assigned type, like "u" for MesonType::UP. + */ + std::string getShortName(); + + //******************************************************** + //*** SETTERS AND GETTERS ******************************** + //******************************************************** + + /** + * Get type being assigned to a declared object of this class. + */ + MesonType::Type getType() const; + + /** + * Assign type to a declared object of this class. + */ + void setType(Type type); + +private: + + /** + * Type associated to a declared object of this class. + */ + MesonType::Type m_type; +}; + +} /* namespace PARTONS */ + +#endif /* MESON_TYPE_H */ diff --git a/include/partons/modules/meson_wave_function/MWFGK19.h b/include/partons/modules/meson_wave_function/MWFGK19.h new file mode 100644 index 0000000000000000000000000000000000000000..9761ffbd8b73b1354b0b34c571670a22842f5827 --- /dev/null +++ b/include/partons/modules/meson_wave_function/MWFGK19.h @@ -0,0 +1,75 @@ +#ifndef MWF_GK19_H +#define MWF_GK19_H + +/** + * @file MWFGK19.h + * @author Pawel Sznajder (NCBJ) + * @date April 22, 2019 + * @version 1.0 + */ + +#include + +#include "../../beans/MesonType.h" +#include "MWFModule.h" + +namespace PARTONS { +class RunningAlphaStrongModule; +} /* namespace PARTONS */ + +namespace PARTONS { + +/** + * @class MWFGK19 + * @brief Abstract class that provides a skeleton to implement a Meson Wave Function (MWF) module. + */ +class MWFGK19: public MWFModule { + +public: + + static const unsigned int classId; ///< ID assigned by BaseObjectRegistry + + /** + * Constructor. + * See BaseObject::BaseObject and ModuleObject::ModuleObject for more details. + * + * @param className name of child class. + */ + MWFGK19(const std::string &className); + + /** + * Default destructor. + */ + virtual ~MWFGK19(); + + virtual MWFGK19* clone() const; + + virtual void resolveObjectDependencies(); + + virtual double compute(double tau, const NumA::Vector2D& b, double muF2, + MesonType::Type mesonType, double mesonHelicity); + +protected: + + /** + * Copy constructor. + * @param other Object to be copied. + */ + MWFGK19(const MWFGK19 &other); + +private: + + /** + * Initial factorization scale. + */ + double m_MuF2_ref; + + /** + * Pointer to the running coupling module to be used. + */ + RunningAlphaStrongModule* m_pRunningAlphaStrongModule; +}; + +} /* namespace PARTONS */ + +#endif /* MWF_GK19_H */ diff --git a/include/partons/modules/meson_wave_function/MWFModule.h b/include/partons/modules/meson_wave_function/MWFModule.h new file mode 100644 index 0000000000000000000000000000000000000000..79a2a17bd8f3bb9ec0c75b071f951caf494116f9 --- /dev/null +++ b/include/partons/modules/meson_wave_function/MWFModule.h @@ -0,0 +1,70 @@ +#ifndef MWF_MODULE_H +#define MWF_MODULE_H + +/** + * @file MWFModule.h + * @author Pawel Sznajder (NCBJ) + * @date April 22, 2019 + * @version 1.0 + */ + +#include + +#include "../../beans/MesonType.h" +#include "../../ModuleObject.h" + +namespace NumA { +class Vector2D; +} /* namespace NumA */ + +namespace PARTONS { + +/** + * @class MWFModule + * @brief Abstract class that provides a skeleton to implement a Meson Wave Function (MWF) module. + */ +class MWFModule: public ModuleObject { + +public: + + /** + * Constructor. + * See BaseObject::BaseObject and ModuleObject::ModuleObject for more details. + * + * @param className name of child class. + */ + MWFModule(const std::string &className); + + /** + * Default destructor. + */ + virtual ~MWFModule(); + + virtual MWFModule* clone() const = 0; + + /** + * Evaluate meson wave function. + * @param tau Fractional longitudinal momentum carried by X quark. + * @param b Impact parameter (in XX) + * @param muF2 Factorization scale squared (in \f$GeV^{2}\f$) + * @param mesonType Meson type. + * @param mesonHelicity Meson helicity. + */ + virtual double compute(double tau, const NumA::Vector2D& b, double muF2, + MesonType::Type mesonType, double mesonHelicity) = 0; + +protected: + + /** + * Copy constructor. + * @param other Object to be copied. + */ + MWFModule(const MWFModule &other); + + virtual void initModule(); + virtual void isModuleWellConfigured(); +}; + +} /* namespace PARTONS */ + +#endif /* MWF_MODULE_H */ diff --git a/src/partons/ModuleObjectFactory.cpp b/src/partons/ModuleObjectFactory.cpp index 62481ce034b721b549ed6c6699fbe56d66786a61..f13531cc23597a93106dcdb77275c6d77ed2f669 100644 --- a/src/partons/ModuleObjectFactory.cpp +++ b/src/partons/ModuleObjectFactory.cpp @@ -11,6 +11,7 @@ #include "../../include/partons/modules/evolution/gpd/GPDEvolutionModule.h" #include "../../include/partons/modules/gpd_border_function/GPDBorderFunctionModule.h" #include "../../include/partons/modules/gpd_subtraction_constant/GPDSubtractionConstantModule.h" +#include "../../include/partons/modules/meson_wave_function/MWFModule.h" #include "../../include/partons/modules/observable/Observable.h" #include "../../include/partons/modules/overlap/IncompleteGPDModule.h" #include "../../include/partons/modules/process/DVCS/DVCSProcessModule.h" @@ -21,7 +22,6 @@ namespace PARTONS { - ModuleObjectFactory::ModuleObjectFactory(BaseObjectFactory* pBaseObjectFactory) : BaseObject("ModuleObjectFactory"), m_pBaseObjectFactory( pBaseObjectFactory) { @@ -209,7 +209,8 @@ ProcessModule* ModuleObjectFactory::newProcessModule( return static_cast(newModuleObject(className)); } -DVCSProcessModule* ModuleObjectFactory::newDVCSProcessModule(unsigned int classId) { +DVCSProcessModule* ModuleObjectFactory::newDVCSProcessModule( + unsigned int classId) { return static_cast(newModuleObject(classId)); } @@ -235,17 +236,27 @@ ActiveFlavorsThresholdsModule* ModuleObjectFactory::newActiveFlavorsThresholdsMo ActiveFlavorsThresholdsModule* ModuleObjectFactory::newActiveFlavorsThresholdsModule( const std::string &className) { - return static_cast(newModuleObject(className)); + return static_cast(newModuleObject( + className)); } ScalesModule* ModuleObjectFactory::newScalesModule(unsigned int classId) { return static_cast(newModuleObject(classId)); } -ScalesModule* ModuleObjectFactory::newScalesModule(const std::string &className) { +ScalesModule* ModuleObjectFactory::newScalesModule( + const std::string &className) { return static_cast(newModuleObject(className)); } +MWFModule* ModuleObjectFactory::newMWFModule(unsigned int classId) { + return static_cast(newModuleObject(classId)); +} + +MWFModule* ModuleObjectFactory::newMWFModule(const std::string &className) { + return static_cast(newModuleObject(className)); +} + XiConverterModule* ModuleObjectFactory::newXiConverterModule( unsigned int classId) { return static_cast(newModuleObject(classId)); diff --git a/src/partons/beans/MesonType.cpp b/src/partons/beans/MesonType.cpp new file mode 100644 index 0000000000000000000000000000000000000000..6fc07831b8f182a76a205124afcd5e8e0e08e58c --- /dev/null +++ b/src/partons/beans/MesonType.cpp @@ -0,0 +1,61 @@ +#include "../../../include/partons/beans/MesonType.h" + +namespace PARTONS { + +MesonType::MesonType() : + m_type(MesonType::UNDEFINED) { +} + +MesonType::MesonType(const MesonType& other) { + m_type = other.m_type; +} + +MesonType::MesonType(Type type) : + m_type(type) { +} +MesonType::operator MesonType::Type() const { + return m_type; +} + +std::string MesonType::toString() const { + + switch (m_type) { + + case RHO0: + return "RHO0"; + break; + case PHI: + return "PHI"; + break; + default: + return "UNDEFINED"; + } +} + +std::string MesonType::getShortName() { + + switch (m_type) { + + case RHO0: + return "rho0"; + break; + case PHI: + return "phi"; + break; + default: + return "UNDEFINED"; + } +} + +MesonType::Type MesonType::getType() const { + return m_type; +} + +void MesonType::setType(Type type) { + m_type = type; +} + +MesonType::~MesonType() { +} + +} /* namespace PARTONS */ diff --git a/src/partons/modules/meson_wave_function/MWFGK19.cpp b/src/partons/modules/meson_wave_function/MWFGK19.cpp new file mode 100644 index 0000000000000000000000000000000000000000..b97f246851e22632f312d08765eba74ebeddcb7a --- /dev/null +++ b/src/partons/modules/meson_wave_function/MWFGK19.cpp @@ -0,0 +1,244 @@ +#include "../../../../include/partons/modules/meson_wave_function/MWFGK19.h" + +#include +#include +#include +#include +#include + +#include "../../../../include/partons/BaseObjectRegistry.h" +#include "../../../../include/partons/modules/running_alpha_strong/RunningAlphaStrongStandard.h" +#include "../../../../include/partons/ModuleObjectFactory.h" +#include "../../../../include/partons/Partons.h" + +namespace PARTONS { + +const unsigned int MWFGK19::classId = + BaseObjectRegistry::getInstance()->registerBaseObject( + new MWFGK19("MWFGK19")); + +MWFGK19::MWFGK19(const std::string &className) : + MWFModule(className), m_pRunningAlphaStrongModule(0), m_MuF2_ref(1.) { +} + +MWFGK19::MWFGK19(const MWFGK19 &other) : + MWFModule(other), m_MuF2_ref(other.m_MuF2_ref) { + + if (other.m_pRunningAlphaStrongModule != 0) { + m_pRunningAlphaStrongModule = + (other.m_pRunningAlphaStrongModule)->clone(); + } else { + m_pRunningAlphaStrongModule = 0; + } +} + +MWFGK19::~MWFGK19() { + + if (m_pRunningAlphaStrongModule != 0) { + delete m_pRunningAlphaStrongModule; + m_pRunningAlphaStrongModule = 0; + } +} + +MWFGK19* MWFGK19::clone() const { + return new MWFGK19(*this); +} + +void MWFGK19::resolveObjectDependencies() { + + //run for mother class + MWFModule::resolveObjectDependencies(); + + //set pointer + m_pRunningAlphaStrongModule = + Partons::getInstance()->getModuleObjectFactory()->newRunningAlphaStrongModule( + RunningAlphaStrongStandard::classId); +} + +double modulus(const NumA::Vector2D& b){ + + return sqrt(pow(b.getX(),2.0)+pow(b.getY(),2.0)); + +} + +double GegenbauerN2(double lambda, double x){ + + return -lambda+2*lambda*(1+lambda)*pow(x,2.0); + +} + +//double sudakov(double tau, const NumA::Vector2D& b, double Q2){ +// double Nf=3.0; +// double beta0 = 11 - 2./3.*Nf; +// double LambdaQCD = 0.22; +// double bhat = -log(modulus(b)*LambdaQCD); +// double qhat = log(tau*sqrt(Q2)/(sqrt(2)*LambdaQCD)); +// +// double maximumArray[3]; +// maximumArray[0] = tau*sqrt(Q2); +// maximumArray[1] = (1-tau)*sqrt(Q2); +// maximumArray[2] = 1./modulus(b); +// +// double muR = maximumArray[0]; +// +// for(int i=1;i<3;i++){ +// if(muR < maximumArray[i]) +// muR = maximumArray[i]; +// } +// +// double sudakov; +// +// if(modulus(b) > sqrt(2)/(tau*sqrt(Q2))){ +// sudakov = 8./(3.*beta0)*(qhat*log(qhat/bhat) - qhat + bhat); +// } +// else{ +// sudakov = 0.; +// } +// +// return sudakov; +//} + +//double Sudakov(double tau, const NumA::Vector2D& b, double Q2){ +// +// double Nf=3.0; +// double beta0 = 11 - 2./3.*Nf; +// double LambdaQCD = 0.22; +// double bhat = -log(modulus(b)*LambdaQCD); +// +// double maximumArray[3]; +// maximumArray[0] = tau*sqrt(Q2); +// maximumArray[1] = (1-tau)*sqrt(Q2); +// maximumArray[2] = 1./modulus(b); +// +// double muR = maximumArray[0]; +// +// for(int i=1;i<3;i++){ +// if(muR < maximumArray[i]) +// muR = maximumArray[i]; +// } +// +// +// double Sudakov = sudakov(tau, b, Q2) + sudakov(1 - tau, b, Q2) - 4/beta0*log(log(muR/LambdaQCD)/bhat); +// +// return Sudakov; +// +//} + +double Heaviside(double x){ + if(x < 0.0) + return 0.; + else + return 1.; +} + +//double CrossSection(double tau,double b){ +// double Nc = 3.; +// double x; +// double W = 75.; +// double Q2 = 4.0; +// double Mp = 0.938; +// double Mv = 0.0; +// double CScoefficient = 1./(16.*M_PI*(pow(W,2.) - pow(Mp,2.))*sqrt(pow(W,4.) + pow(Q2,2.) + pow(Mp,4.) + 2.*pow(W,2.)*Q2 - 2.*pow(W,2.)*pow(Mp,2.) + 2.*Q2*pow(Mp,2.))); +// double chargeU = 2./3.; +// double chargeD = -1./3.; +// double chargeS = -1./3.; +// double CPhiStrange = 1.0; +// +// double CF = (Nc^2 - 1.)/(2.*Nc); +// double xBj = Q2/(pow(W,2.) + Q2 - pow(Mp,2.)); +// double xi = xBj/(2. - xBj)*(1. + pow(Mv,2.)/Q2); +// double gsl_sf_bessel_K0 (sqrt(-(1. - tau)/(2.*xi)*(xi - x))*b*sqrt(Q2)); +// +// return 0.0; +// +//} + +// double PropagatorQuarkLong[x_, xi_, tau_, Q2_, b_] := CF*sqrt(2./Nc)*sqrt(Q2)/xi +// *(1./(2.*M_PI)*gsl_sf_bessel_K0 (sqrt(-(1. - tau)/(2.*xi)*(xi - x))*b*sqrt(Q2))*Heaviside(xi - x) +// + I/4*HankelH1[0, Sqrt[-(1 - tau)/(2*xi)*(x + xi)]*b*Sqrt[Q2]] * HeavisideTheta[x + xi]) +// - CF*Sqrt[2/Nc]*Sqrt[Q2]/xi +// *(1/(2*Pi)*BesselK[0, Sqrt[(1 - tau)/(2*xi)*(-xi - x)]*b*Sqrt[Q2]]*HeavisideTheta[-xi - x] +// + I/4*HankelH1[0, Sqrt[(1 - tau)/(2*xi)*(x - xi)]*b*Sqrt[Q2]]*HeavisideTheta[x - xi]); + + +double MWFGK19::compute(double tau, const NumA::Vector2D& b, double muF2, + MesonType::Type mesonType, double mesonHelicity) { + + double WF = 0.0; + + double Nc = 3.0; + double LambdaQCD = 0.22; + + double gamma2L = 50./81.; + double gamma2T = 40./81.; + + double rhof0L = 209.0; + double rhof0T = 167.0; + double phif0L = 221.0; + double phif0T = 177.0; + + + double rhofL = rhof0L; + double rhofT = rhof0T*pow(m_pRunningAlphaStrongModule->compute(muF2)/m_pRunningAlphaStrongModule->compute(m_MuF2_ref),4./27.); + double phifL = phif0L; + double phifT = phif0T*pow(m_pRunningAlphaStrongModule->compute(muF2)/m_pRunningAlphaStrongModule->compute(m_MuF2_ref),4./27.); + + double rhoaL = 0.75; + double rhoaT = 1.0; + double phiaL = 0.7; + double phiaT = 0.95; + + double rhoB20L = 0.0; + double rhoB20T = 0.1; + double phiB20L = 0.0; + double phiB20T = 0.1; + + double rhoB2L = rhoB20L*pow(m_pRunningAlphaStrongModule->compute(muF2)/m_pRunningAlphaStrongModule->compute(m_MuF2_ref),gamma2L); + double rhoB2T = rhoB20T*pow(m_pRunningAlphaStrongModule->compute(muF2)/m_pRunningAlphaStrongModule->compute(m_MuF2_ref),gamma2T); + double phiB2L = phiB20L*pow(m_pRunningAlphaStrongModule->compute(muF2)/m_pRunningAlphaStrongModule->compute(m_MuF2_ref),gamma2L); + double phiB2T = phiB20T*pow(m_pRunningAlphaStrongModule->compute(muF2)/m_pRunningAlphaStrongModule->compute(m_MuF2_ref),gamma2T); + + if(mesonType == MesonType::RHO0 && mesonHelicity == 1.){ + + debug(__func__, "some deb"); + info(__func__, "some info"); + warn(__func__, "some warn"); + + WF = 2*M_PI*sqrt(2*Nc) * rhofL*(tau*(1 - tau)) * (1 + rhoB2L * GegenbauerN2(3./2., (2*tau - 1))) * + exp(-(tau*(1 - tau))*pow(modulus(b),2.0)/(4*pow(rhoaL,2.0))); + + + } + else if(mesonType == MesonType::RHO0 && mesonHelicity == 2.){ + + WF = 2*M_PI*sqrt(2*Nc) * rhofT*(tau*(1 - tau)) * (1 + rhoB2T * GegenbauerN2(3./2., (2*tau - 1))) * + exp(-(tau*(1 - tau))*pow(modulus(b),2.0)/(4*pow(rhoaT,2.0))); + + } + else if(mesonType == MesonType::PHI && mesonHelicity == 1.){ + + WF = 2*M_PI*sqrt(2*Nc) * phifL*(tau*(1 - tau)) * (1 + phiB2L * GegenbauerN2(3./2., (2*tau - 1))) * + exp(-(tau*(1 - tau))*pow(modulus(b),2.0)/(4*pow(phiaL,2.0))); + + } + else if(mesonType == MesonType::PHI && mesonHelicity == 2.){ + + WF = 2*M_PI*sqrt(2*Nc) * phifT*(tau*(1 - tau)) * (1 + phiB2T * GegenbauerN2(3./2., (2*tau - 1))) * + exp(-(tau*(1 - tau))*pow(modulus(b),2.0)/(4*pow(phiaT,2.0))); + + } + else{ + throw ElemUtils::CustomException(getClassName(), __func__, ElemUtils::Formatter() << "XXX for meson " << MesonType(mesonType).toString() << " not implemented" ); + } + + + + return gsl_sf_bessel_K0 (-4.); + + //return 2*M_PI* sqrt(2.0 * Nc)* rhofL *(tau*(1-tau))*exp(-1.0*tau*(1.0-tau)*(pow(b.getX(),2.0))/(4*pow(rhoaL,2.0))); + + //return 2*M_PI* sqrt(2.0 * Nc)* rhofL *(tau*(1-tau))*exp(-1.0*tau*(1.0-tau)*(bX^2+bY^2)/(4*pow(rhoaL,2.0))); + //return m_pRunningAlphaStrongModule->compute(muF2) * gsl_sf_bessel_J0(0.1); +} + +} /* namespace PARTONS */ diff --git a/src/partons/modules/meson_wave_function/MWFModule.cpp b/src/partons/modules/meson_wave_function/MWFModule.cpp new file mode 100644 index 0000000000000000000000000000000000000000..890390e21e4d34e260a0daffab6a1dbc23561599 --- /dev/null +++ b/src/partons/modules/meson_wave_function/MWFModule.cpp @@ -0,0 +1,22 @@ +#include "../../../../include/partons/modules/meson_wave_function/MWFModule.h" + +namespace PARTONS { + +MWFModule::MWFModule(const std::string &className) : + ModuleObject(className) { +} + +MWFModule::MWFModule(const MWFModule &other) : + ModuleObject(other) { +} + +MWFModule::~MWFModule() { +} + +void MWFModule::initModule() { +} + +void MWFModule::isModuleWellConfigured() { +} + +} /* namespace PARTONS */