#include<iostream>
#include<cmath>
#include<fstream>
#include<string>
#include<TCanvas.h>
#include <TFile.h>
#include<TApplication.h>
#include<TGraphErrors.h>
#include<TH1D.h>
#include<TF1.h>
#include<TStyle.h>
#include<TAxis.h>
#include<TLatex.h>
#include<TMath.h>
#include<TTree.h>
#include<cstdlib>
#include<stdlib.h>
#include<TRandom.h>
#include "TH2F.h"
#include "TCanvas.h"
#include "TPad.h"
#include<cstring>
#include<libgen.h>   //serve per basename
#include<sstream>  // serve per stringstream
#include<iomanip>
#include <stdio.h>

using namespace std;

const int ADC = 8;
const int Ndet = 6;
string Dir = "/Root/";
TRandom *r = new TRandom();

double SIGMA(double E, int j) {
	double A[ADC] = {0, 9.60227, 18.4655, 11.6821, 24.0247, 14.0132, 6.0811, 0};
	double B[ADC] = {1, 1.11395, 0.78216, 1.14578, 0.44079, 1.12809, 1.33874, 1};
	double C[ADC] = {0, 0.00786, 0.01097, 0.01077, 0.01095, 0.00927, 0.00577, 0};
	double sigma = A[j] + B[j]*sqrt(E) + C[j]*E;
	return sigma;
}
  void anaprin(){
  	
	string RUN = "1024";
	string FileIn =  "run"+RUN+"_coinc.root";;
	string FileWithoutExtension = FileIn.substr(0, FileIn.rfind("."));
	string FileOut = FileWithoutExtension+".hist";
	FileIn = Dir + FileIn;
	double Emin = 0;
	double Emax = 9000;
	double dE = 1;
	double Ethr = 10;
	int Nbin = (int)((Emax-Emin)/dE);
  	TFile *RootIn= TFile::Open(FileIn.c_str(),"READ");
  	TTree *hIn = (TTree*) (RootIn->Get("T_coinc"));
  	TFile *RootOut = new TFile(FileOut.c_str(),"RECREATE");

	vector<Float_t> *E = new vector <Float_t>();
	vector<Long64_t> *T = new vector <Long64_t>();
	int mult;
	
	E=0;
	T=0;

	hIn->SetBranchAddress("Energy", &E);
	hIn->SetBranchAddress("Timestamp", &T);
	hIn->SetBranchAddress("Multiplicity", &mult);

	TH1D *Hraw[Ndet];
	TH1D *H[Ndet];

	string Name[Ndet], NameRaw[Ndet];
	ostringstream ss[ADC];
	for (int j=1; j<=Ndet; j++) {
		ss[j] << j;
		Name[j] = "BGO"+ss[j].str();
		NameRaw[j] = "BGO"+ss[j].str()+"raw";
	}
  	for (int j=1; j<=Ndet; j++) {
		Hraw[j] = new TH1D(NameRaw[j].c_str(), NameRaw[j].c_str(), Nbin, Emin, Emax);
		H[j] = new TH1D(Name[j].c_str(), Name[j].c_str(), Nbin, Emin, Emax);
  	}
	TH1D *HANDraw = new TH1D("ANDraw", "ANDraw", Nbin, Emin, Emax);
	TH1D *HAND = new TH1D("AND", "AND", Nbin, Emin, Emax);
	long int NEvents = hIn->GetEntries();
	double Stat = (double)NEvents;
	cout << "Number of Events in " << FileIn <<" : " << NEvents << endl;
	long int Step = (long int)(NEvents/10);
	double ROImin = 7000;
	double ROImax = 8200; 
	double sigma;
	for(long int i=0; i<NEvents; i++) {
		if (i%Step==0) cout << (int)(i*1e2/NEvents) << "% events processed.." << endl;
		hIn->GetEntry(i);
		double Esum = 0;
		double ES = 0;
		if (E->size()>0) {
			for (int j=0; j<Ndet; j++) {
				if (E->at(j)<0) E->at(j)=0; 
				Esum = Esum + E->at(j);
			}
			if (Esum>Ethr) HANDraw->Fill(Esum);
			for (int j=1; j<=Ndet; j++) {
				if (E->at(j)>Ethr) Hraw[j]->Fill(E->at(j));
				double E1 = ROImin + r->Gaus(0, SIGMA(ROImin,j));
				double E2 = ROImax + r->Gaus(0, SIGMA(ROImax,j));		
				if(Esum>E1 && Esum<E2 && E->at(j)>0) {
					H[j]->Fill(E->at(j));
					ES = ES + E->at(j);
				}
			}
			if (ES>0) HAND->Fill(ES);
		}
	}
	TH1D *HrawOR = (TH1D*)Hraw[1]->Clone();
	HrawOR->SetName("ORraw");
	HrawOR->SetTitle("ORraw");
	for (int j=2; j<=Ndet; j++) HrawOR->Add(Hraw[j]);
	TH1D *HOR = (TH1D*)H[1]->Clone();
	HOR->SetName("OR");
	HOR->SetTitle("OR");
	for (int j=2; j<=Ndet; j++) HOR->Add(H[j]);
	RootIn->Close();
	RootOut->Write();
	RootOut->Close();
	cout<<endl;
	cout<<"Results on file: "<<FileOut<<endl;
	cout<<endl;
	delete E;
	delete T;
}