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

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

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>
#include "TRandom3.h"
#include "JFit/JGradient.hh"
#include "JTools/JElement.hh"
#include "JTools/JQuantile.hh"
#include "JMath/JGauss.hh"
#include "JMath/JConstants.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::JGradient algorithm.

Author
mdejong

Definition in file JGradientFitToGauss.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 121 of file JGradientFitToGauss.cc.

122{
123 using namespace std;
124 using namespace JPP;
125
126 int numberOfEvents;
128 JGauss precision;
129 ULong_t seed;
130 int debug;
131
132 try {
133
134 JParser<> zap("Program to test JGradient algorithm.");
135
136 zap['n'] = make_field(numberOfEvents) = 0;
137 zap['@'] = make_field(gauss) = JGauss(0.0, 1.0, 1000.0, 10.0);
138 zap['e'] = make_field(precision) = JGauss(0.05, 0.05, 25.0, 25.0);
139 zap['S'] = make_field(seed) = 0;
140 zap['d'] = make_field(debug) = 3;
141
142 zap(argc, argv);
143 }
144 catch(const exception& error) {
145 FATAL(error.what() << endl);
146 }
147
148 gRandom->SetSeed(seed);
149
150 if (numberOfEvents == 0) {
151
152 buffer.push_back(element_type(-3.00000, 16.00000));
153 buffer.push_back(element_type(-2.50000, 22.00000));
154 buffer.push_back(element_type(-2.00000, 70.00000));
155 buffer.push_back(element_type(-1.50000, 152.00000));
156 buffer.push_back(element_type(-1.00000, 261.00000));
157 buffer.push_back(element_type(-0.50000, 337.00000));
158 buffer.push_back(element_type( 0.00000, 378.00000));
159 buffer.push_back(element_type( 0.50000, 360.00000));
160 buffer.push_back(element_type( 1.00000, 253.00000));
161 buffer.push_back(element_type( 1.50000, 129.00000));
162 buffer.push_back(element_type( 2.00000, 72.00000));
163 buffer.push_back(element_type( 2.50000, 22.00000));
164 buffer.push_back(element_type( 3.00000, 11.00000));
165
166 JGradient gradient(1000000, 0, 1.0e-4, debug);
167
168 // start values
169
170 fit = JGauss(0.5, 0.5, 700.0, 0.0);
171
172 gradient.push_back(JModifier_t("mean", new JGaussEditor(fit, JGauss(1.0, 0.0, 0.0, 0.0)), 1.0e-2));
173 gradient.push_back(JModifier_t("sigma", new JGaussEditor(fit, JGauss(0.0, 1.0, 0.0, 0.0)), 1.0e-2));
174 gradient.push_back(JModifier_t("signal", new JGaussEditor(fit, JGauss(0.0, 0.0, 1.0, 0.0)), 5.0e-0));
175 gradient.push_back(JModifier_t("background", new JGaussEditor(fit, JGauss(0.0, 0.0, 0.0, 1.0)), 5.0e-1));
176
177 gradient(getChi2);
178
179 } else {
180
181 JQuantile Q[] = { JQuantile("mean "),
182 JQuantile("sigma "),
183 JQuantile("signal "),
184 JQuantile("background") };
185
186 const size_t nx = 21;
187 const double xmin = -5.0;
188 const double xmax = +5.0;
189
190 for (int i = 0; i != numberOfEvents; ++i) {
191
192 STATUS("event: " << setw(10) << i << '\r'); DEBUG(endl);
193
194 buffer.clear();
195
196 for (double x = xmin, dx = (xmax - xmin) / (nx - 1); x < xmax + 0.5*dx; x += dx) {
197
198 const double value = gauss(x);
199
200 buffer.push_back(element_type(x, gRandom->Poisson(value)));
201 }
202
203 JGradient gradient(1000000, 0, 1.0e-4, debug);
204
205 // start values
206
207 fit = JGauss(0.5, 0.5, 700.0, 0.0);
208
209 gradient.push_back(JModifier_t("mean", new JGaussEditor(fit, JGauss(1.0, 0.0, 0.0, 0.0)), 1.0e-2));
210 gradient.push_back(JModifier_t("sigma", new JGaussEditor(fit, JGauss(0.0, 1.0, 0.0, 0.0)), 1.0e-2));
211 gradient.push_back(JModifier_t("signal", new JGaussEditor(fit, JGauss(0.0, 0.0, 1.0, 0.0)), 5.0e-0));
212 gradient.push_back(JModifier_t("background", new JGaussEditor(fit, JGauss(0.0, 0.0, 0.0, 1.0)), 5.0e-1));
213
214 const double chi2 = gradient(getChi2);
215
216 DEBUG("Final value " << fit << endl);
217 DEBUG("Chi2 " << chi2 << endl);
218
219 Q[0].put(fit.mean - gauss.mean);
220 Q[1].put(fit.sigma - gauss.sigma);
221 Q[2].put(fit.signal - gauss.signal);
222 Q[3].put(fit.background - gauss.background);
223 }
224
225 for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
226 NOTICE((i == 0 ? longprint : shortprint) << Q[i]);
227 }
228
229 for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
230 ASSERT(fabs(Q[i].getMean()) < Q[i].getSTDev());
231 }
232
233 ASSERT(Q[0].getSTDev() < precision.mean);
234 ASSERT(Q[1].getSTDev() < precision.sigma);
235 ASSERT(Q[2].getSTDev() < precision.signal);
236 ASSERT(Q[3].getSTDev() < precision.background);
237 }
238
239 return 0;
240}
double getMean(vector< double > &v)
get mean of vector content
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
Utility class to parse command line options.
Definition JParser.hh:1698
const double xmax
const double xmin
double getChi2(const double P)
Get chi2 corresponding to given probability.
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).
Conjugate gradient fit.
Definition JGradient.hh:75
Auxiliary data structure for editable parameter.
Definition JGradient.hh:49
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