#include "TFile.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TH2.h"
#include "TGaxis.h"
#include "TSystem.h"
#include "TTree.h"
#include "TLegend.h"
#include "TPaveStats.h"
#include "TMath.h"
#include <iostream>

void emit()
{
		char Secondary[50];
			sprintf(Secondary, "B1Secondaryparticles");
	TFile *fin = TFile::Open("C:/B1.root");
	 	fin->ls ();
		TGaxis::SetMaxDigits(3);
			TTree *ts=0;
			fin->GetObject(Secondary,ts);
			gStyle->SetOptFit();
			gStyle->SetOptStat(111110);
			ts->SetLineColor(kBlue);
			
		TCanvas *c37 = new TCanvas("c37","Plot",1280,1024);
		TH2F *htempdirxsp = new TH2F("htempdirxsp", "", 500, -2., 2.,500, 1.5705,1.571);
		 TString varhemitxxpsp = TString::Format("SecondaryParticleAng.x : PositionDirection.x >> htempdirxsp");
         ts->Draw(varhemitxxpsp);
 			ts->GetHistogram()->SetTitle("Primary e^{#pm} emittance ");
 			ts->SetScanField(0);
  			htempdirxsp->GetXaxis()->SetTitle("x (mm)"); 
			htempdirxsp->GetYaxis()->SetTitle("x' (rad)");
   			htempdirxsp->SetFillStyle(1001);
   			htempdirxsp->SetFillColorAlpha(kBlue, 0.35);
	   		htempdirxsp->Draw("COLZ");
	   		gPad->Modified();
			gPad->Update(); // make sure it's really (re)drawn
   			TLegend* legdirxsp = new TLegend(0.25, 0.7, .35, .75);
   			legdirxsp->SetHeader("Legend");
  			legdirxsp->SetNColumns(1);
   			legdirxsp->AddEntry(htempdirxsp, "MC", "l");
   			legdirxsp->Draw(); 
   			c37->Update();
   			gPad->Modified();
			gPad->Update();
			//----------------EMITTANCE CALCULATION-----------------//
			Double_t emittance=0.;
			Double_t x=0;
			Double_t x2=0;	
			Double_t xp=0;
			Double_t xp2=0;
			Double_t xxp=0;
			Double_t sigmax2=0;
			Double_t sigmaxp2=0;
			Double_t sigmaxsigmaxp=0;
			int jentry;
			int nentries = ts->GetEntries();
			ts->SetBranchAddress("PositionDirection.x", &x);
			ts->SetBranchAddress("SecondaryParticleAng.x", &xp);
		    for(jentry=0; jentry<nentries;jentry++) {
				ts->GetEntry(jentry);
				x+=x;
				x2+=x*x;
				xp+=xp;
				xp2+=xp*xp;
				xxp+=x*xp;
			}
			x*=1./nentries;
			x2*=1./nentries;
			xp*=1./nentries;
			xp2*=1./nentries;
			xxp*=1./nentries;
			sigmax2=x2-x*x;
			sigmaxp2=xp2-xp*xp;
			sigmaxsigmaxp=xxp-x*xp;
			std::cout << "N = " << nentries << "\t x = " << x << "\t x2 = " << 2 << " \t xp = " << xp <<  " \t xp2= " << xp2 << " \t xxp= " <<  xxp << std::endl;
			emittance=sigmax2*sigmaxp2-sigmaxsigmaxp;
			emittance=TMath::Sqrt(emittance);
			emittance*=1000000;
			//----------------EMITTANCE CALCULATION-----------------//
			TPaveStats *statsdirxsp = (TPaveStats*)htempdirxsp->GetListOfFunctions()->FindObject("stats");
			statsdirxsp->SetTextColor(kBlue);
			statsdirxsp->SetX1NDC(0.80); statsdirxsp->SetX2NDC(0.98);
			statsdirxsp->SetY1NDC(0.77); statsdirxsp->SetY2NDC(0.92);
			statsdirxsp->AddText(TString::Format("emittance = %g", emittance));
			statsdirxsp->DrawClone();
			gPad->Update();
			c37->Print("C:/simulazione/emittancexxp.pdf");
			delete c37;	
		
}