Problem in monitor_nD banana geometry

[mads-bertelsen]

When Monitor_nD is used with the banana geometry, intersect is still done with a cylinder that contains the top and bottom. In the following code the intersect variable contains the return value from cylinder_intersect, this return variable will be 1 if only intersections with the curved part are found and larger if intersections with the top or bottom is found.

Intersections with the top and bottom needs to be removed and the following code attempts to remove them, but does not seem to work in all cases.

https://github.com/McStasMcXtrace/McCode/blob/62c86860a926668edf12e61fa4760e660152fbdd/mcstas-comps/monitors/Monitor_nD.comp#L479-L511

Cases that can create problems:

1. The intersection routine has some error which is taken into account on line 505, but not in the checks on line 491 and 492.
2. Ray starts outside and intersects with both top and bottom plate. Both the checks on line 491 and 492 are true, and the t0
3. Error on intersection calculation is larger than FLT_EPSILON and in the direction of the inside of the cylinder, then the point is not removed even if it is not located at the edge.

As a temporary measure I have used a larger value than FLT_EPSILON and included this in the checks on line 491 and 492, then the problematic data from the top and bottom disappears. This does remove a tiny fraction of the data from the desired part of the detector.

A proper fix could be a dedicated banana_intersect function that simply is a copy of the cylinder_intersect function without the flat parts, removing the need for the flaky filter code afterwards.

The Monitor_nD code is the same for 2.X and 3.X branches and affect both.

1318 could be used to test. Relevant test cases are beam hitting top and bottom and/or source being outside the cylinder.

[mads-bertelsen]

The issue is further complicated by the use of gravity as PROP_DT will propagate with gravitation, where the temporary point check previously used linear propagation.

[mads-bertelsen]

I have chosen another temporary fix as some edge cases still had too large numerical uncertainty and required relatively large floating point error buffers to fix. Instead I am using the return value from cylinder_intersect to filter the events directly, as it contains information on which face the ray intersected. The code is pasted below, can make a pull request on it if this is the desired direction, though there are a lot of cases to consider.

if (intersect)
  {
    if ((abs(Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_CYLIND)
     || (abs(Vars.Flag_Shape) == DEFS.SHAPE_BOX) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_BANANA)
     || (geometry && strlen(geometry) && strcmp(geometry,"0") && strcmp(geometry, "NULL")) )
    {
      /* check if we have to remove the top/bottom with BANANA shape */
      if (abs(Vars.Flag_Shape) == DEFS.SHAPE_BANANA) {

        if (intersect == 1) { // Entered and left through sides
            if (t0 < 0 && t1 > 0) {
              t0 = t;  /* neutron was already inside ! */
            }
            if (t1 < 0 && t0 > 0) { /* neutron exit before entering !! */
              t1 = t;
            }/* t0 is now time of incoming intersection with the detection area */
            if ((Vars.Flag_Shape < 0) && (t1 > 0)) {
              PROP_DT(t1); /* t1 outgoing beam */
            } else {
              PROP_DT(t0); /* t0 incoming beam */
            }

        } else if (intersect == 3 || intersect == 5) { // Entered from top or bottom
            if ((Vars.Flag_Shape < 0) && (t1 > 0)) {
              PROP_DT(t1); /* t1 outgoing beam */
            } else {
              intersect=0;
              Flag_Restore=1;
            }
        } else if (intersect == 9 || intersect == 17) { // Left from top or bottom
            if ((Vars.Flag_Shape < 0) && (t1 > 0)) {
              intersect=0;
              Flag_Restore=1;
            } else {
              PROP_DT(t0); /* t0 incoming beam */
            }
        } else if (intersect == 13 || intersect == 19) { // Went through top/bottom on entry and exit
            intersect=0;
            Flag_Restore=1;
        }

      } else {
        if (t0 < 0 && t1 > 0)
          t0 = t;  /* neutron was already inside ! */
        if (t1 < 0 && t0 > 0) /* neutron exit before entering !! */
          t1 = t;
        /* t0 is now time of incoming intersection with the detection area */
        if ((Vars.Flag_Shape < 0) && (t1 > 0))
          PROP_DT(t1); /* t1 outgoing beam */
        else
          PROP_DT(t0); /* t0 incoming beam */
        /* Final test if we are on lid / bottom of banana/sphere */
        if (abs(Vars.Flag_Shape) == DEFS.SHAPE_BANANA || abs(Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) {
          if (Vars.Cylinder_Height && fabs(y) >= Vars.Cylinder_Height/2 - FLT_EPSILON) {
            intersect=0;
            Flag_Restore=1;
          }
        }
      }
    }
  }

[willend]

@mads-bertelsen I think the latter approach looks nice enough, and a PR is of course most welcome.

Wrt. the more general issue of the intersection-routines e.g. cylinder_intersect vs. PROP_DT in presence of gravity, the effect is there, but for the problem in question typically of quite small order. (Detectors of dimension ~a few m.)

One way to solve this in more general terms would be to create parabolic intersection routines for every geometry (could become algorithmically challenging) or at least within free space potentially something along the lines of

1. Use the normal intersection routine to approximate distance to the object / find the relevant "entry-plane"
2. Solve eq of motion for intersection time with that plane for a better estimate
3. Use PROP_DT as usual

... or some sort of iterative, piecewise solution with linear propagation (which of course easily becomes computationally heavy).

[mads-bertelsen]

We can certainly work on a more general way to handle gravity that would eliminate these types of logic issues that are really easy to miss.

The gravity issue in this case just needs to be larger than FLT_EPSILON to become a problem, as the ray would fail to reach the top of the cylinder and it would then be counted.

Will create the PR for both versions today and we can continue discussion on the new algorithm there.

[willend]

@mads-bertelsen looking at the status of the PR's mentioned, I believe we should be "ready for launch" on this issue also?