#include "TF1.h"
#include "TMath.h"
#include "TGraphErrors.h"
#include "Fit/Fitter.h"
#include "Fit/BinData.h"
#include "Fit/Chi2FCN.h"
#include "TH1.h"
#include "TList.h"
#include "Math/WrappedMultiTF1.h"
#include "HFitInterface.h"
#include "TCanvas.h"
#include "TStyle.h"

// definition of shared parameter

int iparS[4] = {0,1,2,3};
int iparR[4] = {0,1,2,3};

// Create the GlobalCHi2 structure

struct GlobalChi2 {
   GlobalChi2(  ROOT::Math::IMultiGenFunction & f1,
                ROOT::Math::IMultiGenFunction & f2) :
      fChi2_1(&f1), fChi2_2(&f2) {}

   double operator() (const double *par) const {
      double p1[4];
      for (int i = 0; i < 4; ++i) p1[i] = par[iparS[i] ];

      double p2[4];
      for (int i = 0; i < 4; ++i) p2[i] = par[iparR[i] ];

      return (*fChi2_1)(p1) + (*fChi2_2)(p2);
   }

   const  ROOT::Math::IMultiGenFunction * fChi2_1;
   const  ROOT::Math::IMultiGenFunction * fChi2_2;
};

//define fitfunctios

Double_t sigma_base(Double_t * s, Double_t * par)
{
double alpha_1 = par[0];
double a_1 = par[1];
double alpha_2 = 0.703;
double a_2 = par[2];
double alpha_3 = 0.435;
double a_3 = par[3];
double EA_1=a_1*TMath::Cos(TMath::Pi()*alpha_1/2)*TMath::Power(s[0]*s[0],alpha_1);
double MA_1=-a_1*TMath::Sin(TMath::Pi()*alpha_1/2)*TMath::Power(s[0]*s[0],alpha_1);
double EA_2=a_2*TMath::Cos(TMath::Pi()*alpha_2/2)*TMath::Power(s[0]*s[0],alpha_2);
double MA_2=-a_2*TMath::Sin(TMath::Pi()*alpha_2/2)*TMath::Power(s[0]*s[0],alpha_2);
double EA_3=a_3*TMath::Sin(TMath::Pi()*alpha_3/2)*TMath::Power(s[0]*s[0],alpha_3);
double MA_3=a_3*TMath::Cos(TMath::Pi()*alpha_3/2)*TMath::Power(s[0]*s[0],alpha_3);
double EA=EA_1+EA_2-EA_3;
double MA=MA_1+MA_2-MA_3;
return (1.556*TMath::Pi()/(s[0]*s[0]))*MA;
};

Double_t rho_base(Double_t * s, Double_t * par)
{
double alpha_1 = par[0];
double a_1 = par[1];
double alpha_2 = 0.703;
double a_2 = par[2];
double alpha_3 = 0.435;
double a_3 = par[3];
double EA_1=a_1*TMath::Cos(TMath::Pi()*alpha_1/2)*TMath::Power(s[0]*s[0],alpha_1);
double MA_1=-a_1*TMath::Sin(TMath::Pi()*alpha_1/2)*TMath::Power(s[0]*s[0],alpha_1);
double EA_2=a_2*TMath::Cos(TMath::Pi()*alpha_2/2)*TMath::Power(s[0]*s[0],alpha_2);
double MA_2=-a_2*TMath::Sin(TMath::Pi()*alpha_2/2)*TMath::Power(s[0]*s[0],alpha_2);
double EA_3=a_3*TMath::Sin(TMath::Pi()*alpha_3/2)*TMath::Power(s[0]*s[0],alpha_3);
double MA_3=a_3*TMath::Cos(TMath::Pi()*alpha_3/2)*TMath::Power(s[0]*s[0],alpha_3);
double EA=EA_1+EA_2-EA_3;
double MA=MA_1+MA_2-MA_3;
return EA/MA;
};



void tr(){

TGraphErrors *dS = new TGraphErrors("datS.dat","%lg %lg %lg %lg");

TGraphErrors *dR = new TGraphErrors("datR.dat","%lg %lg %lg %lg");
   
TF1 *sigma = new TF1("sigma",sigma_base,1,100000);
sigma->SetParameters(1.1, -1.5, -4, 4);

TF1 *rho = new TF1("rho",rho_base,1,100000);
rho->SetParameters(1.1, -1.5, -4, 4);


   // perform now global fit


     
   ROOT::Math::WrappedMultiTF1 wsigma(*sigma,1);
   ROOT::Math::WrappedMultiTF1 wrho(*rho,1);
 

   ROOT::Fit::DataOptions opt;
  
   ROOT::Fit::BinData::ErrorType kValueError;  

   ROOT::Fit::DataRange rangeS;
   rangeS.SetRange(1,100000);
   ROOT::Fit::BinData dataS(opt,rangeS);
   ROOT::Fit::FillData(dataS, dS);

   ROOT::Fit::DataRange rangeR;
   rangeR.SetRange(1,100000);
   ROOT::Fit::BinData dataR(opt,rangeR);
   ROOT::Fit::FillData(dataR, dR); 

   ROOT::Fit::Chi2Function chi2_sigma(dataS, wsigma);
   ROOT::Fit::Chi2Function chi2_rho(dataR, wrho);

   GlobalChi2 globalChi2(chi2_sigma, chi2_rho);

   ROOT::Fit::Fitter fitter;


   const int Npar = 4;
   double par0[Npar] = {1.1, -1.5, -4, 4};
   fitter.Config().SetParamsSettings(4,par0);


   fitter.Config().MinimizerOptions().SetPrintLevel(0);
   fitter.Config().SetMinimizer("Minuit","Migrad");

   // fit FCN function directly
   fitter.FitFCN(4,globalChi2,0,dataS.Size()+dataR.Size(),true);
   ROOT::Fit::FitResult result = fitter.Result();
   result.Print(std::cout);

}