#include "Riostream.h"
#include "TH1.h"
#include "TF1.h"

//R(t) Fitting Function
Double_t Fit1(Double_t *x, Double_t *par) {
           Double_t fit_val = par[0]*exp(-par[1]*x[0])+par[2]; //<====
           return fit_val;
}
void trio() {
  TFile *f = new TFile("basic.root", "read");
  TH1F *h1= (TH1F*)f->Get("h");
  h1->GetXaxis()->SetTitle("Energy (keV)");
  h1->GetYaxis()->SetTitle("Counts");

  h1->GetXaxis()->SetRangeUser(0,8191);
  h1->GetYaxis()->SetRangeUser(0,300);

  TH1* RAW=(TH1*) h1->Clone();
  

  RAW->SetDirectory(0);
  RAW->SetName("total_spectra");
// create a TF1 with the range from 10 to 1900 and 4 parameters
  TF1 *RT = new TF1("RT",Fit1,0,8191,3);
  RT->SetParameter(0,50.0); //<====
  RT->SetParameter(1,0.0004);
  RT->SetParameter(2,-3.0);
  RT->Draw();
  h1->Draw();
  h1->Fit("RT","R");
// Writing fuction fitted data in an histogram
    TH1F *h2 = new TH1F("h2","background",5212,-0.5,8190.5);
    h2->SetDirectory(0);
   for (int i=0;i<=8192;i++)  //<======
 {
   float x=h2->GetBinCenter(i);
   float y=RT->Eval(x);
   h2->SetStats(0);
   h2->Fill(x,y); //<======
 
   h2->GetXaxis()->CenterTitle();
   h2->GetYaxis()->CenterTitle();
   h2->GetXaxis()->SetTitle("Energy (keV)");
   h2->GetYaxis()->SetTitle("Counts");
   h2->SetLineWidth(2);
   h2->SetLineColor(2);
}
TString resultfile="back.root";
TFile* fresult=TFile::Open( resultfile.Data(), "RECREATE" );
   //fresult->cd();
   h2->Write();
   h2->SetStats(0);
   h2->GetXaxis()->CenterTitle();
   h2->GetYaxis()->CenterTitle();
   h2->GetXaxis()->SetTitle("Energy (keV)");
   h2->GetYaxis()->SetTitle("Counts");
   h2->GetXaxis()->SetRangeUser(0,8191);
   h2->GetYaxis()->SetRangeUser(0,300);

cout <<"Wrote result to "<<resultfile.Data() << endl;
delete fresult;

// difference of the two histograms
TH1* diff=(TH1*) RAW->Clone();
diff->SetDirectory(0);
diff->SetName("back_subtracted ");
diff->GetXaxis()->SetRangeUser(0,8191);
//diff->SetMaximum(3000.);//along Y
diff->GetYaxis()->SetRangeUser(0,300);
diff->Add(h2,-1.0);
diff->SetStats(0);
diff->GetXaxis()->CenterTitle();
diff->GetYaxis()->CenterTitle();
diff->GetXaxis()->SetTitle("Energy (keV)");
diff->GetYaxis()->SetTitle("Counts");

diff->GetXaxis()->CenterTitle();
diff->GetYaxis()->CenterTitle();
std::cout << diff->Integral(diff->FindFixBin(410.), diff->FindFixBin(422.)) << std::endl;
std::cout << diff->Integral(diff->FindFixBin(831.), diff->FindFixBin(843.)) << std::endl;
std::cout << diff->Integral(diff->FindFixBin(1172.), diff->FindFixBin(1184.)) << std::endl;
std::cout << diff->Integral(diff->FindFixBin(1512.), diff->FindFixBin(1526.)) << std::endl;
std::cout << diff->Integral(diff->FindFixBin(2656.), diff->FindFixBin(2673.)) << std::endl;
std::cout << diff->Integral(diff->FindFixBin(2958.), diff->FindFixBin(2976.)) << std::endl;
std::cout << diff->Integral(diff->FindFixBin(3288.), diff->FindFixBin(3307.)) << std::endl;
std::cout << diff->Integral(diff->FindFixBin(3794.), diff->FindFixBin(3814.)) << std::endl;
std::cout << diff->Integral(diff->FindFixBin(4804.), diff->FindFixBin(4826.)) << std::endl;
TString resultfile1="back_subtracted.root";
TFile* fresult1=TFile::Open( resultfile1.Data(), "RECREATE" );
diff->Write();
diff->Draw();
cout <<"Wrote result to "<<resultfile1.Data() << endl;
delete fresult1;
}