#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <sstream>
#include <string>
#define NDIM 100

void fitDataFromFile()
{

  std::string fileNameString = "datiFLin.txt";
  ifstream inputFile; 

  inputFile.open(fileNameString.c_str());
  if (!inputFile.good())
    {
      cout<<"Error while opening file <"<<fileNameString<<">"<<endl;
      return;
    }
  else cout<<"File  <"<<fileNameString<<"> open in reading mode"<<endl;

  int Nlines=0;
  double a1, a2, a3, a4;
  double xmeas[NDIM] = {0};
  double ymeas[NDIM] = {0};
  double exmeas[NDIM] = {0};
  double eymeas[NDIM] = {0};
  std::string aLine;
  while (! inputFile.eof() ) 
    {
      getline(inputFile, aLine);
      if (aLine.size()==0)
	{
	  //cout<<"Empty string; skipping"<<endl; /// handling EOF emty line 
	  continue;
	}
      //std::cout<<Nlines<< " aline = <"<<aLine<<">"<<endl;
      istringstream aLineStream(aLine);
      aLineStream >> a1 >> a2 >> a3 >>a4;
      if (Nlines==NDIM)
	{
	  cout<<"More than "<<NDIM<<" lines of data found in the file; stopping here all subsequent data will be ignored"<<endl;
	  break;
	}
      cout<<"Reading from file  x = "<<a1<<" y = "<<a2<<" ex = "<<a3<<" ey = "<<a4<<endl;
      xmeas[Nlines] = a1;
      ymeas[Nlines] = a2;
      exmeas[Nlines] = a3;
      eymeas[Nlines] = a4;
      ++Nlines;
    }
  if (Nlines<=NDIM) cout<<"End of file or maximum size of data arrays reached; N = "<<Nlines<<" lines of data available"<<endl;
  //for (int j=0; j<NDIM; ++j) cout<<" j, ey "<<j<<" "<<eymeas[j]<<endl;
  inputFile.close();


  TF1 * fun= new TF1("fun","[0]+[1]*x",0,20);
  fun->SetParameters(1.,4.);
  /// now you can manipulate teh data held in the arrays 
  TGraphErrors *gr= new TGraphErrors(Nlines,xmeas,ymeas,exmeas,eymeas);
  gr->Fit("fun");

  TVirtualFitter * fitter = TVirtualFitter::GetFitter();

  double p1=fitter->GetParameter(1);
  double p0=fitter->GetParameter(0);
  // estraggo e stampo la matrice di covarianza
  double * covMatrix = fitter->GetCovarianceMatrix();
  printf("\n p1=%f p0=%f\n",p1,p0);
  printf("\n c11=%f   c12=%f ",covMatrix[0],covMatrix[1]);
  printf("\n c21=%f   c22=%f\n ",covMatrix[2],covMatrix[3]);	
  printf("\n p0+/-errp0 = %f +/- %f",p0,sqrt(covMatrix[0]));
  printf("\n p1+/-errp1 = %f +/- %f",p1,sqrt(covMatrix[3]));
  double chi2 = fun->GetChisquare();
  double pchi2 = fun->GetProb();
   printf("\n\n chi^2=%f Prob(chi^2)=%f\n\n",chi2,pchi2);
  
  
}