Created attachment 63 [details]
escaping photon

It seems that G4TessellatedSolid::DistanceToIn does not check "on surface" condition.
http://www-geant4.kek.jp/lxr-dev/source/geometry/solids/specific/src/G4TessellatedSolid.cc#L672

If I modified the code as follows,

      if (distFromSurface > 0.5*kCarTolerance && dist >= 0.0 && dist < minDist)
      {
        minDist  = dist;
      }
      if(-0.5*kCarTolerance <= dist && dist <= 0.5*kCarTolerance){
        return 0.0;
      }

photons are reflected at boundaries. But after a few reflection, they escape the BGO block again. Before they escape, "StepTooSmall" process is repeated several times. I attache an image which shows this situation.

Created attachment 64 [details]
verbose output of optical process

When a photon escapes from BGO block as shown in the attached image (id=64), std::sqrt(sqrDist) is larger than 0.5*kCarTolerance (see G4TriangularFacet::Intersect), though it must be small enough because the photon is at a surface. Therefore, the distance to next boundary is set to kInfinity.

Any idea?

I'm currently investigating this and modifying G4TessellatedSolid, G4TriangularFacet and G4VFacet.

*** Problem 1105 has been marked as a duplicate of this problem. ***

(In reply to comment #6)
> *** Problem 1105 has been marked as a duplicate of this problem. ***

Confirmed.

I've investigated this problem a few weeks ago and was not able to find any report on bugzilla. I have not expected a new on either.

In bug 1105 I reported a "quick fix". Here is the complete list of changes to the sources of Geant 4.9.3:

** Changes to G4TriangularFacet.cc:

Line 544:  
change "" distance = -0.5*kCarTolerance; ""  
to "" distance = 0.0; ""

Insert after Line 697 :
"" distance=std::fabs(distance); ""

**In G4TessellatedSolid.cc:

Line 672:

change the following line:

"" if (distFromSurface > 0.5*kCarTolerance && dist >= 0.0 && dist < minDist)if (distFromSurface > 0.5*kCarTolerance && dist >= 0.0 && dist < minDist) ""

to 

"" if (distFromSurface > -0.5*kCarTolerance && dist >= 0.0 && dist < minDist) ""


I did a quick DIRC simulation with this fix and the results look fine (qualitative). Photons propagate through the radiator like a charm. Anyhow, I have not done a quantitative analysis yet. I have to do finish other code first.

Hope this will help until an official solution is available.

There appear to be three errors which need to be treated:

(1) The treatment of rays entering the G4TessellatedSolid for points starting at the surface is incorrect (as stated by the bug-reporters).

(2) In G4TriangularFacet the distance to the surface is within +/- kCarTolerance it is set to -0.5 kCarTolerance when really it should be 0.0 (just confirmed by Oliver).

(3) The calculation of sqrDist in G4TriangularFacet::Distance(p) can suffer rounding errors when p is very close to the trianglar facet.  in some cases it's better to use (D+s*E[0]+t*E[1]).mag2() instead, although this too may suffer rounding errors!!

Before I implement and distribute a bug correction, I need to confirm something offline with Gabriele and would appreciate a copy of the application (or a simplified version) which Oliver is using.  Is this possible Oliver or is the application too difficult to port?

Pete,

There is a modified version of exampleN06 posted at:

http://hypernews.slac.stanford.edu/HyperNews/geant4/get/opticalphotons/309/1/1.html

which also shows the problem.

Peter

Hello Pete.

Thank you for looking at this problem. I can also help you in debugging with my simulation code which is focused in only optical photon propagation in BGO scintillation blocks which is created from CAD files. If you already have modified versions of G4TessellatedSolid and G4TriangularFacet, and it is possible to distribute it by private e-mails, I will test them.

(In reply to comment #9)

> Before I implement and distribute a bug correction, I need to confirm something
> offline with Gabriele and would appreciate a copy of the application (or a
> simplified version) which Oliver is using.  Is this possible Oliver or is the
> application too difficult to port?

My main application would be a little bit overkill I guess, but I could send you a minimal example which propagates a single photon (ore more) in an G4 extruded solid (in fact a tessellated) octagon by total internal reflection. This volume can be replaced by a polyhedra to compare the results. 

Bug (3) occurs when I decrease the photons angle relative to the face normal, therefore it was not apparent to me in my previous testcase. I can reproduce this problem, and for me it looks like G4ThreeVector v = Distance(p) in G4TriangularFacet::Distance(const G4ThreeVector,const G4double,const G4bool) has rounding errors and therefore the projection on the face normal (dir) indicates "wrong direction". In that case dist will remain equal to kInfinity.

I think this could be the main problem, because rounding errors in sqrDist seem to occur all the time - even on "working reflections".

Ok, after a deeper look into the code I think the following change in G4TriangularFacet::Distance(const G4ThreeVector&) could help until a more performant solution is available:

G4ThreeVector v = D + s*E[0] + t*E[1];
G4double vmag2 = v.mag2();
if(sqrDist>vmag2) sqrDist = vmag2;
return v;

I don't like the additional mag2, but with this change I got my second testcase working.

(In reply to comment #13)
> Ok, after a deeper look into the code I think the following change in
> G4TriangularFacet::Distance(const G4ThreeVector&) could help until a more
> performant solution is available:
> G4ThreeVector v = D + s*E[0] + t*E[1];
> G4double vmag2 = v.mag2();
> if(sqrDist>vmag2) sqrDist = vmag2;
> return v;
> I don't like the additional mag2, but with this change I got my second testcase
> working.

This is the approach I've been testing.  I agree that the calculation from v.mag2() is less desirable, but not sure there's a simple way around it.

> This is the approach I've been testing.  I agree that the calculation from
> v.mag2() is less desirable, but not sure there's a simple way around it.

I assume that only some of the sqrDist computations contribute to the problem, eg.:

s       = s / det;
t       = t / det;
sqrDist = s*(a*s + b*t + 2.0*d) + t*(b*s + c*t + 2.0*e) + f; 

in "region 0". If this is correct (I've not testet enough constellations), it would be enough to replace them with ""G4ThreeVector v = D + s*E[0] + t*E[1]; sqrDist=v.mag2(); return v; ""

In fact, changing the above sqrDist computation in region 0 is enough to avoid the bug in my second testcase - but this may change at different incident angles or positions. And replacing this computation by a mag2() should be even faster, right? 

Anyway, this would need more testing that you might have done already.

1) Replacing only the computation in "region 0" works fine for me (tested multiple incident angles, but always larger than 48 deg !!). 

2) I found another problem (4) that causes this bug: 
The quads Distance(const G4ThreeVector&, const G4double, const G4bool) has problems with rounding errors in dist. I replaced this code with the same logic as in G4TriangularFacet what solved all observed problems. Again, this has to be tested with higher statistics.

I've been trying to solve the problem by making the equations which determine sqrDist more accurate.  In general terms, the roundoff error results from the nature of double-precision numbers A, B and C which are such that:

sqrDist = A + B + C

but when sqrDist is very small:

(A + B) + C != A + (B + C) i.e. commutativity breaks down

I've implemented a fast method of adding the numbers in the right order to always maximise accuracy ... but this still doesn't resolve the problem because of the limiting accuracy of double defining A, B and C.

So we're stuck with calculating sqrDist from the square of the magnitude of the vector D + s*E[0] + t*E[1].

I'll upload the source code later this morning.

(In reply to comment #17)
> I'll upload the source code later this morning.

Hello Pete,

Thank you so much for using your time for tracking this problem. Is the new code available in public or for the developer team only? I would love to try the new one for my purpose.

Akira

Created attachment 66 [details]
vis output of modified exampleN06 using optPhoton.mac

Created attachment 67 [details]
vis output of MinimumTest3

Created attachment 68 [details]
Modified G4TessellatedSolid.cc to deal with optical photons

Created attachment 69 [details]
Modified G4TriangularFacet.cc to overcome bugs with optical photons

Created attachment 70 [details]
Modified ExN06DetectorConstruction.cc with correction for -Z face definition

I've just attached updated G4TessellatedSolids.cc and G4TriangularFacet.cc files updated as discussed in comments #9 and #17.  Thanks also to Akira and Oliver for comments and suggested corrections.  It seems to work OK with MinimumTest (see graphic), I stumbled a bit with the updated N06 example Peter pointed me to on the SLAC web-site - I needed to update the ExN06DetectorConstruction from Akira (-Z side facets pointing in wrong direction) to get total internal reflection on all surfaces until photon absorption (see graphic & attached C++ file).

I note that in the detailed tracking for the optical photon, there are many instances of "StepTooSmall" as mentioned by Akira.  However, I'm assuming this is a natural artifact of the reflection process where the photon doesn't move and only its direction is changed in the step.

I'll get this uploaded to the G4 repository this afternoon, hopefully, with Fan Lei's help.

(In reply to comment #24)

Nice to see an official patch. One question: 

I'm currently using a self patched version that seems to work fine - but to get rid of some rounding errors, I had to change the Quad's code too. Do you use Quads or Triangles in your test geometry? 

> I note that in the detailed tracking for the optical photon, there are many
> instances of "StepTooSmall" as mentioned by Akira.  However, I'm assuming this
> is a natural artifact of the reflection process where the photon doesn't move
> and only its direction is changed in the step.

Exactly. As far as I understand the boundary process, this "step to small" HAS to occur on every reflection because the navigator assigns the reflected photon to the neighbour volume. In the next step the photon is still located on the boundary and reenters the first volume, so the resulting steplength is 0 in a perfect world, or "too small" in floats.

(In reply to comment #25)

Thanks for the comment, Oliver.  I did recall your ealier comment on  G4QuadrangularFacet and was uncertain whether this really needed updating as well.  Therefore during lunchtime I did re-implement the test geometry from Akira using Quads rather than triangular facets and repeated the run.  For me this worked fine without any modification to G4QuadrangularFacet.  Furthermore, I cannot see why Quads needs to be changed because .... it cheats!  Internally, a Quad is treated as two triangular facets.  However, do you think I'm missing something, Oliver?

(In reply to comment #26)
> ... For me
> this worked fine without any modification to G4QuadrangularFacet.  Furthermore,
> I cannot see why Quads needs to be changed because .... it cheats!  Internally,
> a Quad is treated as two triangular facets.  

I would say: yes and no. The Quad wraps two triangles, but not all functions are exactly equal to the implementation in G4TriangularFacet.

Distance(const G4ThreeVector&, const G4double, const bool) behaves not exactly like the function for triangles: The Quad one uses kDirTolerance, while the triangle uses kCarTolerance and sqrt(sqrDist) which is not available in the Quad. 

When I tested my fix, I still observed a few photon tracks that left the volume. All of them were reflected at Quads. Triangles worked fine. After some debugging I decided to go dirty, make Quad a friend of the triangle and use sqrDist and kCarTolerance like its done in G4TriangularFacet::Distance( , , )   

After these changes, everything worked fine. Therefore I think there might be a problem with that function. 

> However, do you think I'm missing
> something, Oliver?

I'm not sure. I don't see why your changes to the Distance computation should solve this specific problem.

However, I try to find some time to check these issues again with a clean geant tree and your changes.

Created attachment 71 [details]
lost photon after reflection at a Quad

The last steps at the bottom of the file show that the photon is not reentering the radiator volume after total internal reflection.

I used the files uploaded by Pete with a clean geant4 tree.

I'll try to find out why this happens. My previous solution was based on an "educated guess" and is not usable in a mainstream release. I think this has to be understood.

The dirTolerance was slightly too small in the given case (dir larger than 1E-14 but less than 1E-13, so same order of magnitude).

Is there any reason why dirTolerance is fixed to 1E-14 ? Maybe its just too small and therefore causing problems in case of precision errors. I'll try 1E-13.

And why are quads and triangles using different distance implementations, isn't it reasonable to use the better/faster/more precise one for both?

"dirTolerance = 1E-13" resolved all issues I have observed, but maybe its still too small in other situations.

(In reply to comment #29)
> The dirTolerance was slightly too small in the given case (dir larger than
> 1E-14 but less than 1E-13, so same order of magnitude).
> Is there any reason why dirTolerance is fixed to 1E-14 ? Maybe its just too
> small and therefore causing problems in case of precision errors. I'll try
> 1E-13.

It's set to 1E-14 to be close to the tolerance for correct computation of double precision numbers.  I believe it was originally introduced to determine the transition from 3D triangle-line intersection to 2D triangle-line intersection.  However, I agree that some tuning may be necessary based on more extensive experience.

> And why are quads and triangles using different distance implementations, isn't
> it reasonable to use the better/faster/more precise one for both?

Quad uses the fact that it's defined as two separate triangles, and calculates the distance for each, but then taking the minimum.  I developed G4TessellatedSolids etc purely as a prototype and purely for as an internal project for simulating microelectronics geometries.  At that time I did want to develop G4QuadrangularFacet in a similar manner to G4TriangularFacet, but for a prototype it was ultimately quicker set the quad to be two adjacent and co-planar triangles.  You're absolutely correct that, given time, one should have a dedicated an efficient algorithm to treat distances for quads.  I guess I should say, "watch this space..."

The fix is now available in release 9.4-beta.