#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 emitcont()
{
	ofstream results;
	results.open("C:/simulazione/emittancexxp_control_c=1.txt",  ios::out);
 //	results.open("C:/simulazione/emittanceyyp.txt",  ios::out);
		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);
		TH2F *htempdirxsp = new TH2F("htempdirxsp", "", 500, -0.02, 0.02,500, -0.0005,0.0005);
	//	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;
			double c=1;
			int jentry;
			int nentries = ts->GetEntries();
			double xv[10];
			double xVv[10];
			double xpv[10];
			double k=1.e12;
			int n=0;
			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(jentry=0; jentry<nentries;jentry++) {
		   // for(jentry=0; jentry<20;jentry++) {
				ts->GetEntry(jentry);
			//	if(xp+xpv[0] <= abs(c*xp)  && x+(k*xv[0]-k*xVv[0]) <= abs(c*x) ){
				if(xp+xpv[0] <= abs(c*xp) ){
				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];
				std::cout << "jentry = " << jentry << "\t xp0 = " << xpv[0] << "\t xp=" << xp<<  std::endl;
				n++;
				}
			}
			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;
			emittance=sigmax2*sigmaxp2-sigmaxsigmaxp;
			emittance=TMath::Sqrt(emittance);
			emittance/=1e6;
			std::cout << "N = " << nentries << "\t n = " << n <<  "\t c = " << c << "\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 <<  "N = " << nentries << "\t n = " << n <<  "\t c = " << c << "\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;
			//----------------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_c=1.pdf");
		//	c37->Print("C:/simulazione/emittanceyyp.pdf");
			results.close(); //close the file
			delete c37;	
		
}