/*****************************************************************************
*
* McXtrace, X-ray tracing package
*         Copyright, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*         Synchrotron SOLEIL, Saint-Aubin, France
*
* Component: Fluorescence
*
* %Identification
* Written by: E. Farhi
* Date:       March 2022
* Origin:     Synchrotron SOLEIL
*
* Sample model handling absorption, fluorescence, Compton and Rayleigh scattering.
*
* %Description
* Sample that models many photon-matter interactions:
* - absorption (photon excites an electron and creates a hole)
* - fluorescence (excited electrons emit light while falling into lower states)
* - Compton scattering  (inelastic, incoherent)
* - Rayleigh scattering (elastic,   coherent)
*
* An important option to enhance statistics is to set 'p_interact' to, say,
* 30 percent (0.3) in order to force a fraction of the beam to scatter. This
* will result on a larger number of scattered events, retaining intensity.
*
* The 'material' specification is given as a chemical formulae, e.g. "LaB6". It
* may also be given as a file name (CIF/LAU/LAZ/FullProf format) in which case 
* the formulae is guessed (but may be approximative).
*
* By setting the 'order' to 1, the absorption along the scattered path is handled.
* A higher 'order' will handle multiple scattering events, and final absorption.
* For instance, a value order>=2 handles e.g. fluorescence iterative cascades 
* in the material. Leaving 'order=0' handles the single scattering only.
*
* <b>Sample shape:</b>
* Sample shape may be a cylinder, a sphere, a box or any other shape
*   box/plate:       xwidth x yheight x zdepth (thickness=0)
*   hollow box/plate:xwidth x yheight x zdepth and thickness>0
*   cylinder:        radius x yheight (thickness=0)
*   hollow cylinder: radius x yheight and thickness>0
*   sphere:          radius (yheight=0 thickness=0)
*   hollow sphere:   radius and thickness>0 (yheight=0)
*   any shape:       geometry=OFF file
*
*   The complex geometry option handles any closed non-convex polyhedra.
*   It computes the intersection points of the photon ray with the object
*   transparently, so that it can be used like a regular sample object.
*   It supports the OFF, PLY and NOFF file format but not COFF (colored faces).
*   Such files may be generated from XYZ data using:
*     qhull < coordinates.xyz Qx Qv Tv o > geomview.off
*   or
*     powercrust coordinates.xyz
*   and viewed with geomview or java -jar jroff.jar (see below).
*   The default size of the object depends of the OFF file data, but its
*   bounding box may be resized using xwidth,yheight and zdepth.
*
* <b>Concentric components:</b>
* This component has the ability to contain other components when used in
* hollow cylinder geometry (namely sample environment, e.g. cryostat and
* furnace structure). Such component 'shells' should be split into input and
* output side surrounding the 'inside' components. First part must then use
* 'concentric=1' flag to enter the inside part. The component itself must be
* repeated to mark the end of the concentric zone. The number of concentric
* shells and number of components inside is not limited.
*
* COMPONENT F_in = Fluorescence(material="Al", concentric=1, ...)
* AT (0,0,0) RELATIVE sample_position
*
* COMPONENT something_inside ... // e.g. the sample itself or other materials
*
* COMPONENT F_out = COPY(F_in)(concentric=0)
* AT (0,0,0) RELATIVE sample_position
*
* The computation is made via the XRayLib (apt install libxrl-dev).
*
* Example: Fluorescence(material="LaB6",
*  xwidth=0.001,yheight=0.001,zdepth=0.0001, p_interact=0.99,
*  target_index=1, focus_xw=0.0005, focus_yh=0.0005)
*
* This sample component can advantageously benefit from the SPLIT feature, e.g.
* SPLIT COMPONENT pow = Fluorescence(...)
*
* %Parameters
* material:  [str]    Chemical formulae, e.g. "LaB6", "Pb2SnO4". If may also be a CIF/LAZ/LAU file.
* weight:    [g/mol]  Atomic/molecular weight of material.
* density:   [g/cm^3] Density of material. V_rho=density/weight/1e24*N_A.
* rho:       [AA-3]   Density of scattering elements (nb atoms/unit cell V_0).
* packing_factor: [1] How dense is the material compared to bulk 0-1.
* radius:    [m]      Outer radius of sample in (x,z) plane. cylinder/sphere.
* xwidth:    [m]      width for a box sample shape.
* yheight:   [m]      Height of sample in vertical direction for box/cylinder shapes.
* zdepth:    [m]      depth for a box sample shape.
* thickness: [m]      Thickness of hollow sample Negative value extends the hollow volume outside of the box/cylinder.
* concentric: [1]     Indicate that this component has a hollow geometry and may contain other components. It should then be duplicated after the inside part (only for box, cylinder, sphere) [1]
* geometry:  [str]    Name of an Object File Format (OFF) or PLY file for complex geometry. The OFF/PLY file may be generated from XYZ coordinates using qhull/powercrust.
* p_interact: [1]     Force a given fraction of the beam to scatter, keeping intensity right, to enhance small signals (-1 unactivate).
* focus_xw:  [m]      Horiz. dimension of a rectangular area.
* focus_yh:  [m]      Vert.  dimension of a rectangular area.
* focus_aw:  [deg]    Horiz. angular dimension of a rectangular area.
* focus_ah:  [deg]    Vert.  angular dimension of a rectangular area.
* focus_r:   [m]      Radius of disk containing target. Use 0 for full space.
* target_index: [1]   Relative index of component to focus at, e.g. next is +1.
* target_x:  [m]      Position of target to focus at, along X.
* target_y:  [m]      Position of target to focus at, along Y.
* target_z:  [m]      Position of target to focus at, along Z.
* flag_compton: [1]   When 0, the Compton  scattering is ignored.
* flag_rayleigh:[1]   When 0, the Rayleigh scattering is ignored.
* flag_lorentzian:[1] When 1, the line shapes are assumed to be Lorentzian, else Gaussian
* order:     [1]      Limit multiple fluorescence up to given order. Last iteration is absorption only.
*
* CALCULATED PARAMETERS:
* type: scattering event type 0=fluorescence, 1=Rayleigh, 2=Compton, 3=transmit
*
* %Link
* The XRayLib https://github.com/tschoonj/xraylib
* %Link
* Fluorescence https://en.wikipedia.org/wiki/Fluorescence
* %Link
* Rayleigh https://en.wikipedia.org/wiki/Rayleigh_scattering
* %Link
* Compton https://en.wikipedia.org/wiki/Compton_scattering
* %Link
* X-ray absorption edges http://skuld.bmsc.washington.edu/scatter/AS_periodic.html
* %Link
* X-ray fluorescence spectra http://www.xrfresearch.com/xrf-spectra/
* %Link
* X-ray edges and fluo lines https://physics.nist.gov/PhysRefData/XrayTrans/Html/search.html
*
* %E
***********************************************************/

