clear message closing #92

src/Classic/Solvers/CollimatorPhysics.cpp

@@ -714,23 +714,11 @@ void  CollimatorPhysics::Rot(double &px, double &pz, double &x, double &z, doubl
     double Psixz;
     double pxz;
     // Calculate the angle between the px and pz momenta to change from beam coordinate to lab coordinate
-
-    /*
-    if (px>=0 && pz>=0) Psixz = atan(px/pz);
-    else if (px>0 && pz<0)
-	Psixz = atan(px/pz) + Physics::pi;
-    else if (px<0 && pz>0)
-	Psixz = atan(px/pz) + 2*Physics::pi;
-    else
-	Psixz = atan(px/pz) + Physics::pi;
-    */
-
     Psixz = atan2(px,pz);
 
     // Apply the rotation about the random angle thetacou & change from beam
     // coordinate system to the lab coordinate system using Psixz (2 dimensions)
     pxz = sqrt(px*px + pz*pz);
-
     if(coord==1) {
     	x = x + deltas * px/normP + xplane*cos(Psixz);
     	z = z - xplane * sin(Psixz);
@@ -745,6 +733,13 @@ void  CollimatorPhysics::Rot(double &px, double &pz, double &x, double &z, doubl
 
 }
 
+Vector_t ArbitraryRotation(Vector_t &W, Vector_t &Rorg, double Theta) {
+  double C=cos(Theta);
+  double S=sin(Theta);
+  W = W / sqrt(dot(W,W));
+  return Rorg * C + cross(W,Rorg) * S + W * dot(W,Rorg) * (1.0-C);
+} 
+
 /// Coulomb Scattering: Including Multiple Coulomb Scattering and large angle Rutherford Scattering.
 /// Using the distribution given in Classical Electrodynamics, by J. D. Jackson.
 //--------------------------------------------------------------------------
@@ -777,17 +772,6 @@ void  CollimatorPhysics::CoulombScat(Vector_t &R, Vector_t &P, double &deltat) {
     if(collshape_m == "CCollimator")
         R = R * 1000.0;
 
-    double P2 = gsl_rng_uniform(rGen_m);
-    if ((P2 < 0.0047) && enableRutherfordScattering_m) {
-        double P3 = gsl_rng_uniform(rGen_m);
-        double thetaru = 2.5 * sqrt(1 / P3) * sqrt(2.0) * theta0;
-        double P4 = gsl_rng_uniform(rGen_m);
-        if(P4 > 0.5)
-            thetaru = -thetaru;
-	coord = 0; // no change in coordinates but one in momenta-direction
-	Rot(P(0),P(2),R(0),R(2), xplane, normP, thetaru, deltas, coord);
-    }
-
     // y-direction: See Physical Review, "Multiple Scattering"
     z1 = gsl_ran_gaussian(rGen_m,1.0);
     z2 = gsl_ran_gaussian(rGen_m,1.0);
@@ -807,15 +791,23 @@ void  CollimatorPhysics::CoulombScat(Vector_t &R, Vector_t &P, double &deltat) {
     if(collshape_m == "CCollimator")
         R = R * 1000.0;
 
-    P2 = gsl_rng_uniform(rGen_m);
+    // Rutherford-scattering
+    double P2 = gsl_rng_uniform(rGen_m);
     if ((P2 < 0.0047) && enableRutherfordScattering_m) {
         double P3 = gsl_rng_uniform(rGen_m);
-        double thetaru = 2.5 * sqrt(1 / P3) * sqrt(2.0) * theta0;
-        double P4 = gsl_rng_uniform(rGen_m);
-        if(P4 > 0.5)
-            thetaru = -thetaru;
-	coord = 0; // no change in coordinates but one in momenta-direction
-       	Rot(P(1),P(2),R(1),R(2), yplane, normP, thetaru, deltas, coord);
+        //double thetaru = 2.5 * sqrt(1 / P3) * sqrt(2.0) * theta0;
+        double thetaru = 2.5 * sqrt(1 / P3) * 2.0 * theta0;
+        double phiru = 2.0 * M_PI * gsl_rng_uniform(rGen_m);
+	double th0=atan2(sqrt(P(0)*P(0)+P(1)*P(1)),fabs(P(2)));
+	Vector_t W,X;
+	X(0)=cos(phiru)*sin(thetaru);
+	X(1)=sin(phiru)*sin(thetaru);
+	X(2)=cos(thetaru);
+	X*=sqrt(dot(P,P));
+	W(0)=-P(1);
+	W(1)= P(0);
+	W(2)= 0.0;
+	P = ArbitraryRotation(W, X,th0);
     }
 }