Jpp - the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JSimplexFitToGauss.cc File Reference

Program to test JFIT::JSimplex algorithm. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>
#include "TROOT.h"
#include "TFile.h"
#include "TF1.h"
#include "TH1D.h"
#include "TMath.h"
#include "JFit/JSimplex.hh"
#include "JTools/JElement.hh"
#include "JTools/JQuantile.hh"
#include "JROOT/JRootToolkit.hh"
#include "JMath/JGauss.hh"
#include "JMath/JConstants.hh"
#include "Jeep/JTimer.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Program to test JFIT::JSimplex algorithm.

Author
mdejong

Definition in file JSimplexFitToGauss.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 57 of file JSimplexFitToGauss.cc.

58 {
59  using namespace std;
60  using namespace JPP;
61 
62  string outputFile;
63  int numberOfEvents;
64  JGauss gauss;
65  JGauss precision;
66  int debug;
67 
68  try {
69 
70  JParser<> zap("Program to test JSimplex algorithm.");
71 
72  zap['o'] = make_field(outputFile) = "";
73  zap['n'] = make_field(numberOfEvents) = 100;
74  zap['@'] = make_field(gauss) = JGauss(0.0, 1.0, 1000.0, 100.0);
75  zap['e'] = make_field(precision) = JGauss(0.05, 0.05, 25.0, 25.0);
76  zap['d'] = make_field(debug) = 3;
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 
98  JQuantile Q[] = { JQuantile("mean "),
99  JQuantile("sigma "),
100  JQuantile("signal "),
101  JQuantile("background") };
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 
108  JTimer timer;
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 
119  vector<JElement_t> data;
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 
126  JSimplex<JGauss> fit;
127 
128  fit.step.resize(4);
129 
130  fit.step[0] = JGauss(0.1, 0.0, 0.0, 0.0);
131  fit.step[1] = JGauss(0.0, 0.1, 0.0, 0.0);
132  fit.step[2] = JGauss(0.0, 0.0, gauss.signal * 0.01, 0.0);
133  fit.step[3] = JGauss(0.0, 0.0, 0.0, gauss.background * 0.01);
134 
135  fit.value = JGauss(h0.GetMean(),
136  h0.GetRMS(),
137  h0.GetEntries() * (xmax - xmin) / nx - h0.GetMinimum(),
138  h0.GetMinimum());
139 
140  DEBUG("Start value " << fit.value << endl);
141 
142  for (size_t i = 0; i != fit.step.size(); ++i) {
143  DEBUG("Step size " << i << ' ' << fit.step[i] << endl);
144  }
145 
146  timer.start();
147 
148  const double chi2 = fit(g1, data.begin(), data.end());
149 
150  timer.stop();
151 
152  DEBUG("Final value " << fit.value << endl);
153  DEBUG("Chi2 " << chi2 << endl);
154 
155  const double Y[] = { fit.value.mean - gauss.mean,
156  fit.value.sigma - gauss.sigma,
157  fit.value.signal * nx / (xmax - xmin) - gauss.signal,
158  fit.value.background * nx - gauss.background };
159 
160  for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
161  Q[i].put (Y[i]);
162  H[i].Fill(Y[i]);
163  }
164  }
165 
166  for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
167  NOTICE((i == 0 ? longprint : shortprint) << Q[i]);
168  }
169 
170  if (debug >= notice_t) {
171  timer.print(cout, true, micro_t);
172  }
173 
174  if (outputFile != "") {
175 
176  TFile out(outputFile.c_str(), "recreate");
177 
178  for (int i = 0; i != sizeof(H)/sizeof(H[0]); ++i) {
179  out << H[i];
180  }
181 
182  out.Write();
183  out.Close();
184  }
185 
186  for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
187  ASSERT(fabs(Q[i].getMean()) < Q[i].getSTDev());
188  }
189 
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);
194 
195  return 0;
196 }
Utility class to parse command line options.
Definition: JParser.hh:1500
double mean
Definition: JGauss.hh:161
micro
Definition: JScale.hh:29
Auxiliary data structure for running average, standard deviation and quantiles.
Definition: JQuantile.hh:43
#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
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
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
Auxiliary class for CPU timing and usage.
Definition: JTimer.hh:32
#define NOTICE(A)
Definition: JMessage.hh:64
int debug
debug level
Definition: JSirene.cc:63
std::ostream & shortprint(std::ostream &out)
Set short printing.
Definition: JManip.hh:144
JModel_t value
Definition: JSimplex.hh:240
std::ostream & longprint(std::ostream &out)
Set long printing.
Definition: JManip.hh:172
#define FATAL(A)
Definition: JMessage.hh:67
double signal
Definition: JGauss.hh:163
Simple fit method based on Powell&#39;s algorithm, see reference: Numerical Recipes in C++...
Definition: JSimplex.hh:42
double background
Definition: JGauss.hh:164
double sigma
Definition: JGauss.hh:162
double gauss(const double x, const double sigma)
Gauss function (normalised to 1 at x = 0).
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
std::vector< JModel_t > step
Definition: JSimplex.hh:241
Double_t g1(const Double_t x)
Function.
Definition: JQuantiles.cc:25