// 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(); }