#include "Fit/BinData.h"
#include "Fit/Chi2FCN.h"
#include "Math/WrappedMultiTF1.h"
#include "HFitInterface.h"
#include "TCanvas.h"
#include "TStyle.h"
#include <vector>
#include <iterator>
#include <boost/shared_ptr.hpp>
#include <boost/make_shared.hpp>
#include "TFile.h"
#include "TGraphErrors.h"
#include "TF1.h"
#include "Experiment.h"

using namespace std;

int main() { 
   TFile f("graphs.root");
  char Tfilename[50];
  int exp_n = 25;
  double T_max = 10.0;
  double T_min = 0.0;
  int n_min = 1000;
  int n_runs = 1001;
  
   
   //set the parameters to be passed!!!
 double par[34], par_init[34];
 
 // evolution parameters omega T_2
 par[0] = 0.4;        par_init[0] = 0.4;
 par[1] = 0.05;       par_init[1] = 0.05;
 
 //measurement noise eps_1 eps_2
 par[2] = 0.025;      par_init[2] = 0.0;
 par[3] = 0.075;      par_init[3] = 0.0;
 
 //x_state
  par[4] = 0.96;     par_init[4] = 1.0;
  par[5] = 0.0;      par_init[5] = 0.0;
  par[6] = 0.17;     par_init[6] = 0.0;
 
 //y-state
  par[7] = 0.09;     par_init[7] = 0.0;
  par[8] = 0.96;     par_init[8] = 1.0;
  par[9] = 0.09;     par_init[9] = 0.0;
 
 //z-state
  par[10] = 0.0;     par_init[10] = 0.0;
  par[11] = 0.0;     par_init[11] = 0.0;
  par[12] = 1.0;     par_init[12] = 1.0;
 
 //-x-state
  par[13] = -0.96;   par_init[13] = -1.0;
  par[14] = 0.09;    par_init[14] = 0.0;
  par[15] = 0.09;    par_init[15] = 0.0;
 
 //-z-state
  par[16] = 0.09;    par_init[16] = 0.0;
  par[17] = 0.09;    par_init[17] = 0.0;
  par[18] = -0.96;   par_init[18] = -1.0;
 
 //x-measurement 
  par[19] = 0.96;    par_init[19] = 1.0;
  par[20] = 0.08;    par_init[20] = 0.0;
  par[21] = 0.17;    par_init[21] = 0.0;
 
 //y-measurement
  par[22] = 0.08;    par_init[22] = 0.0;
  par[23] = 0.96;    par_init[23] = 1.0;
  par[24] = 0.17;    par_init[24] = 0.0;
 
 //z-measurement
  par[25] = 0.08;    par_init[25] = 0.0;
  par[26] = 0.96;    par_init[26] = 1.0;
  par[27] = 0.17;    par_init[27] = 0.0;
 
 //-x-measurement
  par[28] = -0.96;   par_init[28] = -1.0;
  par[29] = 0.09;    par_init[29] = 0.0;
  par[30] = 0.09;    par_init[30] = 0.0;
 
 //-z-measurement
  par[31] = 0.09;    par_init[31] = 0.0;
  par[32] = 0.09;    par_init[32] = 0.0;
  par[33] = -0.96;   par_init[33] = -1.0;
   
   // set parameters array.
   /*
  par[0]~par[2] ~ r_x, r_y,r_z
  par[3]~par[7] ~ R_x, R_y, R_z, eps_1, eps_2
  par[8]~par[9] ~ omega, T_2
  */
  vector< vector< int > > par_i(exp_n,vector<int>(10));
  for(int i=0; i<exp_n; i++){
	par_i[i][6] = 2;
	par_i[i][7] = 3;
	par_i[i][8] = 0;
	par_i[i][9] = 1;   
  }
  
  for(int i=0; i<3; i++){
	par_i[0][i] = par_i[1][i] = par_i[2][i] = par_i[3][i] = par_i[4][i] = i+4;
	par_i[5][i] = par_i[6][i] = par_i[7][i] = par_i[8][i] = par_i[9][i] = i+7;
	par_i[10][i] = par_i[11][i] = par_i[12][i] = par_i[13][i] = par_i[14][i] = i+10;
	par_i[15][i] = par_i[16][i] = par_i[17][i] = par_i[18][i] = par_i[19][i] = i+13;
	par_i[20][i] = par_i[21][i] = par_i[22][i] = par_i[23][i] = par_i[24][i] = i+16;
	
	par_i[0][i+3] = par_i[5][i+3] = par_i[10][i+3] = par_i[15][i+3] = par_i[20][i+3] = i+19;
	par_i[1][i+3] = par_i[6][i+3] = par_i[11][i+3] = par_i[16][i+3] = par_i[21][i+3] = i+22;
	par_i[2][i+3] = par_i[7][i+3] = par_i[12][i+3] = par_i[17][i+3] = par_i[22][i+3] = i+25;
	par_i[3][i+3] = par_i[8][i+3] = par_i[13][i+3] = par_i[18][i+3] = par_i[23][i+3] = i+28;
	par_i[4][i+3] = par_i[9][i+3] = par_i[14][i+3] = par_i[19][i+3] = par_i[24][i+3] = i+31;
  }
  
  vector< boost::shared_ptr < TGraphErrors > > gr_ptr(exp_n);
  vector< boost::shared_ptr < TF1 > > my_func_ptr;
  vector< boost::shared_ptr < ModelA > > exp_data_ptr;
  
  
  for(int exp_i=0; exp_i<exp_n; exp_i++){   
	 exp_data_ptr.push_back( boost::make_shared < ModelA > (par_i[exp_i]));
	 my_func_ptr.push_back(boost::make_shared<TF1>("f1", exp_data_ptr[exp_i], &ModelA::expectation_value, 0, T_max, 34));
  }
  
  for(int i_runs=n_min; i_runs<n_runs; i_runs++){
	for(int exp_i=0;exp_i<exp_n;exp_i++){ 
	  sprintf(Tfilename, "gr%d_%d", exp_i,i_runs);
	  gr_ptr[exp_i].reset((TGraphErrors*)f.Get(Tfilename));
	}  

	ROOT::Fit::DataOptions opt;
	ROOT::Fit::DataRange range;
	
	T_min = 0.0;
	T_max = 10.0;
	
	range.SetRange(T_min,T_max);
	
	vector< ROOT::Fit::BinData > data;
	vector< ROOT::Math::WrappedMultiTF1 > wf;
	vector< boost::shared_ptr < ROOT::Math::IMultiGenFunction > > chi2;
	
	for(int i=0; i<exp_n;i++){
	    
	  data.push_back(ROOT::Fit::BinData(opt,range));
	  ROOT::Fit::FillData(data[i],gr_ptr[i].get());
	  wf.push_back(ROOT::Math::WrappedMultiTF1(*(my_func_ptr[i].get()),1));
	  chi2.push_back(boost::make_shared<ROOT::Fit::Chi2Function>(data[i], wf[i]));
	}
	(*chi2[0])(par);
  }

  

  return 0;
}