Jpp
JTools/JHistogram2D.cc
Go to the documentation of this file.
1 
2 #include <string>
3 #include <iostream>
4 #include <iomanip>
5 
6 #include "TROOT.h"
7 #include "TFile.h"
8 #include "TH2D.h"
9 #include "TMath.h"
10 #include "TRandom.h"
11 
12 #include "JTools/JHistogram1D_t.hh"
14 #include "JTools/JCollection.hh"
15 #include "JTools/JMultiGrid.hh"
16 #include "JTools/JMultiSet.hh"
17 #include "JTools/JQuadrature.hh"
18 #include "JTools/JFunction1D_t.hh"
20 #include "JTools/JToolsToolkit.hh"
21 #include "JTools/JMultiPDF.hh"
22 
23 #include "Jeep/JParser.hh"
24 #include "Jeep/JMessage.hh"
25 
26 
27 namespace {
28  /**
29  * Function g1.
30  *
31  * \param x x value
32  * \return function value
33  */
34  inline Double_t g1(const Double_t x)
35  {
36  return TMath::Gaus(x, 0.0, 1.0, kTRUE);
37  }
38 
39 
40  /**
41  * Integral of g1.
42  *
43  * \param x x value
44  * \return integral value
45  */
46  inline Double_t G1(const Double_t x)
47  {
48  return 0.5 * (1.0 + TMath::Erf(sqrt(0.5)*x));
49  }
50 
51  /**
52  * Function g2.
53  *
54  * \param x x value
55  * \param y y value
56  * \return function value
57  */
58  inline Double_t g2(const Double_t x, const Double_t y)
59  {
60  return g1(x) * g1(y);
61  }
62 
63 
64  /**
65  * Integral of g2.
66  *
67  * \param x x value
68  * \param y y value
69  * \return integral value
70  */
71  inline Double_t G2(const Double_t x, const Double_t y)
72  {
73  return G1(x) * G1(y);
74  }
75 }
76 
77 
78 /**
79  * \file
80  *
81  * Example program to test JTOOLS::JHistogram2D.
82  * \author mdejong
83  */
84 int main(int argc, char **argv)
85 {
86  using namespace std;
87 
88  string outputFile;
89  int numberOfEvents;
90  int numberOfBins;
91  bool quadrature;
92  bool subtract;
93  int debug;
94 
95  try {
96 
97  JParser<> zap("Example program to test 2D histogram.");
98 
99  zap['o'] = make_field(outputFile) = "histogram.root";
100  zap['n'] = make_field(numberOfEvents);
101  zap['N'] = make_field(numberOfBins);
102  zap['Q'] = make_field(quadrature);
103  zap['D'] = make_field(subtract);
104  zap['d'] = make_field(debug) = 3;
105 
106  zap(argc, argv);
107  }
108  catch(const exception &error) {
109  FATAL(error.what() << endl);
110  }
111 
112 
113  if (numberOfEvents <= 0) {
114  FATAL("No events." << endl);
115  }
116 
117  using namespace JPP;
118 
119  const Int_t nx = 1000;
120  const Double_t xmin = -3.0;
121  const Double_t xmax = +3.0;
122 
123  typedef JMultiHistogram<JHistogram1D_t,
124  JMapList<JHistogramMap_t> > JHistogram2D_t;
125 
126  JHistogram2D_t histogram;
127 
128 
129  // abscissa
130 
131  if (!quadrature)
132  configure(histogram, make_grid(numberOfBins + 1, xmin, xmax));
133  else
134  configure(histogram, make_set(JGaussHermite(numberOfBins + 1)));
135 
136 
137  // fill
138 
139  for (int i = 0; i != numberOfEvents; ++i) {
140 
141  const double x = gRandom->Gaus(0.0, 1.0);
142  const double y = gRandom->Gaus(0.0, 1.0);
143 
144  histogram.fill(x, y, 1.0);
145  }
146 
147 
148  histogram.div((double) numberOfEvents);
149 
150 
151  // interpolation based on function values
152 
153  JMultiPDF<JPolint2Function1D_t,
154  JMapList<JPolint2FunctionalMap> > f2(histogram);
155 
156 
157  // interpolation based on integral values
158 
159  accumulate(histogram);
160 
161  JMultiFunction<JPolint3Function1H_t,
162  JMapList<JPolint3FunctionalMapH> > F2;
163 
164  copy(histogram, F2);
165 
166  F2.compile();
167 
168 
169  TFile out(outputFile.c_str(), "recreate");
170 
171  TH2D h0("h0", NULL, nx, xmin, xmax, nx, xmin, xmax);
172  TH2D h1("h1", NULL, nx, xmin, xmax, nx, xmin, xmax);
173  TH2D h2("h2", NULL, nx, xmin, xmax, nx, xmin, xmax);
174 
175 
176  for (int i = 1; i <= h0.GetXaxis()->GetNbins(); ++i) {
177  for (int j = 1; j <= h0.GetYaxis()->GetNbins(); ++j) {
178 
179  const Double_t x = h0.GetXaxis()->GetBinCenter(i);
180  const Double_t y = h0.GetYaxis()->GetBinCenter(j);
181 
182  h0.SetBinContent(i, j, g2(x,y));
183 
184  try {
185  h1.SetBinContent(i, j, F2(x,y).fp.fp);
186  h2.SetBinContent(i, j, f2(x,y));
187  }
188  catch(const exception& error) {
189  //ERROR(error.what() << endl);
190  }
191  }
192  }
193 
194  if (subtract) {
195  h1.Add(&h0, -1.0);
196  h2.Add(&h0, -1.0);
197  }
198 
199  out.Write();
200  out.Close();
201 }
202 
quadrature
bool quadrature
Definition: JResultPDF.cc:23
JCollection.hh
JMessage.hh
JHistogramMap_t.hh
g1
Double_t g1(const Double_t x)
Function.
Definition: JQuantiles.cc:25
JQuadrature.hh
JTOOLS::j
int j
Definition: JPolint.hh:634
JMultiSet.hh
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1493
main
int main(int argc, char **argv)
Definition: JTools/JHistogram2D.cc:84
JMultiPDF.hh
JAANET::copy
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:152
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:37
JFunction1D_t.hh
JLANG::make_set
T make_set(T __begin, T __end, JResult_t std::iterator_traits< T >::value_type::*value, const JComparator_t &comparator)
Method to exclude outliers from already sorted data.
Definition: JVectorize.hh:171
debug
int debug
debug level
Definition: JSirene.cc:59
JToolsToolkit.hh
JFunctionalMap_t.hh
JParser.hh
JTOOLS::JHistogram1D_t
JHistogram1D< JElement2D< double, double >, JCollection > JHistogram1D_t
Type definition of a 1 dimensional histogram based on a JCollection.
Definition: JHistogram1D_t.hh:22
JTOOLS::make_grid
JGrid< JAbscissa_t > make_grid(const int nx, const JAbscissa_t Xmin, const JAbscissa_t Xmax)
Helper method for JGrid.
Definition: JGrid.hh:177
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
std
Definition: jaanetDictionary.h:36
JHistogram1D_t.hh
JMultiGrid.hh
G1
Double_t G1(const Double_t x)
Integral of method g1.
Definition: JQuantiles.cc:37
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JTOOLS::accumulate
void accumulate(T &value, JBool< false > option)
Accumulation of value.
Definition: JToolsToolkit.hh:321
JTOOLS::configure
void configure(const T &value, const JAbstractCollection< JAbscissa_t > &bounds, JBool< false > option)
Configuration of value.
Definition: JToolsToolkit.hh:285