// Overall, it makes the code simpler, more readable, and more stable: less magic buffer-size number etc.
// V2: Add tight gaussians

# include <iostream>
# include <fstream>
# include <math.h>
# include <iomanip>
# include <cmath>
# include <stdlib.h>
# include <cstdlib>
//# include <fstream.h>
# include <string.h>
# include <string>
# include <cstring> 
# include <vector>
//# include <dos.h> //For Sleep() 
#include <typeinfo>

#include "TGraph.h"
#include "TAxis.h"
# include "TROOT.h"
# include "TFile.h"
# include "TTree.h"
# include "TBrowser.h"
# include "TH1.h"
# include "TH2.h"
# include "TH3.h"
# include "TRandom.h"

# include <TLorentzVector.h>

double myEnergy(double peak) {
	double A0 = 1.24703e+01;
	double A1 = 1.23646e+03;
	double energy = (peak - A0)/A1;
	return energy;
}

double mydEnergy(double peak, double FWHM) {
	double A0 = 1.24703e+01;
	double A1 = 1.23646e+03;
	double dA0 = 2.77963e+01;
	double dA1 = 7.81646e+01;	
	double denergy = sqrt(pow(FWHM/A1,2)+pow((-(peak-A0)*dA1)/pow(A1,2),2)
	+pow(-dA0/A1,2));
	return denergy;
}

// Linear background function
Double_t background(Double_t *x, Double_t *par) {
  return par[0] + par[1]*x[0];
}

