Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:33:33

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 }