// CPP includes
#include <cstdlib>
#include <cstdio>
#include <iostream>
#include <cmath>
#include <vector>
#include <string>
#include <iostream>
#include <sstream>
#include <fstream>

// C includes
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

// ROOT includes
// ROOT includes
#include <TSystem.h>
#include <TROOT.h>
#include <TApplication.h>
#include <TRint.h>
#include <TDirectory.h>
#include <TFile.h>
#include <TChain.h>
#include <TTree.h>
#include <TNtuple.h>
#include <TBranch.h>
#include <TH1.h>
#include <TAxis.h>
#include <TCanvas.h>
#include <TColor.h>
#include <TPad.h>
#include <TStyle.h>
#include <TKey.h>
#include <TLegend.h>
#include <TText.h>
#include <TPaveText.h>
#include <TLine.h>
#include <TMarker.h>
#include <TLatex.h>
#include <TMath.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TF2.h>
#include <TString.h>
#include "TObjString.h"
#include "TRandom.h"
#include "TRandom3.h"
#include "TMinuit.h"

using namespace std;



#define VERSION "v1.06082015";

// Fitting function
double fitf_cb(Double_t *x, Double_t *par) 
{
  Double_t arg = 0;
  if (par[2] != 0) arg = (x[0] - par[1])/par[2];
  Double_t fitval = par[0]*TMath::Exp(-2.772589*arg*arg)+par[3];
  return fitval;
}


Bool_t reject;
Double_t fline(Double_t *x, Double_t *par)
{
    if (reject && x[0] > 267 && x[0] < 274) {
      TF1::RejectPoint();
      return 0;
   }
  Double_t arg = 0;
 // arg = (x[0] - par[1])/par[2]/sqrt(2);
  Double_t fitval = par[2]/(1+TMath::Exp((x[0]-par[0])/0.75/par[1]));  // *TMath::Erfc(arg);
  return fitval;
}

// Fitting function
double fitf_bckg(Double_t *x, Double_t *par) 
{
  Double_t arg = 0;
  arg = par[0]+par[1]*x[0]+par[2]*x[0]*x[0]+par[3]*x[0]*x[0]*x[0];
  Double_t fitval = arg;
  return fitval;
}



Double_t fcub(Double_t *x, Double_t *par)
{
    if (reject && x[0] > 174&& x[0] < 190) {
      TF1::RejectPoint();
      return 0;
   }
  Double_t fitval = par[0]+par[1]*x[0]+par[2]*x[0]*x[0] + par[3]*x[0]*x[0]*x[0];
  return fitval;
}


Double_t fcub2(Double_t *x, Double_t *par)
{
  Double_t fitval = par[0]+par[1]*x[0]+par[2]*x[0]*x[0] + par[3]*x[0]*x[0]*x[0];
  return fitval;
}


Double_t mygaus(Double_t *x, Double_t *par)
{
   Double_t arg = 0;
   if (par[2] != 0) 
   {  
     arg = (x[0] - par[1])/par[2];
     Double_t fitval = par[0]/ (sqrt(2.0*TMath::Pi())*par[2])*exp(-0.5*arg*arg);
     return fitval;
   }
   else
     return 0;
}


Double_t mygaus2(Double_t *x, Double_t *par)
{
   Double_t arg = 0;
   if (par[2] != 0) 
   {  
     arg = (x[0] - par[1])/par[2];
     Double_t fitval = par[0]/ (sqrt(2.0*TMath::Pi())*par[2])*exp(-0.5*arg*arg);
     return fitval;
   }
   else
     return 0;
}


// Sum of background and 2 gaussian functions
Double_t fitFunction(Double_t *x, Double_t *par) {
  return fcub2(x,par) + mygaus(x,&par[4]) + mygaus2(x,&par[7]);
}




