//ROOT stuff
#include <iostream>
#include "TLine.h"
#include "TAttMarker.h"
#include "TLegend.h"
#include "TGraph.h"
#include "TGraphErrors.h"
//#include "TMultiGraph.h"
#include "TCanvas.h"
#include "TLatex.h"
#include "TDirectory.h"
#include "TAxis.h"
#include "TPaveText.h"
//#include "TString.h"
#include "TFile.h"
#include "TH1F.h"
#include "TF1.h"
//////////////////////////////
//#include <iomanip>      // std::setprecision
//#include <sstream>      // std::istringstream or ostringstream
#include <string>       // std::string
//#include <vector>
#include <stdio.h>
#include <fstream>
#include <math.h>
#include <stdlib.h>
#include <cstdlib>
#include <cstdio>
#include <cmath>


using namespace std;
void cernsample(){


  const char *filename = "/home/ic00025/Desktop/NumuneSpectrum/SourceData/60Codata.txt";
  TGraph *g1 = new TGraph(filename, "%lg %*lg %lg"); // X axis = channels
  g1->SetTitle("g1;channel;counts");
  TGraph *g2 = new TGraph(filename, "%*lg %lg %lg"); // X axis = energies
  g2->SetTitle("g2;energy;counts");
  if ((g1->GetN() < 2) || (g2->GetN() < 2)) return; // just a precaution

  TF1 *fg = new TF1("fg", "gaus");
  Int_t ch_min = 360, ch_max = 460; // "channel numbers"
  TF1 *f1 = new TF1("f1", "gaus(0) + pol1(3)", ch_min, ch_max);
  f1->SetParameters(30000, 430, 10, 5000, 0);
  g1->Fit(f1, "WEMR");
  cout << "gaus peak integral (mean -+ 5 sigma) = "<< fg->Integral(f1->GetParameter(1) - 5.0 * f1->GetParameter(2), f1->GetParameter(1) + 5.0 * f1->GetParameter(2),  f1->GetParameters() ) << endl;
 
  
  TF1 *f2 = new TF1("f2", "gaus(0) + pol1(3)", g2->GetX()[ch_min], g2->GetX()[ch_max]);
  f2->SetParameters(30000, 1180, 25, 5000, 0);
  g2->Fit(f2, "WEMR");
  double s, sc;
  s = fg->Integral(f2->GetParameter(1) - 5.0 * f2->GetParameter(2), f2->GetParameter(1) + 5.0 * f2->GetParameter(2), f2->GetParameters());
  // note: sc = s / "average X bin width"
  sc = s  / ((g2->GetX()[ch_max] - g2->GetX()[ch_min]) / (ch_max - ch_min));
  cout << "gaus peak integral (mean -+ 5 sigma) = " << s<< " ( corrected = " << sc << " )" <<endl;
 
  TCanvas *c = new TCanvas("c", "c");
  c->Divide(1, 2);
  c->cd(1);
  g1->Draw("AP");
  c->cd(2);
  g2->Draw("AP");
  c->cd(0);
}