Jpp test-rotations-old
the software that should make you happy
Loading...
Searching...
No Matches
JSimplexFitArray.cc
Go to the documentation of this file.
1#include <string>
2#include <iostream>
3#include <iomanip>
4#include <vector>
5#include <cmath>
6
7#include "TROOT.h"
8#include "TFile.h"
9#include "TH1D.h"
10#include "TRandom3.h"
11#include "TMath.h"
12
13#include "JFit/JSimplex.hh"
14#include "JTools/JArray.hh"
16#include "JTools/JGrid.hh"
17#include "JTools/JQuantile.hh"
18#include "JROOT/JRootToolkit.hh"
19
20#include "Jeep/JPrint.hh"
21#include "Jeep/JParser.hh"
22#include "Jeep/JMessage.hh"
23
24namespace {
25
26 using JTOOLS::JArray;
27 using TMath::PoissonI;
28
29 static const int NUMBER_OF_DIMENSIONS = 3;
30
31
32 /**
33 * Type definition of a point in ND-space.
34 */
35 typedef JArray<NUMBER_OF_DIMENSIONS, double> point_type;
36
37
38 /**
39 * Model.
40 */
41 struct model_type :
42 public JArray<1 + 2 * NUMBER_OF_DIMENSIONS, double>
43 {
44 /**
45 * Default constructor.
46 */
47 model_type() :
48 JArray<1 + 2 * NUMBER_OF_DIMENSIONS, double>()
49 {}
50
51
52 /**
53 * Constructor.
54 *
55 * \param value first value
56 * \param args remaining values
57 */
58 template<class ...Args>
59 model_type(argument_type value, const Args& ...args) :
60 JArray<1 + 2 * NUMBER_OF_DIMENSIONS, double>(value, args...)
61 {}
62
63
64 /**
65 * Model function.
66 *
67 * \param point point
68 * \return value
69 */
70 inline double operator()(const point_type& point) const
71 {
72 double value = (*this)[0];
73
74 for (int i = 0; i != NUMBER_OF_DIMENSIONS;++i) {
75 value *= gauss(point[i] - (*this)[2*i + 1], (*this)[2*i + 2]);
76 }
77
78 return value;
79 }
80
81
82 /**
83 * Gauss function (normalised to 1 at x = 0).
84 *
85 * \param x x
86 * \param sigma sigma
87 * \return function value
88 */
89 static inline double gauss(const double x, const double sigma)
90 {
91 const double u = x / sigma;
92
93 if (fabs(u) < 10.0)
94 return exp(-0.5*u*u);
95 else
96 return 0.0;
97 }
98 };
99
100
101 /**
102 * Hit.
103 */
104 struct hit_type :
105 public point_type
106 {
107 /**
108 * Constructor.
109 *
110 * \param point point
111 * \param value value
112 * \param sigma sigma
113 */
114 hit_type(const point_type& point, const double value, const double sigma) :
115 point_type(point),
116 value(value),
117 sigma(sigma)
118 {}
119
120 double value;
121 double sigma;
122 };
123
124
125 /**
126 * Fit function based on classical chi2.
127 *
128 * \param model model
129 * \param hit hit
130 * \return chi2
131 */
132 inline double chi2(const model_type& model, const hit_type& hit)
133 {
134 const double u = (model(hit) - hit.value) / hit.sigma;
135
136 return u*u;
137 }
138
139
140 /**
141 * Fit function based on logarithm of likelihood.
142 *
143 * \param model model
144 * \param hit hit
145 * \return chi2
146 */
147 inline double likelihood(const model_type& model, const hit_type& hit)
148 {
149 return -log(PoissonI(hit.value, model(hit)));
150 }
151}
152
153
154/**
155 * \file
156 *
157 * Program to test JFIT::JSimplex algorithm.
158 * \author mdejong
159 */
160int main(int argc, char **argv)
161{
162 using namespace std;
163 using namespace JPP;
164
166
167 string outputFile;
168 int numberOfEvents;
169 array_type parameters;
170 ULong_t seed;
171 int debug;
172
173 try {
174
175 JParser<> zap("Program to test JSimplex algorithm.");
176
177 zap['o'] = make_field(outputFile) = "simplex.root";
178 zap['n'] = make_field(numberOfEvents);
179 zap['M'] = make_field(parameters, "model parameters") = array_type(100.0, 1.0);
180 zap['S'] = make_field(seed) = 0;
181 zap['d'] = make_field(debug) = 2;
182
183 zap(argc, argv);
184 }
185 catch(const exception& error) {
186 FATAL(error.what() << endl);
187 }
188
189
190 gRandom->SetSeed(seed);
191
192
193 const double top = parameters[0];
194 const double sigma = parameters[1];
195
196 model_type model;
197
198 model[0] = top;
199
200 for (int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
201 model[i*2 + 1] = gRandom->Uniform(-sigma, +sigma);
202 model[i*2 + 2] = gRandom->Gaus(sigma, 0.2*sigma);
203 }
204
205 const JGrid<double> grid(11, -3.0*sigma, +3.0*sigma);
206
207
208 JSimplex<model_type> simplex;
209
212
213 double (*fit)(const model_type&, const hit_type&) = likelihood; // fit function
214
215
216 vector<TH1D*> H1;
217
218 TH1D h0("chi2/NDF", NULL, 200, 0.0, 10.0);
219
220 H1.push_back(new TH1D(MAKE_CSTRING("p" << 0), NULL, 101, -20.0, +20.0));
221
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));
224 }
225
226
227 for (int counter = 0; counter != numberOfEvents; ++counter) {
228
229 STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);
230
231
232 // generate data
233
234 vector<hit_type> data;
235
237
239
240 const double value = gRandom->Poisson(model(p1));
241 const double sigma = (value > 1.0 ? sqrt(value) : 0.7);
242
243 data.push_back(hit_type(p1, value, sigma));
244 }
245
246
247 // start values
248
249 JQuantile Q[NUMBER_OF_DIMENSIONS];
250
251 double value = 0.0;
252
253 for (const auto& hit : data) {
254
255 for (int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
256 Q[i].put(hit[i], hit.value);
257 }
258
259 if (hit.value > value) {
260
261 simplex.value[0] = hit.value;
262
263 for (int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
264 simplex.value[2*i + 1] = hit[i];
265 }
266
267 value = hit.value;
268 }
269 }
270
271 for (int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
272 simplex.value[2*i + 2] = 1.0 * Q[i].getSTDev();
273 }
274
275
276 // steps
277
278 simplex.step.resize(2*NUMBER_OF_DIMENSIONS + 1);
279
280 model_type step;
281
282 simplex.step[0] = model_type();
283 simplex.step[0][0] = 1.0;
284
285 for (int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
286
287 simplex.step[2*i + 1] = model_type();
288 simplex.step[2*i + 1][2*i + 1] = 0.1;
289
290 simplex.step[2*i + 2] = model_type();
291 simplex.step[2*i + 2][2*i + 2] = 0.1;
292 }
293
294 DEBUG("start values " << simplex.value << endl);
295
296 for (size_t i = 0; i != simplex.step.size(); ++i) {
297 DEBUG("step: " << setw(2) << i << ' ' << simplex.step[i] << endl);
298 }
299
300
301 // fit
302
303 const double chi2 = simplex(fit, data.begin(), data.end());
304 const int ndf = data.size() - simplex.step.size();
305
306
307 h0.Fill(chi2 / ndf);
308
309 for (size_t i = 0; i != model.size(); ++i) {
310 H1[i]->Fill(model[i] - simplex.value[i]);
311 }
312
313 DEBUG("chi2 / NDF" << FIXED(12,3) << chi2 << " / " << ndf << endl);
314
315 DEBUG("final values " << simplex.value << endl);
316 DEBUG("model values " << model << endl);
317
318 if (debug >= debug_t) {
319
320 for (const auto& hit : data) {
321
322 cout << "hit: ";
323
324 for (int i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
325 cout << FIXED(7,3) << hit[i] << ' ';
326 }
327
328 const double value = simplex.value(hit);
329
330 cout << FIXED(7,3) << hit.value << ' '
331 << FIXED(7,3) << value << ' '
332 << FIXED(7,3) << (hit.value - value) / hit.sigma << endl;
333 }
334 }
335
336 }
337 STATUS(endl);
338
339
340 if (outputFile != "") {
341
342 TFile out(outputFile.c_str(), "recreate");
343
344 out << h0;
345
346 for (size_t i = 0; i != sizeof(H1)/sizeof(H1[0]); ++i) {
347 out << *(H1[i]);
348 }
349
350 out.Write();
351 out.Close();
352 }
353}
string outputFile
TPaveText * p1
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define STATUS(A)
Definition JMessage.hh:63
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
I/O formatting auxiliaries.
#define MAKE_CSTRING(A)
Make C-string.
Definition JPrint.hh:72
Double_t g1(const Double_t x)
Function.
Definition JQuantiles.cc:25
int main(int argc, char **argv)
Simple fit method based on Powell's algorithm, see reference: Numerical Recipes in C++,...
Definition JSimplex.hh:44
JModel_t value
Definition JSimplex.hh:240
std::vector< JModel_t > step
Definition JSimplex.hh:241
Utility class to parse command line options.
Definition JParser.hh:1698
One dimensional array of template objects with fixed length.
Definition JArray.hh:43
const double sigma[]
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.
Definition JManip.hh:448
Auxiliary data structure for return type of make methods.
Definition JVectorize.hh:28
ND array iterator.
Simple data structure for an abstract collection of equidistant abscissa values.
Definition JGrid.hh:40
Auxiliary data structure for running average, standard deviation and quantiles.
Definition JQuantile.hh:46
double getSTDev() const
Get standard deviation.
Definition JQuantile.hh:281
void put(const double x, const double w=1.0)
Put value.
Definition JQuantile.hh:133