//
// Author: E.Gorini
// last modification 28/1/2014
//
//   to execute:
//   [0] .x GenFit.C("Brewster.dat"); 
// or, to not plot pulls
//   [0] .x GenFit.C("Brewster.dat",false)
//
#include <TH1>
#include <TF1>
#include <TLegend>
#include <TCanvas>
#include <Riostream>
#include <vector>

void GenFit(string file, bool makepulls=true){
  // reinitialization
  gROOT->Reset();
  bool makestat=true;
  
  // here open the data file
  ifstream in;
  in.open(file.c_str());
// uscire se file non esiste. 
  if(in.is_open()==0){
     in.close();
     gROOT->Reset();
    cout << "File " << file <<" do not exist, exit";
    break;
  } 
  // reads plot title and X and Y labels
  string Title,Xtitle,Ytitle;
  getline(in,Title);
//  in >> Title; 
  getline(in,Xtitle); 
  getline(in,Ytitle);
  cout << "Plot Title is: \""<< Title <<"\", ";
  cout << "x axis title is: \""<< Xtitle <<"\", ";
  cout << "y axis title is: \""<< Ytitle <<"\".\n";
  // reads function
  string Function;
  getline(in,Function);
  cout << "Function to fit is: "<< Function << "\n";
  // read parameter number and names
  int numpar; 
  in >> numpar;
  vector<string> parnames;
  for (int np=0;np<numpar;np++){
    string sbuf;
    in >> sbuf;
    parnames.push_back(sbuf);
  }
  cout << "Number of Parameters is "<< parnames.size() << " with value: ";
  // read parameter initial values
  TVector pvalue(numpar);
  for (int np=0;np<numpar;np++)in >> pvalue(np);
  for (int np=0;np<numpar;np++){
    cout << parnames[np]<< "="<< pvalue[np] << ", ";  
  }
  cout << "\n";
  // create the vectors, limited at 100 for the points, 20 for the parameters
  // vectors will be anyway resized if necessary.....
  Int_t n = 100, m = 20;
  TVectorF xx(n),yy(n),ex(n),ey(n),pull(n);
  Float_t xxx,yyy,exx,eyy;
  bool flag;
  Int_t nlines=0;
  // now reads data x,y, ex,ey 
  while (in.good()) {
    in >> flag >> xxx >> yyy >> exx >> eyy;
    if(flag){
      xx(nlines)=xxx;
      yy(nlines)=yyy;
      ex(nlines)=exx;
      ey(nlines)=eyy;
      nlines++;
      flag=false;
    } 
  }
  cout << "found " << nlines << " data points, ";
  //resize here TVectorF
  xx.ResizeTo(nlines);
  yy.ResizeTo(nlines);
  ex.ResizeTo(nlines);
  ey.ResizeTo(nlines);
  //   xx.Print();
  in.close();

  if (makestat) {
	gStyle->SetOptStat(1);
	gStyle->SetOptFit(111111);
  }
    
  // create and set canvas style
  TCanvas *c1 = new TCanvas("c1",Title.c_str(),10,10,900,700);
  c1->SetFillColor(19);
  c1->SetGrid(1,1); 
  if(makepulls)c1->Divide(1,2);
  c1->cd(1);  
  
  // plot the measurement 
  gr = new TGraphErrors(xx,yy,ex,ey);
  gr->SetTitle(Title.c_str());
  gr->SetMarkerColor(kBlue);
  gr->SetMarkerStyle(21);
  gr->GetXaxis()->SetTitle(Xtitle.c_str());
  gr->GetYaxis()->SetTitle(Ytitle.c_str());
  gr->Draw("AP");
  c1->Update();
  
  // create the function
  TF1 *fff = new TF1("fff",Function.c_str(),1,1000);
  
  // set parameter names and initial values
  for (int np=0;np<numpar;np++){
    fff->SetParName(np,parnames[np].c_str());
    fff->SetParameter(np,pvalue(np));
  }	
  
  // fit   
  cout << "Start Fitting.....\n"; 
  gr->Fit("fff");
  

	if(makepulls){	
// make plot of pulls
	   c1->cd(2);
	   c1->SetGrid(1,1);
       for (int bin=0;bin<nlines;bin++)pull(bin) = (yy(bin)-fff->Eval(xx(bin),0.0,0.0))/ey(bin);
		re = new TGraph(xx,pull);
		re->SetTitle(Title.c_str());
		re->SetMarkerColor(kBlue);
		re->SetMarkerStyle(21);
		re->SetMaximum(5.0);
		re->SetMinimum(-5.0);
		re->GetXaxis()->SetTitle(Xtitle.c_str());
        re->GetYaxis()->SetTitle("Pulls");
		re->Draw("AP");
		// draw a line 
		tl= new TLine(xx(0),0.0,xx(nlines-1),0.0);
		tl->SetLineColor(kRed);
		tl->SetLineWidth(2);
		tl->SetLineStyle(2);
		tl->Draw();
		c1->Update();
	}	
	
  
  // save the plots in pdf
  
  string TitleOut=Title+".pdf";
  c1->Print(TitleOut.c_str());  
}