forked from cepgen/cepgen
-
Notifications
You must be signed in to change notification settings - Fork 0
/
IntegratorPlain.cpp
74 lines (61 loc) · 2.68 KB
/
IntegratorPlain.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
/*
* CepGen: a central exclusive processes event generator
* Copyright (C) 2013-2023 Laurent Forthomme
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <gsl/gsl_monte_plain.h>
#include "CepGen/Core/Exception.h"
#include "CepGen/Integration/IntegratorGSL.h"
#include "CepGen/Modules/IntegratorFactory.h"
namespace cepgen {
/// Plain integration algorithm randomly sampling points in the phase space
class IntegratorPlain final : public IntegratorGSL {
public:
explicit IntegratorPlain(const ParametersList& params);
static ParametersDescription description();
Value integrate(Integrand&) override;
private:
const int ncvg_;
};
IntegratorPlain::IntegratorPlain(const ParametersList& params)
: IntegratorGSL(params), ncvg_(steer<int>("numFunctionCalls")) {}
Value IntegratorPlain::integrate(Integrand& integrand) {
setIntegrand(integrand);
//--- launch integration
std::unique_ptr<gsl_monte_plain_state, void (*)(gsl_monte_plain_state*)> pln_state(
gsl_monte_plain_alloc(function_->dim), gsl_monte_plain_free);
double result, abserr;
int res = gsl_monte_plain_integrate(function_.get(),
&xlow_[0],
&xhigh_[0],
function_->dim,
ncvg_,
gsl_rng_.get(),
pln_state.get(),
&result,
&abserr);
if (res != GSL_SUCCESS)
throw CG_FATAL("Integrator:integrate") << "Error while performing the integration!\n\t"
<< "GSL error: " << gsl_strerror(res) << ".";
return Value{result, abserr};
}
ParametersDescription IntegratorPlain::description() {
auto desc = IntegratorGSL::description();
desc.setDescription("Plain (trial/error) integrator");
desc.add<int>("numFunctionCalls", 50000);
return desc;
}
} // namespace cepgen
REGISTER_INTEGRATOR("plain", IntegratorPlain);