Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JGandalfFitToGauss.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <vector>
5 #include <cmath>
6 
7 #include "TROOT.h"
8 #include "TFile.h"
9 #include "TF1.h"
10 #include "TH1D.h"
11 #include "TMath.h"
12 
13 #include "JFit/JGandalf.hh"
14 #include "JTools/JElement.hh"
15 #include "JTools/JQuantile.hh"
16 #include "JROOT/JRootToolkit.hh"
17 #include "JMath/JGauss.hh"
18 #include "JMath/JConstants.hh"
19 
20 #include "Jeep/JTimer.hh"
21 #include "Jeep/JPrint.hh"
22 #include "Jeep/JParser.hh"
23 #include "Jeep/JMessage.hh"
24 
25 
26 namespace {
27 
28  using namespace JPP;
29  using TMath::PoissonI;
30 
31  typedef JElement2D<double, double> JElement_t;
32  typedef JGandalf<JGauss>::result_type result_type;
33 
34  /**
35  * Fit function.
36  *
37  * \param gauss gauss
38  * \param point point
39  * \return chi2
40  */
41  inline result_type g1(const JGauss& gauss, const JElement_t& point)
42  {
43  result_type result;
44 
45  const double u = (point.getX() - gauss.mean) / gauss.sigma;
46  const double fs = gauss.signal * exp(-0.5*u*u) / (sqrt(2.0*PI) * gauss.sigma);
47  const double fb = gauss.background;
48  const double f1 = fs + fb;
49 
50  const double p = PoissonI(point.getY(), f1);
51 
52  result.chi2 = -log(p);
53 
54  result.gradient.mean = fs * u / gauss.sigma; // d(f)/d(mean)
55  result.gradient.sigma = fs * u*u / gauss.sigma - fs / gauss.sigma; // d(f)/d(sigma)
56  result.gradient.signal = fs / gauss.signal; // d(f)/d(signal)
57  result.gradient.background = 1.0; // d(f)/d(background)
58 
59  result.gradient *= 1.0 - point.getY()/f1; // d(chi2)/d(f)
60 
61  return result;
62  }
63 }
64 
65 
66 /**
67  * \file
68  *
69  * Program to test JFIT::JGandalf algorithm.
70  * \author mdejong
71  */
72 int main(int argc, char **argv)
73 {
74  using namespace std;
75  using namespace JPP;
76 
77  string outputFile;
78  int numberOfEvents;
79  JGauss gauss;
80  JGauss precision;
81  int debug;
82 
83  try {
84 
85  JParser<> zap("Program to test JGandalf algorithm.");
86 
87  zap['o'] = make_field(outputFile) = "";
88  zap['n'] = make_field(numberOfEvents) = 100;
89  zap['@'] = make_field(gauss) = JGauss(0.0, 1.0, 1000.0, 100.0);
90  zap['e'] = make_field(precision) = JGauss(0.05, 0.05, 25.0, 25.0);
91  zap['d'] = make_field(debug) = 3;
92 
93  zap(argc, argv);
94  }
95  catch(const exception& error) {
96  FATAL(error.what() << endl);
97  }
98 
99  using namespace JFIT;
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 
115  JQuantile Q[] = { JQuantile("mean "),
116  JQuantile("sigma "),
117  JQuantile("signal "),
118  JQuantile("background") };
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 
125  JTimer timer;
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 
138  vector<JElement_t> data;
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 
145  JGandalf<JGauss> fit;
146 
147  fit.parameters.resize(4);
148 
149  fit.parameters[0] = &JGauss::mean;
150  fit.parameters[1] = &JGauss::sigma;
151  fit.parameters[2] = &JGauss::signal;
152  fit.parameters[3] = &JGauss::background;
153 
154  fit.value = JGauss(h0.GetMean(),
155  h0.GetRMS(),
156  h0.GetEntries() * (xmax - xmin) / nx - h0.GetMinimum(),
157  h0.GetMinimum());
158 
159  DEBUG("Start value " << fit.value << endl);
160 
161  timer.start();
162 
163  const double chi2 = fit(g1, data.begin(), data.end());
164 
165  timer.stop();
166 
167  DEBUG("Final value " << fit.value << endl);
168  DEBUG("Chi2 " << chi2 << endl);
169 
170  const double Y[] = { fit.value.mean - gauss.mean,
171  fit.value.sigma - gauss.sigma,
172  fit.value.signal * nx / (xmax - xmin) - gauss.signal,
173  fit.value.background * nx - gauss.background };
174 
175  for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
176  Q[i].put (Y[i]);
177  H[i].Fill(Y[i]);
178  }
179  }
180 
181  for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
182  NOTICE((i == 0 ? longprint : shortprint) << Q[i]);
183  }
184 
185  if (debug >= notice_t) {
186  timer.print(cout, true, micro_t);
187  }
188 
189  if (outputFile != "") {
190 
191  TFile out(outputFile.c_str(), "recreate");
192 
193  for (int i = 0; i != sizeof(H)/sizeof(H[0]); ++i) {
194  out << H[i];
195  }
196 
197  out.Write();
198  out.Close();
199  }
200 
201  for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
202  ASSERT(fabs(Q[i].getMean()) < Q[i].getSTDev());
203  }
204 
205  ASSERT(Q[0].getSTDev() < precision.mean);
206  ASSERT(Q[1].getSTDev() < precision.sigma);
207  ASSERT(Q[2].getSTDev() < precision.signal);
208  ASSERT(Q[3].getSTDev() < precision.background);
209 
210  return 0;
211 }
Utility class to parse command line options.
Definition: JParser.hh:1500
double mean
Definition: JGauss.hh:161
The elements in a collection are sorted according to their abscissa values and a given distance opera...
micro
Definition: JScale.hh:29
Quantile calculator.
Definition: JQuantile.hh:88
#define STATUS(A)
Definition: JMessage.hh:63
notice
Definition: JMessage.hh:32
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
string outputFile
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
Gauss function object.
Definition: JGauss.hh:173
#define ASSERT(A,...)
Assert macro.
Definition: JMessage.hh:90
I/O formatting auxiliaries.
Mathematical constants.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
return result
Definition: JPolint.hh:727
Auxiliary class for CPU timing and usage.
Definition: JTimer.hh:32
#define NOTICE(A)
Definition: JMessage.hh:64
static const double PI
Mathematical constants.
int debug
debug level
Definition: JSirene.cc:63
std::ostream & shortprint(std::ostream &out)
Set short printing.
Definition: JManip.hh:143
General purpose messaging.
std::ostream & longprint(std::ostream &out)
Set long printing.
Definition: JManip.hh:171
#define FATAL(A)
Definition: JMessage.hh:67
double signal
Definition: JGauss.hh:163
Utility class to parse command line options.
double background
Definition: JGauss.hh:164
Data structure for return value of fit function.
Definition: JGandalf.hh:69
double sigma
Definition: JGauss.hh:162
double gauss(const double x, const double sigma)
Gauss function (normalised to 1 at x = 0).
2D Element.
Definition: JElement.hh:46
double u[N+1]
Definition: JPolint.hh:739
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
Double_t g1(const Double_t x)
Function.
Definition: JQuantiles.cc:25
int main(int argc, char *argv[])
Definition: Main.cpp:15