/**********************************************************************
  *                                                                    *
  * Copyright (c) 2003 INFN - Sezione di Napoli                        *
  *                                                                    *
  * For more information (including a list of authors) see             *
  * the README file                                                    *
  *                                                                    *
  * This library is free software; you can redistribute it and/or      *
  * modify it under the terms of the GNU General Public License        *
  * as published by the Free Software Foundation; either version 2     *
  * of the License, or (at your option) any later version.             *
  *                                                                    *
  * This library is distributed in the hope that it will be useful,    *
  * but WITHOUT ANY WARRANTY; without even the implied warranty of     *
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the       *
  * GNU General Public License for more details.                       *
  *                                                                    *
  * You should have received a copy of the GNU General Public License  *
  * along with this library (see file COPYING); if not, write to the   *
  * Free Software Foundation, Inc., 59 Temple Place, Suite 330,        *
  * Boston, MA 02111-1307 USA, or contact the authors.                 *
  *                                                                    *
  **********************************************************************/
#include "StatTools/Pdf/Ext/Defaults.h"
#include "StatTools/Fit/Ext/Defaults.h"
#include "StatTools/Funct/Argus.h"
#include "StatTools/Funct/Gaussian.h"
#include "StatTools/Funct/Parameter.h"
#include "StatTools/Funct/Function.h"
#include "StatTools/Pdf/Argus.h"
#include "StatTools/Pdf/RandomGeneratorPoissonian.h"
#include "StatTools/Pdf/RandomGeneratorGaussian.h"
#include "StatTools/Pdf/RandomGeneratorMixture.h"
#include "StatTools/Pdf/RandomGeneratorSample.h"
#include "StatTools/Fit/Chi2.h"
#include "StatTools/Fit/Chi2Fitter.h"
#include "StatTools/ToyMC/Experiment.h"
#include "TFile.h"
#include "TH1.h"

using namespace std;
using namespace Fit;
using namespace Funct;
using namespace ToyMC;

const double xmin = 5.20;
const double threshold = 5.29;
const double csi = 0.5;
const double mean = 5.28;
const double sigma = 0.002;
const double s_b   = 0.3;
const unsigned int nevents = 100000;

template<typename T> 
struct functor {
  static double f( double * x, double * ) { return (*_t)( *x ); }
  static void set( const T& t ) { _t = & t; }
private:
  static T const * _t;
};

template<typename T> 
T const * functor<T>::_t = 0;


RANDOM_GENERATOR_SAMPLE( Argus, 1024, xmin, threshold )

int main()
{
  Pdf::Argus a( xmin, threshold, csi );
  Pdf::Gaussian g( mean, sigma );  
  typedef Pdf::Mixture< Pdf::Gaussian, Pdf::Argus > Shape;
  Shape shape( g, a, s_b );
  Pdf::Poissonian signal( nevents );
  Experiment< Pdf::Poissonian, Shape > experiment( signal, shape );
  UniformPartition partition( 100, xmin, threshold ); 
  const int bins = partition.bins();
  const double binw = partition.binWidth();

  TFile file( "testfunctfit.root", "RECREATE" );
  TH1F histo( "signal", "argus + gaussian", bins, xmin, threshold );

  SampleErr<double> sample( bins );
  experiment.generate( sample, partition );
  int bin = 1;
  for( SampleErr<double>::const_iterator i = sample.begin(); i != sample.end(); i++, bin++ )
  {
    double n = i->content, dn = i->error;
    histo.SetBinContent( bin, n );
    histo.SetBinError( bin, dn );
  }
  histo.Write();

  Parameter fit_csi( csi );
  Parameter fit_threshold( threshold );
  Parameter fit_mean( mean );
  Parameter fit_sigma( sigma );
  Argus< X, Parameter, Parameter > argus( fit_threshold, fit_csi );
  Gaussian< X, Parameter, Parameter > gaussian( fit_mean, fit_sigma );  

  double area_argus = integral< X >( argus, xmin, threshold ); 
  Parameter fit_nargus( nevents * ( 1 - s_b ) * binw / area_argus );
  Parameter fit_ngaussian( nevents * s_b * binw );
  Function< X > f( fit_nargus * argus + fit_ngaussian * gaussian );
  Chi2< Function< X > > chi2( f, partition );
  Chi2Fitter< Chi2< Function< X > > > fit( chi2 );
  fit.addParameter( "nargus", fit_nargus.ptr() );
  fit.addParameter( "threshold", fit_threshold.ptr() );
  fit.addParameter( "csi", fit_csi.ptr() );
  fit.addParameter( "ngaussian", fit_ngaussian.ptr() );
  fit.addParameter( "mean", fit_mean.ptr() );
  fit.addParameter( "sigma", fit_sigma.ptr() );

  cout << "intial function : " << endl << f << endl;
  cout << "input values: " << endl;
  cout << " 0) nargus:    " << nevents * ( 1 - s_b ) * binw / area_argus << endl;
  cout << " 1) threshold: " << threshold << endl;
  cout << " 2) csi:       " << csi << endl;
  cout << " 3) ngauss:    " << nevents * s_b * binw << endl;
  cout << " 4) mean:      " << mean << endl;
  cout << " 5) sigma:     " << sigma << endl;

  double par[] = { fit_nargus, fit_threshold, fit_csi, fit_ngaussian, fit_mean, fit_sigma };
  double err[] = { 1, 1, 1, 1, 1, 1 };
  fit.fit( par, err, sample );

  cout << "fit function : " << endl << f << endl;
  cout << "fit values: " << endl;
  cout << " 0) nargus:    " << par[ 0 ] << " +/- " << err[ 0 ] << endl;
  cout << " 1) threshold: " << par[ 1 ] << " +/- " << err[ 1 ] << endl;
  cout << " 2) csi:       " << par[ 2 ] << " +/- " << err[ 2 ] << endl;
  cout << " 3) ngauss:    " << par[ 3 ] << " +/- " << err[ 3 ] << endl;
  cout << " 4) mean:      " << par[ 4 ] << " +/- " << err[ 4 ] << endl;
  cout << " 5) sigma:     " << par[ 5 ] << " +/- " << err[ 5 ] << endl;


  functor< Function< X > >::set( f );
  TF1 rootfit( "fit", functor< Function< X > >::f, xmin, threshold );
  rootfit.Write();

  return 0;
}