#include "TMath.h"
#include "TMatrixD.h"
#include "TMatrixDLazy.h"
#include "TVectorD.h"
#include "TDecompLU.h"
#include "TDecompSVD.h"
#include "TMinuit.h"

TH1F *histo_cum = 0;



double pol7_model(Double_t *x, Double_t *par)
{
  return par[0] + par[1]*x + par[2]*pow(x,2) + par[3]*pow(x,3) + par[4]*pow(x,4) + par[5]*pow(x,5) + par[6]*pow(x,6) +par[7]*pow(x,7);

}
  

void fcn(Int_t &npar, Double_t *gin, Double_t &value, Double_t *par, Int_t iflag)
{


  printf("CHECK FCN\n");
  std::cout << histo_cum << std::endl;
  value = 0. ;

  Int_t nPoints = 56;
  Double_t *residual=new Double_t[nPoints];


  

  for (int i = 1; i < 56;i++)
  {
    double x = histo_cum->GetBinCenter(i);

    double measure = histo_cum->GetBinContent(i);

    double func = pol7_model(x,par);

    residual[i] = func - measure ;
  }


  for (Int_t i=0;i<nPoints;i++) 
  {
    for (Int_t j=0;j<nPoints;j++) 
    {
      value += residual[i]*residual[j]*mat[i][j];
    }
  }
  delete[] residual;
}


using namespace std;

void script_matr(char *fname1)
{

  gROOT->SetStyle("Plain");
  gStyle->SetOptTitle(1);
  gStyle->SetOptStat(111011);
  //gStyle->SetOptStat(0);
  gStyle->SetOptFit(1111);


  TFile *f1 = 0;
 
  f1 = new TFile(fname1);

   

  TH1F *histo = (TH1F*)f1->Get("hFts_S1_Ch1_L1");
  


  Int_t entries_before = histo->Integral();

 printf("Number of entries before range cut = %d\n",entries_before);

  //h1->SetAxisRange(0.,325.,"X");
  Int_t entries1 = histo->Integral();

  printf("Number of entries after range cut = %d\n",entries1);
  
  TH1F* histo_cum = histo->Clone("r-t-relation");

  //Int_t entries1 = h1_sub->Integral();

  Int_t entries1 = histo->Integral();

  printf("After cut entries = %d\n", entries1);
  //int n_bins_1 = h1_sub->GetSize();
  int n_bins_h1 = histo->GetSize();
  printf("Amount of bins in h1 = %d\n",n_bins_h1);


  Double_t evNum1=0.;

  for (int i = 1; i <= n_bins_h1; i++)
  {
    printf("evNum1 = %f\n",evNum1);
    evNum1+=histo->GetBinContent(i);   
    histo_cum ->SetBinContent(i,evNum1);
    histo_cum ->SetBinError(i,TMath::Sqrt(evNum1));
  }
  histo_cum->Scale((1.48/entries1));

  std::cout << histo_cum << std::endl;

  
 


  TMatrixD mat(56,56);

  TMatrixD B (56,56);
 



  for (int i2 = 0; i2 < 56; i2++)
  {
    for (int j2 = i2; j2 < 56; j2++)
    {
      mat[i2][j2] = histo_cum->GetBinContent(i2+1);
      //printf("Matrix : mat[%d][%d] = %f\n",i2,j2,mat[i2][j2]);
    }
  }

  for (int i3 = 1; i3 < 56; i3++)
  {
    for (int j3 = 0; j3 < i3; j3++)
    {
      mat[i3][j3] = histo_cum->GetBinContent(j3+1);
      //printf("Matrix : mat[%d][%d] = %f\n",i2,j2,mat[i2][j2]);
    }
  }




  
  B = mat;

 const Double_t det = mat.Determinant();
  cout << "det: " << det <<endl;
  
  TDecompSVD svd(mat);
  const TVectorD singularVal = svd.GetSig();


  svd.SetTol(1.0e-38);
  singularVal.Print();



  svd.Invert(mat);

 // mat.Print();
 // B.Print();
 // (mat*B).Print();
  //(B*mat).Print();


  printf("CHECK1\n");

  double step[8];
  int ierflg = 0;

  Double_t arglist[20];

  TMinuit *minimization = new TMinuit(8);

  minimization->SetFCN(fcn); 


  for (int i = 0; i < 8; i++)
	  step[i] = 0.0;

  minimization->mnparm(0,	"p0",	0.109,	0.01, 0, 0, ierflg);
	minimization->mnparm(1,	"p1",	0.007026,		0.0, 0, 0, ierflg);
	minimization->mnparm(2,	"p2",	0.0001906,		0.0, 0, 0, ierflg);
	minimization->mnparm(3,	"p3",-3.032e-6,	0.0, 0, 0, ierflg);
	minimization->mnparm(4,	"p4", 1.953e-8,		0.0, 0, 0, ierflg);
	minimization->mnparm(5,	"p5", -6.574e-11,	0.0, 0, 0, ierflg);
	minimization->mnparm(6,	"p6", 1.141e-13,		0.0, 0, 0, ierflg);
	minimization->mnparm(7,	"p7", -8.069e-17,		0.0, 0, 0, ierflg);
	
  

  printf("step1\n");

  arglist[0] = 10000;    // <<< Can't remember what is this.
	arglist[1] = 0.1;      // <<< Can't remember what is this.

  printf("step2\n");
 
	minimization->SetMaxIterations(10000);

   printf("step3\n")	;

	minimization->mnexcm("MINIMIZE", arglist, 1, ierflg);

  printf("step4\n")	;
	
	Double_t vfin[8], vfin_err[8];
	for (int i = 0; i < 8; i++)
  {
	  minimization->GetParameter(i, vfin[i], vfin_err[i]);
    printf("vfin[%d] = %f, vfin_err[%d] = %f\n",i,vfin[i],i, vfin_err[i]);
  }




  histo->SetName("Drift time distribution with constant backroung cut");
  histo->SetTitle("Drift time distribution with constant backroung cut");


  histo_cum->SetName("r-t-relation");
  histo_cum->SetTitle("r-t-relation");
  




   TFile *fk = new TFile("test_cailbration.root","RECREATE");

   
   //TCanvas *c1 = new TCanvas("c1","c1",800,1000);
   
   //c1->SetGrid();
   //c1->Divide(2,1);
   
   //c1->cd(1)->SetGrid();
   //histo->Draw();
   //histo->Write();
  

   //c1->cd(2)->SetGrid();
   //h1_sub->Draw();
   
   //c1->cd(2)->SetGrid();
   histo_cum->Draw();
   histo_cum->Write();
   

   //c1->Write();
       

   

   
 
  fk->Close();

}