31{
34
35 const char* const count_t = "count";
36 const char* const value_t = "value";
37
38 struct {
39 double signal = 100;
40 double background = 5;
41 double center = 0.0;
43 } parameters;
44
45 string inputFile;
50 size_t N;
51 string option;
52 bool writeFits;
53 ULong_t seed;
55
56 try {
57
59
64
65 JParser<> zap(
"Program to test JRootfit algorithm.");
66
74 zap[
'O'] =
make_field(option) = count_t, value_t;
78
79 zap(argc, argv);
80 }
81 catch(const exception& error) {
82 FATAL(error.what() << endl);
83 }
84
85
86 gRandom->SetSeed(seed);
87
89
90 const size_t nx = 20 *
ls;
91 const double xmin = -5.0 *
ls;
92 const double xmax = +5.0 *
ls;
93
94 const size_t ny = 20 *
ls;
95 const double ymin = -5.0 *
ls;
96 const double ymax = +5.0 *
ls;
97
98 const size_t nz = 20 *
ls;
99 const double zmin = -5.0 *
ls;
100 const double zmax = +5.0 *
ls;
101
102 TH3D* h3 = NULL;
103
104 if (inputFile == "") {
105
106 h3 = new TH3D("h3", NULL, nx, xmin, xmax, ny, ymin, ymax, nz, zmin, zmax);
107
108 h3->Sumw2();
109
110 auto f3 = (
JGauss3X<1>(parameters.center, parameters.sigma) *
113 JP0<1>(parameters.signal) +
114 JP0<2>(parameters.background));
115
116 if (N == 0)
118 else
120
121 } else {
122
123 TFile* in = TFile::Open(inputFile.c_str(), "exist");
124
125 in->GetObject("h3", h3);
126 }
127
128
129
130
131 const double x0 = h3->GetMean(1);
132 const double y0 = h3->GetMean(2);
133 const double z0 = h3->GetMean(3);
134 const double xs = h3->GetStdDev(1) * 0.17;
135 const double ys = h3->GetStdDev(2) * 0.17;
136 const double zs = h3->GetStdDev(3) * 0.17;
137 const double signal = h3->GetMaximum();
138 const double background = h3->GetMinimum() + 0.10;
139
140
141
142
147
148 typedef decltype(
f3) function_type;
149
151
152 for (size_t i = 0; i != getNumberOfParameters<function_type>(); ++i) {
153 cout << setw(2) << i << ' '
154 <<
FIXED(12,5) << getParameter(
f3, i) << endl;
155 }
156
157 {
158 const chrono::steady_clock::time_point t0 = chrono::steady_clock::now();
159
160 const auto result = (option == count_t ?
161 (writeFits ? Fit<m_count>(h3,
f3, {}, X, Y, Z) :
Fit<
m_count>(*h3,
f3, {}, X, Y, Z)) :
162 (writeFits ?
Fit<
m_value>(h3,
f3, {}, X, Y, Z) :
Fit<
m_value>(*h3,
f3, {}, X, Y, Z)));
163
164 const chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
165
166 cout <<
"chi2/NDF " <<
FIXED(7,3) <<
result.getChi2() <<
"/" <<
result.getNDF() << endl;
167 cout <<
"Number of iterations " <<
result.numberOfIterations << endl;
168 cout << "Elapsed time [us] " << setw(8) << chrono::duration_cast<chrono::microseconds>(t1 - t0).count() << endl;
169
170 for (
size_t i = 0; i !=
result.getNumberOfParameters(); ++i) {
171 cout << setw(2) << i << ' '
174 }
175 }
176
177
179
180 out << *h3;
181
182 out.Write();
183 out.Close();
184}
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
double f3(const double x, const double y, const double z)
3D function.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Utility class to parse parameter values.
Utility class to parse command line options.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
void FillRandom(TH1 *h1, const T &f1)
Fill 1D histogram according Poisson statistics with expectation values from given 1D function.
JRootfit_t< JF1_t > Fit(const TH1 &h1, const JF1_t &f1, const index_list &ls=index_list(), const range_type &X=range_type())
Global fit fuction.
Auxiliary data structure for floating point format specification.
Make 3D function of x from 1D function.
Make 3D function of y from 1D function.
Make 3D function of z from 1D function.
Termination class for polynomial function.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Data point for counting measurement.
Data point for value with error.
Auxiliary data structure to list files in directory.