122{
125
126 int numberOfEvents;
129 ULong_t seed;
131
132 try {
133
134 JParser<> zap(
"Program to test JGradient algorithm.");
135
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
167
168
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
178
179 } else {
180
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
204
205
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
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) {
227 }
228
229 for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
231 }
232
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.
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
Utility class to parse command line options.
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).
Auxiliary data structure for editable parameter.