58{
61
63 int numberOfEvents;
67
68 try {
69
70 JParser<> zap(
"Program to test JSimplex algorithm.");
71
77
78 zap(argc, argv);
79 }
80 catch(const exception& error) {
81 FATAL(error.what() << endl);
82 }
83
84 ASSERT(numberOfEvents > 0);
85
87
88 TF1 fs("fs", "exp(-0.5 * (x-[0])*(x-[0]) / ([1]*[1]))");
89 TF1 fb("fb", "1.0");
90
91 fs.FixParameter(0,
gauss.mean);
92 fs.FixParameter(1,
gauss.sigma);
93
94 const Int_t nx = 21;
95 const Double_t
xmin = -5.0;
96 const Double_t
xmax = +5.0;
97
102
103 TH1D
H[] = { TH1D(
"ha",
"", 101, -0.1, +0.1),
104 TH1D("hb", "", 101, -0.1, +0.1),
105 TH1D("hc", "", 101, -100.0, +100.0),
106 TH1D("hd", "", 101, -100.0, +100.0) };
107
109
110 for (int i = 0; i != numberOfEvents; ++i) {
111
112 STATUS(
"event: " << setw(10) << i <<
'\r');
DEBUG(endl);
113
114 TH1D h0("h0", NULL, nx, xmin, xmax);
115
116 h0.FillRandom(
"fs", (Int_t)
gauss.signal);
117 h0.FillRandom(
"fb", (Int_t)
gauss.background);
118
120
121 for (Int_t i = 1; i <= h0.GetNbinsX(); ++i) {
122 data.push_back(JElement_t(h0.GetBinCenter (i),
123 h0.GetBinContent(i)));
124 }
125
127
129
134
136 h0.GetRMS(),
137 h0.GetEntries() * (xmax - xmin) / nx - h0.GetMinimum(),
138 h0.GetMinimum());
139
141
142 for (
size_t i = 0; i != fit.
step.size(); ++i) {
143 DEBUG(
"Step size " << i <<
' ' << fit.
step[i] << endl);
144 }
145
147
148 const double chi2 = fit(
g1,
data.begin(),
data.end());
149
151
153 DEBUG(
"Chi2 " << chi2 << endl);
154
155 const double Y[] = { fit.
value.mean -
gauss.mean,
158 fit.
value.background * nx -
gauss.background };
159
160 for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
163 }
164 }
165
166 for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
168 }
169
172 }
173
175
177
178 for (
int i = 0; i !=
sizeof(
H)/
sizeof(
H[0]); ++i) {
180 }
181
182 out.Write();
183 out.Close();
184 }
185
186 for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
188 }
189
194
195 return 0;
196}
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.
Auxiliary class for CPU timing and usage.
void print(std::ostream &out, const JScale_t scale=milli_t) const
Print timer data.
Simple fit method based on Powell's algorithm, see reference: Numerical Recipes in C++,...
std::vector< JModel_t > step
Utility class to parse command line options.
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).