// $Id$
// Author: Peter Barnes 2008-10

#include "iostream"
#include "string"

#include "TF1.h"
#include "TGraphErrors.h"

using namespace std;

const char * src = "DebugFitBug";

//______________________________________________________________________________
TGraphErrors *
DrawFitBug (const char * dataf = "Test/DrawFitBug.data.txt")
{
  // Illustrate odd behavour of fits when drawing w/ or w/o Logx

  TGraphErrors *g1 = ReadGraphData(dataf);
  
  // my eyes must be going bad
  g1->SetMarkerStyle(kFullCircle);
  g1->SetMarkerSize(0.7);
  g1->SetMarkerColor(kBlue);
  g1->SetLineColor(kBlue);

  // make sure we ask for fit statistics, to see that we can get them
  gStyle->SetOptFit(1);     // default: Chisquare/DoF, parameter names, values, errors
  gStyle->SetOptTitle();

  //  Create the graphs and fit them ___
  //  *before* we have a canvas/pad
  
  //  Correct plot:
  TF1 *f1 = NewFitFunction("Fit(\"RN\"), Add");
  g1->SetTitle(f1->GetName());
  g1->Fit(f1, "RN");     // do the fit, but don't attach to graph
  g1->GetListOfFunctions()->Add(f1);  // attach by hand, to avoid plotting bug

  //  Obvious, but incorrect plot
  TF1 *f2 = NewFitFunction("Fit(\"R\")");
  TGraphErrors *g2 = new TGraphErrors(*g1);
  g2->SetTitle(f2->GetName());
  g2->Fit(f2, "R");  // do the fit, attaching to graph

  //  Should fix the error, but doesn't
  TF1 *f3 = NewFitFunction("Fit(\"R\"), Update()");
  TGraphErrors *g3 = new TGraphErrors(*g1);
  g3->SetTitle(f3->GetName());
  g3->Fit(f3, "R");  // do the fit, attaching to graph
  g3->GetFunction(f3->GetName())->Update();


  // Now we make a pad and draw ________
  
  TCanvas *c1 = new TCanvas("DrawFitBug", "DrawFitBug", 1200, 400);
  c1->Divide(3,1);

  c1->cd(1);  gPad->SetLogx();  g1->Draw("AP");
  c1->cd(2);  gPad->SetLogx();  g2->Draw("AP");
  c1->cd(3);  gPad->SetLogx();  g3->Draw("AP");

  return g1;
}

//______________________________________________________________________________
TGraphErrors *
ReadGraphData(const char * dataf)
{
  // Read ascii points and error from dataf

  ifstream in;
  in.open(dataf);

  int index;
  double gate, mom, err;

  // discard 5 header lines
  cout<<src<<": reading data file:  header"<<std::flush;
  string line;
  for (index=0; index<5; index++) {
	getline(in, line);
	//cout<<line<<endl;
  }

  TGraphErrors *g = new TGraphErrors();

  // read graph data
  index = -1;
  while (in >> index >> gate >> mom >> err ) {
	g->SetPoint     (index, gate, mom);
	//g->SetPointError(index,    0, err);
	g->SetPointError(index,    0, 0);
	if ( index && (index % 30) ) cout<<"."<<std::flush;
  }
  cout<<"= "<<++index<<" points."<<endl;

  return g;
}

//______________________________________________________________________________
Double_t
Fit_Y2F(Double_t *x, Double_t *par)
{
  double t   = x[0];
  double R2F = par[0];
  double tau = par[1];
  double R0  = par[2];

  double y2f = R2F * (1 - (1 - exp( -1 * t / tau) ) / ( t / tau)) + R0;

  return y2f;
}

//______________________________________________________________________________
TF1 *
NewFitFunction(const char * name)
{
  TF1 *f = new TF1(name, &Fit_Y2F, 0,0,3);

  const double maxGate = 2E6;

  // domain range
  f->SetRange(1, maxGate);

  // parameter names and ranges
  f->SetParName  (0, "R_{2F}");
  f->SetParLimits(0, -10, 10);  // let's not get carried away
  f->SetParName  (1, "#lambda^{-1} (ns)");
  f->SetParLimits(1, 0.1, maxGate);
  f->SetParName  (2, "R^{0}_{2F}");
  f->SetParLimits(2, -10, 10);

  // parameter seeds
  f->SetParameter(0, 0.1);  
  f->SetParameter(1, maxGate);
  f->SetParameter(2, 0);
  
  // plotting attributes
  f->SetMarkerSize(0);
  f->SetLineWidth(2);
  f->SetLineColor(kGreen);

  return f;
}