#include <fstream>
#include <string>
#include <vector>

void tFit(int tBin,int fitLevel){
  gROOT->Reset();//Reset session memory
  gROOT->SetStyle("Plain");//Make style plain
  gStyle->SetOptStat(0);//Remove the statistics box from histogram
  char fString[120];//Character string to hold file name
  char hId[120];

  //Function prototypes
  int readVec(char *fToGo,vector<double> *vec);




  sprintf(fString,"./all.root");//Copy filename to character string
  sprintf(hId,"hCasStarTval_%d",tBin);




  TFile *hFile1 = new TFile(fString);//Open file named ./all.root
  TH1D *hMass; hFile1->GetObject(hId,hMass);//get histogram



  hMass->SetName("hMass");
  hMass->GetXaxis()->SetTitle("Invariant Mass of #pi^{0} #Xi^{-} GeV/c^{2}");
  hMass->GetYaxis()->SetTitle("Counts");
  hMass->GetXaxis()->SetTitleSize(.04);
  hMass->GetYaxis()->SetTitleSize(.04);
  hMass->Rebin(4);
  
  TF1 *f1 = new TF1("f1","gaus");//Fit function of type gauss
  TF1 *f2 = new TF1("f2","gaus(0)+pol2(3)");//Fit function of type gauss
  // hMass->Fit("f1","R","O");
  f1->SetLineWidth(3);    
  f2->SetLineWidth(3);    

  double fitLow = 1.45;
  double fitHigh = 1.8;

  //READ FIT PARAMETERS
  vector<double> vecPars;
  sprintf(fString,"./fitPars/tFitPars_%d.dat",tBin);
  readVec(fString,&vecPars);
  int nPars = vecPars.size();
  cout<<"file = "<<fString<<endl;
  cout<<"nPars = "<<nPars<<endl;
  int rebin = 2;
  double centerSet = 1.535;
  double sigmaSet = -1;
  if (nPars == 5) {
    fitLow = vecPars.at(0);
    fitHigh = vecPars.at(1);
    rebin = (int) vecPars.at(2);
    centerSet = vecPars.at(3);
    sigmaSet = vecPars.at(4);

    if (centerSet < 0){
      f1->FixParameter(1,-centerSet);
      f2->FixParameter(1,-centerSet);
    } else {
      f1->SetParameter(1,centerSet);
      f2->FixParameter(1,-centerSet);
    }
    if (sigmaSet < 0){
      f1->FixParameter(2,-sigmaSet);
      f2->FixParameter(1,-centerSet);
      cout<<"WHAT"<<endl;
    } else {
      f1->SetParameter(2,sigmaSet);
      f2->FixParameter(1,-centerSet);
    }
  }
  cout<<"fitPars = "<<fitLow<<" , "<<fitHigh
      <<" , "<<rebin<<" , "<<centerSet<<" , "<<sigmaSet<<endl;

  //f1->FixParameter(0,0.1);
  //f1->FixParameter(1,0.1);
  //f1->FixParameter(2,0.1);
  hMass->Rebin(rebin);

  hMass->GetXaxis()->SetRangeUser(1.4,1.8);
  double fitLow1 = 1.5;
  double fitHigh1 = 1.56;
  hMass->Fit("f1","R","",fitLow1,fitHigh1);

  if (fitLevel == 1) return;

  // hCasV2->SetLineWidth(2);
  hMass->SetLineColor(1);
  //hMass->GetYaxis()->SetTitle("Counts");
  //hMass->GetXaxis()->SetTitle("Mass (GeV)");
  hMass->Draw("e");



  double center = f1->GetParameter(1);
  double sigma  = f1->GetParameter(2);
  f2->SetParameter(1,center);
  f2->FixParameter(2,sigma);
  
  hMass->Fit("f2","R","",fitLow,fitHigh);
 
  //Save the gaussian parameters of center and sigma
  //double center = fitFun->GetParameter(1);
  //double width  = fitFun->GetParameter(3)/2;//gamma value
 

  //NOTE the integral function wants bin numbers and not x-position.
  //It is a pain, but we do what we do
  int binLow = hMass->FindBin(center-3*sigma);
  int binHigh = hMass->FindBin(center+3*sigma);

  TH1D *hBK = (TH1D*) hMass->Clone("hBK");

  f2->SetParameter(0,0);
  hBK->Eval(f2);
  hBK->GetXaxis()->SetRangeUser(fitLow,fitHigh);
  hBK->SetLineWidth(2);
  hBK->SetLineColor(2);
  hBK->Draw("csame");
  
  //Integrate
  
  double signalYield = hMass->Integral(binLow,binHigh);
  double backgroundYield = hBK->Integral(binLow,binHigh);
  //double yieldError = sqrt(signalYield+2*backgroundYield);
  cout<<"Yield = "<<signalYield<<" +/- "<<sqrt(signalYield)<<endl;
  cout<<"Background Yield = "<<backgroundYield<<endl;
  //cout<<"Yield Error = "<<yieldError<<endl;
  if(tBin>0)
    {
      char tstring[120];
      sprintf(tstring, "t bin =%d GeV^{2}/c^{4} ",tBin);
      TLatex *t3 = new TLatex();
      t3->SetNDC();
      t3->SetTextSize(.05);
      double xLocationFraction =.1;
      double yLocationFraction =.93;
      t3->DrawLatex(xLocationFraction,yLocationFraction,tstring);

    }

  TFile *hOut = new TFile("myOutFile.root","RECREATE");
  hOut->cd();
  hBK->Write(); 

  hOut->Write();
  return;
}

int readVec(char *fToGo,vector<double> *vec)
{
  string line;
  stringstream os(line);
  string temp;
  ifstream myfile (fToGo);
  //if (myfile.fail() == false){ 
    if (myfile.is_open())
      {
	int j = 0;
	while ( myfile.good() )
	  {
	    int i = 0;
	    getline(myfile,line);
	    stringstream os(line);
	    while (os >> temp) {
	      vec->push_back(atof(temp.c_str()));
	      i++;
	    }
	    j++;
	  }
	myfile.close();
      }
    else cout << "Unable to open file :"<<endl;
    return 0;
    //}
}