// nucleon radius 1.35 fm
//	lead radius 7 
const int MAX_A = 200*208;
const int A = 208;
const double Pb_r = 7;
const double Nu_r = 1.3;
const int nNumbers = 1;
const double pi = TMath::Pi();
double participant[nNumbers];
int mm = 0;
int parti = 0;
double l1[2 * 416];
double l2[2 * 416];

bool checkDist(double x1, double y1, double x2, double y2)
{
	if ((x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2) < 4 * Nu_r*Nu_r)
		return true;
	else return false;
}

void overfunction(double b)
{
	printf("Hello world\n");
	using namespace std;
	const double sx = b;
	const double sy = 0;
	const double sz = 0;
	double nuccord1[3][MAX_A];
	double nuccord2[3][MAX_A];
	bool collision0[A];
	bool collision1[A];
	double l1pointx[2 * A][MAX_A];
	double l1pointy[2 * A][MAX_A];
	double l2pointx[2 * A][MAX_A];
	double l2pointy[2 * A][MAX_A];
	//nucleon1
	int cnn1 = 0;
	while (cnn1 < A)
	{
		double x = (gRandom->Rndm() - 0.5) * 2 * Pb_r;
		double y = (gRandom->Rndm() - 0.5) * 2 * Pb_r;
		double z = (gRandom->Rndm() - 0.5) * 2 * Pb_r;
		if (x*x + y*y +z*z< Pb_r*Pb_r)
		{
			nuccord1[0][cnn1] = x;
			nuccord1[1][cnn1] = y;
			nuccord1[2][cnn1] = z;
			collision0[cnn1] = false;
			cnn1++;
		}
	}
	//nucleon2
	int cnn2 = 0;
	while (cnn2 < A)
	{
		double x = (gRandom->Rndm() - 0.5) * 2 * Pb_r;
		double y = (gRandom->Rndm() - 0.5) * 2 * Pb_r;
		double z = (gRandom->Rndm() - 0.5) * 2 * Pb_r;
		if (x*x + y*y +z*z < Pb_r*Pb_r)
		{
			nuccord2[0][cnn2] = x + sx;
			nuccord2[1][cnn2] = y + sy;
			nuccord2[2][cnn2] = z + sz;
			collision1[cnn2] = false;
			cnn2++;
		}
	}

	double theta;
	for (int i = 0; i < A; i++)
		for (int j = 0; j < A; j++)
		{
			if (checkDist(nuccord1[0][i], nuccord1[1][i], nuccord2[0][j], nuccord2[1][j]))
			{
				collision0[i] = true;
				collision1[j] = true;
			}
		}
	double m;
	for (int j = 0; j < A; j++)
	{
		if (collision0[j])
		{
			for (int i=0; i < 2; i++)
			{

				theta = (gRandom->Rndm() - 0.50001)*pi;
				m = tan(theta);
				printf("%lf\n", m);
				double yinte1 = nuccord1[1][j];
				l1pointx[parti][i] = (-m*yinte1 + sqrt(m*m*Pb_r*Pb_r + Pb_r*Pb_r - yinte1*yinte1)) / (m*m + 1);
				l1pointy[parti][i] = m*l1pointx[parti][i] + yinte1;
				l2pointx[parti][i] = ((b - yinte1*m) - sqrt((b - yinte1*m)*(b - yinte1*m) - (m*m + 1)*(yinte1*yinte1 + b*b - Pb_r*Pb_r))) / (m*m + 1);
				l2pointy[parti][i] = m*l1pointx[parti][i] + yinte1;
				printf("%lf\n", m);
			}
			parti++;
		}
	}
	for (int j = 0; j < A; j++)
	{
		if (collision1[j])
		{
			for (int i=0; i < 2; i++)
			{

				theta = (gRandom->Rndm() - 0.50001)*pi;
				double m= tan(theta);
				double yinte2 = nuccord2[1][j];
				l1pointx[parti][i] = ((-m*yinte2 + sqrt(m*m*Pb_r*Pb_r + Pb_r*Pb_r - yinte2*yinte2)) / (m*m + 1));
				l1pointy[parti][i] = m*l1pointx[parti][i] + yinte2;
				l2pointx[parti][i] = ((b - yinte2*m) - sqrt((b - yinte2*m)*(b - yinte2*m) - (m*m + 1)*(yinte2*yinte2 + b*b - Pb_r*Pb_r))) / (m*m + 1);
				l2pointy[parti][i] = m*l1pointx[parti][i] + yinte2;
				printf("%lf\n", m);
			}
			parti++;
		}
	}

	for (int i=0; i < 2; i++)
	{
		for (int j=0; j < parti; j++)
		{
			l1[mm++] = sqrt(l1pointx[j][i] * l1pointx[j][i] + l1pointy[j][i] * l1pointy[j][i]);
		}
	}
	for (int i=0; i < 2; i++)
	{
		for (int j=0; j < parti; j++)
		{
			l2[mm++] = sqrt(l2pointx[j][i] * l2pointx[j][i] + l2pointy[j][i] * l2pointy[j][i]);
		}
	}
}
void l1l2()
{
	auto h0 = new TH2D("Gr", "1", 100, 0, 410, 100, 0, 1.0);
	for (int i = 0; i < nNumbers; i++)
	{
		double rate;
		double b = gRandom->Rndm() * 2 * Pb_r;
		printf("%lf\n", b);
		overfunction(b);
		printf("%lf\n", b);
		printf("%d", parti);
		for (int i = 0; i < 2 * parti; i++)
		{
			h0->Fill(l1[i], l2[i]);
			/*	if (((ecc[i] - 0.53) < 0.01) && ((ecc[i] - 0.53) > -0.01))
				{
					printf("%lf\n", participant[i]);
				}

			*/
			printf("%lf\n", l1[i]);
			printf("%lf\n", l2[i]);

		}
	}
	h0->SetMarkerStyle(2);
	h0->SetMarkerSize(1);
	h0->SetMarkerColor(60);
	h0->SetTitle("Glauber MC Simulation");
	h0->Draw();
}