#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>

void emit()
{
	ofstream results;
	results.open("C:/simulazione/emittancexxp_control_N=3.txt",  ios::out);
 //	results.open("C:/simulazione/emittanceyyp.txt",  ios::out);
		char Secondary[50];
		int n=1;
			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);
		//TH2F *htempdirxsp = new TH2F("htempdirxsp", "", 500, -0.02, 0.02,500, -0.0005,0.0005);
		TH2F *htempdirxsp = new TH2F("htempdirxsp", "", 50, -2, 2,50, -2,2);
	//	TString varhemitxxpsp = TString::Format("SecondaryParticleAng.x: PositionDirection.x >> htempdirxsp");
		 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();
			//----------------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();
			double xv[10];
			double xVv[10];
			double xpv[10];
			double k=1e12;
			int N=3;
			ts->SetBranchAddress("SecondaryParticleVert.x", &xVv[0]);
			ts->SetBranchAddress("PositionDirection.x", &xv[0]);
			ts->SetBranchAddress("SecondaryParticleAng.x", &xpv[0]);
			double previous_x = k*xv[0]-k*xVv[0];
			double previous_x2= (k*xv[0]-k*xVv[0])*(k*xv[0]-k*xVv[0]);
			double previous_xp=xpv[0];
			double previous_xp2=xpv[0]*xpv[0];
			double previous_xxp=(k*xv[0]-k*xVv[0])*xpv[0];
			double sum_x = previous_x;
			double sum_x2 = previous_x2;
			double sum_xp = previous_xp;
			double sum_xp2 = previous_xp2;
			double sum_xxp = previous_xxp;
		   for (int i = 1; i <nentries; i++) {
   			if (( k*xv[i]-k*xVv[i])<previous_x/N  && ((k*xv[i]-k*xVv[i])*(k*xv[i]-k*xVv[i]))<previous_x2/N && xpv[i]<previous_xp/N && (xpv[i]*xpv[i])<previous_xp2/N && ((k*xv[i]-k*xVv[i])*xpv[i])<previous_xxp/N) {
      		 previous_x = k*xv[i]-k*xVv[i];
			 previous_x2= (k*xv[i]-k*xVv[i])*(k*xv[i]-k*xVv[i]);
			 previous_xp=xpv[i];
			 previous_xp2=xpv[i]*xpv[i];
			previous_xxp=(k*xv[i]-k*xVv[i])*xpv[i];
      		 sum_x +=previous_x;
      		 sum_x2 +=previous_x2;
      		 sum_xp +=previous_xp;
      		 sum_xp2 +=previous_xp2;
      		 sum_xxp +=previous_xxp;
      			n++;
   				} 
   			}   			
			sum_x*=1./n;
			sum_x2*=1./n;
			sum_xp*=1./n;
			sum_xp2*=1./n;
			sum_xxp*=1./n;
			sigmax2=sum_x2-sum_x*sum_x;
			sigmaxp2=sum_xp2-sum_xp*sum_xp;
			sigmaxsigmaxp=sum_xxp-sum_x*sum_xp;
			emittance=sigmax2*sigmaxp2-sigmaxsigmaxp;
			emittance=TMath::Sqrt(emittance);
			emittance/=1e6;
			std::cout << "N = " << n << "\t x = " << sum_x << "\t x2 = " << sum_x2 << " \t xp = " << sum_xp <<  " \t xp2= " << sum_xp2 << " \t xxp= " <<  sum_xxp << std::endl;
			std::cout << "sigmax2 = " << sigmax2 << "\t sigmaxp2 = " << sigmaxp2 << "\t sigmaxsigmaxp = " << sigmaxsigmaxp << " \t emittance = " << emittance  << std::endl;
		    results << "N = " << n << "\t x = " << sum_x << "\t x2 = " << sum_x2 << " \t xp = " << sum_xp <<  " \t xp2= " << sum_xp2 << " \t xxp= " <<  sum_xxp << endl;
		    results << "sigmax2 = " << sigmax2 << "\t sigmaxp2 = " << sigmaxp2 << "\t sigmaxsigmaxp = " << sigmaxsigmaxp << " \t emittance = " << emittance  << endl;
			//----------------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 (nm*rad) = %g", emittance));
			statsdirxsp->DrawClone();
			gPad->Update();
			c37->Print("C:/simulazione/emittancexxp_control_N=3.pdf");
		//	c37->Print("C:/simulazione/emittanceyyp.pdf");
			results.close(); //close the file
			delete c37;	
		
}