using namespace std;
int main(){

  // controls
  const char *inputFileName = "Europium_152.Spe";
  bool WantToFit = true; // switch for fitting
  bool WantLog = true;
	
  const char *GraphTitle = "Europium";
  const char *GraphXLabel = "Channel Number";
  const char *GraphYLabel = "Counts";
	
  const char *FitFunc =  "gaus";//"[0]*x + [1]"; 
  const char *DrawOption = "AP" ; // see https://root.cern.ch/doc/master/classTGraphPainter.html
  int xmin = 0;
  int xmax = 0;

  // create the coordinate arrays
  vector<double> energy;
  vector<double> counts;

  ifstream inFile(inputFileName);

  if(inFile.is_open()){
    cout<<"Input File was opened successfully"<<endl;
  }

  //Main loop filling arrays from file//
  int tempEnergy = 0;
  double liney;
  while (inFile>>liney) {
    energy.emplace_back(tempEnergy);
    counts.emplace_back(liney);
	tempEnergy++;
	cout << "Energy: " << tempEnergy << " Counts: " << liney << endl;
  }
  cout << "Read " << energy.size() << " lines\n";
  xmax = energy.size() + 1;

  cout << "Creating canvas" << endl;
  // create the canvas
  TCanvas *c1 = new TCanvas("c1",GraphTitle,200,10,700,500);
  c1->SetGrid();
  
  if (WantLog){
	c1->SetLogy();
  }
  
  // create the TGraphErrors and draw it
  TGraphErrors *gr = new TGraphErrors(energy.size(),&energy[0],&counts[0]);
  
  // change title and axes titles here
  gr->SetTitle(GraphTitle);
  gr->GetXaxis()->SetTitle(GraphXLabel);
  gr->GetYaxis()->SetTitle(GraphYLabel);
  
  gr->SetMarkerColor(4);
  gr->SetMarkerStyle(21);
  gr->SetMarkerSize(0.2);
  
  // to set axis limits
  // TAxis *axis = gr->GetXaxis();
  // axis->SetLimits(x[0],x[5]);
  
  gr->Draw(DrawOption); 
   
  // Define Fit Function
  if (WantToFit){
    TF1 *f1 = new TF1("f1", FitFunc, 37, 70);
	TF1 *f2 = new TF1("f2", FitFunc, 93, 109);	
	TF1 *f3 = new TF1("f3", FitFunc, 113, 128);	
	TF1 *f4 = new TF1("f4", FitFunc, 143, 176);
	TF1 *f5 = new TF1("f5", FitFunc, 198, 242);
	TF1 *f6 = new TF1("f6", FitFunc, 290, 335);
	TF1 *f7 = new TF1("f7", FitFunc, 406, 465);
	TF1 *f8 = new TF1("f8", FitFunc, 948, 979);
	TF1 *f9 = new TF1("f9", "[18]*x + [19]", 35, 980);
	TF1 *total = new TF1("ftotal","gaus(0)+gaus(3)+gaus(6)+gaus(9)+gaus(12)+gaus(15)+gaus(18)+gaus(21)+[24]*x + [25]",35,980);
	
	f1->SetParNames("peak 1","position 1","RMS width 1");
	f2->SetParNames("peak 2","position 2","RMS width 2");
	f3->SetParNames("peak 3","position 3","RMS width 3");
	f4->SetParNames("peak 4","position 4","RMS width 4");
	f5->SetParNames("peak 5","position 5","RMS width 5");
	f6->SetParNames("peak 6","position 6","RMS width 6");
	f7->SetParNames("peak 7","position 7","RMS width 7");
	f8->SetParNames("peak 8","position 8","RMS width 8");
	f9->SetParNames("slope","constant");
	
	f9->SetParameters(-2,823);

	Double_t par[26];	
	
	// Fit each function and add it to the list of functions	
    gr->Fit(f1, "R");
	gr->Fit(f2, "R+");
	gr->Fit(f3, "R+");
	gr->Fit(f4, "R+");
	gr->Fit(f5, "R+");
	gr->Fit(f6, "R+");
	gr->Fit(f7, "R+");
	gr->Fit(f8, "R+");
	gr->Fit(f9, "R+");
    
   // Get the parameters from the fit
   f1->GetParameters(&par[0]);
   f2->GetParameters(&par[3]);
   f3->GetParameters(&par[6]);
   f4->GetParameters(&par[9]);
   f5->GetParameters(&par[12]);
   f6->GetParameters(&par[15]);
   f7->GetParameters(&par[18]);
   f8->GetParameters(&par[21]);
   f9->GetParameters(&par[24]);
   f9->SetParameters(-2,823);
	
   // The total is the sum of the three, each has 3 parameters
   total->SetParameters(par);
   gr->Fit(total,"R+");
   
    cout << "Fitting..." << endl;
	
	double peak1 = f1->GetParameter(0);
	double peakCh1 = f1->GetParameter(1);
	double sigma1 = f1->GetParameter(2);
	double FWHM1 = sigma1*(sqrt(8 * (log (2))));
	cout << "peak 1: " << peak1 << endl;
	cout << "FWHM for peak 1: " << FWHM1 << endl;
	cout << "energy for peak 1: " << myEnergy(peakCh1) << "  MeV with absolute error " << 
	mydEnergy(peak1, FWHM1) << endl;
	
	double peak2 = f2->GetParameter(0);
	double peakCh2 = f2->GetParameter(1);
	double sigma2 = f2->GetParameter(2);
	double FWHM2 = sigma2*(sqrt(8 * (log (2))));
	cout << "peak 2: " << peak2 << endl;
	cout << "FWHM for peak 2: " << FWHM2 << endl;
	cout << "energy for peak 2: " << myEnergy(peakCh2) << " MeV with absolute error " << 
	mydEnergy(peak2, FWHM2) << endl;
	
	double peak3 = f3->GetParameter(0);
	double peakCh3 = f3->GetParameter(1);
	double sigma3 = f3->GetParameter(2);
	double FWHM3 = sigma3*(sqrt(8 * (log (2))));
	cout << "peak 3: " << peak3 << endl;
	cout << "FWHM for peak 3: " << FWHM3 << endl;
	cout << "energy for peak 3: " << myEnergy(peakCh3) <<  " MeV with absolute error " << 
	mydEnergy(peak3, FWHM3) << endl;
	
	double peak4 = f4->GetParameter(0);
	double peakCh4 = f4->GetParameter(1);
	double sigma4 = f4->GetParameter(2);
	double FWHM4 = sigma4*(sqrt(8 * (log (2))));
	cout << "peak 4: " << peak4 << endl;
	cout << "FWHM for peak 4: " << FWHM4 << endl;
	cout << "energy for peak 4: " << myEnergy(peakCh4) << " MeV with absolute error " << 
	mydEnergy(peak4, FWHM4) << endl;
	
	double peak5 = f5->GetParameter(0);
	double peakCh5 = f5->GetParameter(1);
	double sigma5 = f5->GetParameter(2);
	double FWHM5 = sigma5*(sqrt(8 * (log (2))));
	cout << "peak 5: " << peak5 << endl;
	cout << "FWHM for peak 5: " << FWHM5 << endl;
	cout << "energy for peak 5: " << myEnergy(peakCh5) << " MeV with absolute error " << 
	mydEnergy(peak5, FWHM5) << endl;
	
	double peak6 = f6->GetParameter(0);
	double peakCh6 = f6->GetParameter(1);
	double sigma6 = f6->GetParameter(2);
	double FWHM6 = sigma6*(sqrt(8 * (log (2))));
	cout << "peak 6: " << peak6 << endl;
	cout << "FWHM for peak 6: " << FWHM6 << endl;
	cout << "energy for peak 6: " << myEnergy(peakCh6) << " MeV with absolute error " << 
	mydEnergy(peak6, FWHM6) << endl;
	
	double peak7 = f7->GetParameter(0);
	double peakCh7 = f7->GetParameter(1);
	double sigma7 = f7->GetParameter(2);
	double FWHM7 = sigma7*(sqrt(8 * (log (2))));
	cout << "peak 7: " << peak7 << endl;
	cout << "FWHM for peak 7: " << FWHM7 << endl;
	cout << "energy for peak 7: " << myEnergy(peakCh7) << " MeV with absolute error " << 
	mydEnergy(peak7, FWHM7) << endl; 
	
	double peak8 = f8->GetParameter(0);
	double peakCh8 = f8->GetParameter(1);
	double sigma8 = f8->GetParameter(2);
	double FWHM8 = sigma8*(sqrt(8 * (log (2))));
	cout << "peak 8: " << peak8 << endl;
	cout << "FWHM for peak 8: " << FWHM8 << endl;
	cout << "energy for peak 8: " << myEnergy(peakCh8) << " MeV with absolute error " << 
	mydEnergy(peak8, FWHM8) << endl; 
  }


  c1->Update();
   
  return 0;
}