#include <TTree.h>
#include <TFile.h>
#include <TF1.h>
#include <TH1D.h>
#include <TRandom.h>
#include <TMath.h>
#include <TNamed.h>
#include <TString.h>
#include <TCanvas.h>

#include <string>
#include <iostream>
#include <iomanip>
#include <fstream>

using namespace std;

void test(const char *input = "data.txt")
{
		
    gRandom->SetSeed(0);
	
	char line[100];
	
	int nenergy=29;
	double *xval=NULL;		//energy of the protons
	double *yval=NULL;		//integral flux
	double *xbin=NULL;		//array for variable bin size
	
	xval = new double[nenergy];
	yval = new double[nenergy];
	xbin = new double[nenergy+1];
	
	//Read protons data from file
	ifstream protfile(input);
    protfile.getline(line, 100); cout << line << endl;
    protfile.getline(line, 100); cout << line << endl;
    for (Int_t i=0;i<nenergy;i++) 
	{
		protfile >> xval[i] >> yval[i]; 
		cout << xval[i] << "    " << yval[i] << endl;
	}	
	protfile.close();
	
	double minf=yval[nenergy-1];
	double maxf=yval[0];
	double mine=xval[0];
	double maxe=xval[nenergy-1];
	
	for (Int_t i=0;i<nenergy;i++) {xbin[i]=xval[i];}	
    xbin[nenergy]=2*xbin[nenergy-1]-xbin[nenergy-2];
	
	TH1D *hfl = new TH1D("hfl","hfl",nenergy,xbin);
	for (Int_t i=0;i<nenergy;i++) {hfl->SetBinContent(i+1,yval[i]);}
	
    // To test that the random number generation reproduces the spectrum	
	TH1D *hrand = new TH1D("hrand","hrand",nenergy,xbin);	
	TH1D *hrand2 = new TH1D("hrand2","hrand2",nenergy,mine,maxe);	
	for (Int_t i=0;i<10000;i++) {hrand->Fill(hfl->GetRandom()); hrand2->Fill(hfl->GetRandom());}
	hrand->Scale(hfl->GetBinContent(1)/hrand->GetBinContent(1));
	hrand2->Scale(hfl->GetBinContent(1)/hrand2->GetBinContent(1)); 
	TCanvas *c1 = new TCanvas("c1","c1",400,400);
	c1->cd(1);
	gPad->SetLogx(); gPad->SetLogy();
	hfl->Draw("H");
	hrand->Draw("Hsame");
	hrand->SetLineColor(kRed);
	hrand2->Draw("Hsame");
	hrand2->SetLineColor(kBlue);
	c1->Modified();
	
	delete [] xval;		//energy of the protons
	delete [] yval;		//integral flux 
	delete [] xbin;				//array for variable bin size
	
	
}