Jpp test-rotations-old
the software that should make you happy
Loading...
Searching...
No Matches
JSimplexFitToGauss.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 "TF1.h"
10#include "TH1D.h"
11#include "TMath.h"
12
13#include "JFit/JSimplex.hh"
14#include "JTools/JElement.hh"
15#include "JTools/JQuantile.hh"
16#include "JROOT/JRootToolkit.hh"
17#include "JMath/JGauss.hh"
18#include "JMath/JConstants.hh"
19
20#include "Jeep/JTimer.hh"
21#include "Jeep/JPrint.hh"
22#include "Jeep/JParser.hh"
23#include "Jeep/JMessage.hh"
24
25namespace {
26
27 using namespace JPP;
28 using TMath::PoissonI;
29
30 typedef JElement2D<double, double> JElement_t;
31
32 /**
33 * Fit function.
34 *
35 * \param gauss gauss
36 * \param point point
37 * \return chi2
38 */
39 inline double g1(const JGauss& gauss, const JElement_t& point)
40 {
41 const double u = (point.getX() - gauss.mean) / gauss.sigma;
42 const double fs = gauss.signal * exp(-0.5*u*u) / (sqrt(2.0*PI) * gauss.sigma);
43 const double fb = gauss.background;
44 const double f1 = fs + fb;
45
46 return -log(PoissonI(point.getY(), f1));
47 }
48}
49
50
51/**
52 * \file
53 *
54 * Program to test JFIT::JSimplex algorithm.
55 * \author mdejong
56 */
57int main(int argc, char **argv)
58{
59 using namespace std;
60 using namespace JPP;
61
62 string outputFile;
63 int numberOfEvents;
65 JGauss precision;
66 int debug;
67
68 try {
69
70 JParser<> zap("Program to test JSimplex algorithm.");
71
72 zap['o'] = make_field(outputFile) = "";
73 zap['n'] = make_field(numberOfEvents) = 100;
74 zap['@'] = make_field(gauss) = JGauss(0.0, 1.0, 1000.0, 100.0);
75 zap['e'] = make_field(precision) = JGauss(0.05, 0.05, 25.0, 25.0);
76 zap['d'] = make_field(debug) = 3;
77
78 zap(argc, argv);
79 }
80 catch(const exception& error) {
81 FATAL(error.what() << endl);
82 }
83
84 ASSERT(numberOfEvents > 0);
85
87
88 TF1 fs("fs", "exp(-0.5 * (x-[0])*(x-[0]) / ([1]*[1]))");
89 TF1 fb("fb", "1.0");
90
91 fs.FixParameter(0, gauss.mean);
92 fs.FixParameter(1, gauss.sigma);
93
94 const Int_t nx = 21;
95 const Double_t xmin = -5.0;
96 const Double_t xmax = +5.0;
97
98 JQuantile Q[] = { JQuantile("mean "),
99 JQuantile("sigma "),
100 JQuantile("signal "),
101 JQuantile("background") };
102
103 TH1D H[] = { TH1D("ha", "", 101, -0.1, +0.1),
104 TH1D("hb", "", 101, -0.1, +0.1),
105 TH1D("hc", "", 101, -100.0, +100.0),
106 TH1D("hd", "", 101, -100.0, +100.0) };
107
108 JTimer timer;
109
110 for (int i = 0; i != numberOfEvents; ++i) {
111
112 STATUS("event: " << setw(10) << i << '\r'); DEBUG(endl);
113
114 TH1D h0("h0", NULL, nx, xmin, xmax);
115
116 h0.FillRandom("fs", (Int_t) gauss.signal);
117 h0.FillRandom("fb", (Int_t) gauss.background);
118
120
121 for (Int_t i = 1; i <= h0.GetNbinsX(); ++i) {
122 data.push_back(JElement_t(h0.GetBinCenter (i),
123 h0.GetBinContent(i)));
124 }
125
127
128 fit.step.resize(4);
129
130 fit.step[0] = JGauss(0.1, 0.0, 0.0, 0.0);
131 fit.step[1] = JGauss(0.0, 0.1, 0.0, 0.0);
132 fit.step[2] = JGauss(0.0, 0.0, gauss.signal * 0.01, 0.0);
133 fit.step[3] = JGauss(0.0, 0.0, 0.0, gauss.background * 0.01);
134
135 fit.value = JGauss(h0.GetMean(),
136 h0.GetRMS(),
137 h0.GetEntries() * (xmax - xmin) / nx - h0.GetMinimum(),
138 h0.GetMinimum());
139
140 DEBUG("Start value " << fit.value << endl);
141
142 for (size_t i = 0; i != fit.step.size(); ++i) {
143 DEBUG("Step size " << i << ' ' << fit.step[i] << endl);
144 }
145
146 timer.start();
147
148 const double chi2 = fit(g1, data.begin(), data.end());
149
150 timer.stop();
151
152 DEBUG("Final value " << fit.value << endl);
153 DEBUG("Chi2 " << chi2 << endl);
154
155 const double Y[] = { fit.value.mean - gauss.mean,
156 fit.value.sigma - gauss.sigma,
157 fit.value.signal * nx / (xmax - xmin) - gauss.signal,
158 fit.value.background * nx - gauss.background };
159
160 for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
161 Q[i].put (Y[i]);
162 H[i].Fill(Y[i]);
163 }
164 }
165
166 for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
167 NOTICE((i == 0 ? longprint : shortprint) << Q[i]);
168 }
169
170 if (debug >= notice_t) {
171 timer.print(cout, true, micro_t);
172 }
173
174 if (outputFile != "") {
175
176 TFile out(outputFile.c_str(), "recreate");
177
178 for (int i = 0; i != sizeof(H)/sizeof(H[0]); ++i) {
179 out << H[i];
180 }
181
182 out.Write();
183 out.Close();
184 }
185
186 for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
187 ASSERT(fabs(Q[i].getMean()) < Q[i].getSTDev());
188 }
189
190 ASSERT(Q[0].getSTDev() < precision.mean);
191 ASSERT(Q[1].getSTDev() < precision.sigma);
192 ASSERT(Q[2].getSTDev() < precision.signal);
193 ASSERT(Q[3].getSTDev() < precision.background);
194
195 return 0;
196}
double getMean(vector< double > &v)
get mean of vector content
string outputFile
The elements in a collection are sorted according to their abscissa values and a given distance opera...
std::ostream & longprint(std::ostream &out)
Set long printing.
Definition JManip.hh:172
std::ostream & shortprint(std::ostream &out)
Set short printing.
Definition JManip.hh:144
Mathematical constants.
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define STATUS(A)
Definition JMessage.hh:63
#define ASSERT(A,...)
Assert macro.
Definition JMessage.hh:90
#define NOTICE(A)
Definition JMessage.hh:64
#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.
Double_t g1(const Double_t x)
Function.
Definition JQuantiles.cc:25
int main(int argc, char **argv)
Auxiliary class for CPU timing and usage.
Definition JTimer.hh:33
void print(std::ostream &out, const JScale_t scale=milli_t) const
Print timer data.
Definition JTimer.hh:172
void stop()
Stop timer.
Definition JTimer.hh:127
void start()
Start timer.
Definition JTimer.hh:106
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
const JPolynome f1(1.0, 2.0, 3.0)
Function.
@ micro_t
micro
Definition JScale.hh:29
@ notice_t
notice
Definition JMessage.hh:32
double gauss(const double x, const double sigma)
Gauss function (normalised to 1 at x = 0).
static const double H
Planck constant [eV s].
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
double background
Definition JGauss.hh:164
double signal
Definition JGauss.hh:163
Gauss function object.
Definition JMathlib.hh:1589
double sigma
sigma
Definition JMathlib.hh:1665
2D Element.
Definition JPolint.hh:1131
Auxiliary data structure for running average, standard deviation and quantiles.
Definition JQuantile.hh:46
void put(const double x, const double w=1.0)
Put value.
Definition JQuantile.hh:133