/// Post: http://root.cern.ch/phpBB2/viewtopic.php?t=10213

/// Doc: http://root.cern.ch/drupal/content/function-integration
///      http://www.gnu.org/software/gsl/manual/html_node/Numerical-Integration-Introduction.html
///      http://project-mathlibs.web.cern.ch/project-mathlibs/sw/html/classROOT_1_1Math_1_1IntegratorOneDim.html

/// !!! This script has to be compiled (>>.L zzztest.c+)
/// !!! Don't pay attention to all the error messages during execution. 
///      - SetAbsTolerance() is not implemented for GaussIntegrator and not used in GaussLegendreIntegrator
///      - IntegratorOneDim::Integral (const vector<double> &pts) does not work for Integrator/kNONADAPTIVE & Integrator/kADAPTIVE whatever the rule used


#include <TF1.h>
#include <Math/Integrator.h>
#include <Math/WrappedTF1.h>
#include <Math/GSLIntegrator.h>
#include <Math/GaussIntegrator.h>
#include <Math/GaussLegendreIntegrator.h>
#include <Math/Functor.h>
#include <Math/OneDimFunctionAdapter.h>

#include <vector>
#include <functional>


Double_t F_rect ( const Double_t *x, const Double_t *par )
{
	if ( x[0]<par[0] || x[0]>par[1] ) { return (0.) ; } else { return (1.) ; }
}

/////

