#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <iostream> 
#include<float.h> 
#include "TRandom.h"
#include "TMath.h"
#include "Math/SpecFunc.h"
#include "TTree.h"
#include "TFile.h"
#include "TNtupleD.h"
#include "TH1.h"
#include "TVirtualFitter.h"
#include <iostream>
#include <fstream>
#include "Riostream.h"
#include "TROOT.h"
#include "TTree.h"
#include "TGraph.h"
#include "TMinuit.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TLine.h"
#include "TPie.h"
#include "TPavesText.h"
#include "TGraphErrors.h"
#include "TStyle.h"
#include "TSystem.h"
#include "TAxis.h"
#include "TH1.h"
#include "TLatex.h"
#include <stdio.h>
#include "math.h"
#include "TMatrixT.h"
#include "TH1F.h"
#include "TF1.h"
#include "TMath.h"
#include "Math/IFunction.h"
//#include "T.h"

#define DEBUG 0

//These data are made global, so that they can be accessed by
//Minuit functions..
double x[500];
double y[500];
double dx[500];
double dy[500];

typedef struct {
  double egmean;
  double Wmean;
  double ctmean;
  double crosslambda;
  double errorlambda;
  double crosssigma;
  double errorsigma;
} DataSet;

double eps1=1e-6;
double eps=1e-10;	// A very small number
const int lmax=9;   // order of polynomials

//legendre polynomial
Double_t myfitLegendre(Double_t *x,Double_t *par){
  Float_t xx =x[0];
  Double_t f = par[0];
  if (xx >=-1 && xx<=1){
    for ( int i=0 ; i<=lmax ; i++ ) {
      f += par[i]*ROOT::Math::assoc_legendre(i,0,xx);
    }
  }
  return f;
}


