#include "TFile.h"
#include "TH1.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TAxis.h"
   
TF1 *Gks, *Gkf, *eLst, *fInt;

double fint(double *x, double *par) {
   double tp = x[0];
   double alpha_k = par[0];
   double gamma   = par[1];
   double lambda_k= par[2];
   double tau_c   = par[3];
   //double econst  = par[4];
   //double Lambdas = par[5];
   double t       = par[6];
    Gks->SetParameters(alpha_k,gamma);
    Gkf->SetParameters(lambda_k,tau_c);
   return lambda_k*Gkf->Eval(tp)*eLst->Eval(tp)*Gks->Eval(t-tp);
}
double gkcalc(double *x, double *par) {
   double t = x[0];
   double econst  = par[4];
   double Lambdas = par[5];
   fInt->SetParameters(par[0],par[1],par[2],par[3],econst,Lambdas,t);
    eLst->SetParameters(econst,Lambdas);
   return eLst->Eval(t)*Gkf->Eval(t) +fInt->Integral(0,t);
}

void gks2() {
   double alpha_k = 0.221243;
   double gamma   = 0.4;
   double lambda_k= 0.172;
   double tau_c   = 0.5;
   double Lambdas = 0.0345;
   double econst  = 1;
    //Calculate Gks
    //parameter[0] = alpha_k
    //parameter[1] = gamma
    Gks = new TF1("Gks","[0]+(1-[0])*exp(-[1]*x)",0,15);
    Gks->SetParameters(alpha_k,gamma);
    //Calculate Gkf
    //parameter[0] = lambda_k
    //parameter[1] = tau_c
    Gkf = new TF1("Gkf","(1-[0]*[1])*exp(-[0]*x)",0,15);
    Gkf->SetParameters(lambda_k,tau_c);
    //Calculate expLamba*t
    //parameter[0] = econst = 1    
    //parameter[1] = Lambda*
    eLst = new TF1("eLst","[0]*exp(-[1]*x)",0,15);
    eLst->SetParameters(econst,Lambdas);
    //calculate the function to be integrated
    fInt = new TF1("fInt",fint,0,15,7);
    fInt->SetParameters(alpha_k,gamma,lambda_k,tau_c,econst,Lambdas,0);
    //Combine Gks and Gkf to for Gkcalc
    //TF1 *Gkcalc = new TF1("Gkcalc","Gkf*eLst*[0]*[1]",0,15);  // This is the bit I can't get!!!!!!!!
    TF1 *Gkcalc = new TF1("Gkcalc",gkcalc,0,15,6);  // This is the bit I can't get!!!!!!!!
    Gkcalc->SetParNames("alpha_k","gamma","lambda_k","tau_c","econst","Lambdas");
    Gkcalc->SetParameters(alpha_k,gamma,lambda_k,tau_c,econst,Lambdas);
    //Create a canvas to draw Gkcalc
    //TCanvas *Gkcalc1 = new TCanvas("Gkcalc1","Gkcalc1",10,10,1400,1200);
    //new TCanvas("Gkcalc1","Gkcalc1",10,10,700,600);
    //open histogram file and fit histogram with Gkcalc
    TFile *file = TFile::Open("MCGcalc.root");
    TCanvas *canvas = (TCanvas*)file->Get("MCGcalc");
    canvas->Draw();
    TH1F *h = (TH1F*)canvas->FindObject("G2values");
    h->Fit(Gkcalc,"w");
  
}