110{
113
115 int numberOfEvents;
119
120 try {
121
122 JParser<> zap(
"Program to test JGandalf algorithm.");
123
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
157 };
158
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
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
202
203 const double chi2 = fit(
g1,
data.begin(),
data.end());
204
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) {
218 }
219 }
220
221 for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
223 }
224
227 }
228
230
232
233 for (
int i = 0; i !=
sizeof(
H)/
sizeof(
H[0]); ++i) {
235 }
236
237 out.Write();
238 out.Close();
239 }
240
241 for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
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
std::ostream & longprint(std::ostream &out)
Set long printing.
std::ostream & shortprint(std::ostream &out)
Set short printing.
#define DEBUG(A)
Message macros.
#define ASSERT(A,...)
Assert macro.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Double_t g1(const Double_t x)
Function.
Place holder for custom implementation.
Auxiliary class for CPU timing and usage.
void print(std::ostream &out, const JScale_t scale=milli_t) const
Print timer data.
Utility class to parse command line options.
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).