// example of a double integral using GSL integrator
// must be run in compiled mode
// .x exampleDoubleIntegral.C+

#include "Math/GSLIntegrator.h"
#include "Math/WrappedFunction.h"
#include "Math/Functor.h"
#include <cmath>
#include <iostream>


using namespace ROOT::Math;

// starting integrand function

double f (double z, double x, double y) { 
  return z*z + x*x + y*y;
} 


// first integrand function (to be integrated in y) 

struct IntegrandFunctionY { 

  IntegrandFunctionY(double z, double x) : fZ(z), fX(x) {}

  // function of y and use x as data member
  double operator() (double y) const { 
    return f(fZ, fX,y); 
  }

  double fZ, fX;
};






// second integrand function (to be integrated in x) 

class IntegrandFunctionX { 

public: 
  IntegrandFunctionX(double z) :
    fZ(z),
    fIntegrator( new GSLIntegrator( 1.E-9,1.E-9,1000) )
  {}

  ~IntegrandFunctionX() { delete fIntegrator; }
 
  // function of y and use x as data member
  double operator() (double x) const { 

    IntegrandFunctionY fy(fZ,x);
    // wrapped the function object fFuncY in the required interface by GSLIntegrator
    // use reference
    WrappedFunction<IntegrandFunctionY > wf(fy); 
    double a = - std::sqrt(1. - x*x);
    double b = - a;

    return fIntegrator->Integral(wf,a,b); 
  }
private: 

  //avoid copying 
  IntegrandFunctionX(const IntegrandFunctionX &) {}
  IntegrandFunctionX & operator=(const IntegrandFunctionX &) {return *this;}
  double fZ; 
  mutable GSLIntegrator  * fIntegrator;
  
};


// Third integrand function (to be integrated in z) 

class IntegrandFunctionZ { 

public: 
  IntegrandFunctionZ() :
    fIntegrator( new GSLIntegrator( 1.E-9,1.E-9,1000) )
  {}

  ~IntegrandFunctionZ() { delete fIntegrator; }
 
  // function of y and use x as data member
  double operator() (double z) const { 

    IntegrandFunctionX fx(z);
    // wrapped the function object fFuncY in the required interface by GSLIntegrator
    // use reference
    WrappedFunction<IntegrandFunctionX > wf(fx); 
    double a = - 1;
    double b =  -a;

    return fIntegrator->Integral(wf,a,b); 
  }
private: 

  //avoid copying 
  IntegrandFunctionZ(const IntegrandFunctionZ &) {}
  IntegrandFunctionZ & operator=(const IntegrandFunctionZ &) {return *this;}
   
  mutable GSLIntegrator  * fIntegrator;
};












void exampleTripleIntegral() { 

  // calculate the integral in x
  GSLIntegrator  ig(1.E-9,1.E-9,1000); 
  IntegrandFunctionZ fz;
  //   Functor1D func( new IntegrandFunctionX() );
  // use  reference as type to avoid copy
  WrappedFunction<const IntegrandFunctionZ &> func(fz);
  double val = ig.Integral(func, 0.,1.);
  std::cout << "integral is " << val << std::endl;
}

int main() { 
  exampleTripleIntegral();
}