Coulomb scattering and rot updated in cuda version

src/CUDA/CudaCollimatorPhysics.cu

 #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;
-
-  /*
-  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);
-
-  if(coord==1) {
-    x = x + deltas * px/normP + xplane*cos(Psixz);
-    z = z - xplane * sin(Psixz);
+  // 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);
+
+  // 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) {
+    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