From 24604042e7ec69e2fabd865cc318d092345d1303 Mon Sep 17 00:00:00 2001 From: Uldis Locans Date: Tue, 19 Sep 2017 13:36:42 +0200 Subject: [PATCH] Coulomb scattering and rot updated in cuda version --- src/CUDA/CudaCollimatorPhysics.cu | 63 +++++++++++++------------------ 1 file changed, 27 insertions(+), 36 deletions(-) diff --git a/src/CUDA/CudaCollimatorPhysics.cu b/src/CUDA/CudaCollimatorPhysics.cu index 2c600d5..8522901 100644 --- a/src/CUDA/CudaCollimatorPhysics.cu +++ b/src/CUDA/CudaCollimatorPhysics.cu @@ -1,13 +1,11 @@ #include "CudaCollimatorPhysics.cuh" //constants used in OPAL -//#define M_P 0.93827231e+00 #define M_P 0.93827204e+00 #define C 299792458.0 #define PI 3.14159265358979323846 #define AVO 6.022e23 #define R_E 2.81794092e-15 -//#define eM_E 0.51099906e-03 #define eM_E 0.51099892e-03 #define Z_P 1 #define K 4.0*PI*AVO*R_E*R_E*eM_E*1e7 @@ -137,38 +135,29 @@ __device__ inline void energyLoss(double &Eng, bool &pdead, curandState &state, * For details: see J. Beringer et al. (Particle Data Group), Phys. Rev. D 86, 010001 (2012), * "Passage of particles through matter" */ -__device__ inline void Rot(double &px, double &pz, double &x, double &z, double &xplane, - double &normP, double &thetacou, double &deltas, int coord, +__device__ inline void Rot(double &px, double &pz, double &x, double &z, double &plane, + double &betaGamma, double &thetacou, double &deltas, int coord, double *par) { - double Psixz; - double pxz; + // Calculate the angle between the px and pz momenta to change from beam coordinate to lab coordinate + const double Psi = atan2(px, pz); + const double pxz = sqrt(px*px + pz*pz); + const double cosPsi = cos(Psi); + const double sinPsi = sin(Psi); + const double cosTheta = cos(thetacou); + const double sinTheta = sin(thetacou); - /* - if (px>=0 && pz>=0) - Psixz = atan(px/pz); - else if (px>0 && pz<0) - Psixz = atan(px/pz) + PI; - else if (px<0 && pz>0) - Psixz = atan(px/pz) + 2*PI; - else - Psixz = atan(px/pz) + PI; - */ - Psixz = atan2(px, pz); - pxz = sqrt(px*px + pz*pz); + // Apply the rotation about the random angle thetacou & change from beam + // coordinate system to the lab coordinate system using Psixz (2 dimensions) + x += deltas * px / betaGamma + plane * cosPsi; + z -= plane * sinPsi; - if(coord==1) { - x = x + deltas * px/normP + xplane*cos(Psixz); - z = z - xplane * sin(Psixz); + if (coord == 1) { + z += deltas * pz / betaGamma; } - if(coord==2) { - x = x + deltas * px/normP + xplane*cos(Psixz); - z = z - xplane * sin(Psixz) + deltas * pz / normP; - } - - px = pxz*cos(Psixz)*sin(thetacou) + pxz*sin(Psixz)*cos(thetacou); - pz = -pxz*sin(Psixz)*sin(thetacou) + pxz*cos(Psixz)*cos(thetacou); + px = pxz * (cosPsi * sinTheta + sinPsi * cosTheta); + pz = pxz * (-sinPsi * sinTheta + cosPsi * cosTheta); } /** @@ -182,11 +171,12 @@ __device__ inline void coulombScat(double3 &R, double3 &P, curandState &state, d double Eng = sqrt(dot(P, P) + 1.0) * M_P - M_P; double gamma = (Eng + M_P) / M_P; - double normP = sqrt(dot(P, P)); double beta = sqrt(1.0 - 1.0 / (gamma * gamma)); + double betaGamma = sqrt(dot(P, P)); double deltas = par[DT_M] * beta * C; + double mass = M_P * 1e9; // in eV - double theta0 = 13.6e6 / (beta * normP * M_P * 1e9) * + double theta0 = 13.6e6 / (beta * betaGamma * mass) * Z_P * sqrt(deltas / par[X0_M]) * (1.0 + 0.038 * log(deltas / par[X0_M])); // x-direction: See Physical Review, "Multiple Scattering" @@ -201,8 +191,9 @@ __device__ inline void coulombScat(double3 &R, double3 &P, curandState &state, d } //__syncthreads(); - double xplane = z1 * deltas * theta0 / sqrt(12.0) + z2 * deltas * theta0 / 2.0; - Rot(P.x, P.z, R.x, R.z, xplane, normP, thetacou, deltas, 1, par); + //double xplane = z1 * deltas * theta0 / sqrt(12.0) + z2 * deltas * theta0 / 2.0; + double xplane = 0.5 * deltas * theta0 * (z1 / sqrt(3.0) + z2); + Rot(P.x, P.z, R.x, R.z, xplane, betaGamma, thetacou, deltas, 0, par); // y-direction: See Physical Review, "Multiple Scattering" z1 = curand_normal_double(&state);//gsl_ran_gaussian(rGen_m,1.0); @@ -215,10 +206,9 @@ __device__ inline void coulombScat(double3 &R, double3 &P, curandState &state, d thetacou = z2 * theta0; } - //__syncthreads(); - - double yplane = z1 * deltas * theta0 / sqrt(12.0) + z2 * deltas * theta0 / 2.0; - Rot(P.y,P.z,R.y,R.z, yplane, normP, thetacou, deltas, 2, par); + //double yplane = z1 * deltas * theta0 / sqrt(12.0) + z2 * deltas * theta0 / 2.0; + double yplane = 0.5 * deltas * theta0 * (z1 / sqrt(3.0) + z2); + Rot(P.y,P.z,R.y,R.z, yplane, betaGamma, thetacou, deltas, 1, par); double P2 = curand_uniform_double(&state);//gsl_rng_uniform(rGen_m); if( (P2 < 0.0047) && enableRutherfordScattering) { @@ -305,6 +295,7 @@ __global__ void kernelCollimatorPhysics(T *data, double *par, curandState *state P.x = P.x * ptot / sq; P.y = P.y * ptot / sq; P.z = P.z * ptot / sq; + coulombScat(R[tid], P, s, p, enableRutherfordScattering); //update particle momentum