#include <TCanvas.h>
#include <TF1.h>
#include <TFile.h>
#include <TGraphErrors.h>
#include <TH1.h>

using namespace std;

void Fit(TH1F* hcharge, string Energy);
TFile* ofile;
int rebin=1;
int RunNumber = 247 ;

//////////////////////////////////////////////////////
void FitCharge(){

  //TFile* ifile = new TFile(Form("Output/histo_charge_%d.root",RunNumber),"read");
  TFile* ifile = new TFile(Form("Output/histo_charge_merged.root"),"read");

  ofile = new TFile("Output/histo_charge_fitted.root","recreate");

  TH1F* hall = (TH1F*)ifile->FindObjectAny("hcharge_all");
  Fit(hall,"AllEnergy");
  ofile->Close();
}

//////////////////////////////////////////////////////
void Fit(TH1F* hcharge,string Energy){
  TGraphErrors* gsigma = new TGraphErrors();
  TGraphErrors* gsigmaZ = new TGraphErrors();
  hcharge->Rebin(rebin);
  hcharge->Draw();

  int Zmin = 31;
  int Zmax = 66;


  int Z = Zmin;

  int NumberOfZ = Zmax - Zmin;
  Double_t para[3*NumberOfZ];

  vector<TF1*> g(NumberOfZ);
  vector<double> Integral(NumberOfZ);
  vector<double> Integral_err(NumberOfZ);
  double total_integral = 0;
  vector<double> Yield(NumberOfZ);

  for(int i=0; i<NumberOfZ; i++){
    g[i] = new TF1(Form("g%i",i),"gaus",Z-0.3,Z+0.6);   
    if(i==0){
      g[i]->SetParameter(0,100);
      g[i]->SetParameter(1,Zmin);
      g[i]->SetParameter(2,0.3);
    }
    else if(i>0){
      g[i]->SetParameter(0,g[i-1]->GetParameter(0));
      g[i]->SetParameter(1,Z);
      g[i]->SetParameter(2,g[i-1]->GetParameter(2));
    }
 
    g[i]->SetParLimits(1,Z-0.2,Z+0.2);
    g[i]->SetParLimits(2,0.22,0.4);

   
    hcharge->Fit(g[i],"qBR");

    //g[i]->Draw("lsame");
    g[i]->GetParameters(&para[3*i]);

    Z++;
  }


  TString total_func = "gaus(0)";
  for(int i=1; i<NumberOfZ; i++){
    TString gaus_func = Form("gaus(%i)",3*i);
    total_func += "+" + gaus_func;
  }

  TF1* total = new TF1("total",total_func,Zmin-0.3, Zmin+NumberOfZ -0.7);
  total->SetParameters(para);
  Z = Zmin;
  for(int i=0; i<NumberOfZ; i++){
    total->SetParameter(3*i,g[i]->GetParameter(0));
    total->SetParameter(3*i+1,g[i]->GetParameter(1));
    total->SetParameter(3*i+2,g[i]->GetParameter(2));
    total->SetParLimits(3*i+1,Z-0.2,Z+0.2);
    total->SetParLimits(3*i+2,g[i]->GetParameter(2)*0.5,g[i]->GetParameter(2)*1.3);
    Z++;
  }
  total->SetLineColor(4);
  total->SetNpx(1000);

  hcharge->Fit(total,"RBq");

  Z = Zmin;
  double MeanSigmaZ = 0 ;
  //TF1* one_gauss = new TF1("one_gaus","gaus");
  for(int i=0; i<NumberOfZ; i++){
      
    double Amplitude = total->GetParameter(3*i);
    double mean = total->GetParameter(3*i+1);
    double sigma = total->GetParameter(3*i+2);
    double Amplitude_err = total->GetParError(3*i);
    double mean_err = total->GetParError(3*i+1);
    double sigma_err = total->GetParError(3*i+2);
 
    /*double Amplitude = g[i]->GetParameter(0);
    double mean = g[i]->GetParameter(1);
    double sigma = g[i]->GetParameter(2);
    double Amplitude_err = g[i]->GetParError(0);
    double mean_err = g[i]->GetParError(1);
    double sigma_err = g[i]->GetParError(2);*/
    
    //one_gauss->SetParameter(0,Amplitude);
    //one_gauss->SetParameter(1,mean);
    //one_gauss->SetParameter(2,sigma);
    //one_gauss->SetParError(0,Amplitude_err);
    //one_gauss->SetParError(1,mean_err);
    //one_gauss->SetParError(2,sigma_err);

    cout << "************************" << endl;
    cout << "Z=  " << Z << endl;
    //cout << "Amplitude= " << Amplitude << "+/-" << Amplitude_err << endl;
    //cout << "Mean= " << mean  << endl;
    cout << "Sigma= " << sigma << "+/-" << sigma_err << endl;
    cout << "Sigma/Z =  " << 2.35*sigma /Z*100. << endl;
    
    //if (2.35*sigma/Z*100 > 1 && 2.35*sigma/Z*100 < 3){MeanSigmaZ += 2.35*sigma/Z*100;}
    
    if (i==0 || i==NumberOfZ){continue;}
    MeanSigmaZ += 2.35*sigma/Z*100;

    gsigma->SetPoint(i-1,Z,sigma);
    gsigma->SetPointError(i-1,0,sigma_err);
    
    gsigmaZ->SetPoint(i-1,Z,2.35*sigma/Z*100);

    Integral[i] = (12.5/rebin)*Amplitude*sigma*sqrt(2*TMath::Pi());
    Integral_err[i] = (12.5/rebin)*sqrt(2*TMath::Pi()*(pow(sigma*Amplitude_err,2) + pow(Amplitude*sigma_err,2)));
    //Integral_err[i] = sqrt(Integral[i]);
    

    total_integral += Integral[i];
    Z++;
  }

  MeanSigmaZ = MeanSigmaZ/NumberOfZ;
  cout << "********************" << endl << "Mean Sigma/Z " << MeanSigmaZ << endl;
  double histo_integral = hcharge->Integral();
  hcharge->Write();
  total->Write();

  string output_filename = "Charge_Yield_" + Energy + ".dat";
  ofstream ofile;
  ofile.open(output_filename.c_str());

  // Yield calculation //
  Z = Zmin;
  for(int i=0; i<NumberOfZ; i++){
    //Yield[i] = Integral[i]/total_integral*200;
    Yield[i] = Integral[i]/histo_integral*200;
    //double yield_err = Integral_err[i]/total_integral*200;
    double yield_err = Integral_err[i]/histo_integral*200;
    ofile << Z << " " << Yield[i] << " " << yield_err << endl;
    Z++;
  }
  ofile.close();

  gsigma->SetMarkerStyle(8);
  gsigma->Write();

  TCanvas* csig = new TCanvas("csig","csig",800,800);
  csig->cd();
  gsigma->Draw("ap");

  TCanvas* csigZ = new TCanvas("csigZ","csigZ",800,800);
  gsigmaZ->SetMarkerStyle(8);
  gsigmaZ->GetYaxis()->SetTitle("2.35*#sigma/Z*100");
  gsigmaZ->GetXaxis()->SetTitle("Z");
  csigZ->cd();
  gsigmaZ->Draw("ap");
  total->Draw("lsame");


}