DEFINE COMPONENT Fluorescence
SETTING PARAMETERS(
  string geometry=0,
  radius=0, thickness=0,
  xwidth=0, yheight=0, zdepth=0,
  int concentric=0,
  string material="LaB6", packing_factor=0, rho=0, density=0, weight=0, 
  p_interact=0, 
  target_x = 0, target_y = 0, target_z = 0, focus_r = 0,
  focus_xw=0, focus_yh=0, focus_aw=0, focus_ah=0, int target_index=0,
  int flag_compton=1, int flag_rayleigh=1, int flag_lorentzian=1, int order=1)

/* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */
DEPENDENCY " @XRLFLAGS@ -DUSE_OFF "
NOACC

/* ========================================================================== */

SHARE %{
  #ifndef XRAYLIB_LINES_MAX
  #define XRAYLIB_LINES_MAX 383
  #define FLUORESCENCE 0  // Fluo
  #define RAYLEIGH     1  // Coherent
  #define COMPTON      2  // Incoherent
  #define TRANSMISSION 3
  #include <xraylib/xraylib.h>
  #endif
  %include "read_table-lib"
  %include "interoff-lib" // for OFF/PLY geometry
  
/* See XrayLib code (c) T. Schoonjans
 * /usr/include/xraylib/xraylib-parser.h      for compoundData
 * /usr/include/xraylib/xraylib.h             for XS
 */

/* inspired from:
 * https://github.com/golosio/xrmc src/photon/photon.cpp (c) Bruno Golosio    
 */
  
/* XRMC_CrossSections: Compute interaction cross sections in [barn/atom]
 * Return total cross section, given Z and E0:
 *   total_xs = XRMC_CrossSections(Z, E0, xs[3]);
 */
double XRMC_CrossSections(int Z, double E0, double *xs) {
  int    i_line;
  
  if (xs == NULL) return 0;
  
  // loop on possible fluorescence lines
  for (i_line=0; i_line<XRAYLIB_LINES_MAX; i_line++) { 
    // cumulative sum of the line cross sections
    xs[FLUORESCENCE] += CSb_FluorLine(Z, -i_line, E0, NULL); /* XRayLib */
  }

  // coherent and incoherent cross sections
  xs[RAYLEIGH] = CSb_Rayl( Z, E0, NULL);
  xs[COMPTON]  = CSb_Compt(Z, E0, NULL);
  
  // total interaction cross section, should converge to CSb_Total(Z, E0, NULL)
  return xs[FLUORESCENCE] + xs[RAYLEIGH] + xs[COMPTON];
} // XRMC_CrossSections

/* XRMC_SelectFromDistribution: select a random element from a distribution
 *   index = XRMC_SelectFromDistribution(cum_sum[N], N);
 * index is returned within 0 and N-1
 * The x_arr must be a continuously increasing cumulated sum, which last element is the max
 */
int XRMC_SelectFromDistribution(double x_arr[], int N)
{
  double x=rand01()*x_arr[N-1];
  if (x<x_arr[0]) {    // x is smaller than lower limit
    return 0;
  }
  if (x>=x_arr[N-1]) { // x is greater or equal to upper limit
    return N-1;
  }
  int id=0, iu=N-1; // lower and upper index of the subarray to search
  while (iu-id>1) { // search until the size of the subarray to search is >1
    int im = (id + iu)/2; // use the midpoint for equal partition
    // decide which subarray to search
    if (x>=x_arr[im]) id=im; // change min index to search upper subarray
    else iu=im; // change max index to search lower subarray
  }

  return id;
} // XRMC_SelectFromDistribution

/* XRMC_SelectInteraction: select interaction type Fluo/Compton/Rayleigh
 * Return the interaction type from a random choice within cross sections 'xs'
 *   type = XRMC_SelectInteraction(xs[3]);
 * 'xs' is computed with XRMC_CrossSections.
 * type is one of FLUORESCENCE | RAYLEIGH | COMPTON
 */
int XRMC_SelectInteraction(double *xs)
{
  double sum_xs, cum_xs[4];
  int    i;
  
  cum_xs[0]=sum_xs=0;
  for (i=0; i< 3; i++) {
    sum_xs += xs[i];
    cum_xs[i+1]= sum_xs;
  }
  return XRMC_SelectFromDistribution(cum_xs, 4);
} // XRMC_SelectInteraction

/* XRMC_SelectFluorescenceEnergy: select outgoing fluo photon energy, when incoming with 'E0'
 *   Ef = XRMC_SelectFluorescenceEnergy(Z, E0, &dE);
 */
double XRMC_SelectFluorescenceEnergy(int Z, double E0, double *dE)
{
  int i_line;
  double sum_xs, cum_xs_lines[XRAYLIB_LINES_MAX+1];

  // compute cumulated XS for all fluo lines
  cum_xs_lines[0] = sum_xs = 0;
  for (i_line=0; i_line<XRAYLIB_LINES_MAX; i_line++) { // loop on fluorescent lines
    double xs = CSb_FluorLine(Z, -i_line, E0, NULL); /* XRayLib */
    // when a line is inactive: E=xs=0
    sum_xs += xs;
    cum_xs_lines[i_line+1] = sum_xs; // cumulative sum of their cross sections
  }
  // select randomly one of these lines
  i_line = XRMC_SelectFromDistribution(cum_xs_lines, XRAYLIB_LINES_MAX); // extract a line
  // get the K shell line width as approximation of fluorescence line width
  if (dE) *dE = AtomicLevelWidth(Z, K_SHELL, NULL); // keV

  return LineEnergy(Z, -i_line, NULL); // fluorescent line energy
} // XRMC_SelectFluorescenceEnergy

  // Function removing spaces from string
  char * removeSpacesFromStr(char *string)
  {
      // non_space_count to keep the frequency of non space characters
      int non_space_count = 0;
   
      //Traverse a string and if it is non space character then, place it at index non_space_count
      for (int i = 0; string[i] != '\0'; i++)
      {
          if (isalpha(string[i]) || isdigit(string[i]))
          {
              string[non_space_count] = string[i];
              non_space_count++; //non_space_count incremented
          }    
      }
      
      //Finally placing final character at the string end
      string[non_space_count] = '\0';
      return string;
  }

  // ok = fluo_get_material(material, formula)
  // extracts material atoms from file header
  // the result is concatenated into 'formula'
  int fluo_get_material(char *filename, char *formula) {
  
    int  ret = 0;
    char Line[65535];
    int flag_found_cif=0;
    FILE *file = Open_File(filename, "r", NULL);
    if (!file)
      exit(fprintf(stderr, "%s: ERROR: can not open file %s\n", 
      __FILE__, filename));
    
    // Read the file, and search tokens in rows
    formula[0]=0;
    while (fgets(Line, sizeof(Line), file) != NULL) {
      char *token = NULL;
      int   flag_exit=0;
      char *first_non_space=NULL;
      char *next_non_space =NULL;
      // CIF: _chemical_formula_structural 'chemical_formulae'
      // CIF: _chemical_formula_sum 'chemical_formulae'
      // LAZ/LAU: # ATOM <at> <trailing>
      // LAZ/LAU: # Atom <at> <trailing>
      // LAZ/LAU: # TITLE <at> <at> ... [ trailing...]
      // CFL: Title <chemical_formulae>
      // CFL: Atom <at> <trailing>
      
      // search for CIF token
      // single line: search " '\'" delimiter after CIF token, marks reading the formula
      if (!strncasecmp(Line, "_chemical_formula_structural", 28) 
       || !strncasecmp(Line, "_chemical_formula_sum", 21) 
       || !strncasecmp(Line, "_chemical_formula_moiety", 24)
       || flag_found_cif) {
        if  (flag_found_cif) { flag_found_cif=0; /* can not span on more that 2 lines */} 
        else flag_found_cif=1;
        // search for delimiter after the CIF token
        char *first_space_after_token=strpbrk(Line, " \'\n");
        if (first_space_after_token) {
          // search for the characters that may compose the formula
          first_non_space = strpbrk(first_space_after_token, "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789()[]");
          if (strchr(first_space_after_token,'?')) {
            // we have an invalid/unknown formula in the CIF. Skip CIF token/line
            flag_found_cif=0;
            continue;
          }
        }
        if (first_non_space) 
          next_non_space  = strpbrk(first_non_space+1, "\'\n\r"); // position of formula end
        if (next_non_space) { 
          flag_exit = 1;
          token = first_non_space;
        }
      } else if (!strncasecmp(Line, "# TITLE", 7) && strchr(Line, '['))
        token = Line+7;
      else if (!strncasecmp(Line, "# Atom", 6))
        token = Line+6;
      else if (!strncasecmp(Line, "Atom", 4))
        token = Line+4;
      else if (!strncasecmp(Line, "Title", 5))
        token = Line+7;
      if (!token) continue;
      if (!strncasecmp(Line, "# TITLE", 7)) {
        first_non_space = Line+7;
        next_non_space  = strchr(Line+7, '[');
        if (next_non_space) flag_exit = 1;
      }
      if (!first_non_space) first_non_space = strtok(token, " \t\n");
      if (!first_non_space) continue;
      if (!next_non_space)  next_non_space  = strtok(NULL,  " \t\n"); // end of formulae
      if (!next_non_space)  next_non_space = Line+strlen(Line);
      // remove spaces
      strncat(formula, first_non_space, next_non_space-first_non_space);
      ret++;
      if (flag_exit) break;
    }
    fclose(file);
    removeSpacesFromStr(formula);
    return(ret);
  } // fluo_get_material
%}

