#include <stdio.h>
#include <TMath.h>
#include <TCanvas.h>
#include <iostream>
#include <TROOT.h>
#include <TChain.h>
#include <TObject.h>
#include <TRandom.h>
#include <TFile.h>
#include <math.h>
#include <TF1Convolution.h>
#include <TF1.h>
#include <TH1F.h>
#include <TGraph.h>
#include <TStopwatch.h>


void alpha1()
{
  //gStyle->SetOptFit(1);
  gROOT->ForceStyle();
  
  //----Channel-to-Energy Calibration----
  
  //Fill x values with the approximate channel of the peaks
  //Fill y values with the energy of the peaks (keV)
  double x1 = 1376; double y1 = 5155.0;//Pu 239
  double x2 = 1729; double y2 = 5486.0;//Am 241
  double x3 = 2068; double y3 = 5805.0;//Cm 244
  double n = 3.0;//Size of the array

  //Linear Regression calculation, no need to modify.
  double sumX = x1+x2+x3; double sumY = y1+y2+y3;
  double sumXY = (x1*y1)+(x2*y2)+(x3*y3);
  double sumXsq = TMath::Power(x1,2)+TMath::Power(x2,2)+TMath::Power(x3,2);
  double m = ((n*sumXY)-(sumX)*(sumY))/((n*sumXsq)-TMath::Power(sumX,2));//slope
  double b = ((sumY)*(sumXsq)-(sumX)*(sumXY))/((n*sumXsq)-TMath::Power(sumX,2));//offset
  
  //----Fit information----

  //Am241 85.2%
  float P = 5486.0;//Energy of the peak
  float dP = 25.0;
  double Pi = P - dP;
  double Pf = P + dP;
  //AM241 12.8%
  float P1 = 5443.0;
  float dP1 = 15.0;
  double Pi1 = P1 - dP1;
  double Pf1 = P1 + dP1;
  //Am241 1.4%
  float P3 = 5388.0;
  float dP3 = 15.0;
  double Pi3 = P3 - dP3;
  double Pf3 = P3 + dP3;

  
  //Gaussian Functions
  TF1 *g1 = new TF1("m1","gaus",Pi,Pf);
  TF1 *g2 = new TF1("m2","gaus",Pi1,Pf1);
  TF1 *g3 = new TF1("m3","gaus",Pi3,Pf3);
  //g1->SetParameters(87350,5486,7.5);

  
  //Lorentzian Function
  TF1 *l1 = new TF1("l1","(1*[2]/(TMath::Pi()*2))*((TMath::Power([0],2))/(TMath::Power(x-[1],2)+TMath::Power([0],2)))",Pi,Pf);
  //l1->SetParameters(0,dP,1000);

  //Pseudo-voigt
  TF1Convolution *f_conv = new TF1Convolution("m1","l1",Pi,Pf,true);
  f_conv->SetRange(Pi,Pf);
  f_conv->SetNofPointsFFT(1000);
  //TF1   *f = new TF1("f",*f_conv, Pi, Pf, f_conv->GetNpar());

  //Actual Voigt
  TF1 *f = new TF1("f","[0] * TMath::Voigt(x, [1], [2], 4)",Pi,Pf);

  
  
  //Global Fit
  //TF1 *total = new TF1("mstotal", "gaus(0)+gaus(3)+gaus(6)",5200,5750);


  //----Creating the Histograms----
  //The histograms are calibrated by including the linear regression
  auto h1=new TH1D("h1",
		  "Linear Scale;Energy (keV);counts",
		      4096, //number of bins
		      m*1+b, //x-axis min
		  m*4097+b);//x-axis max
  auto h2=new TH1D("h2",
		  "Log Scale;Energy (keV);counts",
		      4096, //number of bins
		      m*1+b, //x-axis min
		  m*4097+b);//x-axis max

  //----Filling the Histograms----
  std::ifstream inp; double x;
  std::string answer;
  std::cout << "Open default file? (y or n) ";
  std::cin >> answer;
  if (answer == "y") {
    inp.open("Am241pos6.TKA");
  }
  else if (answer == "n") { 
    std::string filename;
    std::cout << "Enter file name: ";
    std::cin >> filename;
    inp.open(filename);
  }
  else {
    std::cout << "Please type yes(y) or no(n)" << std::endl;
  }
 
  

  for (int i=0; i<4096; i++) {
    
   inp >> x;
   h1->SetBinContent(i,x);
   h2->SetBinContent(i,x);
  }

 
  //----Making the canvas for the histograms----
 Double_t par[12];
  
 TCanvas  *c1 = new TCanvas("c1","Am241 Alpha Detection", 1000,800);

 c1->cd(1);
 
 h1->GetYaxis()->SetRangeUser(0,100000);
 h1->GetXaxis()->SetRangeUser(5300,5560);
 g1->SetLineColor(kGreen); g1->SetLineStyle(0);
 //g2->SetLineColor(kMagenta); g2->SetLineStyle(0);
 //g3->SetLineColor(kOrange); g3->SetLineStyle(0);


 h1->Fit(g1,"RWM");
 //h1->Fit(g2,"NMWR+");
 //h1->Fit(g3,"NMWR+");
 h1->Fit(l1,"MWR+");
 //h1->Fit(f,"RWM+");

 
 //g1->GetParameters(&par[0]);
 //g2->GetParameters(&par[3]);
 //g3->GetParameters(&par[6]);
 //l1->GetParameters(&par[9]);
 //l1->GetParameters(&par[10]);
 
 
 //total->SetLineColor(kRed); 
 //total->SetParameters(par);
 //total->SetParameter(9,5);
 //total->SetParameter(10,0);
 //total->SetParameter(11,0);
 //total->SetParameter(12,0);
 //total->SetParameter(13,0);
 //h1->Fit(total,"MWR+");
 
 


 h1->Draw(); 
 
 
 //----Full width at half maximum----
 //fwhm of a gaussian distribution is = sigma*sqrt(8*ln(2))
}