/**********************************************************************
  *                                                                    *
  * 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/Fit/Chi2.h"
#include "StatTools/Fit/Chi2Fitter.h"
#include "StatTools/Pdf/RandomGeneratorPoissonian.h"
#include "StatTools/Pdf/RandomGeneratorTriangular.h"
#include "StatTools/ToyMC/Experiment.h"
#include "StatTools/Funct/Line.h"
#include "StatTools/Funct/Parabola.h"
#include "TFile.h"
#include "TH1.h"
#include <vector>
#include <algorithm>
#include <iterator>
#include <functional>
#include <numeric>
using namespace std;
using namespace Funct;
using namespace Fit;
using namespace Pdf;
using namespace ToyMC;

int main()
{
  Triangular triangle;
  RandomGenerator<Triangular> r( triangle );

  Poissonian signal( 1000 );
  Experiment<Poissonian, Triangular> experiment( signal, triangle );
  UniformPartition partition( 10, 0.0, 1.0 ); 
  const int bins = partition.bins();

  Line line( 0, 1 );
  Chi2<Line> chi2line( line, partition );
  Chi2Fitter<Chi2<Line> > fitter1( chi2line );
  fitter1.addParameter( "a", &line.a );
  fitter1.addParameter( "b", &line.b );

  Parabola para( 0, 1, 0 );
  Chi2<Parabola> chi2para( para, partition );
  Chi2Fitter<Chi2<Parabola> > fitter2( chi2para );
  fitter2.addParameter( "a", &para.a );
  fitter2.addParameter( "b", &para.b );
  fitter2.addParameter( "c", &para.c );

  TFile file( "testchi2fitter.root", "RECREATE" );
  TH1F pull1( "pulla_line", "pull intercept", 50, -5, 5 );
  TH1F pull2( "pullb_line", "pull slope", 50, -5, 5 );

  TH1F pull3( "pulla_parab", "pull intercept", 50, -5, 5 );
  TH1F pull4( "pullb_parab", "pull slope", 50, -5, 5 );
  TH1F pull5( "pullc_parab", "pull quadratic", 50, -5, 5 );

  SampleErr<double> histo( bins );
  for ( int i = 0; i < 1000; i++ )
    {
      double norm = 2 * experiment.generate( histo, partition ) / bins;
      histo.scale ( 1 / norm );
      double par1[] = { 0, 1 }, err1[] = { 1, 1 };
      fitter1.fit( par1, err1, histo.begin(), histo.end() );
      pull1.Fill( ( par1[ 0 ] - 0 ) / err1[ 0 ] );
      pull2.Fill( ( par1[ 1 ] - 1 ) / err1[ 1 ] );

      double par2[] = { 0, 1, 0 }, err2[] = { 1, 1, 1 };
      fitter2.fit( par1, err2, histo.begin(), histo.end() );
      pull3.Fill( ( par2[ 0 ] - 0 ) / err2[ 0 ] );
      pull4.Fill( ( par2[ 1 ] - 1 ) / err2[ 1 ] );
      pull5.Fill( ( par2[ 2 ] - 0 ) / err2[ 2 ] );
    }
  pull1.Write();
  pull2.Write();
  pull3.Write();
  pull4.Write();
  pull5.Write();
  return 0;
}