#include "TVirtualFitter.h"

gROOT->Reset();

// Decay A function
Double_t decayA(Double_t *x, Double_t *par) {
   Double_t T_halfA = par[1];
   if (T_halfA !=0.) Double_t lambdaA = 0.69314718/T_halfA; // lambda = ln2/T1_2
   Double_t Na0       = par[0];
   Double_t xx      = x[0];
   return Na0*lambdaA*TMath::Exp( -lambdaA*(x[0]));
}

int histgen( void )
{
  // generate data for a test fit    
  Double_t Tstart  = 1.;
  Double_t Tend    = 2880722.;
  
  TF1 *fitFcn =new TF1("decayA",decayA,Tstart,Tend,2);
  fitFcn->SetNpx(500);
  fitFcn->SetLineWidth(2);
  fitFcn->SetLineColor(kRed);
  fitFcn->SetParNames("const","T_{1/2}");
  // Initial # of Counts, Half-Life in seconds
  fitFcn->SetParameters(1.8882258e+12,232945.88);

  hTMod = new TH1F("MC Toy Data",  "Decay data ",127,Tstart,Tend);
  hTMod->SetLineColor(kBlue);
  hTMod->FillRandom("decayA",100000000);
  hTMod->GetListOfFunctions()->Add(fitFcn);

  TFile *f = new TFile("Au198_ToyMC.root","RECREATE");

  Float_t time,Count,errCount;

  TTree *tree = new TTree("data","data from ascii file");
  tree->Branch(    "time", &time,    "time/F");
  tree->Branch(   "Count", &Count,   "Count/F");
  tree->Branch("errCount", &errCount,"errCount/F");
  
  Int_t nbins = hTMod->GetXaxis()->GetNbins();

  for (int bin=0; bin<nbins; bin++) {
     time    = hTMod->GetXaxis()->GetBinCenter(bin);
     Count   = hTMod->GetBinContent(bin);
    errCount = hTMod->GetBinError(bin);
//    printf("%6.f %6.2f %6.2f \n",time,Count,errCount);
    if (bin != 0) tree->Fill(); 
  }

  hTMod->SetName("result");
  hTMod->Write();

  f->Write();
  f->Close();

  return nbins;
}


