#include "utilities.h"
//#include "macro01.h"

TF1* fitSCS(TH1D * hSCS) {
  // creates the fitting functions consisting out of a double gaussian and a tail
  // and fits it to the given histogram
  // also defines the fitting range
  Double_t start = 0.02335;     // 0.02335 
  Double_t end = 0.0236;
  hSCS->GetXaxis()->SetRangeUser(0.023, 0.024);

  TF1 * fFit = new TF1("fFit", "gaus(0) + gaus(3) + [6] * (exp([7]*(x-[8])) + [9])/(exp((x-[8])/[10])+[8])", start, end); 
  // gaus(0) is a substitute for [0]*exp(-0.5*((x-[1])/[2])**2) and (0) means start numbering parameters at 0
  fFit->SetParName(0, "constant1");
  fFit->SetParName(1, "mean1");
  fFit->SetParName(2, "sigma1");
  fFit->SetParName(3, "constant2");
  fFit->SetParName(4, "mean2");
  fFit->SetParName(5, "sigma2");
  fFit->SetParName(6, "m");
  fFit->SetParName(7, "f");
  fFit->SetParName(8, "l");
  fFit->SetParName(9, "d");
  fFit->SetParName(10, "t");

  // provide some help to our dear fitter
  fFit->SetParameter(0, 350);       // height?
  fFit->SetParameter(1, 0.02345);   // mean
  fFit->SetParameter(2, 0.00005);   // sigma
  fFit->SetParameter(3, 400);       // height?
  fFit->SetParameter(4, 0.02354);   // mean
  fFit->SetParameter(5, 0.00005);   // sigma
  fFit->SetParameter(6, 9);         // m
  fFit->SetParameter(7, 1);         // f
  fFit->SetParameter(8, 3);         // l 
  fFit->SetParameter(9, 1);         // d
  fFit->SetParameter(10,2);         // t
  
  //fFit->SetParameter(0, hSCS->GetMaximum());

  hSCS->Draw();
  hSCS->Fit(fFit, "R  Q");   // use function range, don't plot, append to the list of fitted functions, quiet mode


  return fFit;
}

void macro02() {
  // access the ROOT file, the generated plots and histograms and create a new one for the results
  TFile * analysis = TFile::Open("analysis/analysis.root", "READ");
  TDirectory * dirSCS = (TDirectory *) analysis->Get("dirSCS");
  TFile * results = TFile::Open("analysis/results.root", "RECREATE");

  // acess the A/E histograms in the SCS region and fit a double gausian plus tail on them
  TDirectory * dirSCS_fitted = results->mkdir("dirSCS_fitted", "dirSCS_fitted");
  Int_t nSCSregions = SCS_start.size();   // 27
  std::vector<TH1D*> hSCS;
  std::vector<TF1*> fFit;
  std::vector<Double_t> lMeans, uMeans, lMeanErrors, uMeanErrors;
  
  for(Int_t i=0; i<nSCSregions; i++) {
    std::string histoName = "hSCS[" + std::to_string(i) + "]";
    std::string funcName = "fFit" + std::to_string(i);
    hSCS.push_back( (TH1D*) dirSCS->Get(histoName.c_str()) );       // place the next histogram at the end of the vector; therefore the highest entergy histogram is located at the end of the vector after the for loop
    
    TCanvas * cSCS = new TCanvas( ((std::string)"cSCS_" + std::to_string(i)).c_str() , "SCS region");
    fFit.push_back( fitSCS(hSCS[i]) );                              // do the fitting
    cSCS->SaveAs( ((std::string)"cSCS_" + std::to_string(i) + (std::string)".pdf").c_str(), "pdf");
    delete cSCS;
    fFit[i]->SetName(funcName.c_str());

    Double_t lower_mean = fFit[i]->GetParameter(1);                 // and now extract the means of the gausians
    Double_t lower_mean_error = fFit[i]->GetParError(1);
    Double_t upper_mean = fFit[i]->GetParameter(4);
    Double_t upper_mean_error = fFit[i]->GetParError(4);
    //printf("i = %d: lower_mean_error = %f, upper_mean_error = %f\n", i, lower_mean_error, upper_mean_error);
    if (lower_mean >= upper_mean) {printf("Oh no! The lower mean is larger or equal than the upper mean.\n");}
    lMeans.push_back(lower_mean);
    lMeanErrors.push_back(lower_mean_error);
    uMeans.push_back(upper_mean);
    uMeanErrors.push_back(upper_mean_error);
  }

  // create the A/E position vs energy plot both for the lower and the upper peaks
  TGraphErrors * glMeans = getGraphWithYErrors(SCS_start, lMeans, lMeanErrors);
  TGraphErrors * guMeans = getGraphWithYErrors(SCS_start, uMeans, uMeanErrors);
  glMeans->SetName("glMeans");
  guMeans->SetName("guMeans");
  glMeans->SetTitle("mean positions of the lower energy gaus");
  guMeans->SetTitle("mean positions of the higher energy gaus");
  glMeans->GetXaxis()->SetTitle("energy (keV)");
  guMeans->GetXaxis()->SetTitle("energy (keV)");
  glMeans->GetYaxis()->SetTitle("A/E (a.u.)");
  guMeans->GetYaxis()->SetTitle("A/E (a.u.)");
  glMeans->GetYaxis()->SetRangeUser(0.02343, 0.02346);
  guMeans->GetYaxis()->SetRangeUser(0.023523, 0.023542);
  glMeans->SetMarkerStyle(8);
  guMeans->SetMarkerStyle(8);
  printf("%zu, %zu and %zu\n", SCS_start.size(), lMeans.size(), uMeans.size());
 
  // fit the plots with a first oder polynomial
  glMeans->Fit("pol1", "Q");
  guMeans->Fit("pol1", "Q");

  glMeans->Write();
  guMeans->Write();
  //dirSCS_fitted->Write();
  for (Int_t i=0; i<hSCS.size(); i++) {
    hSCS[i]->SetName( ((std::string)"hSCS_" + std::to_string(i)).c_str() );
    hSCS[i]->Write();
  }

  delete results;
}