#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>
#include <vector>
#include <fstream>
#include "TRandom.h"
void emiter()
{
	ofstream results;
	ofstream resultsthf;
 	results.open("emittanceerxxp.txt",  ios::out);
 	//results.open("C:/simulazione/emittanceeryyp.txt",  ios::out);
 	resultsthf.open("emittanceerxxp_2.txt",  ios::out);
 	//resultsthf.open("C:/simulazione/emittanceeryyp_2.txt",  ios::out);
		char Secondary[50];
			sprintf(Secondary, "B1Secondaryparticles");
	TFile *fin = TFile::Open("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);
		TH2F *htempdirxsp = new TH2F("htempdirxsp", "", 500, 0., 0.,500, 0.,0.);
		 TString varhemitxxpsp = TString::Format("SecondaryParticleAng.x : (PositionDirection.x-SecondaryParticleVert.x) >> htempdirxsp");
		 // TString varhemitxxpsp = TString::Format("SecondaryParticleAng.y : PositionDirection.y-SecondaryParticleVert.y >> htempdirxsp");
		 TH2::StatOverflows(true);
         ts->Draw(varhemitxxpsp);
 			//ts->GetHistogram()->SetTitle("Primary e^{+} emittance");
 			ts->GetHistogram()->SetTitle("Secondary #mu^{#pm} emittance");
 			ts->SetScanField(0);
  			htempdirxsp->GetXaxis()->SetTitle("x (mm)"); 
			htempdirxsp->GetYaxis()->SetTitle("x' (rad)");
		//	htempdirxsp->GetXaxis()->SetTitle("y (mm)"); 
		//	htempdirxsp->GetYaxis()->SetTitle("y' (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.15, 0.6, .25, .65);
   			legdirxsp->SetHeader("Legend");
  			legdirxsp->SetNColumns(1);
   			legdirxsp->AddEntry(htempdirxsp, "MC", "l");
   			legdirxsp->Draw(); 
   			c37->Update();
   			gPad->Modified();
			gPad->Update();
			 TCanvas *c0 = new TCanvas("c0","c0",1280,1024);
            TH1F *hscoring = new TH1F("hscoring","Released energy",100,0.,0.);
			//----------------EMITTANCE CALCULATION-----------------//
			int i;
			int Nextraction=100;
			for(i=0; i<Nextraction;i++) {
			Double_t emittance=0.;
			Double_t emittance2=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 nentries = ts->GetEntries();
			double xv[10];
			double xVv[10];
			double xpv[10];
			int k=1000000000;
			int entry;
			ts->SetBranchAddress("SecondaryParticleVert.x", &xVv[0]);
			ts->SetBranchAddress("PositionDirection.x", &xv[0]);
			ts->SetBranchAddress("SecondaryParticleAng.x", &xpv[0]);
		     //ts->SetBranchAddress("SecondaryParticleVert.y", &xVv[0]);
			//ts->SetBranchAddress("PositionDirection.y", &xv[0]);
			//ts->SetBranchAddress("SecondaryParticleAng.y", &xpv[0]);
		   for(int jentry=0; jentry<nentries;jentry++) {
			entry = gRandom->Integer(ts->GetEntries());
				ts->GetEntry(entry);
				x+=k*xv[0]-k*xVv[0];
				x2+=(k*xv[0]-k*xVv[0])*(k*xv[0]-k*xVv[0]);
				xp+=xpv[0];
				xp2+=xpv[0]*xpv[0];
				xxp+=(k*xv[0]-k*xVv[0])*xpv[0];
			/*	x+=k*xv[0];
				x2+=k*xv[0]*k*xv[0];
				xp+=xpv[0];
				xp2+=xpv[0]*xpv[0];
				xxp+=k*xv[0]*xpv[0];*/
			}
			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;
			emittance2=sigmax2*sigmaxp2-sigmaxsigmaxp;
			emittance=TMath::Sqrt(emittance2);
			emittance/=1000;
		    std::cout << "Extraction = " << i << "\t x = " << x << "\t x2 = " << x2 << " \t xp = " << xp <<  " \t xp2= " << xp2 << " \t xxp= " <<  xxp << std::endl;
		    std::cout << "sigmax2 = " << sigmax2 << "\t sigmaxp2 = " << sigmaxp2 << "\t sigmaxsigmaxp = " << sigmaxsigmaxp << " \t emittance = " << emittance  << std::endl;
		    results << "Extraction = " << i << "\t x = " << x << "\t x2 = " << x2 << " \t xp = " << xp <<  " \t xp2= " << xp2 << " \t xxp= " <<  xxp << endl;
		    results << "sigmax2 = " << sigmax2 << "\t sigmaxp2 = " << sigmaxp2 << "\t sigmaxsigmaxp = " << sigmaxsigmaxp << " \t emittance = " << emittance  << endl;
		  	if (emittance2>=0){
		    resultsthf << emittance <<  endl;
		    	hscoring->Fill(emittance);
	}
}
			//----------------EMITTANCE CALCULATION-----------------//
			results.close(); //close the file
			resultsthf.close(); //close the file
            hscoring->SetTitle("Emittance");
   			hscoring->GetXaxis()->SetTitle("Emittance (nm*rad)");
   			hscoring->GetYaxis()->SetTitle("Counts");
   			hscoring->Draw();
   				gPad->Modified();
			gPad->Update(); // make sure it's really (re)drawn
   			TPaveStats *statshscoring = (TPaveStats*)hscoring->GetListOfFunctions()->FindObject("stats");
			statshscoring->SetTextColor(kBlue);
			statshscoring->SetX1NDC(0.80); statshscoring->SetX2NDC(0.98);
			statshscoring->SetY1NDC(0.77); statshscoring->SetY2NDC(0.92);
			c0->Update();
   			gPad->Modified();
			gPad->Update();
			c0->Print("emittanceer_xxp.pdf");
			delete c0;	
			delete c37;	
}