#include <TROOT.h>
#include <TObject.h>
#include <TTree.h>
#include <TFile.h>
#include <TBranch.h>
#include <TH1.h>
#include <TSQLServer.h>
#include <TSQLResult.h>
#include <TSQLRow.h>
#include <TVector3.h>
#include <TRandom.h>


#include <math.h>
#include <stdio.h>
#include <stdlib.h> 

TRandom *pRand;
double cangle;
TVector3 vg[110];
TSQLServer *pServer;
TSQLResult *pResult;
TSQLRow    *pRow;


void GVector(int idet, TVector3& xgeo);
void Simulate(int n);

//TROOT simple("simple", "Copy data in the root file");

int main(int argc, char **argv)
{
	
	if(gROOT->LoadClass("TMySQLServer","libTmysql.dll") < 0)
	{
		printf("TMySQLServer was not loaded \n");
		return -1;
	}
	
	double teta,fi,rteta,rfi,angle;
	double rr = TMath::Pi() / 180.; // convertion into radian
	int    idet,i;
	char   SQL[255];

	//== Detectors
	angle  = 7.5 * rr;						  //== Max angle on the gamma detector
	cangle = 1. - cos(angle);

	pServer = TSQLServer::Connect("mysql://localhost/Ternary"
		,"daniel","XXXXXX");
	if(!pServer)
	{
		printf("Server was not connected\n");
		return -2;
	}
	
	//== Direction on gamma detectors
	sprintf(SQL,"SELECT * FROM gdposition");
	pResult = pServer->Query(SQL);
	while(pRow = pResult->Next())
	{
		sscanf(pRow->GetField(0),"%i",  &idet);
		sscanf(pRow->GetField(1),"%le", &teta);
		sscanf(pRow->GetField(2),"%le", &fi);

		i = idet - 1;

		rteta  = rr*teta;
		rfi    = rr*fi;

		vg[i].SetX(1.);
		vg[i].SetTheta(rteta);
		vg[i].SetPhi(rfi);
		delete pRow;
	}
	delete pResult;
	pServer->Close();

	pRand = new TRandom();
	Simulate(100000);
	return 0;
} //-------------------------------------------------------------------

void Simulate(int n)
{
	//=================================================================
	//== Simulate angular distributions between all detector pairs
	//=================================================================

	int i,i1,i2;
	TVector3 v1,v2;
	TFile Fsave("andist.root","RECREATE");
	if(!Fsave.IsOpen())
	{
		printf("FILE WAS NOT OPEN\n");
		return;
	}

	TH1D d1("d1","",200,-1,1);
	char name[10];
	char title[80];

	//for(i1 = 0; i1 < 110; ++i1)
	{
		//for(i2 = i1+1; i2 < 110; ++i2)
		{
			//d1.Reset();
			//sprintf(name,"d%3.3d_%3.3d",i1+1,i2+1);
			//sprintf(title,"ANGULAR DISTRIBUTION DET%d AND DET%d",i1+1,i2+1);
			//d1.SetName(name);
			//d1.SetTitle(title);

			//printf("DET %3.3d vs DET %3.3d\n",i1+1,i2+1);
			for(i = 0; i < n; ++i)
			{
				GVector(1,v1);
				GVector(2,v2);
				d1.Fill(v1.Unit()*v2.Unit());
			}
			d1.Write();
		}
	}
}



void GVector(int idet, TVector3& xgeo)
{
	//===============================================================
	//== Random vector on gamma detector
	//===============================================================
			
	TVector3 geo;
	double rteta;
	double rfi,cfi;
	double ang;
	int ig;

	rteta = 2.*TMath::Pi()*pRand->Uniform();
	rfi = acos(1. - pRand->Uniform(cangle));

	geo = vg[idet].Orthogonal();
	geo.Rotate(rteta, vg[idet]);

	xgeo = vg[idet];
	xgeo.Rotate(rfi,geo);
}