{
  // Macro loads a file after reconstruction and plots difference between initial direction of particle and angular position of cluster
  
	gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
	basiclibs();
	gSystem->Load("libGeoBase");
	gSystem->Load("libParBase");
	gSystem->Load("libBase");
	gSystem->Load("libMCStack");
	gSystem->Load("libField");
	gSystem->Load("libGen");
	gSystem->Load("libPassive");
	// add the detector libraries you need
	gSystem->Load("libEmc");	 
	gSystem->Load("libGeom.so");
	
	TFile* f = new TFile("full_emc_5-22.root"); //file you want to analyse
	TTree *t=(TTree *) f->Get("cbmsim") ;
	TClonesArray* cluster_array=new TClonesArray("EmcCluster");
	t->SetBranchAddress("EmcCluster",&cluster_array);
	
	TFile* fsim = new TFile("sim_emc_5-22.root"); //file you want to analyse
	TTree *tsim=(TTree *) fsim->Get("cbmsim") ;
	TClonesArray* mctrack_array=new TClonesArray("CbmMCTrack");
	tsim->SetBranchAddress("MCTrack",&mctrack_array);

	TVector3 photon_momentum;

	double cluster_energy;
	double cluster_theta, cluster_phi; //position of the cluster
	double theta, phi; // angular position of the initial particle
	double theta_diff, phi_diff;
	int ndigi, npoint;
	double max_energy=0;
	
	TH1F *h1= new TH1F("h2","Theta difference",100,-5.,5.);
	TH1F *h2= new TH1F("h3","Phi difference",100,-2.,2.);
	TH1F *h3= new TH1F("h1","Cluster energy",100,0.85,1.05);
	
	// Cluster angular position
	// Entrance point is determined by minimal time
		
	// Cluster energy
	for (Int_t j=0; j< t->GetEntriesFast(); j++)
	{
		t->GetEntry(j);
		for (Int_t i=0; i<cluster_array->GetEntriesFast(); i++)
		{
			EmcCluster *cluster=(EmcCluster*)cluster_array->At(i);
			cluster_energy=cluster->energy();
			if ((cluster->numberOfDigis()>1)&&(cluster_energy>0.02))
			h3->Fill(cluster_energy);
		}
	}

	for (Int_t j=0; j< t->GetEntriesFast(); j++)//t->GetEntriesFast()
	{
		t->GetEntry(j);
		tsim->GetEntry(j);
		
		CbmMCTrack *mctrack=(CbmMCTrack *) mctrack_array->At(0);
		photon_momentum=mctrack->GetMomentum();
		theta=photon_momentum.Theta();
		phi=photon_momentum.Phi();
		
	
		// Loop over clusters
		// If we have 1 initial particle and several cluster
		// we can separate cluster from the first interaction by maximum energy
		
		for (Int_t i=0; i<cluster_array->GetEntriesFast(); i++)
		{
			EmcCluster *cluster=(EmcCluster*)cluster_array->At(i);
			cluster_energy=cluster->energy();
			if (cluster_energy>max_energy)
			{
				max_energy=cluster_energy;
				TVector3 cluster_pos=cluster->where();
				cluster_theta=cluster_pos.Theta();
				cluster_phi=cluster_pos.Phi();
			}
						
		}
		max_energy=0;
		
		
		theta_diff=(cluster_theta-theta)*180./TMath::Pi();
		h1->Fill(theta_diff);

		phi_diff=(cluster_phi-phi)*180./TMath::Pi();
		h2->Fill(phi_diff);
		
	}

	TCanvas* c1 = new TCanvas("c1", "", 100, 100, 800, 800); 	
	h3->SetTitle("Cluster energy of 1 GeV photon");
	h3->GetXaxis()->SetTitle("Energy, GeV");
	h3->Draw();

	
	TCanvas* c2 = new TCanvas("c2", "", 100, 100, 800, 800); 	
	h1->SetTitle("Difference between cluster and initial photon theta angle");
	h1->GetXaxis()->SetTitle("#theta_{reco} - #theta_{truth}, degree");
	h1->Draw();

// 	TF1 *f1 = new TF1("f1","gaus(0)",-2.,2.);
// 	Double_t par1[3]={60,0,0.5};
// 	f1->SetParameters(par1);
// 	f1->SetLineColor(2);
// 
// 	h1->Fit("f1","RB");
// 	double mu1=f1->GetParameter(1);
// 	double sigma1=f1->GetParameter(2);

	TCanvas* c3 = new TCanvas("c3", "", 100, 100, 800, 800); 	
	h2->SetTitle("Difference between cluster and initial photon phi angle");
	h2->GetXaxis()->SetTitle("#phi_{reco} - #phi_{truth}, degree");
	h2->Draw();

}