void zzztest ( Double_t xrange=15.,Double_t xshift=0. )
{
	Double_t xmin=-xrange, xmax=xrange ;
	TF1 *f_rect = new TF1 ("f_rect", F_rect, xmin,xmax, 2) ; 
	f_rect->SetParameters(-0.5+xshift,+0.5+xshift) ;
	f_rect->SetNpx(1000) ; f_rect->Draw() ;

	Double_t abstol=1e-12 ; Double_t reltol=1e-12 ; 
	UInt_t size=100000 ;  /// maximum number of sub-intervals (useless below a certain threshold)
	Int_t npoints=-1 ;  /// only for GaussLegendreIntegrator
	Int_t rule=-1 ;  /// only for GSL kADAPTIVE: 1-6 (lower rules are indicated for singular functions while higher for smooth functions to get better accuracies 

	ROOT::Math::WrappedTF1 wf (*f_rect) ;

	ROOT::Math::GaussIntegrator         ig1    ; ig1   .SetFunction(wf) ; ig1   .SetAbsTolerance(abstol) ; ig1   .SetRelTolerance(reltol) ;
	ROOT::Math::GaussLegendreIntegrator ig21e3 ; ig21e3.SetFunction(wf) ; ig21e3.SetAbsTolerance(abstol) ; ig21e3.SetRelTolerance(reltol) ; ig21e3.SetNumberPoints( 1000) ;
	ROOT::Math::GaussLegendreIntegrator ig21e4 ; ig21e4.SetFunction(wf) ; ig21e4.SetAbsTolerance(abstol) ; ig21e4.SetRelTolerance(reltol) ; ig21e4.SetNumberPoints(10000) ;
	ROOT::Math::GSLIntegrator           ig31   ; ig31  .SetFunction(wf) ; ig31  .SetAbsTolerance(abstol) ; ig31  .SetRelTolerance(reltol) ; ig31  .SetIntegrationRule(ROOT::Math::Integration::GKRule(1)) ;  /// kGAUSS15
	ROOT::Math::GSLIntegrator           ig32   ; ig32  .SetFunction(wf) ; ig32  .SetAbsTolerance(abstol) ; ig32  .SetRelTolerance(reltol) ; ig32  .SetIntegrationRule(ROOT::Math::Integration::GKRule(2)) ;  /// kGAUSS21
	ROOT::Math::GSLIntegrator           ig33   ; ig33  .SetFunction(wf) ; ig33  .SetAbsTolerance(abstol) ; ig33  .SetRelTolerance(reltol) ; ig33  .SetIntegrationRule(ROOT::Math::Integration::GKRule(3)) ;  /// kGAUSS31
	ROOT::Math::GSLIntegrator           ig34   ; ig34  .SetFunction(wf) ; ig34  .SetAbsTolerance(abstol) ; ig34  .SetRelTolerance(reltol) ; ig34  .SetIntegrationRule(ROOT::Math::Integration::GKRule(4)) ;  /// kGAUSS41
	ROOT::Math::GSLIntegrator           ig35   ; ig35  .SetFunction(wf) ; ig35  .SetAbsTolerance(abstol) ; ig35  .SetRelTolerance(reltol) ; ig35  .SetIntegrationRule(ROOT::Math::Integration::GKRule(5)) ;  /// kGAUSS51
	ROOT::Math::GSLIntegrator           ig36   ; ig36  .SetFunction(wf) ; ig36  .SetAbsTolerance(abstol) ; ig36  .SetRelTolerance(reltol) ; ig36  .SetIntegrationRule(ROOT::Math::Integration::GKRule(6)) ;  /// kGAUSS61

	Double_t par[] = {f_rect->GetParameter(0), f_rect->GetParameter(1)} ;  /// parameter vector
	ROOT::Math::Functor func ( std::bind2nd ( std::ptr_fun(F_rect) , par ), 1 ) ;
	ROOT::Math::OneDimMultiFunctionAdapter<> wf_ (func) ;

	ROOT::Math::Integrator ig4  (wf_, ROOT::Math::IntegrationOneDim::kADAPTIVESINGULAR, abstol,reltol,size  ) ; 
	ROOT::Math::Integrator ig5  (wf_, ROOT::Math::IntegrationOneDim::kNONADAPTIVE     , abstol,reltol,size  ) ;
	ROOT::Math::Integrator ig61 (wf_, ROOT::Math::IntegrationOneDim::kADAPTIVE        , abstol,reltol,size,1) ;  /// rule number
	ROOT::Math::Integrator ig62 (wf_, ROOT::Math::IntegrationOneDim::kADAPTIVE        , abstol,reltol,size,2) ;
	ROOT::Math::Integrator ig63 (wf_, ROOT::Math::IntegrationOneDim::kADAPTIVE        , abstol,reltol,size,3) ;
	ROOT::Math::Integrator ig64 (wf_, ROOT::Math::IntegrationOneDim::kADAPTIVE        , abstol,reltol,size,4) ;
	ROOT::Math::Integrator ig65 (wf_, ROOT::Math::IntegrationOneDim::kADAPTIVE        , abstol,reltol,size,5) ;
	ROOT::Math::Integrator ig66 (wf_, ROOT::Math::IntegrationOneDim::kADAPTIVE        , abstol,reltol,size,6) ;
	ROOT::Math::Integrator ig7  (wf_, ROOT::Math::IntegrationOneDim::kGAUSS           , abstol,reltol,size  ) ;

	vector<Double_t> range_and_singularities ; 
		range_and_singularities.push_back(xmin) ;
		range_and_singularities.push_back(f_rect->GetParameter(0)) ;
		range_and_singularities.push_back(f_rect->GetParameter(1)) ;
		range_and_singularities.push_back(xmax) ;

	FILE* fout = fopen ("zzztest.res","a") ;
		fprintf (fout,"\n *** x-axis range = ( %+.0f ; %+.0f )  /  rect. range = ( %+.1f ; %+.1f ) *** \n",xmin,xmax,-0.5+xshift,+0.5+xshift) ;
		fprintf (fout,"\nROOT::Math::GaussIntegrator                             : %e",ig1   .Integral(xmin,xmax)) ;  /// range_and_singularities, NOT implemented : http://project-mathlibs.web.cern.ch/project-mathlibs/sw/html/classROOT_1_1Math_1_1GaussIntegrator.html#25a640c405f15558b8d567078df8b4ca
		fprintf (fout,"\nROOT::Math::GaussLegendreIntegrator       , npoints=1e3 : %e",ig21e3.Integral(xmin,xmax)) ;  /// range_and_singularities, NOT implemented : http://project-mathlibs.web.cern.ch/project-mathlibs/sw/html/classROOT_1_1Math_1_1GaussLegendreIntegrator.html#5dac4e5d8c2de14236fd8572d7719df3
		fprintf (fout,"\nROOT::Math::GaussLegendreIntegrator       , npoints=1e4 : %e",ig21e4.Integral(xmin,xmax)) ;  /// range_and_singularities, NOT implemented : http://project-mathlibs.web.cern.ch/project-mathlibs/sw/html/classROOT_1_1Math_1_1GaussLegendreIntegrator.html#5dac4e5d8c2de14236fd8572d7719df3
		fprintf (fout,"\nROOT::Math::GSLIntegrator                 , rule=1      : %e",ig31  .Integral(xmin,xmax)) ;
		fprintf (fout,"\nROOT::Math::GSLIntegrator                 , rule=2      : %e",ig32  .Integral(xmin,xmax)) ;
		fprintf (fout,"\nROOT::Math::GSLIntegrator                 , rule=3      : %e",ig33  .Integral(xmin,xmax)) ;
		fprintf (fout,"\nROOT::Math::GSLIntegrator                 , rule=4      : %e",ig34  .Integral(xmin,xmax)) ;
		fprintf (fout,"\nROOT::Math::GSLIntegrator                 , rule=5      : %e",ig35  .Integral(xmin,xmax)) ;
		fprintf (fout,"\nROOT::Math::GSLIntegrator                 , rule=6      : %e",ig36  .Integral(xmin,xmax)) ;
		fprintf (fout,"\nROOT::Math::GSLIntegrator *               , rule=1      : %e",ig31  .Integral(range_and_singularities)) ;
		fprintf (fout,"\nROOT::Math::GSLIntegrator *               , rule=2      : %e",ig32  .Integral(range_and_singularities)) ;
		fprintf (fout,"\nROOT::Math::GSLIntegrator *               , rule=3      : %e",ig33  .Integral(range_and_singularities)) ;
		fprintf (fout,"\nROOT::Math::GSLIntegrator *               , rule=4      : %e",ig34  .Integral(range_and_singularities)) ;
		fprintf (fout,"\nROOT::Math::GSLIntegrator *               , rule=5      : %e",ig35  .Integral(range_and_singularities)) ;
		fprintf (fout,"\nROOT::Math::GSLIntegrator *               , rule=6      : %e",ig36  .Integral(range_and_singularities)) ;
		fprintf (fout,"\nROOT::Math::Integrator/kADAPTIVESINGULAR                : %e",ig4   .Integral(xmin,xmax)) ;
		fprintf (fout,"\nROOT::Math::Integrator/kNONADAPTIVE                     : %e",ig5   .Integral(xmin,xmax)) ;
		fprintf (fout,"\nROOT::Math::Integrator/kADAPTIVE          , rule=1      : %e",ig61  .Integral(xmin,xmax)) ;
		fprintf (fout,"\nROOT::Math::Integrator/kADAPTIVE          , rule=2      : %e",ig62  .Integral(xmin,xmax)) ;
		fprintf (fout,"\nROOT::Math::Integrator/kADAPTIVE          , rule=3      : %e",ig63  .Integral(xmin,xmax)) ;
		fprintf (fout,"\nROOT::Math::Integrator/kADAPTIVE          , rule=4      : %e",ig64  .Integral(xmin,xmax)) ;
		fprintf (fout,"\nROOT::Math::Integrator/kADAPTIVE          , rule=5      : %e",ig65  .Integral(xmin,xmax)) ;
		fprintf (fout,"\nROOT::Math::Integrator/kADAPTIVE          , rule=6      : %e",ig66  .Integral(xmin,xmax)) ;
		fprintf (fout,"\nROOT::Math::Integrator/kGAUSS                           : %e",ig7   .Integral(xmin,xmax)) ;
		fprintf (fout,"\nROOT::Math::Integrator/kADAPTIVESINGULAR *              : %e",ig4   .Integral(range_and_singularities)) ;
		fprintf (fout,"\nROOT::Math::Integrator/kNONADAPTIVE      *              : %e",ig5   .Integral(range_and_singularities)) ;
		fprintf (fout,"\nROOT::Math::Integrator/kADAPTIVE         *, rule=1      : %e",ig61  .Integral(range_and_singularities)) ;
		fprintf (fout,"\nROOT::Math::Integrator/kADAPTIVE         *, rule=2      : %e",ig62  .Integral(range_and_singularities)) ;
		fprintf (fout,"\nROOT::Math::Integrator/kADAPTIVE         *, rule=3      : %e",ig63  .Integral(range_and_singularities)) ;
		fprintf (fout,"\nROOT::Math::Integrator/kADAPTIVE         *, rule=4      : %e",ig64  .Integral(range_and_singularities)) ;
		fprintf (fout,"\nROOT::Math::Integrator/kADAPTIVE         *, rule=5      : %e",ig65  .Integral(range_and_singularities)) ;
		fprintf (fout,"\nROOT::Math::Integrator/kADAPTIVE         *, rule=6      : %e",ig66  .Integral(range_and_singularities)) ;
		fprintf (fout,"\nROOT::Math::Integrator/kGAUSS            *              : %e",ig7   .Integral(range_and_singularities)) ;
		fprintf (fout,"\nTF1::Integral(xmin,xmax)                                : %e",f_rect->Integral(xmin,xmax)) ;
		fprintf (fout,"\nTF1::Integral(xmin,xmax,(Double_t*)(NULL),1e-18)        : %e",f_rect->Integral(xmin,xmax,(Double_t*)(NULL),1e-18)) ;
		fprintf (fout,"\n") ;
	fclose (fout) ;

	return ;
}