/* ========================================================================== */

DECLARE %{
  struct   compoundData *compound;
  DArray1d cum_massFractions;
  DArray1d cum_xs_fluo;
  DArray1d cum_xs_Compton;
  DArray1d cum_xs_Rayleigh;
  int  shape;
  off_struct offdata;
  int  n_fluo;
  int  n_Compton;
  int  n_Rayleigh;
  double p_fluo;
  double p_Compton;
  double p_Rayleigh;
  
  // cached values for SPLIT
  int    reuse_nb;
  int    reuse_intersect;
  double reuse_ki;
  double reuse_l0;
  double reuse_l1;
  double reuse_l2;
  double reuse_l3;
  double reuse_sigma_barn;
  DArray1d reuse_cum_xs_fluo;
  DArray1d reuse_cum_xs_Compton;
  DArray1d reuse_cum_xs_Rayleigh;
  char *filename;
%}

INITIALIZE %{

  /* energies en [keV], angles in [radians], XRL CSb cross sections are in [barn/atom] */
  double     E0, dE;
  xrl_error *error = NULL;
  int        i;
  
  XRayInit();
  
  shape=-1; /* -1:no shape, 0:cyl, 1:box, 2:sphere, 3:any-shape  */
  if (geometry && strlen(geometry) && strcmp(geometry, "NULL") && strcmp(geometry, "0")) {
    #ifndef USE_OFF
    fprintf(stderr,"Error: You are attempting to use an OFF geometry without -DUSE_OFF. You will need to recompile with that define set!\n");
    exit(-1);
    #else
    if (off_init(geometry, xwidth, yheight, zdepth, 0, &offdata)) {
      shape=3; thickness=0; concentric=0;
    }
    #endif
  }
  else if (xwidth && yheight && zdepth)  shape=1; /* box */
  else if (radius > 0 &&  yheight)       shape=0; /* cylinder */
  else if (radius > 0 && !yheight)       shape=2; /* sphere */

  if (shape < 0)
    exit(fprintf(stderr,"Fluorescence: %s: sample has invalid dimensions.\n"
                        "ERROR       Please check parameter values (xwidth, yheight, zdepth, radius).\n", NAME_CURRENT_COMP));
  
  if (!material || !strlen(material) || !strcmp(material, "NULL") || !strcmp(material, "0")) 
    exit(fprintf(stderr, "Fluorescence: ERROR: %s: Null material specification\n", NAME_CURRENT_COMP));
    
  // test if the material is given as a file
  char path[1024];
  char formula[65536];
  formula[0]='\0';
  FILE *file = Open_File(material, "r", path);
  if (file != NULL) {
    fclose(file);
    // open the material structure file (laz/lau/cif...)
    // search (case sensitive) along lines
    if (!fluo_get_material(material, formula))
      exit(fprintf(stderr, "ERROR: %s: file %s does not contain material formulae.\n", NAME_CURRENT_COMP, material));
    
    fprintf(stderr, "Fluorescence: INFO: %s: found material %s from file %s\n", NAME_CURRENT_COMP, formula, material);
    strcpy(material, formula);
    // CIF: _chemical_formula_structural 'chemical_formulae'
    // CIF: _chemical_formula_sum 'chemical_formulae'
    // LAZ/LAU: # ATOM <at> <trailing>
    // LAZ/LAU: # Atom <at> <trailing>
    // LAZ/LAU: # TITLE <at> <at> ... [ trailing...]
    // CFL: Title <chemical_formulae>
    // CFL: Atom <at> <trailing>
    
  } else filename = NULL;
    
  compound = CompoundParser(material, &error); /* XRayLib */
  if (error != NULL) {
    exit(fprintf(stderr, "ERROR: %s: Invalid material %s: %s\n",
      NAME_CURRENT_COMP, material, error->message));
  }
  xrl_error_free(error);
  
  /* compute total density for raw material and display information ========= */  
  if (weight <= 0) weight = compound->molarMass; /* g/mol */
  MPI_MASTER(
    printf("%s: Material %s mass fractions:\n", 
      NAME_CURRENT_COMP, material);
  )
  
  double mat_density       = 0; /* g/cm3 */
  double sum_massFractions = 0;
  cum_massFractions        = create_darr1d(compound->nElements+1);
  cum_xs_fluo              = create_darr1d(compound->nElements+1);
  cum_xs_Compton           = create_darr1d(compound->nElements+1);
  cum_xs_Rayleigh          = create_darr1d(compound->nElements+1);
  cum_massFractions[0]     = 0;
  
  /* print material information, and check for elements */
  for (i=0; i< compound->nElements; i++) {
    int    Z      = compound->Elements[i];
    error = NULL;
    double Z_dens = ElementDensity(Z, &error);
    if (error != NULL)
      exit(fprintf(stderr, "ERROR: %s: Z=%i %s\n", NAME_CURRENT_COMP, Z, error->message));
    mat_density           += compound->massFractions[i]*Z_dens;
    sum_massFractions     += compound->massFractions[i];
    cum_massFractions[i+1] = sum_massFractions;
    MPI_MASTER(
      printf("  | %6.2g %%: Z=%3i %3s %8.3g [g/mol] %8.3g [g/cm3]\n",
        compound->massFractions[i]*100, Z, AtomicNumberToSymbol(Z,NULL), AtomicWeight(Z, NULL),
        Z_dens);
    )
  }
  xrl_error_free(error);
  if (density <= 0)        density        = mat_density;   /* g/cm3 */
  if (packing_factor <= 0) packing_factor = density/mat_density;
  
  /* molar volume [cm^3/mol] = weight [g/mol] / density [g/cm^3] */
  /* atom density per Angs^3 = [mol/cm^3] * N_Avogadro *(1e-8)^3 */
  if (!rho) rho = density/weight/1e24*NA; // atom density [at/Angs-3]
  MPI_MASTER(
    printf("%s: Material %s M=%g [g/mol] density=%g [g/cm3] rho=%g [at/Angs-3]",
      NAME_CURRENT_COMP, material, weight, density, rho);
    if (fabs(packing_factor-1) > 1e-2)
      printf(" packing_factor=%g", packing_factor);
    printf("\n");
  )
  if (0 < packing_factor && packing_factor < 1) rho *= packing_factor;
  
  /* target for scattering ================================================== */
  if (!target_index && !target_x && !target_y && !target_z) target_index=1;
  if (target_index)
  {
    Coords ToTarget;
    ToTarget = coords_sub(POS_A_COMP_INDEX(INDEX_CURRENT_COMP+target_index),POS_A_CURRENT_COMP);
    ToTarget = rot_apply(ROT_A_CURRENT_COMP, ToTarget);
    coords_get(ToTarget, &target_x, &target_y, &target_z);
  }
  if (!(target_x || target_y || target_z)) {
    MPI_MASTER(
    printf("Fluorescence: %s: The target is not defined. Using 4PI.\n",
      NAME_CURRENT_COMP);
    );
  }
  
  n_fluo = n_Compton = n_Rayleigh = 0;
  p_fluo = p_Compton = p_Rayleigh = 0;
  
  // cached variables set to 0 (for SPLIT)
  reuse_nb=0;
  reuse_ki=reuse_l0=reuse_l3=0;
  reuse_cum_xs_fluo              = create_darr1d(compound->nElements+1);
  reuse_cum_xs_Compton           = create_darr1d(compound->nElements+1);
  reuse_cum_xs_Rayleigh          = create_darr1d(compound->nElements+1);
%}