void run(){
    
    gRandom->SetSeed(0.);
 
    FILE*infile;
    if((infile=fopen("Output.dat","r"))==NULL){
          printf("   ERROR: Cannot open cycles file.\n");
          exit(0);}
    int binNb,x;
    int cont =0;
    double Emin=2.E-05*1.e+6;
    double Emax=2.E-05*1024*1.e+6;
    TH1F*h = new TH1F("h","h",1024,0,1023);
    TH1F*h2 = new TH1F("h2","h2",1024,Emin,Emax );
    while(1)
         {
          fscanf(infile,"%d %d ",&binNb ,&x);
         
          h->Fill(binNb,x);
          h2->SetBinContent(binNb,x);
          cont = cont +1;
         if(feof(infile)!=NULL) break;
         }
       
  TCanvas*c1 = new TCanvas("c1","c1",500,500);
  c1->cd();
  h->Draw();

// Fit background by rejecting peaks
  
  reject = kTRUE;
  TF1 *cub_fit_p3    = new TF1("cub_fit_p3",fcub,158,230,4);
  cub_fit_p3->SetParameter(0,72000);
  cub_fit_p3->SetParameter(1,-1200);
  cub_fit_p3->SetParameter(2,6.76);
  cub_fit_p3->SetParameter(2,-0.012);
  
  cub_fit_p3->SetLineColor(4);
  h->Fit(cub_fit_p3,"R+");
  
 // Draw background 

  TF1 *cub_fit_p3b    = new TF1("cub_fit_p3",fcub2,158,230,4);

  cub_fit_p3b->FixParameter(0,cub_fit_p3->GetParameter(0));
  cub_fit_p3b->FixParameter(1,cub_fit_p3->GetParameter(1));
  cub_fit_p3b->FixParameter(2,cub_fit_p3->GetParameter(2));
  cub_fit_p3b->FixParameter(3,cub_fit_p3->GetParameter(3));
  
  cub_fit_p3b->SetLineColor(3);
  cub_fit_p3b->Draw("same");
  

  
  
  
// Fits one peak  
  TF1 *gaus_fit_p3    = new TF1("gaus_fit_p3","gaus",184,189);
  gaus_fit_p3->SetLineColor(5);
  h->Fit(gaus_fit_p3,"R+");  
  
  double mean_p3 = gaus_fit_p3->GetParameter(1);
  double sigma_p3 = gaus_fit_p3->GetParameter(2);
  
  
  
 // Fits other peak 
  
  TF1 *gaus_fit_p4    = new TF1("gaus_fit_p4","gaus",176,184);
  gaus_fit_p4->SetLineColor(6);
  h->Fit(gaus_fit_p4,"R+");  

  double mean_p4 = gaus_fit_p4->GetParameter(1);
  double sigma_p4 = gaus_fit_p4->GetParameter(2);



// Reconstruct peaks
  
  TH1D*h3= new TH1D("h3","h3",1024,0,1024);
  for(int i= 158;i<230;i++) h3->SetBinContent(i+1, cub_fit_p3b->Eval(i));
  h3->SetLineColor(1);
  h3->Draw("same");


   cout << mean_p3 << " " << sigma_p3<< endl;
   int xx;
   for(int i= 1;i<7000;i++) {
       h3->Fill((int)gRandom->Gaus(mean_p3,sigma_p3/3));
     }
   for(int i= 1;i<1500;i++) {
       h3->Fill((int)gRandom->Gaus(mean_p4,sigma_p4/5));
     }

     h3->Draw("same");
  
 
   
   
// Fit everything   

   TF1 *fitFcn = new TF1("fitFcn",fitFunction,158,230,10);
   fitFcn->SetNpx(500);
   fitFcn->SetLineWidth(4);
   fitFcn->SetLineColor(kBlue);
    
   // first try without starting values for the parameters
   // This defaults to 1 for each param. 
   // this results in an ok fit for the polynomial function
   // however the non-linear part (lorenzian) does not 
   // respond well.
   fitFcn->SetParameter(0,cub_fit_p3->GetParameter(0));
   fitFcn->SetParameter(1,cub_fit_p3->GetParameter(1));
   fitFcn->SetParameter(2,cub_fit_p3->GetParameter(2));
   fitFcn->SetParameter(3,cub_fit_p3->GetParameter(3));
   fitFcn->SetParameter(4,90000);
   fitFcn->SetParameter(5,181);
   fitFcn->SetParameter(6,3);
   fitFcn->SetParameter(7,90000);
   fitFcn->SetParameter(8,186);
   fitFcn->SetParameter(9,2);


   h->Fit("fitFcn","R+");
   cout << fitFcn->GetNDF() << endl;
   cout << fitFcn->GetChisquare() << endl;
   cout << fitFcn->GetProb() << endl;
   

  TCanvas*c2 = new TCanvas("c2","c2",500,500);
  c2->cd();
  h2->Draw("c");
  gPad->SetLogy(1);  
  
}