mcxtrace: fix FluoPowder on arm (1e-307 -> 0)

Author: Emmanuel FARHI

mcxtrace-comps/examples/Tests_samples/Test_SX/Test_SX.instr

@@ -14,13 +14,18 @@
 *
 * Simply a model source illuminating a SX sample.
 * The sample itself is a Mo bulk crystal.
+* The idea is to compare the fluorescence and diffraction patterns:
+* - index=1: use Single_crystal
+* - index=2: use FluoCrystal
+* - index=3: use Fluorescence+Single_crystal in a GROUP
 *
 * %Example: index=1 Detector: psd4pi_I=1.18843e-13
-* %Example: index=2 Detector: psd4pi_I=4.80386e-13
+* %Example: index=2 Detector: psd4pi_I=6.89931e-12
+* %Example: index=3 Detector: psd4pi_I=4.4547e-13
 *
 * %Parameters
 * reflections: [str] List of powder reflections, LAU/CIF format.
-* index:       [1]   Index of the sample component to use. 1=Single_crystal, 2=FluoCrystal
+* index:       [1]   Index of the sample component to use. 1=Single_crystal, 2=FluoCrystal, 3=Fluorescence+Single_crystal in a GROUP
 *
 * %End
 *******************************************************************************/
@@ -58,6 +63,24 @@ EXTEND %{
   else Stype=type;
 %}
 
+COMPONENT FluoG = Fluorescence(
+  radius = .5e-4, yheight = 1e-3, material=reflections)
+WHEN (index == 3)
+AT (0, 0, 0) RELATIVE sample_pos
+GROUP FluSX
+EXTEND %{
+  if (SCATTERED) Stype = type;
+%}
+
+COMPONENT SX = COPY(sample)(sigma_inc=-1)
+WHEN (index == 3)
+AT (0, 0, 0) RELATIVE sample_pos
+GROUP FluSX
+EXTEND %{
+  if (SCATTERED) Stype=DIFFRACTION;
+  else ABSORB;
+%}
+
 COMPONENT psd4pi = PSD_monitor_4PI(
     nx = 180, ny = 180, filename = "psd4pi", radius = 0.1, restore_xray = 1)
 AT (0, 0, 0) RELATIVE sample

mcxtrace-comps/samples/FluoPowder.comp

@@ -947,7 +947,7 @@ do { /* while (intersect) Loop over multiple scattering events */
   if (intersect) { /* scattering event */
     int    i_Z=-1, Z, line=-1;
     double solid_angle;
-    double theta, dsigma, alpha;
+    double theta, dsigma, alpha=0;
     double Ef, dE;
     double kf=k, kf_x, kf_y, kf_z;  // next particle direction
 
@@ -1027,12 +1027,12 @@ do { /* while (intersect) Loop over multiple scattering events */
       { /* relate height of detector to the height on DS cone */
         arg = sin(powder_d_phi*DEG2RAD/2)/sin(2*theta);
         /* If full Debye-Scherrer cone is within powder_d_phi, don't focus */
-        if (arg < -1 || arg > 1) alpha = 0;
+        if (arg < -1 || arg > 1) powder_d_phi = 0;
         /* Otherwise, determine alpha to rotate from scattering plane
          *    into powder_d_phi focusing area */
         else alpha = 2*asin(arg);
       }
-      if (alpha) {
+      if (powder_d_phi) {
         /* Focusing */
         alpha  = fabs(alpha);
         alpha0 = 0.5*randpm1()*alpha;
@@ -1097,13 +1097,13 @@ do { /* while (intersect) Loop over multiple scattering events */
         if (Ex!=0 || Ey!=0 || Ez!=0){
           double EE=sqrt(Ex*Ex+Ey*Ey+Ez*Ez);
           double s=scalar_prod(kf_x,kf_y,kf_z,Ex,Ey,Ez)/EE;
-          p*=(1-s)*(1-s);
-        }else{
+          p *= (1-s)*(1-s);
+        } else {
           /*unpolarized light in - means an effective reduction according to only theta*/
-          p*=(1+cos(theta)*cos(theta))*0.5;
+          p *= (1+cos(theta)*cos(theta))*0.5;
         }
         if (alpha) {
-          p    *= alpha/PI;
+          p *= alpha/PI;
         }
         n_diff++;
         p_diff += p;

mcxtrace-comps/samples/PowderN.comp

@@ -1079,7 +1079,7 @@ TRACE
           }
           if (d_phi_thread) {
             /* Focusing */
-            alpha = fabs(alpha);
+            alpha  = fabs(alpha);
             alpha0 = 0.5*randpm1()*alpha;
             if(focus_flip){
                 alpha0+=M_PI_2;