Jpp  16.0.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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"
15 #include "JTools/JArrayIterator.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 
24 namespace {
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  */
160 int main(int argc, char **argv)
161 {
162  using namespace std;
163  using namespace JPP;
164 
165  typedef JArray<2, double> array_type;
166 
167  string outputFile;
168  int numberOfEvents;
169  array_type parameters;
170  UInt_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 
211  JSimplex<model_type>::MAXIMUM_ITERATIONS = 10000;
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 
236  for (JArrayIterator<NUMBER_OF_DIMENSIONS, double> g1(grid); g1; ++g1) {
237 
238  const JArray<NUMBER_OF_DIMENSIONS, double> p1 = *g1;
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 }
Utility class to parse command line options.
Definition: JParser.hh:1500
Q(UTCMax_s-UTCMin_s)-livetime_s
debug
Definition: JMessage.hh:29
int main(int argc, char *argv[])
Definition: Main.cc:15
TPaveText * p1
#define STATUS(A)
Definition: JMessage.hh:63
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:151
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
string outputFile
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable STRING $argv[2] set_array QUANTILES set_variable FORMULA *[0] exp(-0.5 *(x-[1])*(x-[1])/([2]*[2]))" set_variable MODULE `getModule -a $DETECTOR -L "$STRING 0"` source JAcoustics.sh -- typeset -A TRIPODS get_tripods $WORKDIR/tripod.txt TRIPODS XMEAN
I/O formatting auxiliaries.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
int debug
debug level
Definition: JSirene.cc:63
One dimensional array of template objects with fixed length.
Definition: JArray.hh:40
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
Utility class to parse command line options.
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
double gauss(const double x, const double sigma)
Gauss function (normalised to 1 at x = 0).
double u[N+1]
Definition: JPolint.hh:755
Double_t g1(const Double_t x)
Function.
Definition: JQuantiles.cc:25