int test() {



  gSystem->Load("libMathMore");



  // This object will hold the full set of data...
  DataSet myDataSet[1500];
  
  // Declare some temporary variables so that the data can be read
  // in... 
  double egmean;
  double Wmean;
  double ctmean;
  double crosslambda,errorlambda;
  double crosssigma,errorsigma;
  
  // Open up the file containing all the data...
  ifstream datafile("data.txt"); 
  int index=0;
  // Read this file until the end is reached. 
  if (datafile.is_open()) {
    while (! datafile.eof() ) {
      // Each line is read into individual variables...
      datafile >> egmean >> Wmean >> ctmean
	       >> crosslambda >> errorlambda
	       >> crosssigma >> errorsigma;
      // If the record contains a cross section then put the
      // correct numbers into the appropriate RooRealVar...
      if (crosslambda!=0.0) {
	myDataSet[index].egmean = egmean;
	myDataSet[index].Wmean  = Wmean;
	myDataSet[index].ctmean = ctmean;
	myDataSet[index].crosslambda  = crosslambda;
	myDataSet[index].errorlambda  = errorlambda;
	myDataSet[index].crosssigma  = crosssigma;
	myDataSet[index].errorsigma  = errorsigma;
	index++;
      }
    }
    datafile.close(); //	We're finished with the text file
  }
  // Just a check that the data file was opened correctly... 
  else cout << "Unable to open file"; 
  
  // Print out a summary of the data taken in so far..
  int nDataPoints = index;
  cout << "Found " << nDataPoints << " valid data points" << endl; 
     
  TCanvas *c1 = new TCanvas("c1","c1",700,1200);
  c1->SetFillStyle(0);
  c1->SetFillColor(0);
  c1->SetFrameFillStyle(0);
  gPad->SetRightMargin(0.04);
  gPad->SetLeftMargin(0.2);
  gPad->SetTopMargin(0.06);
  gPad->SetBottomMargin(0.3);
  c1->SetBorderMode(0);
  c1->SetBorderSize(0.1);
  c1->SetFrameBorderMode(8);
  c1->Divide(8,10,0,0);

  double egamma_loop = 0.;
  //  loop for each gamma energy
  const int N_energies = int ((2.950-0.944)/0.025);
  for ( Int_t energy=0 ; energy<= N_energies; ++energy) {
    std::cout << " === Processing energy " << energy << " === " << std::endl;
    if (energy<=14){
      egamma_loop = (Double_t)energy * 0.025 + 0.944;
    }
    if (energy>=15){
      egamma_loop = (Double_t)energy * 0.025 + 0.945;
    }
    if  (energy>=19){
      egamma_loop = (Double_t)energy * 0.025 + 0.946;
    }
    if (energy>=28){
      egamma_loop = (Double_t)energy * 0.025 + 0.947;
    }
    if (energy>=30){
      egamma_loop = (Double_t)energy * 0.025 + 0.948;
    }
    if (energy>=33){
      egamma_loop = (Double_t)energy * 0.025 + 0.949;
    }
    if (energy>=37){
      egamma_loop = (Double_t)energy * 0.025 + 0.950;
    }
    if (energy>=42){
      egamma_loop = (Double_t)energy * 0.025 + 0.951;
    }
    if (energy>=49){
      egamma_loop = (Double_t)energy * 0.025 + 0.952;
    }
    if (energy>=52){
      egamma_loop = (Double_t)energy * 0.025 + 0.953;
    }
    if(energy==57){
      continue;}
    if(energy==58){ 
      continue;}
    if (energy>=59){
      egamma_loop = (Double_t)energy * 0.025 + 0.955;
    }
    if (energy>=70){
      egamma_loop = (Double_t)energy * 0.025 + 0.954;
    }
    if (energy>=73){
      egamma_loop = (Double_t)energy * 0.025 + 0.953;
    }
    if (energy>=75){
      egamma_loop = (Double_t)energy * 0.025 + 0.952;
    }
    if (energy>=77){
      egamma_loop = (Double_t)energy * 0.025 + 0.951;
    }
    if (energy>=79){
      egamma_loop = (Double_t)energy * 0.025 + 0.950;
    }
    
    //  An example of selecting a subset of the data. Here we choose
    //  cross sections for cos(theta) = 0.05. Note the way that the
    //  comparison is done, to avoid problems with rounding errors. 
    int nSelected = 0;
    for ( int i=0 ; i<nDataPoints ; i++ ) {
      if ( fabs(myDataSet[i].egmean - egamma_loop)< eps ) {
	x[nSelected]  = myDataSet[i].ctmean;
	dx[nSelected] = 0.0;
	y[nSelected]  = myDataSet[i].crosslambda;
	dy[nSelected] = myDataSet[i].errorlambda;
	nSelected++;
	cout << i << "\t" 
	     << myDataSet[i].ctmean << "\t"
	     << myDataSet[i].egmean << "\t"
	     << myDataSet[i].crosslambda << "\t"
	     << myDataSet[i].errorlambda << endl;
      }
    }
    
    c1->cd(energy+1);
    // put results of selection into TGraphErrors..
    TGraphErrors *gr = new TGraphErrors(nSelected,x,y,dx,dy);
    gr -> SetName(Form("graph_%i", energy));
    gr->SetTitle("");
    gr->SetMarkerColor(4);
    gr->SetMarkerStyle(20);
    gr->SetMarkerSize(1.5);
    gr->SetFillColor(0);
    gr->SetMarkerColor(kBlue);
    gr->SetLineWidth(1);
    gr->GetXaxis()->SetTitle("cos(#theta^{K+}_{CM})");
    gr->GetXaxis()->SetRangeUser(-1,0.95);
    gr->GetYaxis()->SetTitle("d#sigma/d#Omega[#mub/sr]");
    gr->GetXaxis()->CenterTitle();
    gr->GetYaxis()->CenterTitle();
    gr->GetXaxis()->SetTitleSize(0.15);
    gr->GetXaxis()->SetTitleOffset(0.9);
    gr->GetXaxis()->SetLabelSize(0.13);
    gr->GetXaxis()->SetLabelOffset(0.0005);
    gr->GetYaxis()->SetTitleSize(0.15);
    gr->GetYaxis()->SetTitleOffset(0.8);
    gr->GetYaxis()->SetLabelOffset(-0.0009);
    gr->GetYaxis()->SetLabelSize(0.13);
    gr->Draw("AP");
    
    // Fit Associated legendre Polynomial 
    for ( Int_t nOrder=1 ; nOrder<=lmax ; nOrder++ ) {
      std::cout << "Processing Order " << nOrder << std::endl;
      TF1 *f1 = new TF1(Form("f1_%i", lmax),myfitLegendre,-1,1,nOrder);
      double parmin = -30;
      double parmax = 30;
      for (int jj = 0; jj < nOrder; ++jj) {
	f1->SetParameter(jj, 1.);
	f1->SetParLimits(jj,parmin,parmax);
      }
      gr->Fit(Form("f1_%i", lmax));
      // Get the fitter and extract the covariance matrix from it into a  TMatrix..
      TVirtualFitter::SetPrecision(eps1); 
      TVirtualFitter *fit = TVirtualFitter::GetFitter();
      double my_chi2 = f1 -> GetChisquare();
      std::cout << "Fit chi2: " << my_chi2 << std::endl;
      /*
      if (!fit->GetCovarianceMatrix()) {
	std::cout << "WARNING: Fit did not converge, trying to refit..." << std::endl;
	gr->Fit(Form("f1_%i", lmax));
	my_chi2 = f1 -> GetChisquare();
	std::cout << "Fit chi2: " << my_chi2 << std::endl;
      }
      */
      if (my_chi2 > 0. && fit->GetCovarianceMatrix()) {
	TMatrixT<double>matrix(nOrder,nOrder);
	for (int lines=0 ; lines<nOrder ; lines++ ) {
	  for (int columns=0 ; columns<nOrder ; columns++ ) {
	    std::cout << "Trying to get element (" << lines << ", " << columns << ")" << std::endl; 
	    matrix(lines,columns) = fit->GetCovarianceMatrixElement(lines,columns);
	  } // columns
	} // lines
	matrix.Print();
      } else {
	std::cerr << "ERROR: fit did not converge! for l=" << lmax << endl;
      }
    } // legendre polynoms
  }
  c1 -> Print(Form("%s.eps", c1 -> GetName()));
}