Different hole-ice correction mechanism: Generic propagation through different media

[fiedl]

Current hole-ice-correction mechanism ("Algorithm (a)")

1. Randomize relative scattering and absorption lengths for one simulation step
2. Loop over layers of different ice properties and calculate physical distance to the next interaction point
3. Calculate distance through the hole ice relative to the overall distance between two interaction points
4. Calculate correction for the distance between these two interaction points based on the hole-ice properties relative to the ice properties outside

The current hole-ice-correction mechanism has a number of shortcomings:

• Relative rather than absolute hole-ice properties
• Large errors when jumping over several ice layers with largely differing ice properties
• Errors when several hole-ice cylinders (nested or stacked) are between two interaction points

These issues could be solved by using a different hole-ice-correction mechanism ("Algorithm (b)"):

1. Randomize relative scattering and absorption lengths for one simulation step
2. Determine all medium changes (layers, hole ice) between two interaction points
3. Loop over all medium changes with different ice properties and calculate physical distance to the next interaction point

But this would require a significant rewrite of both hole_ice.c and propagation_kernel.c.cl.

[fiedl]

Back to the roots

When this code snippet replaces the code that is responsible for handling propagating through different media

// PROPAGATION THROUGH DIFFERENT MEDIA 2018: Layers, Cylinders
// -----------------------------------------------------------------------------

// After this code, we need:
// - abs_lens_left
// - distancePropagated

// abs_lens_left is already given.
floating_t sca_step_left = -my_log(RNG_CALL_UNIFORM_OC);
int currentPhotonLayer = min(max(findLayerForGivenZPos(photonPosAndTime.z), 0), MEDIUM_LAYERS-1);
floating_t distancePropagated = sca_step_left * getScatteringLength(currentPhotonLayer, photonDirAndWlen.w);
floating_t distanceToAbsorption = abs_lens_left * getAbsorptionLength(currentPhotonLayer, photonDirAndWlen.w);
{
  if (distanceToAbsorption < distancePropagated)
  {
    distanceToAbsorption = ZERO;
  } else {
    distanceToAbsorption -= distancePropagated;
  }

  abs_lens_left = distanceToAbsorption / getAbsorptionLength(currentPhotonLayer, photonDirAndWlen.w);
}

printf("HOLE ICE 2018 DEBUG\n");
printf("  sca_step_left = %f\n", sca_step_left);
printf("  distancePropagated = %f\n", distancePropagated);
printf("  abs_lens_left = %f\n", abs_lens_left);
printf("  distanceToAbsorption = %f\n", distanceToAbsorption);

i.e. there are no medium changes handled at all, then the propagation looks like this:

[fiedl]

First attempt to a medium-changes loop

// PROPAGATION THROUGH DIFFERENT MEDIA 2018: Layers, Cylinders
// -----------------------------------------------------------------------------

// We know how many scattering and absorption lengths the photon will
// travel in this step. But these lengths are local properties.
// Therefore, we need to loop over all media in range and convert
// these into geometrical distances.

// After this code, we need:
// - abs_lens_left
// - distancePropagated

// abs_lens_left is already given.
floating_t sca_step_left = -my_log(RNG_CALL_UNIFORM_OC);
floating_t distancePropagated = 0;
floating_t distanceToAbsorption = 0;
{
  int number_of_medium_changes = 2;
  floating_t distances_to_medium_changes[MEDIUM_LAYERS] = {0.0, 1.0, 1.5};
  floating_t local_scattering_lengths[MEDIUM_LAYERS] = {100.0, 100.0, 100.0};
  floating_t local_absorption_lengths[MEDIUM_LAYERS] = {100.0, 100.0, 0.0};

  for (int j = 0; (j < number_of_medium_changes) && (sca_step_left > 0); j++) {
    const floating_t max_distance_in_current_medium = distances_to_medium_changes[j+1] - distances_to_medium_changes[j];
    sca_step_left -= my_divide(max_distance_in_current_medium, local_scattering_lengths[j]);
    abs_lens_left -= my_divide(max_distance_in_current_medium, local_absorption_lengths[j]);
    distancePropagated += max_distance_in_current_medium;
    distanceToAbsorption += max_distance_in_current_medium;
    if (sca_step_left < 0) {
      distancePropagated += sca_step_left * max_distance_in_current_medium;
      sca_step_left = 0;
    }
    if (abs_lens_left < 0) {
      distanceToAbsorption += abs_lens_left * max_distance_in_current_medium;
      abs_lens_left = 0;
    }
  }

  // Spend the rest of the budget with the last medium properties.
  distancePropagated += sca_step_left * local_scattering_lengths[number_of_medium_changes];
  distanceToAbsorption += abs_lens_left * local_absorption_lengths[number_of_medium_changes];

  if (distanceToAbsorption < distancePropagated) {
    distancePropagated = distanceToAbsorption;
    distanceToAbsorption = ZERO;
    abs_lens_left = ZERO;
  }
}

