73{
76
78 int numberOfEvents;
82
83 try {
84
85 JParser<> zap(
"Program to test JGandalf algorithm.");
86
92
93 zap(argc, argv);
94 }
95 catch(const exception& error) {
96 FATAL(error.what() << endl);
97 }
98
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
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
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
164
165 const double chi2 = fit(
g1,
data.begin(),
data.end());
166
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,
175 fit.value.background * nx -
gauss.background };
176
177 for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
180 }
181 }
182
183 for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
185 }
186
189 }
190
192
194
195 for (
int i = 0; i !=
sizeof(
H)/
sizeof(
H[0]); ++i) {
197 }
198
199 out.Write();
200 out.Close();
201 }
202
203 for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
205 }
206
211
212 return 0;
213}
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).