update DKS version

CMakeLists.txt

@@ -2,7 +2,7 @@ CMAKE_MINIMUM_REQUIRED (VERSION 3.2)
 PROJECT (DKS)
 SET (DKS_VERSION_MAJOR 1)
 SET (DKS_VERSION_MINOR 1)
-SET (DKS_VERSION_PATCH 0)
+SET (DKS_VERSION_PATCH 1)
 set (DKS_VERSION ${DKS_VERSION_MAJOR}.${DKS_VERSION_MINOR}.${DKS_VERSION_PATCH})
 SET (PACKAGE \"dks\")
 SET (PACKAGE_BUGREPORT \"locans.uldis@psi.ch\")

src/Algorithms/CollimatorPhysics.h

@@ -16,7 +16,8 @@ public:
   
   virtual ~DKSCollimatorPhysics() { }
   
-  virtual int CollimatorPhysics(void *mem_ptr, void *par_ptr, int numpartices) = 0;
+  virtual int CollimatorPhysics(void *mem_ptr, void *par_ptr, int numpartices, 
+				bool enableRutherforScattering = true) = 0;
 
   virtual int CollimatorPhysicsSoA(void *label_ptr, void *localID_ptr, 
 				   void *rx_ptr, void *ry_ptr, void *rz_ptr, 

src/CUDA/CudaCollimatorPhysics.cu

@@ -23,9 +23,10 @@
 #define X0_M 9
 #define I_M 10
 #define DT_M 11
+#define LOWENERGY_THR 12
 
 #define BLOCK_SIZE 128
-#define NUMPAR 12
+#define NUMPAR 13
 
 __device__ inline double dot(double3 &d1, double3 &d2) {
 
@@ -41,6 +42,23 @@ __device__ inline double3 cross(double3 &lhs, double3 &rhs) {
   return tmp;
 }
 
+__device__ inline double3 ArbitraryRotation(double3 &W, double3 &Rorg, double Theta) {
+  double c=cos(Theta);
+  double s=sin(Theta);
+  double dotW = sqrt(dot(W,W));
+  W.x = W.x / dotW;
+  W.y = W.y / dotW;
+  W.z = W.z / dotW;
+
+  double dotWR = dot(W, Rorg) * (1.0 - c);
+  double3 crossW = cross(W, Rorg);
+  double3 tmp;
+  tmp.x = Rorg.x * c + crossW.x * s + W.x * dotWR;
+  tmp.y = Rorg.y * c + crossW.y * s + W.y * dotWR;
+  tmp.z = Rorg.z * c + crossW.z * s + W.z * dotWR;
+  return tmp;
+} 
+
 __device__ inline bool checkHit(double &z, double *par) {
 
   /* check if particle is in the degrader material */
@@ -89,7 +107,7 @@ __device__ inline void energyLoss(double &Eng, bool &pdead, curandState &state,
     Eng = Eng + delta_E / 1E3;
   }
 
-  pdead = ((Eng<1E-4) || (dEdx>0));
+  pdead = ( (Eng < par[LOWENERGY_THR]) || (dEdx > 0) );
 
 }
 
@@ -100,6 +118,7 @@ __device__ inline void Rot(double &px, double &pz, double &x, double &z, double
   double Psixz;
   double pxz;
 
+  /*
   if (px>=0 && pz>=0)
     Psixz = atan(px/pz);
   else if (px>0 && pz<0)
@@ -108,7 +127,8 @@ __device__ inline void Rot(double &px, double &pz, double &x, double &z, double
     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) {
@@ -125,7 +145,9 @@ __device__ inline void Rot(double &px, double &pz, double &x, double &z, double
   pz = -pxz*sin(Psixz)*sin(thetacou) + pxz*cos(Psixz)*cos(thetacou);
 }
 
-__device__ inline void coulombScat(double3 &R, double3 &P, curandState &state, double* par) {
+__device__ inline void coulombScat(double3 &R, double3 &P, curandState &state, double* par,
+				   bool enableRutherfordScattering) 
+{
 
   double Eng = sqrt(dot(P, P) + 1.0) * M_P - M_P;
   double gamma = (Eng + M_P) / M_P;
@@ -148,20 +170,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 P2 = curand_uniform_double(&state);//gsl_rng_uniform(rGen_m);
-  if(P2 < 0.0047) {
-    double P3 = curand_uniform_double(&state);//gsl_rng_uniform(rGen_m);
-    double thetaru = 2.5 * sqrt(1 / P3) * sqrt(2.0) * theta0;
-    double P4 = curand_uniform_double(&state);//gsl_rng_uniform(rGen_m);
-    if(P4 > 0.5)
-      thetaru = -thetaru;
-    Rot(P.x,P.z,R.x,R.z, xplane, normP, thetaru, deltas, 0, par);
-  }
-
   // y-direction: See Physical Review, "Multiple Scattering"
   z1 = curand_normal_double(&state);//gsl_ran_gaussian(rGen_m,1.0);
   z2 = curand_normal_double(&state);//gsl_ran_gaussian(rGen_m,1.0);
@@ -178,14 +189,23 @@ __device__ inline void coulombScat(double3 &R, double3 &P, curandState &state, d
   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);
 
-  P2 = curand_uniform_double(&state);//gsl_rng_uniform(rGen_m);
+  double P2 = curand_uniform_double(&state);//gsl_rng_uniform(rGen_m);
-  if(P2 < 0.0047) {
+  if( (P2 < 0.0047) && enableRutherfordScattering) {
     double P3 = curand_uniform_double(&state);//gsl_rng_uniform(rGen_m);
-    double thetaru = 2.5 * sqrt(1 / P3) * sqrt(2.0) * theta0;
+    //double thetaru = 2.5 * sqrt(1 / P3) * sqrt(2.0) * theta0;
-    double P4 = curand_uniform_double(&state);//gsl_rng_uniform(rGen_m);
+    double thetaru = 2.5 * sqrt(1 / P3) * 2.0 * theta0;
-    if(P4 > 0.5)
+    double phiru = 2.0 * M_PI * curand_uniform_double(&state);
-      thetaru = -thetaru;
+    double th0=atan2(sqrt(P.x*P.x+P.y*P.y),fabs(P.z));
-    Rot(P.y,P.z,R.y,R.z, yplane, normP, thetaru, deltas, 0, par);
+    double3 W,X;
+    
+    double dotP = sqrt(dot(P,P));
+    X.x = cos(phiru)*sin(thetaru) * dotP;
+    X.y = sin(phiru)*sin(thetaru) * dotP;
+    X.z = cos(thetaru) * dotP;
+    W.x = -P.y;
+    W.y = P.x;
+    W.z = 0.0;
+    P = ArbitraryRotation(W, X, th0);
   }
 
 }
@@ -193,7 +213,7 @@ __device__ inline void coulombScat(double3 &R, double3 &P, curandState &state, d
 
 template <typename T>
 __global__ void kernelCollimatorPhysics(T *data, double *par, curandState *state,
-					int numparticles)
+					int numparticles, bool enableRutherfordScattering)
 {
 
   //get global id and thread id
@@ -235,7 +255,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); 
+	coulombScat(R[tid], P, s, p, enableRutherfordScattering); 
 
 	data[idx].Pincol = P;
       } else {
@@ -258,7 +278,8 @@ __global__ void kernelCollimatorPhysics(T *data, double *par, curandState *state
 }
 
 __global__ void kernelCollimatorPhysics2(CUDA_PART2_SMALL data, double *par, 
-					 curandState *state, int numparticles)
+					 curandState *state, int numparticles,
+					 bool enableRutherfordScattering)
 {
 
   //get global id and thread id
@@ -296,7 +317,7 @@ __global__ void kernelCollimatorPhysics2(CUDA_PART2_SMALL data, double *par,
     	P.x = P.x * ptot / sq;
     	P.y = P.y * ptot / sq;
     	P.z = P.z * ptot / sq;
-    	coulombScat(R[tid], P, s, p); 
+    	coulombScat(R[tid], P, s, p, enableRutherfordScattering); 
       
     	data.Pincol[idx] = P;
     } else {
@@ -663,7 +684,8 @@ struct less_then
   }
 };
 
-int CudaCollimatorPhysics::CollimatorPhysics(void *mem_ptr, void *par_ptr, int numparticles)
+int CudaCollimatorPhysics::CollimatorPhysics(void *mem_ptr, void *par_ptr, int numparticles,
+					     bool enableRutherfordScattering)
 {
 
   int threads = BLOCK_SIZE;
@@ -676,7 +698,8 @@ int CudaCollimatorPhysics::CollimatorPhysics(void *mem_ptr, void *par_ptr, int n
   kernelCollimatorPhysics<<<blocks, threads, smem_size>>>((CUDA_PART_SMALL*)mem_ptr, 
 							  (double*)par_ptr,
 							  m_base->cuda_getCurandStates(),
-							  numparticles);
+							  numparticles,
+							  enableRutherfordScattering);
 
   cudaError_t err = cudaGetLastError();
   if (err != cudaSuccess)

src/CUDA/CudaCollimatorPhysics.cuh

@@ -110,7 +110,7 @@ public:
    *
    */
   int CollimatorPhysics(void *mem_ptr, void *par_ptr, 
-			int numpartices);
+			int numpartices, bool enableRutherforScattering = true);
 
   int CollimatorPhysicsSoA(void *label_ptr, void *localID_ptr, 
 			   void *rx_ptr, void *ry_ptr, void *rz_ptr, 

src/DKSOPAL.cpp

@@ -209,18 +209,20 @@ int DKSOPAL::callMultiplyComplexFields(void *mem_ptr1, void *mem_ptr2, int size,
 
 int DKSOPAL::callCollimatorPhysics(void *mem_ptr, void *par_ptr, 
 				   int numparticles, int numparams,
-				   int &numaddback, int &numdead) 
+				   int &numaddback, int &numdead, 
+				   bool enableRutherforScattering) 
 {
 
-  return dkscol->CollimatorPhysics(mem_ptr, par_ptr, numparticles);
+  return dkscol->CollimatorPhysics(mem_ptr, par_ptr, numparticles, enableRutherforScattering);
 
 }
 
 
-int DKSOPAL::callCollimatorPhysics2(void *mem_ptr, void *par_ptr, int numparticles) 
+int DKSOPAL::callCollimatorPhysics2(void *mem_ptr, void *par_ptr, int numparticles,
+				    bool enableRutherforScattering) 
 {
 
-  return dkscol->CollimatorPhysics(mem_ptr, par_ptr, numparticles);
+  return dkscol->CollimatorPhysics(mem_ptr, par_ptr, numparticles, enableRutherforScattering);
   
 }
 

src/DKSOPAL.h

@@ -159,14 +159,16 @@ public:
    */
   int callCollimatorPhysics(void *mem_ptr, void *par_ptr, 
 			    int numparticles, int numparams, 
-			    int &numaddback, int &numdead);
+			    int &numaddback, int &numdead,
+			    bool enableRutherfordScattering = true);
   
   /** 
    * Monte carlo code for the degrader from OPAL classic/5.0/src/Solvers/CollimatorPhysics.cpp on device.
    * For specifics check OPAL docs and CudaCollimatorPhysics class documentation.
    * TODO: opencl and mic implementations.
    */
-  int callCollimatorPhysics2(void *mem_ptr, void *par_ptr, int numparticles);
+  int callCollimatorPhysics2(void *mem_ptr, void *par_ptr, int numparticles, 
+			     bool enableRutherfordScattering = true);
 
   /** 
    * Monte carlo code for the degrader from OPAL classic/5.0/src/Solvers/CollimatorPhysics.cpp on device.

src/MIC/MICCollimatorPhysics.cpp

@@ -368,7 +368,9 @@ void energyLoss(double &Eng, double &dEdx, double *par, double *randv, int ri) {
 
 }
 
-int MICCollimatorPhysics::CollimatorPhysics(void *mem_ptr, void *par_ptr, int numparticles) {
+int MICCollimatorPhysics::CollimatorPhysics(void *mem_ptr, void *par_ptr, int numparticles,
+					    bool enableRutherforScattering) 
+{
 
   //cast device memory pointers to appropriate types
   MIC_PART_SMALL *data = (MIC_PART_SMALL*) mem_ptr;

src/MIC/MICCollimatorPhysics.h

@@ -40,7 +40,8 @@ public:
 
   ~MICCollimatorPhysics() {  };
 
-  int CollimatorPhysics(void *mem_ptr, void *par_ptr, int numparticles);
+  int CollimatorPhysics(void *mem_ptr, void *par_ptr, int numparticles, 
+			bool enableRutherforScattering = true);
 
   int CollimatorPhysicsSoA(void *label_ptr, void *localID_ptr, 
 			   void *rx_ptr, void *ry_ptr, void *rz_ptr, 

src/OpenCL/OpenCLCollimatorPhysics.cpp

@@ -34,7 +34,7 @@ TODO:
 2. boost.compute sort for user defined structure crashes
 */
 int OpenCLCollimatorPhysics::CollimatorPhysics(void *mem_ptr, void *par_ptr, 
-					       int numparticles) 
+					       int numparticles, bool enableRutherforScattering) 
 {
   /*
   //set number of total threads, and number threads per block

src/OpenCL/OpenCLCollimatorPhysics.h

@@ -52,7 +52,8 @@ public:
   }
 
   /* execute degrader code on device */
-  int CollimatorPhysics(void *mem_ptr, void *par_ptr, int numparticles);
+  int CollimatorPhysics(void *mem_ptr, void *par_ptr, int numparticles, 
+			bool enableRutherforScattering = true);
 
   int CollimatorPhysicsSoA(void *label_ptr, void *localID_ptr, 
 			   void *rx_ptr, void *ry_ptr, void *rz_ptr,