Utilities/test/testIntegral.cc

0001 #include <cppunit/extensions/HelperMacros.h>
0002 #include "PhysicsTools/Utilities/interface/Operations.h"
0003 #include "PhysicsTools/Utilities/interface/Integral.h"
0004 #include "PhysicsTools/Utilities/interface/Variables.h"
0005 #include "PhysicsTools/Utilities/interface/Expression.h"
0006 #include "PhysicsTools/Utilities/interface/Parameter.h"
0007
0008 class testIntegral : public CppUnit::TestFixture {
0009   CPPUNIT_TEST_SUITE(testIntegral);
0010   CPPUNIT_TEST(checkAll);
0011   CPPUNIT_TEST_SUITE_END();
0012
0013 public:
0014   void setUp() {}
0015   void tearDown() {}
0016   void checkAll();
0017 };
0018
0019 NUMERICAL_INTEGRAL(X, Expression, TrapezoidIntegrator);
0020
0021 struct gauss {
0022   static const double c;
0023   double operator()(double x) const { return c * exp(-x * x); }
0024 };
0025
0026 struct gauss2 : public gauss {};
0027
0028 struct gauss3 : public gauss {};
0029
0030 struct gauss4 : public gauss {};
0031
0032 struct gaussPrimitive {
0033   double operator()(double x) const { return erf(x); }
0034 };
0035
0036 const double gauss::c = 2. / sqrt(M_PI);
0037
0038 NUMERICAL_FUNCT_INTEGRAL(gauss, TrapezoidIntegrator);
0039
0040 DECLARE_FUNCT_PRIMITIVE(gauss2, gaussPrimitive);
0041
0042 NUMERICAL_FUNCT_INTEGRAL(gauss3, GaussLegendreIntegrator);
0043
0044 NUMERICAL_FUNCT_INTEGRAL(gauss4, GaussIntegrator);
0045
0046 CPPUNIT_TEST_SUITE_REGISTRATION(testIntegral);
0047
0048 void testIntegral::checkAll() {
0049   using namespace funct;
0050   X x;
0051   double epsilon = 1.e-6;
0052   // symbolic primitive exists
0053   CPPUNIT_ASSERT(fabs(integral<X>(x, 0, 1) - 0.5) < epsilon);
0054   CPPUNIT_ASSERT(fabs(integral<X>(x ^ num<2>(), 0, 1) - 1. / 3.) < epsilon);
0055
0056   TrapezoidIntegrator integrator(1000);
0057
0058   // numerical integration
0059   Expression f_x = x;
0060   CPPUNIT_ASSERT(fabs(integral<X>(f_x, 0, 1, integrator) - 0.5) < epsilon);
0061
0062   Expression f_x2 = (x ^ num<2>());
0063   CPPUNIT_ASSERT(fabs(integral<X>(f_x2, 0, 1, integrator) - 1. / 3.) < epsilon);
0064
0065   Parameter c("c", gauss::c);
0066   Expression f = c * exp(-(x ^ num<2>()));
0067   CPPUNIT_ASSERT(fabs(integral<X>(f, 0, 1, integrator) - erf(1)) < epsilon);
0068
0069   // trapezoid integration
0070   gauss g;
0071   CPPUNIT_ASSERT(fabs(integral_f(g, 0, 1, integrator) - erf(1)) < epsilon);
0072
0073   // user-defined integration
0074   gauss2 g2;
0075   CPPUNIT_ASSERT(fabs(integral_f(g2, 0, 1) - erf(1)) < epsilon);
0076
0077   // Gauss-Legendre integration
0078   gauss3 g3;
0079   GaussLegendreIntegrator integrator2(10000, 1.e-5);
0080   CPPUNIT_ASSERT(fabs(integral_f(g3, 0, 1, integrator2) - erf(1)) < epsilon);
0081
0082   // Gauss-Legendre integration
0083   gauss4 g4;
0084   GaussIntegrator integrator3(1.e-6);
0085   CPPUNIT_ASSERT(fabs(integral_f(g4, 0, 1, integrator3) - erf(1)) < epsilon);
0086
0087   // automatic (trivial) integration
0088   Parameter pi("pi", M_PI);
0089   CPPUNIT_ASSERT(fabs(integral_f(pi, 0, 2) - 2 * pi()) < epsilon);
0090 }