27 using TMath::PoissonI;
29 static const int NUMBER_OF_DIMENSIONS = 3;
35 typedef JArray<NUMBER_OF_DIMENSIONS, double> point_type;
42 public JArray<1 + 2 * NUMBER_OF_DIMENSIONS, double>
48 JArray<1 + 2 * NUMBER_OF_DIMENSIONS, double>()
58 template<
class ...Args>
59 model_type(argument_type value,
const Args& ...
args) :
60 JArray<1 + 2 * NUMBER_OF_DIMENSIONS, double>(value,
args...)
70 inline double operator()(
const point_type& point)
const
72 double value = (*this)[0];
74 for (
int i = 0; i != NUMBER_OF_DIMENSIONS;++i) {
75 value *=
gauss(point[i] - (*
this)[2*i + 1], (*
this)[2*i + 2]);
89 static inline double gauss(
const double x,
const double sigma)
114 hit_type(
const point_type& point,
const double value,
const double sigma) :
132 inline double chi2(
const model_type&
model,
const hit_type& hit)
134 const double u = (
model(hit) - hit.value) / hit.sigma;
147 inline double likelihood(
const model_type&
model,
const hit_type& hit)
149 return -log(PoissonI(hit.value,
model(hit)));
160 int main(
int argc,
char **argv)
175 JParser<> zap(
"Program to test JSimplex algorithm.");
185 catch(
const exception& error) {
186 FATAL(error.what() << endl);
190 gRandom->SetSeed(seed);
193 const double top = parameters[0];
194 const double sigma = parameters[1];
200 for (
int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
213 double (*fit)(
const model_type&,
const hit_type&) = likelihood;
218 TH1D h0(
"chi2/NDF", NULL, 200, 0.0, 10.0);
220 H1.push_back(
new TH1D(
MAKE_CSTRING(
"p" << 0), NULL, 101, -20.0, +20.0));
222 for (
size_t i =1; i !=
model.size(); ++i) {
223 H1.push_back(
new TH1D(
MAKE_CSTRING(
"p" << i), NULL, 101, -0.1, +0.1));
227 for (
int counter = 0; counter != numberOfEvents; ++counter) {
229 STATUS(
"event: " << setw(10) << counter <<
'\r');
DEBUG(endl);
238 const JArray<NUMBER_OF_DIMENSIONS, double>
p1 = *
g1;
240 const double value = gRandom->Poisson(
model(
p1));
241 const double sigma = (value > 1.0 ? sqrt(value) : 0.7);
253 for (
const auto& hit :
data) {
255 for (
int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
256 Q[i].
put(hit[i], hit.value);
259 if (hit.value > value) {
261 simplex.
value[0] = hit.value;
263 for (
int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
264 simplex.
value[2*i + 1] = hit[i];
271 for (
int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
278 simplex.
step.resize(2*NUMBER_OF_DIMENSIONS + 1);
282 simplex.
step[0] = model_type();
283 simplex.
step[0][0] = 1.0;
285 for (
int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
287 simplex.
step[2*i + 1] = model_type();
288 simplex.
step[2*i + 1][2*i + 1] = 0.1;
290 simplex.
step[2*i + 2] = model_type();
291 simplex.
step[2*i + 2][2*i + 2] = 0.1;
294 DEBUG(
"start values " << simplex.
value << endl);
296 for (
size_t i = 0; i != simplex.
step.size(); ++i) {
297 DEBUG(
"step: " << setw(2) << i <<
' ' << simplex.
step[i] << endl);
303 const double chi2 = simplex(fit,
data.begin(),
data.end());
304 const int ndf =
data.size() - simplex.
step.size();
309 for (
size_t i = 0; i !=
model.size(); ++i) {
313 DEBUG(
"chi2 / NDF" <<
FIXED(12,3) << chi2 <<
" / " << ndf << endl);
315 DEBUG(
"final values " << simplex.
value << endl);
320 for (
const auto& hit :
data) {
324 for (
int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
325 cout <<
FIXED(7,3) << hit[i] <<
' ';
328 const double value = simplex.
value(hit);
330 cout <<
FIXED(7,3) << hit.value <<
' '
331 <<
FIXED(7,3) << value <<
' '
332 <<
FIXED(7,3) << (hit.value - value) / hit.sigma << endl;
346 for (
size_t i = 0; i !=
sizeof(H1)/
sizeof(H1[0]); ++i) {
General purpose messaging.
#define DEBUG(A)
Message macros.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
I/O formatting auxiliaries.
#define MAKE_CSTRING(A)
Make C-string.
Double_t g1(const Double_t x)
Function.
int main(int argc, char **argv)
Simple fit method based on Powell's algorithm, see reference: Numerical Recipes in C++,...
std::vector< JModel_t > step
Utility class to parse command line options.
double gauss(const double x, const double sigma)
Gauss function (normalised to 1 at x = 0).
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Auxiliary data structure for floating point format specification.
Auxiliary data structure for return type of make methods.