#include<iostream>
#include<iomanip>

#include "TMath.h"
#include "TH1.h"
#include "TF1.h"
#include "TCanvas.h"

TH1F *hTMod;

// Sum of const background and decay functions
Double_t fitFunction(Double_t *x, Double_t *par) {
  Double_t result  = par[0];
  Double_t C       = par[1];
  Double_t T_halfA = par[2];
  Double_t T_halfB = par[3];
  if (T_halfA !=0.) Double_t lambdaA = 0.69314718/T_halfA; // lambda = ln2/T1_2
  if (T_halfB !=0.) Double_t lambdaB = 0.69314718/T_halfB;

  Double_t Tcoll   = 1000.; // collection time = 0.5 Sec 30 mSec per chanel
  Double_t Na0     = 0.;
  
  if (lambdaA !=0) 
     Na0 = C*(1.-TMath::Exp(-lambdaA*Tcoll))/(Tcoll*lambdaA);
  Double_t Nb0     = 0.;
  if (lambdaB !=0 && (lambdaB-lambdaA) !=0) 
     Nb0 = C*( 1.-TMath::Exp(-lambdaB*Tcoll))/(Tcoll*lambdaB ) + 
             ( TMath::Exp(-lambdaB*Tcoll) - TMath::Exp(-lambdaA*Tcoll) )/
             ( Tcoll*(lambdaB-lambdaA) );
  Double_t xx      = x[0];
  Double_t arg     = 0.;
  if ( (lambdaB - lambdaA ) != 0 ) 
     arg = 1. / (lambdaB - lambdaA );
  result += lambdaB*TMath::Exp(-lambdaB*xx)*(Nb0 - lambdaA*Na0*arg) +
            Na0*lambdaA*arg*(2.*lambdaB-lambdaA)*TMath::Exp(-lambdaA*xx);
  return result;
}

void histgen()
{
  // generate data for a test fit
    
  Double_t Tstart  = 200.;
  Double_t Tend    = 5000.+200.;
  
  TF1 *fitFcn =new TF1("fitFcn",fitFunction,Tstart,Tend,4);
  fitFcn->SetNpx(500);
  fitFcn->SetLineWidth(2);
  fitFcn->SetLineColor(kRed);
//   fitFcn->SetParNames("const","C","ThalfA","ThalfB");
  fitFcn->SetParameters(1.,16778222.,450.,925.);
  hTMod = new TH1F("MC Toy Data",  "Decay data",500,Tstart,Tend);
  hTMod->SetLineColor(kBlue);
  hTMod->FillRandom("fitFcn",2426000);
  hTMod->GetListOfFunctions()->Add(fitFcn);
  
  TFile f("toyMC.root","recreate");
  
  hTMod->SetName("result");
  hTMod->Write();
}


void tamuMCtoy_0( void ){

    gROOT->Reset();

    gRandom = new TRandom3();
    
    std::string type = "Minuit";
//     std::string type = "Fumili2";
//     std::string type = "Fumili";
//     std::string type = "Minuit2";
    TVirtualFitter::SetDefaultFitter(type.c_str() );

    gStyle->SetOptFit(1111111);
    gStyle->SetOptStat(1111111);
    gSystem->Load("libMathCore");

    histgen();
    
    TFile *f = new TFile("toyMC.root");
        
    Double_t Tstart  = 200.;
    Double_t Tend    = 5000. + 200.;

    hTMod     = (TH1F*)f->Get("result");
        
    TCanvas *tamuMCtoy_0 = new TCanvas("tamuMCtoy_0","tamuMCtoy_0",10,10,1000,900);

    // create a TF1 with the range from 90 to 15000 and 4 parameters
    // first 2 chennels are not considered in the fit!!!
    TF1 *fitFcn  = new TF1("fitFcn",fitFunction,Tstart,Tend,4);
    TF1 *fitFcn0  = new TF1("fitFcn0",fitFunction,Tstart,Tend,4);
    fitFcn0->SetNpx(500);
    fitFcn0->SetLineWidth(2);
    fitFcn0->SetLineColor(kMagenta);

    fitFcn->SetNpx(500);
    fitFcn->SetLineWidth(2);
    fitFcn->SetLineColor(kRed);

    // Sets initial parameter nvalues
    fitFcn0->SetParameters(1.,1.2e+4,430.,920.);
    // Sets parameters range
    fitFcn0->SetParLimits(0, 1.e-1,50.); // const
    fitFcn0->SetParLimits(1,1.e+4,5.e+7);// C
    fitFcn0->SetParLimits(2, 400.,550.); // ThalfA
    fitFcn0->SetParLimits(3, 800.,950.); // ThalfB

    // Sets initial  parameter names
    fitFcn0->SetParNames("const","C","ThalfA","ThalfB");
      
    hTMod->Fit("fitFcn0");
    
    cout << endl;
    cout << endl;

    fitFcn0->FixParameter(0,1.);
    fitFcn0->FixParameter(3,925.); // half life of 925.0(0) mSec

    hTMod->Fit("fitFcn0","VLE");
    hTMod->SetLineColor(kBlue);
    tamuMCtoy_0->Modified();        
    tamuMCtoy_0->SetLogy();
    
    // Get parameter from the fit    
    Double_t cnst    = fitFcn0->GetParameter(0);
    Double_t C       = fitFcn0->GetParameter(1);
    Double_t lambdaA = fitFcn0->GetParameter(2);
    Double_t lambdaB = fitFcn0->GetParameter(3);

    // Draw fit    
    fitFcn->SetParameter(0,1.);
    fitFcn->SetParameter(1,16778222.);
    fitFcn->SetParameter(2,450.);
    fitFcn->SetParameter(3,925.);
    fitFcn->Draw("same");    
    tamuMCtoy_0->Modified();        
        
    cout << setprecision(8)<< "\n--->   const: "<< setw(10) << cnst << setw(4)<<"+/-" << setw(12)<< fitFcn0->GetParError(0) << endl;
    cout << setprecision(8)<< "--->       C: "<< setw(10)<< C << setw(4)<<"+/-" << setw(12)<< fitFcn0->GetParError(1) << endl;
    cout << setprecision(8)<< "---> lambdaA: "<< setw(10)<<lambdaA <<  setw(4)<<"+/-" << setw(12)<< fitFcn0->GetParError(2) << endl;
    cout << setprecision(8)<< "---> lambdaB: "<< setw(10)<<lambdaB <<  setw(4)<<"+/-" << setw(12)<< fitFcn0->GetParError(3) << endl;
    cout << setprecision(8)<< "\n--->    chi2: "<< setw(10)<<fitFcn0->GetChisquare() << endl;
    cout << endl;
    
//     hTMod->Draw("same");
    hTMod->Draw("e");
    tamuMCtoy_0->Update();
    f->ls();
    cout << endl;
    fitFcn->Print();
    cout << endl;
    cout << endl;
    fitFcn0->Print();
}