Jpp  master_rocky-40-g5f0272dcd
the software that should make you happy
JRootfitToGauss3D.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <vector>
5 #include <type_traits>
6 #include <chrono>
7 
8 #include "TROOT.h"
9 #include "TFile.h"
10 #include "TH3D.h"
11 
12 #include "JROOT/JRootToolkit.hh"
13 #include "JROOT/JRootTestkit.hh"
14 #include "JROOT/JRootfit.hh"
15 
16 #include "JMath/JMathlib3D.hh"
17 
18 #include "JLang/JManip.hh"
19 
20 #include "Jeep/JProperties.hh"
21 #include "Jeep/JParser.hh"
22 
23 
24 /**
25  * \file
26  *
27  * Program to test JRootfit algorithm.
28  * \author mdejong
29  */
30 int main(int argc, char **argv)
31 {
32  using namespace std;
33  using namespace JPP;
34 
35  const char* const count_t = "count";
36  const char* const value_t = "value";
37 
38  struct {
39  double signal = 100;
40  double background = 5;
41  double center = 0.0;
42  double sigma = 1.0;
43  } parameters;
44 
45  string inputFile;
46  string outputFile;
47  range_type X;
48  range_type Y;
49  range_type Z;
50  size_t N;
51  string option;
52  bool writeFits;
53  UInt_t seed;
54  int debug;
55 
56  try {
57 
58  JProperties properties;
59 
60  properties.insert(gmake_property(parameters.signal));
61  properties.insert(gmake_property(parameters.background));
62  properties.insert(gmake_property(parameters.center));
63  properties.insert(gmake_property(parameters.sigma));
64 
65  JParser<> zap("Program to test JRootfit algorithm.");
66 
67  zap['f'] = make_field(inputFile) = "";
68  zap['o'] = make_field(outputFile);
69  zap['@'] = make_field(properties) = JPARSER::initialised();
70  zap['x'] = make_field(X) = range_type();
71  zap['y'] = make_field(Y) = range_type();
72  zap['z'] = make_field(Z) = range_type();
73  zap['n'] = make_field(N) = 0;
74  zap['O'] = make_field(option) = count_t, value_t;
75  zap['w'] = make_field(writeFits);
76  zap['S'] = make_field(seed) = 0;
77  zap['d'] = make_field(debug) = 1;
78 
79  zap(argc, argv);
80  }
81  catch(const exception& error) {
82  FATAL(error.what() << endl);
83  }
84 
85 
86  gRandom->SetSeed(seed);
87 
88  const size_t ls = 2;
89 
90  const size_t nx = 20 * ls;
91  const double xmin = -5.0 * ls;
92  const double xmax = +5.0 * ls;
93 
94  const size_t ny = 20 * ls;
95  const double ymin = -5.0 * ls;
96  const double ymax = +5.0 * ls;
97 
98  const size_t nz = 20 * ls;
99  const double zmin = -5.0 * ls;
100  const double zmax = +5.0 * ls;
101 
102  TH3D* h3 = NULL;
103 
104  if (inputFile == "") {
105 
106  h3 = new TH3D("h3", NULL, nx, xmin, xmax, ny, ymin, ymax, nz, zmin, zmax);
107 
108  h3->Sumw2();
109 
110  auto f3 = (JGauss3X<1>(parameters.center, parameters.sigma) *
111  JGauss3Y<2>(parameters.center, parameters.sigma) *
112  JGauss3Z<3>(parameters.center, parameters.sigma) *
113  JP0<1>(parameters.signal) +
114  JP0<2>(parameters.background));
115 
116  if (N == 0)
117  FillRandom(h3, f3);
118  else
119  FillRandom(h3, f3, N);
120 
121  } else {
122 
123  TFile* in = TFile::Open(inputFile.c_str(), "exist");
124 
125  in->GetObject("h3", h3);
126  }
127 
128 
129  // determine start values from data
130 
131  const double x0 = h3->GetMean(1);
132  const double y0 = h3->GetMean(2);
133  const double z0 = h3->GetMean(3);
134  const double xs = h3->GetStdDev(1) * 0.17;
135  const double ys = h3->GetStdDev(2) * 0.17;
136  const double zs = h3->GetStdDev(3) * 0.17;
137  const double signal = h3->GetMaximum();
138  const double background = h3->GetMinimum() + 0.10;
139 
140 
141  // fit
142 
143  auto f3 = (JGauss3X<1>(x0, xs) *
144  JGauss3Y<2>(y0, ys) *
145  JGauss3Z<3>(z0, zs) *
146  JP0<1>(signal)) + JP0<2>(background);
147 
148  typedef decltype(f3) function_type;
149 
151 
152  for (size_t i = 0; i != getNumberOfParameters<function_type>(); ++i) {
153  cout << setw(2) << i << ' '
154  << FIXED(12,5) << getParameter(f3, i) << endl;
155  }
156 
157  {
158  const chrono::steady_clock::time_point t0 = chrono::steady_clock::now();
159 
160  const auto result = (option == count_t ?
161  (writeFits ? Fit<m_count>(h3, f3, {}, X, Y, Z) : Fit<m_count>(*h3, f3, {}, X, Y, Z)) :
162  (writeFits ? Fit<m_value>(h3, f3, {}, X, Y, Z) : Fit<m_value>(*h3, f3, {}, X, Y, Z)));
163 
164  const chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
165 
166  cout << "chi2/NDF " << FIXED(7,3) << result.getChi2() << "/" << result.getNDF() << endl;
167  cout << "Number of iterations " << result.numberOfIterations << endl;
168  cout << "Elapsed time [us] " << setw(8) << chrono::duration_cast<chrono::microseconds>(t1 - t0).count() << endl;
169 
170  for (size_t i = 0; i != result.getNumberOfParameters(); ++i) {
171  cout << setw(2) << i << ' '
172  << FIXED(12,5) << result.getValue(i) << " +/- "
173  << FIXED(12,5) << result.getError(i) << endl;
174  }
175  }
176 
177 
178  TFile out(outputFile.c_str(), "recreate");
179 
180  out << *h3;
181 
182  out.Write();
183  out.Close();
184 }
string outputFile
I/O manipulators.
Functional algebra in 3D.
#define FATAL(A)
Definition: JMessage.hh:67
int debug
debug level
Definition: JSirene.cc:69
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2142
double f3(const double x, const double y, const double z)
3D function.
Definition: JPolynome3D.cc:23
Utility class to parse parameter values.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
int main(int argc, char **argv)
Utility class to parse parameter values.
Definition: JProperties.hh:501
Utility class to parse command line options.
Definition: JParser.hh:1698
ROOT Fit.
Definition: JRootfit.hh:940
const double sigma[]
Definition: JQuadrature.cc:74
const double xmax
Definition: JQuadrature.cc:24
const double xmin
Definition: JQuadrature.cc:23
int getParameter(const std::string &text)
Get parameter number from text string.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
void FillRandom(TH1 *h1, const T &f1)
Fill 1D histogram according Poisson statistics with expectation values from given 1D function.
Definition: JRootTestkit.hh:33
JRange< int > range_type
Definition: JSTDTypes.hh:14
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:448
Make 3D function of x from 1D function.
Definition: JMathlib3D.hh:25
Make 3D function of y from 1D function.
Definition: JMathlib3D.hh:97
Make 3D function of z from 1D function.
Definition: JMathlib3D.hh:169
Termination class for polynomial function.
Definition: JMathlib.hh:1487
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:68
Auxiliary data structure to list files in directory.
Definition: JFilesystem.hh:20