Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JTrigonometric1D.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <sstream>
4 #include <iomanip>
5 
6 #include "TRandom3.h"
7 
8 #include "JLang/JException.hh"
10 #include "JTools/JFunction1D_t.hh"
11 #include "JTools/JGrid.hh"
12 #include "JTools/JQuantile.hh"
13 #include "JTools/JConstants.hh"
14 
15 #include "Jeep/JParser.hh"
16 #include "Jeep/JMessage.hh"
17 
18 namespace {
19 
20  /**
21  * Test method for interpolation with derivatives.
22  *
23  * Accuracy of nth derivative of N degree polynome
24  * is equal to n-jth derivative of N-j degree polynome.
25  */
26  template<int N>
27  int do_main(int argc, char **argv)
28  {
29  using namespace std;
30  using namespace JPP;
31 
32  unsigned int numberOfEvents;
33  unsigned int numberOfBins;
34  int debug;
35 
36  try {
37 
38  JParser<> zap("Example program to test 1D interpolation of a polynome.");
39 
40  zap['n'] = make_field(numberOfEvents);
41  zap['N'] = make_field(numberOfBins) = 21;
42  zap['d'] = make_field(debug) = 3;
43 
44  zap(argc, argv);
45  }
46  catch(const exception &error) {
47  FATAL(error.what() << endl);
48  }
49 
50  if (N < 1) {
51  FATAL("Number of degrees " << N << " < 1." << endl);
52  }
53 
54 
55  using namespace JPP;
56 
57  // Analytical function and its derivatives.
58 
59  vector<JTrigonometric> f1(1, JTrigonometric(sin));
60 
61  for (int i = 1; i != N + 1; ++i) {
62  f1.push_back(f1.rbegin()->getDerivative());
63  }
64 
65 
66  typedef JResultPolynome<N, double> JResult_t;
67  typedef JPolintFunction1D<N, JElement2D<double, double>, JCollection, JResult_t> JFunction1D_t;
68 
69  JFunction1D_t polint;
70 
71  const double xmin = 0.0;
72  const double xmax = PI;
73 
74  //const double dx = ((N+1)/2) * (xmax - xmin) / (double) (numberOfBins - 2*((N+1)/2));
75  const double dx = 0.0 * PI;
76 
77  const JGrid<double> grid(numberOfBins + 1, xmin - dx, xmax + dx);
78 
79  polint.configure(grid, f1[0]);
80 
81  polint.compile();
82 
83  polint.setExceptionHandler(new typename JFunction1D_t::JDefaultResult(JMATH::zero));
84 
85 
86  if (numberOfEvents != 0) {
87 
88  JQuantile Q[N+1];
89 
90  for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
91 
92  ostringstream os;
93 
94  os << "f^" << i;
95 
96  Q[i].setTitle(os.str());
97  }
98 
99  for (unsigned int i = 0; i != numberOfEvents; ++i) {
100 
101  const double x = gRandom->Uniform(xmin, xmax);
102  const JResult_t result = polint(x);
103 
104  for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
105  Q[i].put(f1[i](x) - result.y[i]);
106  }
107  }
108 
109  cout << "\nInterpolation with " << N << " degrees polynomial:" << endl;
110 
111  for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
112  Q[i].print(cout, i == 0);
113  }
114  }
115 
116  return 0;
117  }
118 }
119 
120 /**
121  * \file
122  *
123  * Example program to test 1D interpolation of a trigonometric function.
124  * \author mdejong
125  */
126 int main(int argc, char **argv)
127 {
128  if (do_main<2>(argc, argv) != 0) { return 1; }
129  if (do_main<3>(argc, argv) != 0) { return 1; }
130  if (do_main<4>(argc, argv) != 0) { return 1; }
131  if (do_main<5>(argc, argv) != 0) { return 1; }
132 }
Utility class to parse command line options.
Definition: JParser.hh:1410
Exceptions.
static const double PI
Constants.
Definition: JConstants.hh:20
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:94
Constants.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1836
int debug
debug level
Definition: JSirene.cc:59
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:65
Utility class to parse command line options.
int main(int argc, char *argv[])
Definition: Main.cpp:15