#include "TH1.h"
#include "TSpectrum.h"
#include "TFile.h"
#include "TGraph.h"
#include "TF1.h"
#include "TMath.h"

#include <string>
#include <cmath>
#include <vector>
#include <iostream>
#include <map>
#include <sstream>


constexpr int eblow = 10; //bins
constexpr int rebinFac = 1; //
constexpr int maxnpeaks = 10; //max peaks to look for
constexpr double sigmaPeak = 4.0;
constexpr double thrPeak = 0.05;

constexpr double ene60Co1 = 1173.2;// keV
constexpr double ene60Co2 = 1332.5;// keV
constexpr double ene137Cs = 661.66;// keV

const std::vector<std::string> histoNames_v = { "EAn01",
						"EAn02",
						"EAn03",
						"EAn04",
						"EAn05",
						"EAn06",
						"EAn07",
						"EAn08",
						"EAn11",
						"EAn12",
						"EAn13",
						"EAn14",
						"EAn15",
						"EAn16",
						"EAn17",
						"EAn18",
						"EDy01",
						"EDy02",
						"EDy03",
						"EDy04",
						"EDy05",
						"EDy06",
						"EDy07",
						"EDy08",
						"EDy11",
						"EDy12",
						"EDy13",
						"EDy14",
						"EDy15",
						"EDy16",
						"EDy17",
						"EDy18"
						};

//const std::vector<std::string> histoNames_v = { "EDy18"	};


double ff (double *x,double*p){return p[0]*std::exp(-(*x - p[1])*(*x - p[1])/(2*p[2]*p[2]))+p[3];};
double ff2 (double *x,double*p){return p[0]*std::exp(-(*x - p[1])*(*x - p[1])/(2*p[2]*p[2]))+p[3]   + p[4]*ROOT::Math::erfc( (*x- p[1])/( sqrt(2)*p[2])) ;};



auto myTF1 = new TF1("myGaus", ff2 , 0 ,1 , 5 );

std::vector<TH1*> h_v;
std::vector<TH1*> hFit_v;

void getSigmaCal(std::string ifilename )
{
  TFile ifile( ifilename.c_str(), "read" );
  
  TH1 * h = nullptr;
  for( auto & histoName_it : histoNames_v )
  {
    ifile.GetObject( histoName_it.c_str(), h );
    if( h )
    {
      h->GetXaxis()->SetRangeUser(550,750);
      myTF1->SetParameters(1000,662,30,200,100);
      h->Fit("myGaus","");
      double mean(h->GetFunction("myGaus")->GetParameters()[1]);
      double sigma(h->GetFunction("myGaus")->GetParameters()[2]);
      
      std::cout << h->GetName() << "\t" << mean << "\t" << sigma << "\t" << 235*sigma/mean << "%\t" <<  std::endl;
      
      h_v.push_back( (TH1*) h->Clone( std::string( std::string(h->GetName()) + "_c" ).c_str() ) );
      h_v.back()->SetDirectory(0);
      hFit_v.push_back( (TH1*)h->GetFunction("myGaus")->CreateHistogram()->Clone( std::string( std::string(h->GetName()) + "_cFit" ).c_str() ) );
      hFit_v.back()->SetDirectory(0);
    }
    
  }
  ifile.Close();
  
  TFile ofile("calFit5.root","recreate");
  for(auto & vit : h_v )
    vit->Write();
//   for(auto & vit : hFit_v )
//     vit->Write();
  
  
  ofile.Close();
  
  return;
}