void Au198_4( void ) {
   //Draw a graph with error bars
   // To see the output of this macro:
   // 1) root [0] .x Au198_4.cxx
   // Author: V.V. Golovko @ 9 Aug 2007
   // modified: by V.V. Golovko @16 Aug 2007
   // modified: by V.V. Golovko @17 Oct 2007
   // modified: by V.V. Golovko @23 Oct 2007
   // to fit Half-Life data for 198-Au
 
//     std::string type = "Minuit";
//     std::string type = "Fumili2";
//     std::string type = "Fumili";
    std::string type = "Minuit2";
   TVirtualFitter::SetDefaultFitter(type.c_str() );
    
   TVirtualFitter * minuit = TVirtualFitter::Fitter(0,2);

   gStyle->SetOptFit(0);
   gStyle->SetOptStat(0);

// generate data from Toy Monte Carlo
// "fit function" as a analitical function    
   Int_t nlines = histgen();

   printf("\nnlines= %i \n\n", nlines);

   Float_t time,Count,errCount;
   TFile *f = new TFile("Au198_ToyMC.root");
   TTree *tree = (TTree*)f->Get("data");

   TBranch *b_time     = tree->GetBranch("time");
   TBranch *b_Count    = tree->GetBranch("Count");
   TBranch *b_errCount = tree->GetBranch("errCount");
   tree->SetBranchAddress("time",         &time, &b_time);
   tree->SetBranchAddress("Count",       &Count, &b_Count);
   tree->SetBranchAddress("errCount", &errCount, &b_errCount);

   Int_t nevent = (Int_t)tree->GetEntries();
   printf("\nnevent= %i \n\n", nevent);

   Float_t *x  = new Float_t[nevent];
   Float_t *ex = new Float_t[nevent];
   Float_t *y  = new Float_t[nevent];
   Float_t *ey = new Float_t[nevent];
                            
   for (int bin=0; bin<nevent; bin++) {
         b_time->GetEntry(bin);
        b_Count->GetEntry(bin);
     b_errCount->GetEntry(bin);
        x[bin] = time;
       ex[bin] = 0.0;
        y[bin] = Count;
       ey[bin] = errCount;
    if (bin < 5) printf("%8.f %8.2f %8.2f \n",x[bin],y[bin],ey[bin]);
   }

   TCanvas *c1 = new TCanvas("198Au","198Au",200,10,1000,900);

// ------------>Primitives in pad: 198Au_1
   TPad *Au_1 = new TPad("198Au_1", "198Au_1",0.01,0.31,0.99,0.99);
   Au_1->Draw();
   Au_1->cd();
   Au_1->Range(-396099,1.40112,3.56489e+06,7.38995);
   Au_1->SetLogy();

   TGraphErrors *gr = new TGraphErrors(nevent,x,y,ex,ey);

   gr->SetTitle("");
   gr->SetMarkerColor(4);
   gr->SetMarkerSize(0.5);
   gr->SetMarkerStyle(20);
   gr->GetYaxis()->SetTitle("Counts");
   gr->GetXaxis()->SetTitle("T_{1/2}, s");
   
   gr->Draw("ap"); 

   // initiate fit function
   TF1 *fitFcn = new TF1("fitFcn",decayA,0.,28807222.,2);
   // Sets parameter names
   fitFcn->SetParNames("const","T_{1/2}");

   // Sets initial parameter nvalues
   fitFcn->SetParameters(7.e+12, 8.e+5);
   // Sets parameters range
   fitFcn->SetParLimits(0, 1.e+12,9.9e+12); // const 1.88822e+12
   fitFcn->SetParLimits(1, 1.e+5, 9.9e+5);  // Thalf

   cout << endl;
   cout << endl;
   cout << "-----> " << endl;

   gr->Fit(fitFcn,"V+","eprL"); // do the fit

   //get result
   Double_t minParams[4];
   Double_t parErrors[4];
    
   // get the fit parameters from the fit
   for (int i=0; i<2; i++)
   {
     minParams[i] = minuit->GetParameter(i);
     parErrors[i] = minuit->GetParError(i);
   }

   double chi2, edm, errdef;
   int nvpar, nparx;
   
   minuit->GetStats(chi2,edm,errdef,nvpar,nparx);

   // setting parametrs to the fit function for drawing
   fitFcn->SetParameters(minParams);
   fitFcn->SetParErrors(parErrors);
   fitFcn->SetChisquare(chi2);
   fitFcn->SetNDF(gr->GetXaxis()->GetNbins()-nvpar);

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

   c1->Update();
   c1->Modify();

   // get the fit parameters from minuit
   Double_t par,parerr;
   par        = fitFcn->GetParameter(1);
   parerr     = fitFcn->GetParError(1);
   Int_t  ndf = fitFcn->GetNDF();

   printf("ndf  = %f, \n", ndf);
   printf("Chi^2= %f, \n", chi2);
   printf("number of var. parameters = %i, \n",nvpar);

   gr->Draw("AP");
   fitFcn->Draw("same");

   // Get parameter from the fit    
   Double_t cnstP;
   if ( ndf != 0.) {
     if (chi2/ndf <= 1.) cnstP = par + parerr;
       else cnstP = par + parerr*sqrt(chi2/ndf);
   }
   
   // Text Labels
   if ( ndf != 0.) {
     if (chi2/ndf <= 1.) tex = new TLatex(1.e6,6.e5,Form("#bar{T_{1/2}}:  %6.4f^{ }#pm^{ }%6.1g^{ }d",par/86400.,parerr));
       else tex = new TLatex(1.e6,6.e5,Form("#bar{T_{1/2}}:  %6.4f^{ }#pm^{ }%6.1g^{ }d",par/86400.,parerr*sqrt(chi2/ndf)/86400.));
   }
   tex->SetTextSize(0.06);
   tex->SetLineWidth(2);
   tex->SetTextColor(kBlue);
   tex->SetTextAngle(0);
   tex->Draw();
      
   // Get parameter from the fit    
   Double_t cnstM;
   if ( ndf != 0.) {
     if (chi2/ndf <= 1.) cnstM = par - parerr;
       else cnstM = par - parerr*sqrt(chi2/ndf);
   }

   if ( ndf != 0.) {
     if (chi2/ndf <= 1.) printf("\ncnst: %g Up: %g Down: %g Error: %g\n\n", par, cnstP, cnstM, parerr);
       else printf("\ncnst: %g Up: %g Down: %g Error corrected: %g\n\n", par, cnstP, cnstM, parerr*sqrt(chi2/ndf));
   }
   
   if ( ndf != 0.) {
     if (chi2/ndf <= 1.) printf("\nIn days: cnst: %6.4f Error: %6.4f\n\n", fitFcn->GetParameter(1)/86400., fitFcn->GetParError(1)/86400.);
       else printf("\nIn days: cnst: %6.4f Error corrected: %6.4f\n\n", 
                   fitFcn->GetParameter(1)/86400., fitFcn->GetParError(1)*sqrt(chi2/ndf)/86400.);
   }

   Float_t *res   = new Float_t[nevent];
   Float_t *resEr = new Float_t[nevent];
   for (int bin=0; bin<nevent; bin++) {
     res[bin]   = (y[bin] - fitFcn->Eval(x[bin]))/ey[bin];
     resEr[bin] = (y[bin] - fitFcn->Eval(x[bin]))/sqrt(fitFcn->Eval(x[bin]) * y[bin]);
     if (bin < 5) printf("x=%8.0f, y=%8.2f, ey=%8.4f sigma=%8.4f res=%8.4f, resEr=%8.4f\n",
                          x[bin],y[bin],ey[bin],sqrt(fitFcn->Eval(x[bin])),res[bin],resEr[bin]);
   }

// ------------>Primitives in pad: 198Au_2
   c1->cd();
   TPad *Au_2 = new TPad("198Au_2", "198Au_1",0.01,0.31,0.99,0.99);
   //TPad *Au_2 = new TPad("198Au_2", "198Au_2",0.01,0.01,0.99,0.29);
   Au_2->SetFrameFillStyle(4000); //<====== make pad transparent
   Au_2->Draw();
   Au_2->cd();
   Au_2->Range(-396099,-8.04516,3.56489e+06,4.90887);
   Au_2->SetBorderSize(2);
   Au_2->SetFrameFillColor(0);
      
   TGraphErrors *gr1 = new TGraphErrors(nevent,x,res,0,resEr);

   gr1->SetTitle(0);
   gr1->SetMarkerColor(kBlue);
   gr1->SetMarkerSize(0.5);
   gr1->SetMarkerStyle(20);
   gr1->GetYaxis()->SetTitle("#frac{Counts - fit}{#sqrt{fit}}");
   gr1->GetXaxis()->SetTitle("T_{1/2}, s");
   gr1->Draw("APLY+"); //<====== note option Y+

   TF1 *func = new TF1("func","pol0",0.,28807222.);
   func->SetLineColor(kMagenta);
   func->SetLineWidth(2);
   func->SetParameter(0,0.);
   func->Draw("same");
   c1->Update();
   c1->Modify();


   Double_t ymean = TMath::Mean(nevent,gr1->GetY());
   printf("ymean=%g\n",ymean);
   

   c1->Print("Au198Toy_MC_HL.ps");

   FILE *his_out = fopen("DataHis_Au198.dat","wt");

   if(!his_out)
   {
     printf("Error open file: DataHis_Au198.dat \n");
     exit(1);
   }

   for (int bin=0; bin<nevent; bin++) {
    fprintf(his_out,"%9.f %9.f %9.f \n",x[bin],y[bin],ey[bin]);
   }    
   fclose(his_out);
}