#include "Riostream.h"
#include "TFile.h"
#include "TH1.h"
#include "TNtuple.h"
#include "TF1.h"
#include "TSpectrum.h"
void fithisto() 
{
ifstream in;
in.open(Form("eu_gu_1.dat"));  


Int_t npeaks = 5, nbins = 8055;
Double_t source[nbins];
TFile*f= new TFile("basic.root","RECREATE");
TH1F *h1 = new TH1F("h1","x distribution",8055.,128.5,8183.5);
TH1F *back = (TH1F*) f->Get("back1");
TNtuple*ntuple= new TNtuple("ntuple","data from asciifile","x:y");

  TGraph *g = new TGraph("eu_gu_1.dat");
  // g->Sort(); // just a precaution
  Int_t n = g->GetN();//get the no. of points on the TGraph *g
  Double_t *x = new Double_t[(n + 1)];
  x[0] = (3.0 * g->GetX()[0] - g->GetX()[1]) / 2.0;
  x[n] = (3.0 * g->GetX()[(n - 1)] - g->GetX()[(n - 2)]) / 2.0;
  for (Int_t i = 1; i < n; i++)
    { x[i] = (g->GetX()[(i - 1)] + g->GetX()[i]) / 2.0; }
  TH1D *h = new TH1D("h", "spectrum;x-value(Channels);#counts", n, x);
  h->FillN(n, g->GetX(), g->GetY());
  h->Sumw2(kFALSE); // restore proper errors
  delete [] x; // no longer needed
  delete g; // no longer needed
  TF1 *f1 = new TF1("f1", "gaus(0) + pol1(3)", 400., 435.); // linear background
  f1->SetParameters(2800., 415., 2., 0., 0.);
  h->Fit("f1", "BR+");
  TF1 *f2 = new TF1("f2", "gaus(0) + pol1(3)", 830., 846.); // linear background
  f2->SetParameters(1000., 838., 2., 0., 0.);
  h->Fit("f2", "BR+");
  TF1 *f3 = new TF1("f3", "gaus(0) + pol1(3)", 1160., 1200.);
  f3->SetParameters(2700., 1178., 2., 0., 0.);
  h->Fit("f3", "BR+");
  TF1 *f4 = new TF1("f4", "gaus(0) + pol1(3)", 1512., 1526.); // clean peak, no background
  f4->SetParameters(280., 1519., 3., 0., 0.);
  h->Fit("f4", "BR+");
  TF1 *f5 = new TF1("f5", "gaus(0) + pol1(3)", 2654., 2674.); // clean peak, no background
  f5->SetParameters(580., 2664., 3., 0., 0.);
  h->Fit("f5", "BR+");
  TF1 *f6 = new TF1("f6", "gaus(0) + pol1(3)", 2954., 2978.); // clean peak, no background
  f6->SetParameters(200., 2966., 3., 0., 0.);
  h->Fit("f6", "BR+");
  TF1 *f7 = new TF1("f7", "gaus(0) + pol1(3)", 3280., 3316.); // clean peak, no background
  f7->SetParameters(540., 3296., 3., 0., 0.);
  h->Fit("f7", "BR+");
  TF1 *f8 = new TF1("f8", "gaus(0) + pol1(3)", 3790., 3820.); // clean peak, no background
  f8->SetParameters(420., 3805., 3., 0., 0.);
  h->Fit("f8", "BR+");
  TF1 *f9 = new TF1("f9", "gaus", 4795., 4835.); // clean peak, no background
  f9->SetParameters(500., 4815., 3.);
  h->Fit("f9", "BR+"); // e.g. "BR+" or "LBR+"


//Area  
  std::cout << h->Integral(h->FindFixBin(400.), h->FindFixBin(435.)) << std::endl;
  std::cout << h->Integral(h->FindFixBin(830.), h->FindFixBin(846.)) << std::endl;
  std::cout << h->Integral(h->FindFixBin(1160.), h->FindFixBin(1200.)) << std::endl;
  std::cout << h->Integral(h->FindFixBin(1512.), h->FindFixBin(1526.)) << std::endl;
  std::cout << h->Integral(h->FindFixBin(2654.), h->FindFixBin(2674.)) << std::endl;
  std::cout << h->Integral(h->FindFixBin(2954.), h->FindFixBin(2978.)) << std::endl;
  std::cout << h->Integral(h->FindFixBin(3280.), h->FindFixBin(3316.)) << std::endl;
  std::cout << h->Integral(h->FindFixBin(3790.), h->FindFixBin(3820.)) << std::endl;
  std::cout << h->Integral(h->FindFixBin(4795.), h->FindFixBin(4835.)) << std::endl;

   
  TSpectrum *s = new TSpectrum(2*npeaks);
  Int_t nfound = s->Search(h,2,"",0.06);


  
  for (int i = 0; i < nbins; i++) source[i] = back->GetBinContent(i + 1);
  // Estimate the background
  s->Background(source,nbins,6,TSpectrum::kBackIncreasingWindow,
                 TSpectrum::kBackOrder2,kFALSE,
                 TSpectrum::kBackSmoothing3,kFALSE);
  // Draw the estimated background
   for (int i = 0; i < nbins; i++) h->SetBinContent(i + 1,source[i]);
   h->SetLineColor(kBlack);
   h->Draw("SAME L");

   


in.close();
f->Write();
}