Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
JGandalfFitToGauss.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 "JFit/JGandalf.hh"
#include "JTools/JElement.hh"
#include "JTools/JQuantile.hh"
#include "JROOT/JRootToolkit.hh"
#include "JMath/JGauss.hh"
#include "JMath/JConstants.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 JGandalfFitToGauss.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 72 of file JGandalfFitToGauss.cc.

73{
74 using namespace std;
75 using namespace JPP;
76
77 string outputFile;
78 int numberOfEvents;
80 JGauss precision;
81 int debug;
82
83 try {
84
85 JParser<> zap("Program to test JGandalf algorithm.");
86
87 zap['o'] = make_field(outputFile) = "";
88 zap['n'] = make_field(numberOfEvents) = 100;
89 zap['@'] = make_field(gauss) = JGauss(0.0, 1.0, 1000.0, 100.0);
90 zap['e'] = make_field(precision) = JGauss(0.05, 0.05, 25.0, 25.0);
91 zap['d'] = make_field(debug) = 3;
92
93 zap(argc, argv);
94 }
95 catch(const exception& error) {
96 FATAL(error.what() << endl);
97 }
98
99 using namespace JFIT;
100
101 ASSERT(numberOfEvents > 0);
102
104
105 TF1 fs("fs", "exp(-0.5 * (x-[0])*(x-[0]) / ([1]*[1]))");
106 TF1 fb("fb", "1.0");
107
108 fs.FixParameter(0, gauss.mean);
109 fs.FixParameter(1, gauss.sigma);
110
111 const Int_t nx = 21;
112 const Double_t xmin = -5.0;
113 const Double_t xmax = +5.0;
114
115 JQuantile Q[] = { JQuantile("mean "),
116 JQuantile("sigma "),
117 JQuantile("signal "),
118 JQuantile("background") };
119
120 TH1D H[] = { TH1D("ha", "", 101, -0.1, +0.1),
121 TH1D("hb", "", 101, -0.1, +0.1),
122 TH1D("hc", "", 101, -100.0, +100.0),
123 TH1D("hd", "", 101, -100.0, +100.0) };
124
125 JTimer timer;
126
127 for (int i = 0; i != numberOfEvents; ++i) {
128
129 STATUS("event: " << setw(10) << i << '\r'); DEBUG(endl);
130
131 TH1D h0("h0", NULL, nx, xmin, xmax);
132
133 h0.Sumw2();
134
135 h0.FillRandom("fs", (Int_t) gauss.signal);
136 h0.FillRandom("fb", (Int_t) gauss.background);
137
139
140 for (Int_t i = 1; i <= h0.GetNbinsX(); ++i) {
141 data.push_back(JElement_t(h0.GetBinCenter (i),
142 h0.GetBinContent(i)));
143 }
144
146
148
149 fit.parameters.resize(4);
150
151 fit.parameters[0] = &JGauss::mean;
152 fit.parameters[1] = &JGauss::sigma;
153 fit.parameters[2] = &JGauss::signal;
154 fit.parameters[3] = &JGauss::background;
155
156 fit.value = JGauss(h0.GetMean(),
157 h0.GetRMS(),
158 h0.GetEntries() * (xmax - xmin) / nx - h0.GetMinimum(),
159 h0.GetMinimum());
160
161 DEBUG("Start value " << fit.value << endl);
162
163 timer.start();
164
165 const double chi2 = fit(g1, data.begin(), data.end());
166
167 timer.stop();
168
169 DEBUG("Final value " << fit.value << endl);
170 DEBUG("Chi2 " << chi2 << endl);
171
172 const double Y[] = { fit.value.mean - gauss.mean,
173 fit.value.sigma - gauss.sigma,
174 fit.value.signal * nx / (xmax - xmin) - gauss.signal,
175 fit.value.background * nx - gauss.background };
176
177 for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
178 Q[i].put (Y[i]);
179 H[i].Fill(Y[i]);
180 }
181 }
182
183 for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
184 NOTICE((i == 0 ? longprint : shortprint) << Q[i]);
185 }
186
187 if (debug >= notice_t) {
188 timer.print(cout, true, micro_t);
189 }
190
191 if (outputFile != "") {
192
193 TFile out(outputFile.c_str(), "recreate");
194
195 for (int i = 0; i != sizeof(H)/sizeof(H[0]); ++i) {
196 out << H[i];
197 }
198
199 out.Write();
200 out.Close();
201 }
202
203 for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
204 ASSERT(fabs(Q[i].getMean()) < Q[i].getSTDev());
205 }
206
207 ASSERT(Q[0].getSTDev() < precision.mean);
208 ASSERT(Q[1].getSTDev() < precision.sigma);
209 ASSERT(Q[2].getSTDev() < precision.signal);
210 ASSERT(Q[3].getSTDev() < precision.background);
211
212 return 0;
213}
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
const double xmax
const double xmin
@ micro_t
micro
Definition JScale.hh:29
@ notice_t
notice
Definition JMessage.hh:32
Auxiliary classes and methods for linear and iterative data regression.
Definition JEnergy.hh:15
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
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