Jpp 20.0.0-rc.2
the software that should make you happy
Loading...
Searching...
No Matches
JGandalfFitAlaMinuitToGauss.cc File Reference

Program to test JFIT::JGandalf algorithm. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>
#include "TROOT.h"
#include "TFile.h"
#include "TF1.h"
#include "TH1D.h"
#include "TMath.h"
#include "JMath/JModel.hh"
#include "JMath/JConstants.hh"
#include "JFit/JGandalf.hh"
#include "JTools/JElement.hh"
#include "JTools/JStats.hh"
#include "JROOT/JRootToolkit.hh"
#include "Jeep/JTimer.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Program to test JFIT::JGandalf algorithm.

Author
mdejong

Definition in file JGandalfFitAlaMinuitToGauss.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 109 of file JGandalfFitAlaMinuitToGauss.cc.

110{
111 using namespace std;
112 using namespace JPP;
113
114 string outputFile;
115 int numberOfEvents;
117 JGauss precision;
118 int debug;
119
120 try {
121
122 JParser<> zap("Program to test JGandalf algorithm.");
123
124 zap['o'] = make_field(outputFile) = "";
125 zap['n'] = make_field(numberOfEvents) = 100;
126 zap['@'] = make_field(gauss) = JGauss(0.0, 1.0, 1000.0, 100.0);
127 zap['e'] = make_field(precision) = JGauss(0.05, 0.05, 25.0, 25.0);
128 zap['d'] = make_field(debug) = 3;
129
130 zap(argc, argv);
131 }
132 catch(const exception& error) {
133 FATAL(error.what() << endl);
134 }
135
136 using namespace JFIT;
137
138 ASSERT(numberOfEvents > 0);
139
141
142 TF1 fs("fs", "exp(-0.5 * (x-[0])*(x-[0]) / ([1]*[1]))");
143 TF1 fb("fb", "1.0");
144
145 fs.FixParameter(0, gauss[0]);
146 fs.FixParameter(1, gauss[1]);
147
148 const Int_t nx = 21;
149 const Double_t xmin = -5.0;
150 const Double_t xmax = +5.0;
151
152 JStats Q[] = {
153 JStats("mean "),
154 JStats("sigma "),
155 JStats("signal "),
156 JStats("background")
157 };
158
159 TH1D H[] = {
160 TH1D("ha", "", 101, -0.1, +0.1),
161 TH1D("hb", "", 101, -0.1, +0.1),
162 TH1D("hc", "", 101, -100.0, +100.0),
163 TH1D("hd", "", 101, -100.0, +100.0)
164 };
165
166 JTimer timer;
167
168 for (int i = 0; i != numberOfEvents; ++i) {
169
170 STATUS("event: " << setw(10) << i << '\r'); DEBUG(endl);
171
172 TH1D h0("h0", NULL, nx, xmin, xmax);
173
174 h0.Sumw2();
175
176 h0.FillRandom("fs", (Int_t) gauss[2]);
177 h0.FillRandom("fb", (Int_t) gauss[3]);
178
180
181 for (Int_t i = 1; i <= h0.GetNbinsX(); ++i) {
182 data.push_back(JElement_t(h0.GetBinCenter (i),
183 h0.GetBinContent(i)));
184 }
185
187
188 fit.parameters.resize(4);
189
190 for (size_t i = 0; i != 4; ++i) {
191 fit.parameters[i] = i;
192 }
193
194 fit.value = JGauss(h0.GetMean(),
195 h0.GetRMS(),
196 h0.GetEntries() * (xmax - xmin) / nx - h0.GetMinimum(),
197 h0.GetMinimum());
198
199 DEBUG("Start value " << fit.value << endl);
200
201 timer.start();
202
203 const double chi2 = fit(g1, data.begin(), data.end());
204
205 timer.stop();
206
207 DEBUG("Final value " << fit.value << endl);
208 DEBUG("Chi2 " << chi2 << endl);
209
210 const double Y[] = { fit.value[0] - gauss[0],
211 fit.value[1] - gauss[1],
212 fit.value[2] * nx / (xmax - xmin) - gauss[2],
213 fit.value[3] * nx - gauss[3] };
214
215 for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
216 Q[i].put (Y[i]);
217 H[i].Fill(Y[i]);
218 }
219 }
220
221 for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
222 NOTICE((i == 0 ? longprint : shortprint) << Q[i]);
223 }
224
225 if (debug >= notice_t) {
226 timer.print(cout, true, micro_t);
227 }
228
229 if (outputFile != "") {
230
231 TFile out(outputFile.c_str(), "recreate");
232
233 for (int i = 0; i != sizeof(H)/sizeof(H[0]); ++i) {
234 out << H[i];
235 }
236
237 out.Write();
238 out.Close();
239 }
240
241 for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
242 ASSERT(fabs(Q[i].getMean()) < Q[i].getSTDev());
243 }
244
245 ASSERT(Q[0].getSTDev() < precision[0]);
246 ASSERT(Q[1].getSTDev() < precision[1]);
247 ASSERT(Q[2].getSTDev() < precision[2]);
248 ASSERT(Q[3].getSTDev() < precision[3]);
249
250 return 0;
251}
double getMean(vector< double > &v)
get mean of vector content
string outputFile
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
#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
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
Double_t g1(const Double_t x)
Function.
Definition JQuantiles.cc:25
Place holder for custom implementation.
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
Utility class to parse command line options.
Definition JParser.hh:1698
@ micro_t
micro
Definition JScale.hh:29
@ notice_t
notice
Definition JMessage.hh:32
Auxiliary classes and methods for linear and iterative data regression.
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).
Gauss function object.
Definition JMathlib.hh:1589
Auxiliary data structure for running average and standard deviation.
Definition JStats.hh:44
void put(const double x, const double w=1.0)
Put value.
Definition JStats.hh:119