MellinMomentFromGPD.cpp 5.46 KB
Newer Older
1
2
#include "../../include/mellin_moment/MellinMomentFromGPD.h"

3
#include <ElementaryUtils/logger/CustomException.h>
4
5
6
7
8
#include <ElementaryUtils/parameters/Parameters.h>
#include <ElementaryUtils/string_utils/Formatter.h>
#include <NumA/functor/one_dimension/Functor1D.h>
#include <NumA/integration/one_dimension/Integrator1D.h>
#include <NumA/integration/one_dimension/IntegratorType1D.h>
9
#include <partons/beans/automation/BaseObjectData.h>
10
11
12
13
14
15
#include <partons/beans/gpd/GPDKinematic.h>
#include <partons/beans/parton_distribution/GluonDistribution.h>
#include <partons/beans/parton_distribution/PartonDistribution.h>
#include <partons/beans/parton_distribution/QuarkDistribution.h>
#include <partons/BaseObjectRegistry.h>
#include <partons/ModuleObjectFactory.h>
16
#include <partons/Partons.h>
17
#include <cmath>
18
#include <map>
19
20
21
22
23
24
25
26
27
28
29
#include <string>


namespace PARTONS {

const unsigned int MellinMomentFromGPD::classId =
		BaseObjectRegistry::getInstance()->registerBaseObject(
				new MellinMomentFromGPD("MellinMomentFromGPD"));

const std::string MellinMomentFromGPD::MELLIN_MOMENT_MODULE_CLASS_NAME =
		"MellinMomentFromGPD";
30
31
const std::string MellinMomentFromGPD::MELLIN_MOMENT_MODULE_GPD_MODEL =
		"GPD_MODULE_GPD_MODEL";
32
33
34
35
36
37
38

MellinMomentFromGPD::MellinMomentFromGPD(const std::string &className) :
		MellinMomentModule(className) {
	m_pGPDModel = 0;
	m_pint = NumA::Integrator1D::newIntegrationFunctor(this,
			&MellinMomentFromGPD::integrant);

39
	initModule();
40
41
42
43
44
}

MellinMomentFromGPD::MellinMomentFromGPD(const MellinMomentFromGPD& other) :
		MellinMomentModule(other), MathIntegratorModule(other) {

45
	m_pGPDModel = other.m_pGPDModel;
46
47
48
49

	m_pint = NumA::Integrator1D::newIntegrationFunctor(this,
			&MellinMomentFromGPD::integrant);

50
51
	initModule();

52
53
54
55
56
57
58
}

MellinMomentFromGPD* MellinMomentFromGPD::clone() const {
	return new MellinMomentFromGPD(*this);
}

void MellinMomentFromGPD::isModuleWellConfigured() {
59
60
61
62
63
64
	MellinMomentModule::isModuleWellConfigured();

	if (m_pGPDModel == 0)
	throw ElemUtils::CustomException(getClassName(), __func__,
	                ElemUtils::Formatter()
	                        << "GPDModule not provided.");
65
}
66

67
void MellinMomentFromGPD::initModule() {
68
69
70
71
72
	MellinMomentModule::initModule();

	NumA::IntegratorType1D::Type integratorType = NumA::IntegratorType1D::DEXP;

	setIntegrator(integratorType);
73
74
75
76
77
78
79
}

void MellinMomentFromGPD::configure(const ElemUtils::Parameters &parameters) {
	MathIntegratorModule::configureIntegrator(parameters);
	MellinMomentModule::configure(parameters);
}

80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
void MellinMomentFromGPD::prepareSubModules(
            const std::map<std::string, BaseObjectData>& subModulesData)
{

    ModuleObject::prepareSubModules(subModulesData);

    std::map<std::string, BaseObjectData>::const_iterator it;

    if (m_pGPDModel != 0) setGPDModule(0);

    it = subModulesData.find(
            GPDModule::GPD_MODULE_CLASS_NAME);

    if (it != subModulesData.end()) {

       	m_pGPDModel = Partons::getInstance()->getModuleObjectFactory()->newGPDModule(
                            (it->second).getModuleClassName());

		info(__func__,
				ElemUtils::Formatter() << "GPD Module is set to: "
						<< m_pGPDModel->getClassName());

		m_pGPDModel->configure((it->second).getParameters());
		m_pGPDModel->prepareSubModules(
                    (it->second).getSubModules());
    } else
    {
    	throw ElemUtils::CustomException(getClassName(), __func__,
    	                ElemUtils::Formatter()
    	                        << "GPDModule not provided.");
	}
}

113
114
115
116
117
118
119
120
double MellinMomentFromGPD::compute(MellinMomentKinematic mKinematic) ///< Compute when everything is set.
		{
	std::vector<double> parameters;
	parameters.push_back(mKinematic.getXi());
	parameters.push_back(mKinematic.getT());
	parameters.push_back(mKinematic.getMuF2());
	parameters.push_back(mKinematic.getMuR2());

121
122
	isModuleWellConfigured();

123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
	return integrate(m_pint, -1.0, 1.0, parameters);

}

double MellinMomentFromGPD::integrant(double x, std::vector<double> par) ///< Integrand
		{

	// Create a GPDKinematic(x, xi, t, MuF2, MuR2) to compute
	GPDKinematic gpdKinematic(x, par[0], par[1], par[2], par[3]);

	if (m_gluon) {
		return pow(x, m_n - 2)
				* m_pGPDModel->compute(gpdKinematic, m_type).getGluonDistribution().getGluonDistribution();
	} else {
		return pow(x, m_n - 1)
				* m_pGPDModel->compute(gpdKinematic, m_type).getQuarkDistribution(
						m_flavor).getQuarkDistribution();
	}
}

143
List<QuarkFlavor> MellinMomentFromGPD::getListOfAvailableQuarkFlavor(
144
145
146
		MellinMomentKinematic mKinematic, const GPDType &gpdType) {

	List<QuarkFlavor> result;
147
	List<GPDType> gpdTypeList = getListOfAvailableGPDType();
148
149
	bool isGPDType = false;

150

151
152
153
154
155
156
	for (unsigned int i = 0; i != gpdTypeList.size(); i++) {
		if (gpdTypeList[i].getType() == gpdType.getType())
			isGPDType = true;
	}

	if (isGPDType) {
157
158
159
		GPDKinematic kinematic(1., mKinematic.getXi(), mKinematic.getT(),
				mKinematic.getMuF2(), mKinematic.getMuR2());

160
161
		List<QuarkDistribution> list =
				m_pGPDModel->compute(kinematic, gpdType).getListOfQuarkDistribution();
162

163
164
165
166
167
168
169
170
171
		for (unsigned int i = 0; i != list.size(); i++) {
			QuarkFlavor flavor = list[i].getQuarkFlavor();
			result.add(flavor);
		}
	}

	return result;
}

172
173
174
175
List<GPDType> MellinMomentFromGPD::getListOfAvailableGPDType(){
	return m_pGPDModel -> getListOfAvailableGPDTypeForComputation();
}

176
177
178
179
180
181
182
183
184
185
186
GPDModule* MellinMomentFromGPD::getGPDModule() const {
	return m_pGPDModel;
}

void MellinMomentFromGPD::setGPDModule(GPDModule* pGPDModel) {
	m_pModuleObjectFactory->updateModulePointerReference(m_pGPDModel,
			pGPDModel);
	m_pGPDModel = pGPDModel;
}

} /* namespace PARTONS */