In a first toy scenario, configure each photon for instant absorption after 1.5 metres.

[fiedl]

First hole-ice attempt

Using this new mechanism, add hole-ice cylinders as medium boundaries: https://github.com/fiedl/clsim/commit/fdeb328bf6bf7b8a1cb236ce71382247a7eb3a73

Doesn't look too bad, but has no ice layers, no ice tilt and no anisotropy yet:

[fiedl]

Nested cylinders

With the new mechanism, nesting cylinders is possible. https://github.com/fiedl/clsim/commit/7a93d4dd20dd35ae11cd8a4bf4a330817e292894 adds hard-coded cylinder ice properties for testing:

An outer cylinder is configured for a small scattering length causing curly paths within the outer cylinder. An inner cylinder, which has a 2cm smaller radius, is configured for instant absorption.

/CC @thoglu Nested cylinders are working well.

[fiedl]

First angular-acceptance test

Test with cable (instant absorption) and bubble column (scattering length: 5mm), direct detection and the new hole-ice code. Configuration: https://github.com/fiedl/hole-ice-study/commit/8bff8d91520c0ae614ffb55f8a2251beb206c6a1.

[2018-03-09 23:20:38] fiedl@wgs16 /afs/ifh.de/group/amanda/scratch/fiedl/hole-ice-study/scripts/AngularAcceptance master ⚡
▶ qsub -l gpu -l tmpdir_size=10G -l s_rt=8:00:00 -l h_rss=2G -m ae batch-job.sh --distance=1.0 --number-of-photons=1e5 --number-of-runs=2 --number-of-parallel-runs=2 --plane-wave

[2018-03-10 00:39:44] fiedl@fiedl-mbp ~/hole-ice-study/scripts/AngularAcceptance master ⚡
▶ lib/plot.rb ~/z/hole-ice-study/scripts/AngularAcceptance/cluster-results
▶ open ~/z/hole-ice-study/scripts/AngularAcceptance/cluster-results/tmp/plot_with_reference.png

[fiedl]

Another test run with 30cm hole-ice radius. Looks even worse.

[fiedl]

When applying parameters close to the previous ones, e.g. this plane-wave example, then the results look good.

old algorithm:

new algorithm:

[fiedl]

Bringing back the ice layers

Now, handling the ice layers needs to be baked into the new code for medium changes.

First experiment, starting photons from below right and shooting them onto an ice layer, which is configured for instant absorption. Scattering and absorption lengths in the other layers is set to 1km.

In this second experiment, the top ice layer isn't set for instant absorption but for its regular properties, whereas the the other ice layers are configured for an interaction length of 1km.

[fiedl]

The following simulations are started with the same conditions and the same randomizing seed.

[2018-03-10 17:50:09] fiedl@fiedl-mbp ~/hole-ice-study/scripts/FiringRange master ⚡
▶ be ruby run.rb --distance=1.0 --number-of-runs=1 --number-of-parallel-runs=1 --save-photon-paths --cpu --angle=45 --number-of-photons=100 --plane-wave

With ice layers turned off:

With ice layers turned on:

As the starting conditions are the same, these shapes look almost the same. But there are several differences in details like little shifts or shortenings:

[fiedl]

Hole ice scattering experiment. Layers are turned on in both runs.

With 100 photons:

Without hole ice:

With hole ice:

With only 2 photons:

Without hole ice:

With hole ice:

[fiedl]

Ice tilt and ice anisotropy have been removed during this rewrite and need to be brought back. But I'll treat this as separate issue https://github.com/fiedl/hole-ice-study/issues/48 because this will be done at a later timer.

Apart from this, the rewrite is done.