#include "TFile.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TH2.h"
#include "TGaxis.h"
#include "TSystem.h"
#include "TTree.h"

void simlemmared()
{
		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 *c45 = new TCanvas("c45","emittance",1280,1024);
		 TH2F *hemitxxpsp = new TH2F("hemitxxpsp", "", 500, 0., 0.,5000, 0.,0.);
         TString varhemitxxpsp = TString::Format("SecondaryParticleAng.x : PositionDirection.x >> hemitxxpsp");
         ts->Draw(varhemitxxpsp);
         gPad->Modified();
         gPad->Update(); // make sure it's really (re)drawn
         hemitxxpsp->SetTitle("Emittance x-x'");
         hemitxxpsp->GetXaxis()->SetTitle("x");
         hemitxxpsp->GetYaxis()->SetTitle("x'");
         hemitxxpsp->Draw();
         hemitxxpsp->SetName("Emittance");
         gPad->Modified();
         gPad->Update(); // make sure it's really (re)drawn
			float x;	
			float xp;
			float xm;
			float xpm;
			float xxpmN;
			float xxpm;
			int jentry;
			int nentries = ts->GetEntries();
			xm= hemitxxpsp->GetMean(1);
			xpm= hemitxxpsp->GetMean(2);
			ts->SetBranchAdress("PositionDirection.x",PositionDirection.x,"x/I");
			ts->SetBranchAdress("SecondaryParticleAng.x",SecondaryParticleAng.x,"xp/I");
		    for(jentry=0; jentry<nentries;jentry++) {
			ts->GetEntry(jentry);
				xxpmN+= (x-xm)*(xp-xpm);
			}
			xxpm=xxpmN/nentries;
}