#include "TFile.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TGaxis.h"
#include "TRandom.h"
#include "TLegend.h"
#include "TPaveStats.h"
#include "TGraph.h"
#include "TSystem.h"
#include "TTree.h"
#include "TTreePlayer.h"
#include "TF1.h"
#include "TCut.h"
#include "TPaletteAxis.h"
#include "TFitResult.h"
#include "TMath.h"


void muons()
{

	float fillcolor=0.35;
	gSystem->mkdir(TString::Format("simulazione" ), kTRUE);
	gStyle->SetOptFit();
	gStyle->SetFitFormat("8.6g");
	gStyle->SetOptStat(111110);
		ofstream results;
		results.open("simulazione/Results-energy.txt",  ios::out);
	TFile *fin = TFile::Open("B1/B1-build/B1.root");
	 	fin->ls ();
		TGaxis::SetMaxDigits(3);
 		if (fin == 0) {
	    	printf("Error: cannot open the file!\n");
	   	}
   		else {
			
			TTree *te=0;
			fin->GetObject("B1Exitparticles",te);
			gStyle->SetOptFit();
			gStyle->SetOptStat(111110);
			te->SetLineColor(kBlue);
			TH1::StatOverflows(true);
			
		 TCanvas *c01 = new TCanvas("c01","c01",1280,1024);
			te->Draw(Form("ExitParticleKinEn >> htemp(380, 20000., 23800.)"));
 			te->GetHistogram()->SetTitle("");
 		//	TPaveText *t=new TPaveText(0.75,0.6,0.8,0.55,"brNDC");
 			te->SetScanField(0);
 			int nentries = te->GetEntries();
  	 		TH1F *htempkinensp= (TH1F*)gPad->GetPrimitive("htemp"); 
  			htempkinensp->GetXaxis()->SetTitle("Kinetic Energy (MeV)"); 
			htempkinensp->GetYaxis()->SetTitle("Counts");
			htempkinensp->GetYaxis()->SetTitleSize(40);
  			htempkinensp->GetYaxis()->SetTitleFont(43);
   			htempkinensp->GetYaxis()->SetTitleOffset(1.1);
   			htempkinensp->GetYaxis()->SetLabelFont(43); 
   			htempkinensp->GetYaxis()->SetLabelSize(40);
   			htempkinensp->GetXaxis()->SetTitleSize(40);
   			htempkinensp->GetXaxis()->SetTitleFont(43);
   			htempkinensp->GetXaxis()->SetTitleOffset(1.1);
   			htempkinensp->GetXaxis()->SetLabelFont(43); 
   			htempkinensp->GetXaxis()->SetLabelSize(40);
   			htempkinensp->SetFillStyle(1001);
   			htempkinensp->SetFillColorAlpha(kBlue, fillcolor);
   			htempkinensp->SetStats(0);
   			double mean= htempkinensp->GetMean();
   			double stdv = htempkinensp->GetRMS();
	   		htempkinensp->Draw();
	   		htempkinensp->SetName("");
	   		gPad->Modified();
			gPad->Update(); // make sure it's really (re)drawn
   			c01->Update();
   			c01->SetLogy();
			gPad->Modified();
			gPad->Update(); // make sure it's really (re)drawn
   			c01->Update();
   			/*TLegend* legkinensp = new TLegend(0.25, 0.7, .35, .75);
   			legkinensp->SetHeader("Legend");
  			legkinensp->SetNColumns(1);
   			legkinensp->AddEntry(htempkinensp, "MC", "l");
   			legkinensp->Draw();*/ 
   			c01->Update();
   			gPad->Modified();
			gPad->Update();
			/*TPaveStats *statskinensp = (TPaveStats*)htempkinensp->GetListOfFunctions()->FindObject("stats");
			statskinensp->SetTextColor(kBlue);
			statskinensp->SetX1NDC(0.80); statskinensp->SetX2NDC(0.98);
			statskinensp->SetY1NDC(0.83); statskinensp->SetY2NDC(0.98);
			//statskinensp->AddText(TString::Format("Integral	= %g", integral));
			statskinensp->DrawClone();*/
			gPad->Update();
			gPad->Update();
			int binmax = htempkinensp->GetMaximumBin(); 
			double xMax = htempkinensp->GetXaxis()->GetBinCenter(binmax);
			int binmin = htempkinensp->GetMinimumBin(); 
			double xMin = htempkinensp->GetXaxis()->GetBinCenter(binmin);
			c01->Print("simulazione/spectrum_muons_Be30.png");
			results << "**************************************** "<< endl;
   			results << "Threshold result "<< endl;
   	 		results << "Entries= " << nentries <<endl;
   	 		results << "Mean= " << mean <<endl;
   	 		results << "Standard dev= " << stdv <<endl;
   	 		results << "Maximum= " << xMax <<endl;
   	 		results << "Minimum= " << xMin <<endl;
			delete c01;		
			
			
			
	}
}