Skip to content
GitLab
Commit 2f622853 authored by Pawel Sznajder's avatar Pawel Sznajder
Browse files

add new observable modules

parent ae3d1037
Hide whitespace changes
Inline Side-by-side
include/partons/modules/observable/DVMP/cross_section/DVMPCrossSectionTotal.h 0 → 100644
View file @ 2f622853
#ifndef DVMPCROSSSECTIONTOTAL_H
#define DVMPCROSSSECTIONTOTAL_H
/**
* @file DVMPCrossSectionTotal.h
* @author Pawel Sznajder (NCBJ)
* @date May 17, 2020
* @version 1.0
*/
#include <ElementaryUtils/parameters/Parameters.h>
#include <stddef.h>
#include <string>
#include <utility>
#include "../../../../beans/gpd/GPDType.h"
#include "../../../../beans/List.h"
#include "../../../../beans/MesonType.h"
#include "../../../../utils/type/PhysicalType.h"
#include "DVMPCrossSectionUUUMinusPhiIntegrated.h"
namespace PARTONS {
/**
* @class DVMPCrossSectionTotal
* @brief Unpolarized cross-section for electro-production integrated over \f$xB\f$, \f$Q^{2}\f$, \f$t\f$ and angles.
*
* Default ranges: \f$1E-4 < xB < 0.7\f$, \f$-1 < t < 0\f$, \f$1 < Q^{2} < 1E3\f$ and \f$0 < y < 1\f$. Full integration over angles.
*
* Unit: \f$\mathrm{nbarn}\f$.
*/
class DVMPCrossSectionTotal: public DVMPCrossSectionUUUMinusPhiIntegrated {
public:
static const std::string DVMP_CROSSSECTION_TOTAL_RANGEXb; ///< String used to set integration xB range via XML scenario.
static const std::string DVMP_CROSSSECTION_TOTAL_RANGET; ///< String used to set integration t range via XML scenario.
static const std::string DVMP_CROSSSECTION_TOTAL_RANGEQ2; ///< String used to set integration Q2 range via XML scenario.
static const std::string DVMP_CROSSSECTION_TOTAL_RANGEY; ///< String used to set integration y range via XML scenario.
static const std::string DVMP_CROSSSECTION_TOTAL_RANGENu; ///< String used to set integration nu range via XML scenario.
static const std::string DVMP_CROSSSECTION_TOTAL_N0; ///< String used to set number of MC integration iterations per cycle via XML scenario.
static const std::string DVMP_CROSSSECTION_TOTAL_N1; ///< String used to set number of MC integration cycles via XML scenario.
static const std::string DVMP_CROSSSECTION_TOTAL_VIRTUALPHOTOPRODUCTION; ///< String used to set if virtual photo-production should be evaluated.
/**
* Unique ID to automatically register the class in the registry.
*/
static const unsigned int classId;
/**
* Function for the integration.
*/
static double DVMPCrossSectionTotalFunction(double* kin, size_t dim,
void* par);
/**
* Constructor.
* @param className Name of class.
*/
DVMPCrossSectionTotal(const std::string &className);
/**
* Destructor.
*/
virtual ~DVMPCrossSectionTotal();
virtual DVMPCrossSectionTotal* clone() const;
virtual void configure(const ElemUtils::Parameters &parameters);
size_t getNI0() const;
void setNI0(size_t nI0);
size_t getNI1() const;
void setNI1(size_t nI1);
const std::pair<double, double>& getRangeQ2() const;
void setRangeQ2(const std::pair<double, double>& rangeQ2);
const std::pair<double, double>& getRangeT() const;
void setRangeT(const std::pair<double, double>& rangeT);
const std::pair<double, double>& getRangexB() const;
void setRangexB(const std::pair<double, double>& rangexB);
const std::pair<double, double>& getRangeY() const;
void setRangeY(const std::pair<double, double>& rangeY);
const std::pair<double, double>& getRangeNu() const;
void setRangeNu(const std::pair<double, double>& rangeNu);
bool isEvaluatePhotoProduction() const;
void setEvaluatePhotoProduction(bool evaluatePhotoProduction);
protected:
/**
* Copy constructor.
* @param other Object to be copied.
*/
DVMPCrossSectionTotal(const DVMPCrossSectionTotal &other);
virtual PhysicalType<double> computeObservable(
const DVMPObservableKinematic& kinematic,
const List<GPDType>& gpdType);
private:
/**
* Print change of ranges.
*/
void printChangeOfRange(const std::string& func, const std::string& name,
const std::pair<double, double>& oldValues,
const std::pair<double, double>& newValues) const;
/**
* Parse range.
*/
std::pair<double, double> parseRange(const std::string& str) const;
size_t m_nI0; ///< Number of iteration in single cycle.
size_t m_nI1; ///< Number of cycles.
std::pair<double, double> m_rangexB; ///< xB integration range.
std::pair<double, double> m_rangeT; ///< t integration range.
std::pair<double, double> m_rangeQ2; ///< Q2 integration range.
std::pair<double, double> m_rangeY; ///< y integration range.
std::pair<double, double> m_rangeNu; ///< nu integration range.
bool m_evaluatePhotoProduction; ///< If true, virtual photo-production will be evaluated.
};
struct DVMPCrossSectionTotalParameters {
DVMPCrossSectionTotal* m_pDVMPCrossSectionTotal; ///< Pointer to DVMPCrossSectionTotal.
double m_E; ///< Beam energy.
MesonType::Type m_mesonType; ///< Meson type.
List<GPDType> m_gpdType; ///< GPD types.
std::pair<double, double> m_yCut; ///< y cut.
std::pair<double, double> m_nuCut; ///< nu cut.
bool m_evaluatePhotoProduction; ///< If true, virtual photo-production will be evaluated.
};
} /* namespace PARTONS */
#endif /* DVMPCROSSSECTIONTOTAL_H */
include/partons/modules/observable/DVMP/cross_section/DVMPCrossSectionUUMinusPi0.h → include/partons/modules/observable/DVMP/cross_section/DVMPCrossSectionUUUMinus.h
View file @ 2f622853
#ifndef DVMPCROSSSECTIONUUMINUSPI0_H
#define DVMPCROSSSECTIONUUMINUSPI0_H
#ifndef DVMPCROSSSECTIONUUUMINUS_H
#define DVMPCROSSSECTIONUUUMINUS_H
/**
* @file DVMPCrossSectionUUMinusPi0.h
* @file DVMPCrossSectionUUUMinus.h
* @author Pawel Sznajder (IPNO)
* @date November 25, 2016
* @version 1.0
......@@ -18,12 +18,12 @@
namespace PARTONS {
/**
* @class DVMPCrossSectionUUMinusPi0
* @brief Unpolarized cross-section for negative beam charge - pi0.
* @class DVMPCrossSectionUUUMinus
* @brief Unpolarized cross-section for negative beam charge.
*
* Unit: \f$\mathrm{nbarn}/\mathrm{GeV}^4\f$.
*/
class DVMPCrossSectionUUMinusPi0: public DVMPObservable {
class DVMPCrossSectionUUUMinus: public DVMPObservable {
public:
......@@ -36,14 +36,14 @@ public:
* Constructor.
* @param className Name of class.
*/
DVMPCrossSectionUUMinusPi0(const std::string &className);
DVMPCrossSectionUUUMinus(const std::string &className);
/**
* Destructor.
*/
virtual ~DVMPCrossSectionUUMinusPi0();
virtual ~DVMPCrossSectionUUUMinus();
virtual DVMPCrossSectionUUMinusPi0* clone() const;
virtual DVMPCrossSectionUUUMinus* clone() const;
protected:
......@@ -51,7 +51,7 @@ protected:
* Copy constructor.
* @param other Object to be copied.
*/
DVMPCrossSectionUUMinusPi0(const DVMPCrossSectionUUMinusPi0 &other);
DVMPCrossSectionUUUMinus(const DVMPCrossSectionUUUMinus &other);
virtual PhysicalType<double> computeObservable(
const DVMPObservableKinematic& kinematic,
......@@ -60,4 +60,4 @@ protected:
} /* namespace PARTONS */
#endif /* DVMPCROSSSECTIONUUMINUSPI0_H */
#endif /* DVMPCROSSSECTIONUUUMINUS_H */
include/partons/modules/observable/DVMP/cross_section/DVMPCrossSectionUUUMinusPhiIntegrated.h 0 → 100644
View file @ 2f622853
#ifndef DVMPCROSSSECTIONUUUMINUSPHIINTEGRATED_H
#define DVMPCROSSSECTIONUUUMINUSPHIINTEGRATED_H
/**
* @file DVMPCrossSectionUUUMinusPhiIntegrated.h
* @author Luca COLANERI (IPNO)
* @date April 24, 2017
* @version 1.0
*/
#include <ElementaryUtils/parameters/Parameters.h>
#include <string>
#include <vector>
#include "../../../../beans/gpd/GPDType.h"
#include "../../../../beans/List.h"
#include "../../../../utils/type/PhysicalType.h"
#include "../../../MathIntegratorModule.h"
#include "DVMPCrossSectionUUUMinus.h"
namespace PARTONS {
/**
* @class DVMPCrossSectionUUUMinusPhiIntegrated
* @brief Unpolarized cross-section for electro-production integrated over \f$\phi\f$.
*
* Unit: \f$\mathrm{nbarn}/\mathrm{GeV}^4\f$.
*/
class DVMPCrossSectionUUUMinusPhiIntegrated: public DVMPCrossSectionUUUMinus,
public MathIntegratorModule {
public:
/**
* Unique ID to automatically register the class in the registry.
*/
static const unsigned int classId;
/**
* Constructor.
* @param className Name of class.
*/
DVMPCrossSectionUUUMinusPhiIntegrated(const std::string &className);
/**
* Destructor.
*/
virtual ~DVMPCrossSectionUUUMinusPhiIntegrated();
virtual DVMPCrossSectionUUUMinusPhiIntegrated* clone() const;
virtual void configure(const ElemUtils::Parameters &parameters);
protected:
/**
* Copy constructor.
* @param other Object to be copied.
*/
DVMPCrossSectionUUUMinusPhiIntegrated(
const DVMPCrossSectionUUUMinusPhiIntegrated &other);
virtual PhysicalType<double> computeObservable(
const DVMPObservableKinematic& kinematic,
const List<GPDType>& gpdType);
/**
* Functor to perform the integration.
*/
NumA::FunctionType1D* m_pFunctionToIntegrateObservable;
/**
* Function to be integrated.
*/
virtual double functionToIntegrateObservable(double x,
std::vector<double> params);
/**
* Initialize functors.
*/
void initFunctorsForIntegrations();
};
} /* namespace PARTONS */
#endif /* DVMPCROSSSECTIONUUUMINUSPHIINTEGRATED_H */
include/partons/modules/process/DVMP/DVMPProcessGK06.h
View file @ 2f622853
#ifndef DVMP_PROCESS_VGG99_H
#define DVMP_PROCESS_VGG99_H
#ifndef DVMP_PROCESS_GK06_H
#define DVMP_PROCESS_GK06_H
/**
* @file DVMPProcessGK06.h
* @author Michel Guidal (IPNO)
* @author Pawel Sznajder (IPNO)
* @author Kemal Tezgin
* @date October 21, 2019
* @version 1.0
*/
......@@ -19,15 +19,7 @@ namespace PARTONS {
/**
* @class DVMPProcessGK06
*
* VGG process model for DVMP.
*
* For the reference see:
* - Prog. Part. Nucl. Phys. 47, 401 (2001)
* - Phys. Rev. Lett. 80 5064 (1998).
* - Phys. Rev. D 60, 094017 (1999).
* - Phys. Rev. D 72, 054013 (2005).
*
* Module based on the original code received from M. Guidal as a private communication.
* TODO
*/
class DVMPProcessGK06: public DVMPProcessModule {
......@@ -92,4 +84,4 @@ private:
} /* namespace PARTONS */
#endif /* DVMP_PROCESS_VGG99_H */
#endif /* DVMP_PROCESS_GK06_H */
src/partons/modules/observable/DVCS/cross_section/DVCSCrossSectionTotal.cpp
View file @ 2f622853
......@@ -26,9 +26,9 @@ const std::string DVCSCrossSectionTotal::DVCS_CROSSSECTION_TOTAL_RANGEQ2 =
const std::string DVCSCrossSectionTotal::DVCS_CROSSSECTION_TOTAL_RANGEY =
"DVCSCrossSectionTotal_rangeY";
const std::string DVCSCrossSectionTotal::DVCS_CROSSSECTION_TOTAL_N0 =
"DVCSCrossSectionTotal_rangeN0";
"DVCSCrossSectionTotal_N0";
const std::string DVCSCrossSectionTotal::DVCS_CROSSSECTION_TOTAL_N1 =
"DVCSCrossSectionTotal_rangeN1";
"DVCSCrossSectionTotal_N1";
const unsigned int DVCSCrossSectionTotal::classId =
BaseObjectRegistry::getInstance()->registerBaseObject(
......
src/partons/modules/observable/DVMP/cross_section/DVMPCrossSectionTotal.cpp 0 → 100644
View file @ 2f622853
#include "../../../../../../include/partons/modules/observable/DVMP/cross_section/DVMPCrossSectionTotal.h"
#include <ElementaryUtils/logger/CustomException.h>
#include <ElementaryUtils/parameters/GenericType.h>
#include <ElementaryUtils/string_utils/Formatter.h>
#include <ElementaryUtils/string_utils/StringUtils.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_vegas.h>
#include <gsl/gsl_rng.h>
#include <cmath>
#include <iostream>
#include <vector>
#include "../../../../../../include/partons/beans/observable/DVMP/DVMPObservableKinematic.h"
#include "../../../../../../include/partons/BaseObjectRegistry.h"
#include "../../../../../../include/partons/FundamentalPhysicalConstants.h"
#include "../../../../../../include/partons/modules/process/DVMP/DVMPProcessModule.h"
#include "../../../../../../include/partons/modules/xi_converter/DVMP/DVMPXiConverterModule.h"
#include "../../../../../../include/partons/utils/type/PhysicalUnit.h"
namespace PARTONS {
const std::string DVMPCrossSectionTotal::DVMP_CROSSSECTION_TOTAL_RANGEXb =
"DVMPCrossSectionTotal_rangeXb";
const std::string DVMPCrossSectionTotal::DVMP_CROSSSECTION_TOTAL_RANGET =
"DVMPCrossSectionTotal_rangeT";
const std::string DVMPCrossSectionTotal::DVMP_CROSSSECTION_TOTAL_RANGEQ2 =
"DVMPCrossSectionTotal_rangeQ2";
const std::string DVMPCrossSectionTotal::DVMP_CROSSSECTION_TOTAL_RANGEY =
"DVMPCrossSectionTotal_rangeY";
const std::string DVMPCrossSectionTotal::DVMP_CROSSSECTION_TOTAL_RANGENu =
"DVMPCrossSectionTotal_rangeNu";
const std::string DVMPCrossSectionTotal::DVMP_CROSSSECTION_TOTAL_N0 =
"DVMPCrossSectionTotal_N0";
const std::string DVMPCrossSectionTotal::DVMP_CROSSSECTION_TOTAL_N1 =
"DVMPCrossSectionTotal_N1";
const std::string DVMPCrossSectionTotal::DVMP_CROSSSECTION_TOTAL_VIRTUALPHOTOPRODUCTION =
"DVMPCrossSectionTotal_virtPhotoProd";
const unsigned int DVMPCrossSectionTotal::classId =
BaseObjectRegistry::getInstance()->registerBaseObject(
new DVMPCrossSectionTotal("DVMPCrossSectionTotal"));
DVMPCrossSectionTotal::DVMPCrossSectionTotal(const std::string &className) :
DVMPCrossSectionUUUMinusPhiIntegrated(className), m_nI0(1000), m_nI1(5), m_evaluatePhotoProduction(
false) {
m_rangexB = std::pair<double, double>(1.E-4, 0.7);
m_rangeT = std::pair<double, double>(-1., 0.);
m_rangeQ2 = std::pair<double, double>(1., 1.E3);
m_rangeY = std::pair<double, double>(0., 1.);
m_rangeNu = std::pair<double, double>(0., 1.E3);
}
DVMPCrossSectionTotal::DVMPCrossSectionTotal(const DVMPCrossSectionTotal& other) :
DVMPCrossSectionUUUMinusPhiIntegrated(other), m_nI0(other.m_nI0), m_nI1(
other.m_nI1), m_evaluatePhotoProduction(
other.m_evaluatePhotoProduction) {
m_rangexB = other.m_rangexB;
m_rangeT = other.m_rangeT;
m_rangeQ2 = other.m_rangeQ2;
m_rangeY = other.m_rangeY;
m_rangeNu = other.m_rangeNu;
}
DVMPCrossSectionTotal::~DVMPCrossSectionTotal() {
}
DVMPCrossSectionTotal* DVMPCrossSectionTotal::clone() const {
return new DVMPCrossSectionTotal(*this);
}
void DVMPCrossSectionTotal::configure(const ElemUtils::Parameters &parameters) {
//run for base
DVMPCrossSectionUUUMinusPhiIntegrated::configure(parameters);
//check
if (parameters.isAvailable(DVMP_CROSSSECTION_TOTAL_RANGEXb)) {
setRangexB(parseRange(parameters.getLastAvailable().getString()));
}
if (parameters.isAvailable(DVMP_CROSSSECTION_TOTAL_RANGET)) {
setRangeT(parseRange(parameters.getLastAvailable().getString()));
}
if (parameters.isAvailable(DVMP_CROSSSECTION_TOTAL_RANGEQ2)) {
setRangeQ2(parseRange(parameters.getLastAvailable().getString()));
}
if (parameters.isAvailable(DVMP_CROSSSECTION_TOTAL_RANGEY)) {
setRangeY(parseRange(parameters.getLastAvailable().getString()));
}
if (parameters.isAvailable(DVMP_CROSSSECTION_TOTAL_RANGENu)) {
setRangeNu(parseRange(parameters.getLastAvailable().getString()));
}
if (parameters.isAvailable(DVMP_CROSSSECTION_TOTAL_N0)) {
setNI0(
ElemUtils::GenericType(
parameters.getLastAvailable().getString()).toUInt());
}
if (parameters.isAvailable(DVMP_CROSSSECTION_TOTAL_N1)) {
setNI1(
ElemUtils::GenericType(
parameters.getLastAvailable().getString()).toUInt());
}
if (parameters.isAvailable(
DVMP_CROSSSECTION_TOTAL_VIRTUALPHOTOPRODUCTION)) {
setEvaluatePhotoProduction(
ElemUtils::GenericType(
parameters.getLastAvailable().getString()).toBoolean());
}
}
double DVMPCrossSectionTotal::DVMPCrossSectionTotalFunction(double* kin,
size_t dim, void* par) {
//parameters
DVMPCrossSectionTotalParameters* params =
static_cast<DVMPCrossSectionTotalParameters*>(par);
double xB = kin[0];
double t = kin[1];
double Q2 = kin[2];
double E = params->m_E;
MesonType::Type mesonType = params->m_mesonType;
//check kinematic range
double y = Q2 / (2 * xB * Constant::PROTON_MASS * E);
double nu = y * E;
double xi =
params->m_pDVMPCrossSectionTotal->getProcessModule()->getXiConverterModule()->compute(
DVMPObservableKinematic(xB, t, Q2, E, 0., mesonType)).getValue();
double tmin = -4. * pow(Constant::PROTON_MASS, 2.) * pow(xi, 2.)
/ (1. - pow(xi, 2.));
double tmax = -1.E6;
double xBmin = 2 * Q2 * E / Constant::PROTON_MASS / (4 * E * E - Q2);
if (xB < xBmin || xB > 1.)
return 0.;
if (t > tmin || t < tmax)
return 0.;
if (y < params->m_yCut.first || y > params->m_yCut.second)
return 0.;
if (nu < params->m_nuCut.first || nu > params->m_nuCut.second)
return 0.;
//flux
double flux = 1.;
if (params->m_evaluatePhotoProduction) {
flux *= Constant::FINE_STRUCTURE_CONSTANT * (1. - xB);
flux /= 2 * M_PI * Q2 * y * E;
flux *= pow(y, 2) * (1. - (2 * pow(Constant::ELECTRON_MASS, 2)) / Q2)
+ 2. / (1. + Q2 / pow(nu, 2)) * (1. - y - Q2 / pow(2 * E, 2));
flux *= nu / xB;
}
//evaluate
double result =
params->m_pDVMPCrossSectionTotal->DVMPCrossSectionUUUMinusPhiIntegrated::computeObservable(
DVMPObservableKinematic(xB, t, Q2, E, 0., mesonType),
params->m_gpdType).getValue();
return result / flux;
}
PhysicalType<double> DVMPCrossSectionTotal::computeObservable(
const DVMPObservableKinematic& kinematic,
const List<GPDType>& gpdType) {
//result and error
double res, err;
//parameters
DVMPCrossSectionTotalParameters params;
params.m_pDVMPCrossSectionTotal = this;
params.m_E = kinematic.getE().getValue();
params.m_mesonType = kinematic.getMesonType();
params.m_gpdType = gpdType;
//Q2 max
double maxQ2 = 2 * Constant::PROTON_MASS * kinematic.getE().getValue();
if (m_rangeQ2.second > maxQ2) {
warn(__func__,
ElemUtils::Formatter() << "Upper limit of Q2 "
<< m_rangeQ2.second
<< " grater than maximum value allowed for this energy, "
<< maxQ2 << ", the later used instead");
}
//ranges
// xB, Q2, t
double xl[3] = { m_rangexB.first, m_rangeT.first, m_rangeQ2.first };
double xu[3] = { m_rangexB.second, m_rangeT.second, (
(m_rangeQ2.second > maxQ2) ? (maxQ2) : (m_rangeQ2.second)) };
// y and nu cut
params.m_yCut = m_rangeY;
params.m_nuCut = m_rangeNu;
// virtual photo-production
params.m_evaluatePhotoProduction = m_evaluatePhotoProduction;
//random
const gsl_rng_type* T = gsl_rng_default;
gsl_rng* r = gsl_rng_alloc(T);
//function
gsl_monte_function G = { &DVMPCrossSectionTotalFunction, 3, &params };
//run
gsl_monte_vegas_state *s = gsl_monte_vegas_alloc(3);
for (size_t i = 0; i < m_nI1; i++) {
gsl_monte_vegas_integrate(&G, xl, xu, 3, m_nI0, r, s, &res, &err);
info(__func__,
ElemUtils::Formatter() << "Intermediate result: cycle: " << i
<< " result: " << res << " error: " << err << " [nb]");
}
//free
gsl_monte_vegas_free(s);
gsl_rng_free(r);
//return
return PhysicalType<double>(res, PhysicalUnit::NB);
}
void DVMPCrossSectionTotal::printChangeOfRange(const std::string& func,
const std::string& name, const std::pair<double, double>& oldValues,
const std::pair<double, double>& newValues) const {
info(func,