Chroma Reemission
====================


`chroma/chroma/cuda/photon.h`::

    335     // absorption 
    336     //   #. advance .time and .position to absorption point
    337     //   #. if BULK_REEMIT(CONTINUE) change .direction .polarization .wavelength
    338     //   #. if BULK_ABSORB(BREAK)  .last_hit_triangle -1  
    339     //
    340     //  huh, branch BULK_REEMIT(CONTINUE) does not set .last_hit_triangle -1 ?
    341     //
    342     if (absorption_distance <= scattering_distance) {
    343         if (absorption_distance <= s.distance_to_boundary)
    344         {
    345             p.time += absorption_distance/(SPEED_OF_LIGHT/s.refractive_index1);
    346             p.position += absorption_distance*p.direction;
    347 
    348             float uniform_sample_reemit = curand_uniform(&rng);
    349             if (uniform_sample_reemit < s.reemission_prob)
    350             {
    351                 p.wavelength = sample_cdf(&rng, s.material1->n,
    352                                           s.material1->wavelength0,
    353                                           s.material1->step,
    354                                           s.material1->reemission_cdf);
                        //
                        // generating a wavelength according to 
                        // the source distribution used in construction 
                        // of the reemission_cdf
                        //
    
    355                 p.direction = uniform_sphere(&rng);
    356                 p.polarization = cross(uniform_sphere(&rng), p.direction);
    357                 p.polarization /= norm(p.polarization);
    358                 p.history |= BULK_REEMIT;
    359                 return CONTINUE;
    360             } // photon is reemitted isotropically
    361             else
    362             {
    363                 p.last_hit_triangle = -1;
    364                 p.history |= BULK_ABSORB;
    365                 return BREAK;
    366             } // photon is absorbed in material1
    367         }
    368     }


See env/geant4/geometry/collada/g4daeview/sample_cdf.py 
for review of `sample_cdf` using inverse CDF sampling to generate distributions. 

`chroma/chroma/cuda/random.h`::

     25 // Draw a random sample given a cumulative distribution function
     26 // Assumptions: ncdf >= 2, cdf_y[0] is 0.0, and cdf_y[ncdf-1] is 1.0
     27 __device__ float
     28 sample_cdf(curandState *rng, int ncdf, float *cdf_x, float *cdf_y)
     29 {
     30     return interp(curand_uniform(rng),ncdf,cdf_y,cdf_x);
     31 }
     32 
     33 // Sample from a uniformly-sampled CDF
     34 __device__ float
     35 sample_cdf(curandState *rng, int ncdf, float x0, float delta, float *cdf_y)
     36 {
     37     float u = curand_uniform(rng);

            //
            // u, uniform random number between 0 and 1
            //
     38 
     39     int lower = 0;
     40     int upper = ncdf - 1;
     41 
     42     while(lower < upper-1)
     43     {
     44         int half = (lower + upper) / 2;
     45 
     46         if (u < cdf_y[half])
     47             upper = half;
     48         else
     49             lower = half;
     50     }

            //
            // find bin of cdf_y that straddles u   
            //
            //     * cdf_y needs to range from 0:1 
            //       (ie a cumulative distribution normalized to 1 at RHS)
            //

     51 
     52     float delta_cdf_y = cdf_y[upper] - cdf_y[lower];
     53 

            //  translate bin int into domain float 

     54     return x0 + delta*lower + delta*(u-cdf_y[lower])/delta_cdf_y;
     55 }