28 using TMath::PoissonI;
39 inline double g1(
const JGauss&
gauss,
const JElement_t& point)
41 const double u = (point.getX() - gauss.
mean) / gauss.
sigma;
42 const double fs = gauss.
signal *
exp(-0.5*u*u) / (sqrt(2.0*
PI) * gauss.
sigma);
44 const double f1 = fs + fb;
46 return -log(PoissonI(point.getY(), f1));
57 int main(
int argc,
char **argv)
70 JParser<> zap(
"Program to test JSimplex algorithm.");
80 catch(
const exception& error) {
81 FATAL(error.what() << endl);
84 ASSERT(numberOfEvents > 0);
88 TF1 fs(
"fs",
"exp(-0.5 * (x-[0])*(x-[0]) / ([1]*[1]))");
91 fs.FixParameter(0, gauss.
mean);
92 fs.FixParameter(1, gauss.
sigma);
95 const Double_t xmin = -5.0;
96 const Double_t xmax = +5.0;
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) };
110 for (
int i = 0; i != numberOfEvents; ++i) {
112 STATUS(
"event: " << setw(10) << i <<
'\r');
DEBUG(endl);
114 TH1D h0(
"h0", NULL, nx, xmin, xmax);
116 h0.FillRandom(
"fs", (Int_t) gauss.
signal);
117 h0.FillRandom(
"fb", (Int_t) gauss.
background);
121 for (Int_t i = 1; i <= h0.GetNbinsX(); ++i) {
122 data.push_back(JElement_t(h0.GetBinCenter (i),
123 h0.GetBinContent(i)));
137 h0.GetEntries() * (xmax - xmin) / nx - h0.GetMinimum(),
142 for (
size_t i = 0; i != fit.
step.size(); ++i) {
143 DEBUG(
"Step size " << i <<
' ' << fit.
step[i] << endl);
148 const double chi2 = fit(
g1, data.begin(), data.end());
153 DEBUG(
"Chi2 " << chi2 << endl);
160 for (
int i = 0; i !=
sizeof(Q)/
sizeof(Q[0]); ++i) {
166 for (
int i = 0; i !=
sizeof(Q)/
sizeof(Q[0]); ++i) {
171 timer.print(cout,
true,
micro_t);
178 for (
int i = 0; i !=
sizeof(
H)/
sizeof(
H[0]); ++i) {
186 for (
int i = 0; i !=
sizeof(Q)/
sizeof(Q[0]); ++i) {
190 ASSERT(Q[0].getSTDev() < precision.mean);
191 ASSERT(Q[1].getSTDev() < precision.sigma);
192 ASSERT(Q[2].getSTDev() < precision.signal);
193 ASSERT(Q[3].getSTDev() < precision.background);
Utility class to parse command line options.
int main(int argc, char *argv[])
The elements in a collection are sorted according to their abscissa values and a given distance opera...
static const double H
Planck constant [eV s].
then fatal Wrong number of arguments fi set_variable STRING $argv[1] set_variable DETECTORXY_TXT $WORKDIR $DETECTORXY_TXT tail read X Y CHI2 RMS printf optimum n $X $Y $CHI2 $RMS awk v Y
then set_variable FORMULA *[0] exp(-0.5 *(x-[1])*(x-[1])/([2]*[2]))*exp(-0.5 *(y-[1])*(y-[1])/([2]*[2]))" JF2 -o $WORKDIR/f2.root -F "$FORMULA" -@ "p0
double getMean(vector< double > &v)
get mean of vector content
#define ASSERT(A,...)
Assert macro.
I/O formatting auxiliaries.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Auxiliary class for CPU timing and usage.
static const double PI
Mathematical constants.
std::ostream & shortprint(std::ostream &out)
Set short printing.
General purpose messaging.
std::ostream & longprint(std::ostream &out)
Set long printing.
Simple fit method based on Powell's algorithm, see reference: Numerical Recipes in C++...
Utility class to parse command line options.
double gauss(const double x, const double sigma)
Gauss function (normalised to 1 at x = 0).
#define DEBUG(A)
Message macros.
std::vector< JModel_t > step
Double_t g1(const Double_t x)
Function.