#include <cmath>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <ostream>
#include <stdlib.h>
#include <string>
#include <sstream>
#include <vector>
#include "Riostream.h"
#include "TFile.h"
#include "TFoam.h"
#include "TH1.h"
#include "TMath.h"
#include "TFoamIntegrand.h"
#include "TRandom3.h"
#include "TF1.h"
#include "TCanvas.h"
#include <TGraphErrors.h>

using namespace std;

float  Gaussian_Distribution(float x,float mean,float sigma)
{

  float norm = 1/sqrt(2*TMath::Pi())/sigma;
  float expo = (x-mean)             /sigma;
  float gdf  = norm * exp (- expo*expo /2);
  return gdf ;

}


Int_t main() {

  int nTot=0;
  float g1,g2,dummy;
  float mean_g1=0,mean_g2=0,sigma_g1=0,sigma_g2=0;
  vector<float> Gauss1,Gauss2; 

  const int nbin=15; float bmin=4.60; float step=0.07;
  float vbin[nbin+1];
  float g1_exp[nbin],g2_exp[nbin];
  float g1_err[nbin],g2_err[nbin];
  float vcen[  nbin];
  float verr[  nbin];

  for(int ib=0;ib<=nbin;ib++){vbin[ib]=bmin+     ib *step;}
  for(int ib=0;ib< nbin;ib++){vcen[ib]=bmin+(0.5+ib)*step;}

  TH1F  *g1h = new TH1F("g1h" ,  "Gauss1",nbin,vbin);
  TH1F  *g2h = new TH1F("g2h" ,  "Gauss2",nbin,vbin);

  // open file
  ifstream fileIN;
  string line;
  fileIN.open("GaussianHistogram.txt");

  // read data from file in a vector and sum elements
  while (fileIN) {
      fileIN >> g1>>g2>>dummy;
      if(!fileIN.eof()){
	nTot++;
	Gauss1.push_back(g1); Gauss2.push_back(g2);
	mean_g1+=g1         ; mean_g2+=g2         ;   
	g1h->Fill(g1)       ; g2h->Fill(g2)       ;
    }
  } fileIN.close();

  // calculate mean
  mean_g1/=float(nTot);
  mean_g2/=float(nTot);

  // calculate sigma
  for(int it=0;it<nTot;it++){
    sigma_g1+=(mean_g1-Gauss1[it])*(mean_g1-Gauss1[it]);
    sigma_g2+=(mean_g2-Gauss2[it])*(mean_g2-Gauss2[it]);
  }
  sigma_g1=sqrt(sigma_g1/float(nTot-1));
  sigma_g2=sqrt(sigma_g2/float(nTot-1));

  for(int it=0;it<nbin;it++){
    float vv=vcen[it];
    float f1=Gaussian_Distribution(vv,mean_g1,sigma_g1);
    float f2=Gaussian_Distribution(vv,mean_g2,sigma_g2);
    g1_exp[it]=nTot*step*f1;
    g1_err[it]=sqrt(nTot*(step*f1)*(1.0-step*f1));
    g2_exp[it]=nTot*step*f2;
    g2_err[it]=sqrt(nTot*(step*f2)*(1.0-step*f2));
  }

   //-----------------------------------------
  cout << "***** Start of Demonstration Program  *****\n" << endl;
  cout << " Total Number of Entries " << fixed 
       << setprecision(0) << setw(8)  << nTot
       << endl;
  
  cout << " Mean1 +/- Sig1 " << fixed
	<< setprecision(2) << setw(8) << mean_g1
	<< setprecision(2) << setw(8) << sigma_g1
	<< " \n"
	<< " Mean2 +/- Sig2 " 
	<< setprecision(2) << setw(8) << mean_g2
	<< setprecision(2) << setw(8) << sigma_g2
	<< endl;

  gStyle->SetOptFit(1111);
  TGraphErrors* a1 = new TGraphErrors(nbin,vcen,g1_exp,verr,g1_err);
  TGraphErrors* a2 = new TGraphErrors(nbin,vcen,g2_exp,verr,g2_err);
  a1->SetMarkerStyle(20);a1->SetMarkerColor(kRed);a1->SetLineColor(kRed);
  a2->SetMarkerStyle(20);a2->SetMarkerColor(kRed);a2->SetLineColor(kRed);

   g1h->SetMaximum(25);
   g2h->SetMaximum(25);
   TCanvas *c1 = new TCanvas("c1","c1",1200,100,1200,900);
   c1->Divide(2);
   c1->cd(1);
   g1h->Fit("gaus");
   a1->Draw("sameep");
   c1->cd(2);
   g2h->Draw();
   a2->Draw("sameep");

   //-----------------------------------------
   cout << "***** End of Demonstration Program  *****" << endl;
   //   usleep(10000000);
   return 0;
}