#include <TGraphErrors.h>
#include <TFitter.h>
#include <TF1.h>
#include <iostream.h>
#include "Foption.h"

extern TVirtualFitter *grFitter;
extern TF1* grF1;
extern Foption_t fitOption;

void GraphFitChisquareNoEX(Int_t &npar, Double_t * /*gin*/, Double_t &f,
			   Double_t *u, Int_t /*flag*/);
void testfit()
{
  TGraph *graph;
  Float_t x[30],y[30],dx[30],dy[30];
  int i;
  for(i=0;i<30;i++){x[i]=i*.5;y[i]=3*x[i]*x[i]*+x[i]-3;dx[i]=0.1;dy[i]=0.01*y[i];}
  graph=new TGraphErrors(30,x,y,dx,dy);
  grFitter=TVirtualFitter::Fitter(graph);
  grFitter->SetFCN(GraphFitChisquareNoEX);
  graph->Fit("pol3","U");
}

void GraphFitChisquareNoEX(Int_t &npar, Double_t * /*gin*/, Double_t &f,
                       Double_t *u, Int_t /*flag*/)
{
//*-*-*-*-*-*Minimization function for Graphs using a Chisquare method*-*-*-*-*
//*-*        =========================================================
//
// In case of a TGraphErrors object, ex, the error along x,  is projected
// along the y-direction by calculating the function at the points x-ex and
// x+ex.
//
// The chisquare is computed as the sum of the quantity below at each point:
//
//                     (y - f(x))**2
//         -----------------------------------
//         ey**2 + ((f(x+ex) - f(x-ex))/2)**2
//
// where x and y are the point coordinates

   Double_t cu,eu,ex,ey,eux,fu,fsum,fm,fp;
   Double_t x[1], xx[1];
   Double_t xm,xp;
   Int_t bin, npfits=0;

   TGraph *gr     = (TGraph*)grFitter->GetObjectFit();
   Int_t n        = gr->GetN();
   Double_t *gx   = gr->GetX();
   Double_t *gy   = gr->GetY();
   Double_t fxmin = grF1->GetXmin();
   Double_t fxmax = grF1->GetXmax();
   npar           = grF1->GetNpar();

   grF1->InitArgs(x,u);
   f      = 0;
   // cout<<"Using New Fitting function"<<endl;
   // cout<<"Fitter name "<<grFitter->GetName()<<endl;
   // cout<<"Graph name "<<gr->GetName()<<endl;
   // cout<<"Function Name "<<grF1->GetName()<<endl;
   grF1->Print();
   //   cout<<"u="<<u<<endl;
   //   cout<<"u[0]="<<u[0]<<endl;
   for (bin=0;bin<n;bin++) {
      x[0] = gx[bin];
      if (!grF1->IsInside(x)) continue;
      cu   = gy[bin];
      TF1::RejectPoint(kFALSE);
      cout<<"About to evaluate the function"<<endl;
      fu   = grF1->EvalPar(x,u);
      cout<<"Evaluated function"<<endl;
      if (TF1::RejectedPoint()) continue;
      fsum = (cu-fu);
      npfits++;
      if (fitOption.W1) {
         f += fsum*fsum;
         continue;
      }
      ex  = gr->GetErrorX(bin);
      ey  = gr->GetErrorY(bin);
      if (ex < 0) ex = 0;
      if (ey < 0) ey = 0;
      if (ex >= 0) {
        xm = x[0] - ex; if (xm < fxmin) xm = fxmin;
        xp = x[0] + ex; if (xp > fxmax) xp = fxmax;
        xx[0] = xm; fm = grF1->EvalPar(xx,u);
        xx[0] = xp; fp = grF1->EvalPar(xx,u);
        eux = 0.5*(fp-fm);
      } else
        eux = 0.;
      eu = ey*ey+eux*eux;
      if (eu <= 0) eu = 1;
      f += fsum*fsum/eu;
   }
   grF1->SetNumberFitPoints(npfits);
}