29 using TMath::PoissonI;
60 this->push_back(
mean);
61 this->push_back(
sigma);
78 inline result_type
g1(
const JGauss&
gauss,
const JElement_t& point)
82 const double u = (point.getX() -
gauss[0]) /
gauss[1];
83 const double fs =
gauss[2] * exp(-0.5*u*u) / (sqrt(2.0*PI) *
gauss[1]);
84 const double fb =
gauss[3];
85 const double f1 = fs + fb;
87 const double p = PoissonI(point.getY(), f1);
96 result.gradient *= 1.0 - point.getY()/f1;
122 JParser<> zap(
"Program to test JGandalf algorithm.");
132 catch(
const exception& error) {
133 FATAL(error.what() << endl);
136 using namespace JFIT;
138 ASSERT(numberOfEvents > 0);
142 TF1 fs(
"fs",
"exp(-0.5 * (x-[0])*(x-[0]) / ([1]*[1]))");
145 fs.FixParameter(0,
gauss[0]);
146 fs.FixParameter(1,
gauss[1]);
149 const Double_t xmin = -5.0;
150 const Double_t xmax = +5.0;
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)
168 for (
int i = 0; i != numberOfEvents; ++i) {
170 STATUS(
"event: " << setw(10) << i <<
'\r');
DEBUG(endl);
172 TH1D h0(
"h0", NULL, nx, xmin, xmax);
176 h0.FillRandom(
"fs", (Int_t)
gauss[2]);
177 h0.FillRandom(
"fb", (Int_t)
gauss[3]);
181 for (Int_t i = 1; i <= h0.GetNbinsX(); ++i) {
182 data.push_back(JElement_t(h0.GetBinCenter (i),
183 h0.GetBinContent(i)));
188 fit.parameters.resize(4);
190 for (
size_t i = 0; i != 4; ++i) {
191 fit.parameters[i] = i;
194 fit.value =
JGauss(h0.GetMean(),
196 h0.GetEntries() * (xmax - xmin) / nx - h0.GetMinimum(),
199 DEBUG(
"Start value " << fit.value << endl);
203 const double chi2 = fit(
g1, data.begin(), data.end());
207 DEBUG(
"Final value " << fit.value << endl);
208 DEBUG(
"Chi2 " << chi2 << endl);
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] };
215 for (
int i = 0; i !=
sizeof(Q)/
sizeof(Q[0]); ++i) {
221 for (
int i = 0; i !=
sizeof(Q)/
sizeof(Q[0]); ++i) {
233 for (
int i = 0; i !=
sizeof(
H)/
sizeof(
H[0]); ++i) {
241 for (
int i = 0; i !=
sizeof(Q)/
sizeof(Q[0]); ++i) {
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]);
double getMean(vector< double > &v)
get mean of vector content
The elements in a collection are sorted according to their abscissa values and a given distance opera...
int main(int argc, char **argv)
std::ostream & longprint(std::ostream &out)
Set long printing.
std::ostream & shortprint(std::ostream &out)
Set short printing.
General purpose messaging.
#define DEBUG(A)
Message macros.
#define ASSERT(A,...)
Assert macro.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
I/O formatting auxiliaries.
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).
Data structure for return value of fit function.
JGauss()
Default constructor.