TRACE %{

int    intersect=0;       /* flag to continue/stop */
int    reuse=0;
int    type=-1;
double l0,  l1,  l2,  l3; /* times for intersections */
double dl0, dl1, dl2, dl; /* time intervals */
int    flag_concentric = 0;
int    flag_ishollow   = 0;
double sigma_barn=0, xs[4]; /* cross sections [barn/atom] fluo/Compton/Rayleigh */
double aim_x=0, aim_y=0, aim_z=1;   /* Position of target relative to scattering point */
int    event_counter   = 0;         /* scattering event counter (multiple fluorescence) */
int    force_transmit  = 0;         /* Flag to handle cross-section weighting in case of finite order */

#ifdef OPENACC
#ifdef USE_OFF
off_struct thread_offdata = offdata;
#endif
#else
#define thread_offdata offdata
#endif

double ki_x,ki_y,ki_z,ki,Ei, pi;
double kf_x,kf_y,kf_z,kf,Ef;

/* Store Initial photon state */

ki_x = kx;
ki_y = ky;
ki_z = kz;
ki   = sqrt(kx*kx+ky*ky+kz*kz); //  Angs-1
kf   = ki;
Ei   = K2E*ki; // keV
pi   = p;      // used to test for multiple fluo weighting and order cutoff

do { /* while (intersect) Loop over multiple scattering events */

  // test for a SPLIT event (same particle comes in)
  if (!event_counter && fabs(reuse_ki - ki) < 1e-15) {
    reuse     = 1;
    reuse_nb++;
    // use cached values and skip actual computation
    intersect = reuse_intersect;
    l0        = reuse_l0;
    l1        = reuse_l1;
    l2        = reuse_l2;
    l3        = reuse_l3;
  } else {
    // we have a different event: compute intersection lengths
    
    /* ========================================================================== */
    /*                                   GEOMETRY                                 */
    /* ========================================================================== */

    /* Intersection photon trajectory / sample (sample surface) */
    if (thickness >= 0) {
      if (shape==0)
        intersect=cylinder_intersect(&l0,&l3, x,y,z,kx,ky,kz, radius,yheight);
      else if (shape==1)
        intersect=box_intersect     (&l0,&l3, x,y,z,kx,ky,kz, xwidth,yheight,zdepth);
      else if (shape==2)
        intersect=sphere_intersect  (&l0,&l3, x,y,z,kx,ky,kz, radius);
      #ifdef USE_OFF
      else if (shape == 3)
        intersect=off_x_intersect(&l0, &l3, NULL, NULL, x, y, z, kx,ky,kz, thread_offdata );
      #endif
    } else {
      if (shape==0)
        intersect=cylinder_intersect(&l0,&l3, x,y,z,kx,ky,kz, radius-thickness,
          yheight-2*thickness > 0 ? yheight-2*thickness : yheight);
      else if (shape==1)
        intersect=box_intersect     (&l0,&l3, x,y,z,kx,ky,kz,
          xwidth-2*thickness > 0 ?  xwidth-2*thickness : xwidth,
          yheight-2*thickness > 0 ? yheight-2*thickness : yheight,
          zdepth-2*thickness > 0 ?  zdepth-2*thickness : zdepth);
      else if (shape==2)
        intersect=sphere_intersect  (&l0,&l3, x,y,z,kx,ky,kz, radius-thickness);
      #ifdef USE_OFF
      else if (shape == 3)
        intersect=off_x_intersect(&l0, &l3, NULL, NULL, x, y, z, kx,ky,kz, thread_offdata );
      #endif
    }


    /* Computing the intermediate lengths */
    if (intersect && p_interact >= 0) {
      flag_ishollow = 0;
      if (thickness > 0) {
        if (shape==0 && cylinder_intersect(&l1,&l2, x,y,z,kx,ky,kz, radius-thickness,
          yheight-2*thickness > 0 ? yheight-2*thickness : yheight))
          flag_ishollow=1;
        else if (shape==2 && sphere_intersect   (&l1,&l2, x,y,z,kx,ky,kz, radius-thickness))
          flag_ishollow=1;
        else if (shape==1 && box_intersect(&l1,&l2, x,y,z,kx,ky,kz,
          xwidth-2*thickness > 0 ? xwidth-2*thickness : xwidth,
          yheight-2*thickness > 0 ? yheight-2*thickness : yheight,
          zdepth-2*thickness > 0 ? zdepth-2*thickness : zdepth))
          flag_ishollow=1;
      } else if (thickness<0) {
        if (shape==0 && cylinder_intersect(&l1,&l2, x,y,z,kx,ky,kz, radius,yheight))
          flag_ishollow=1;
        else if (shape==2 && sphere_intersect   (&l1,&l2, x,y,z,kx,ky,kz, radius))
          flag_ishollow=1;
        else if (shape==1 && box_intersect(&l1,&l2, x,y,z,kx,ky,kz, xwidth, yheight, zdepth))
          flag_ishollow=1;
      }
      if (!flag_ishollow) l1 = l2 = l3; /* no empty space inside */
    } /* if intersect */


    // store values for potential next SPLIT
    reuse           = 0;
    reuse_intersect = intersect;
    reuse_l0        = l0;
    reuse_l1        = l1;
    reuse_l2        = l2;
    reuse_l3        = l3;
    reuse_ki        = ki;
  } // if !reuse (SPLIT)

  if (intersect) { /* the photon hits the sample */

    if (l0 > 0) {  /* we are before the sample */
      PROP_DL(l0); /* propagates photon to the entry of the sample */
    } else if (l1 > 0 && l1 > l0) { /* we are inside first part of the sample */
      /* no propagation, stay inside */
    } else if (l2 > 0 && l2 > l1) { /* we are in the hole */
      PROP_DL(l2); /* propagate to inner surface of 2nd part of sample */
    } else if (l3 > 0 && l3 > l2) { /* we are in the 2nd part of sample */
      /* no propagation, stay inside */
    }

    dl0=l1-(l0 > 0 ? l0 : 0); /* Time in first part of hollow/cylinder/box */
    dl1=l2-(l1 > 0 ? l1 : 0); /* Time in hole */
    dl2=l3-(l2 > 0 ? l2 : 0); /* Time in 2nd part of hollow cylinder */

    if (dl0 < 0) dl0 = 0;
    if (dl1 < 0) dl1 = 0;
    if (dl2 < 0) dl2 = 0;

    /* initialize concentric mode */
    if (concentric && !flag_concentric && l0 >= 0
     && shape==0 && thickness) {
      flag_concentric=1;
    }

    if (flag_concentric == 1) {
      dl1=dl2=0; /* force exit when reaching hole/2nd part */
    }

    if (!dl0 && !dl2) {
      intersect = 0; /* the sample was passed entirely */
    }
  } // if intersect (geometry)
    
  /* ========================================================================== */
  /*                             INTERACTION PROCESS                            */
  /* ========================================================================== */

  if (intersect) {
    double my_s;
    int    i_Z,i;
    int    flag=0;
    double d_path, p_trans, p_scatt, mc_trans, mc_scatt;
    
    /* actual fluorescence calculation */
    
    /* compute total scattering cross section for incoming photon energy Ei */
    /* compute each contribution XS */
    xs[FLUORESCENCE]=xs[COMPTON]=xs[RAYLEIGH]=xs[TRANSMISSION]=sigma_barn=0;
    cum_xs_fluo[0] = cum_xs_Compton[0] = cum_xs_Rayleigh[0] = 0;
    if (!reuse) {
      for (i_Z=0; i_Z< compound->nElements; i_Z++) {
        int    Z   = compound->Elements[i_Z];
        double frac= compound->massFractions[i_Z];
        double xs_Z[3];

        // get Fluorescence xs
        XRMC_CrossSections(Z, Ei, xs_Z); // [barn/atom]
        sigma_barn            += frac*CSb_Total(Z, Ei, NULL); // Photo+Compton+Rayleigh
        cum_xs_fluo[i_Z+1]     = cum_xs_fluo[i_Z]    +frac*xs_Z[FLUORESCENCE];
        cum_xs_Compton[i_Z+1]  = cum_xs_Compton[i_Z] +frac*xs_Z[COMPTON];
        cum_xs_Rayleigh[i_Z+1] = cum_xs_Rayleigh[i_Z]+frac*xs_Z[RAYLEIGH];
        for (i=0; i<3; i++) { xs[i] += frac*xs_Z[i]; }
      } // for Z in compound
      
      // store values into cache for SPLIT
      for (i_Z=0; i_Z< compound->nElements; i_Z++) {
        reuse_cum_xs_fluo[i_Z+1]     = cum_xs_fluo[i_Z+1];
        reuse_cum_xs_Compton[i_Z+1]  = cum_xs_Compton[i_Z+1];
        reuse_cum_xs_Rayleigh[i_Z+1] = cum_xs_Rayleigh[i_Z+1];
      }
      reuse_sigma_barn = sigma_barn;
    } else {
      // reuse cached values (SPLIT)
      for (i_Z=0; i_Z< compound->nElements; i_Z++) {
        cum_xs_fluo[i_Z+1]     = reuse_cum_xs_fluo[i_Z+1];
        cum_xs_Compton[i_Z+1]  = reuse_cum_xs_Compton[i_Z+1];
        cum_xs_Rayleigh[i_Z+1] = reuse_cum_xs_Rayleigh[i_Z+1];
        xs[FLUORESCENCE]      += reuse_cum_xs_fluo[i_Z+1];
        xs[COMPTON]           += reuse_cum_xs_Compton[i_Z+1];
        xs[RAYLEIGH]          += reuse_cum_xs_Rayleigh[i_Z+1];
      }
      
      sigma_barn = reuse_sigma_barn;
    }
    
    /* probability to absorb/scatter */
    my_s   = rho*100*sigma_barn; /* mu, 100: convert from barns to fm^2. my_s in [1/m] */
    d_path = ( dl0 +dl2 );  /* total path lenght in sample */
    
    /* Proba of transmission/interaction along length d_path */    
    p_trans = exp(-my_s*d_path); /* probability to not-interact (transmit) */
    //printf("sigma_barn=%g p_trans=%g\n", sigma_barn, p_trans);
    p_scatt = 1 - p_trans;       /* portion of beam which scatters */
    
    /* force a given fraction of the beam to scatter */
    if (p_interact>0 && p_interact<=1) {
      /* we force a portion of the beam to interact */
      /* This is used to improve statistics */
      mc_trans = 1-p_interact;
    } else {
      mc_trans = p_trans; /* 1 - p_scatt */
    }
    mc_scatt = 1 - mc_trans; /* portion of beam to scatter (or force to) */
    if (mc_scatt <= 0) ABSORB;
    
    if (!force_transmit && mc_scatt > 0 && (mc_scatt >= 1 || rand01() < mc_scatt)) { 
      /* we "scatter" with one of the interaction processes */
        
      dl = -log(1 - rand0max((1 - exp(-my_s*d_path)))) / my_s; /* length */

      /* If t0 is in hole, propagate to next part of the hollow cylinder */
      if (dl1 > 0 && dl0 > 0 && dl > dl0) dl += dl1;

      /* photon propagation to the scattering point */
      PROP_DL(dl);
      p *= fabs(p_scatt/mc_scatt); /* account for p_interact, lower than 1 */
      
    } else { // force_transmit
      /* we go through the material without interaction, and exit */
      if (type <0) type = TRANSMISSION; // 3 transmission
      intersect = 0;
      PROP_DL(dl0+dl2);
      /* attenuate beam by portion which is scattered (and left along) */
      p *= p_trans;
      if (p_interact>0 && p_interact<=1) p /= mc_trans;
      break; // end while (intersect)
    }
    
  } /* if intersect (propagate) */

  if (intersect) { /* scattering event */
    int    i_Z, Z;
    double solid_angle;
    double theta, dsigma;
    double Ef, dE;
    
    /* correct for XS total(photo+Compton+Rayleigh) > sum(fluo+Compton+Rayleigh) */
    dsigma = (xs[FLUORESCENCE]+xs[RAYLEIGH]+xs[COMPTON])/sigma_barn;
    if (dsigma < 1) p *= dsigma; // < 1

    /* MC choose process from cross sections 'xs': fluo, Compton, Rayleigh */
    type = XRMC_SelectInteraction(xs);
    
    /* choose Z (element) on associated XS, taking into account mass-fractions */
    switch (type) {
      case FLUORESCENCE:
        i_Z = XRMC_SelectFromDistribution(cum_xs_fluo,     compound->nElements+1);
        break;
      case RAYLEIGH:
        if (!flag_rayleigh) ABSORB;
        i_Z = XRMC_SelectFromDistribution(cum_xs_Rayleigh, compound->nElements+1);
        break;
      case COMPTON:
        if (!flag_compton) ABSORB;
        i_Z = XRMC_SelectFromDistribution(cum_xs_Compton,  compound->nElements+1);
        break;
      default:
        printf("%s: WARNING: process %i unknown. Absorb.\n", NAME_CURRENT_COMP, type);
        ABSORB;
    }
    Z   = compound->Elements[i_Z];
    
    /* select outgoing vector */
    if ((target_x || target_y || target_z)) {
      aim_x = target_x-x;       /* Vector pointing at target (anal./det.) */
      aim_y = target_y-y;
      aim_z = target_z-z;
    }
    if(focus_aw && focus_ah) {
      randvec_target_rect_angular(&kf_x, &kf_y, &kf_z, &solid_angle,
        aim_x, aim_y, aim_z, focus_aw, focus_ah, ROT_A_CURRENT_COMP);
    } else if(focus_xw && focus_yh) {
      randvec_target_rect(&kf_x, &kf_y, &kf_z, &solid_angle,
        aim_x, aim_y, aim_z, focus_xw, focus_yh, ROT_A_CURRENT_COMP);
    } else {
      randvec_target_circle(&kf_x, &kf_y, &kf_z, &solid_angle, aim_x, aim_y, aim_z, focus_r);
    }
    p *= solid_angle/(4*PI); // correct for selected solid-angle
    NORM(kf_x, kf_y, kf_z);  // normalize the outout direction |kf|=1
    
    // determine final energy
    switch (type) {
      case FLUORESCENCE: /* 0 Fluo: choose line */
        n_fluo++;
        p_fluo += p;
        Ef      = XRMC_SelectFluorescenceEnergy(Z, Ei, &dE);      // dE in keV
        if (dE) {
          if (flag_lorentzian) dE  *= tan(PI/2*randpm1()); // Lorentzian distribution
          else                 dE  *= randnorm();          // Gaussian distribution
          Ef = Ef + dE;
        }
        kf      = Ef*E2K;
        break;
        
      case RAYLEIGH:     /* 1 Rayleigh: Coherent, elastic    */
        n_Rayleigh++;
        p_Rayleigh += p;
        theta      = acos(scalar_prod(kf_x,kf_y, kf_z,ki_x, ki_y,ki_z)/ki);
        dsigma     = DCSb_Rayl(Z,  Ei, theta, NULL); // [barn/at/st]
        p         *= 4*PI*dsigma/xs[RAYLEIGH];
        break;
      
      case COMPTON:      /* 2 Compton: Incoherent: choose final energy */
        n_Compton++;
        p_Compton += p;
        theta      = acos(scalar_prod(kf_x,kf_y, kf_z,ki_x, ki_y,ki_z)/ki);
        dsigma     = DCS_Compt(Z,  Ei, theta, NULL); // [barn/at/st]
        kf         = ComptonEnergy(Ei, theta, NULL)*E2K; /* XRayLib */
        p         *= 4*PI*dsigma/xs[COMPTON];
        break;
    }
    Ef = K2E*kf;
    
    kx = kf*kf_x;
    ky = kf*kf_y;
    kz = kf*kf_z;
    SCATTER;
    event_counter++;
    
    /* exit if multiple scattering order has been reached */
    if (!order) break; // skip final absorption
    // stop when order has been reached, or weighting is very low
    if (order && (event_counter >= order || p/pi < 1e-7)) { force_transmit=1; }
 
  } // if intersect (scatter)
} while(intersect); /* end do (intersect) (multiple scattering loop) */

%}

FINALLY %{
  FreeCompoundData(compound);
  if (filename && filename != material)
    unlink(filename);
  printf("%s: scattered intensity: fluo=%g Compton=%g Rayleigh=%g\n", 
    NAME_CURRENT_COMP, p_fluo, p_Compton, p_Rayleigh);
%}

END