| Summary: | fails to find intersection between line and cone | ||
|---|---|---|---|
| Product: | Geant4 | Reporter: | Richard Jones <richard.t.jones> |
| Component: | geometry/solids | Assignee: | Evgueni.Tcherniaev |
| Status: | RESOLVED FIXED | ||
| Severity: | normal | ||
| Priority: | P4 | ||
| Version: | 10.4 | ||
| Hardware: | All | ||
| OS: | All | ||
| Attachments: |
example demonstrating the problem
Requested log file Fixed G4IntersectingCone.cc |
||
Thank you for reporting a problem and suggesting a solution, I'll investigate the problem. Meanwhile could you provide us a simple example where the problem appears? Ideally, an example of problematic parameters for LineHitsCone1() and LineHitsCone2(). Created attachment 536 [details]
example demonstrating the problem
This is example B1 from the G4 distribution that has been slightly tweaked to demonstrate the bug in G4IntersectingCone.
I have attached a tarball of the B1 example, modified to demonstrate this bug. The principal changes I made are the following, although maybe not all of them were needed to show this issue. 1) change Shape1 from G4Cons to G4Polycone, and add a slight conical pitch to the entry face, instead of being perfectly flat. 2) in the single particle gun, shift the y origin so that it is pointed at the middle of the entry face of Shape1 3) give the single particle gun a small off-axis momentum component so it is not perfectly aligned with the polycone axis. Here is how to reproduce the problem. 1) "make" to build exampleB1-bug-2111 2) "exampleB1-bug-2111 | tee log", wait for the interactive prompt 3) /tracking/verbose 2 4) /run/beamOn 10 5) (escape from graphics) 6) exit Next examine the log, looking at what happens to Track ID=1 after the first step. It should hit the entrance face to Shape1 at z=-100mm, but most of the time it never sees Shape1, passes right through like Shape1 was not there. Now do "mv src/G4IntersectingCone.cc-updated src/G4IntersectingCone.cc" followed by "make", then redo the above steps. Now you should see that the primary track does not pass through the plane at z=-100mm without entering Shape1. -Richard Jones Dear Richard, I can not reproduce the problem, in my case the log files are identical. What version of Geant4 do you use? Could you send me problematic p and v for LineHitsCone(p,v,s1,s2), that I check them with current development version. Evgueni Evgueni, The example is based on B1 from v10.02.p02, and compiled against that release, but I guess that doesn't matter as not a lot is happening right now in the CSG geometry code. Did you make sure that the G4IntersectingCone.cc was completely removed from the src directory of my exampleB1-bug-2111 directory before you did the make? -Richard I've just re-done everything from scratch to be sure that G4IntersectingCone.cc is not in the src/ folder. The log files are identical. There was one fix in Geant4 10.02 which could affect the result: 07-Mar-2016 E.Tcherniaev (geom-specific-V10-02-05) - Fix in G4IntersectingCone for smaller precision constant in LineHitsCone1() and LineHitsCone2(). Fixes issue observed in G4GenericPolycone. Addressing problem report #1794. Evgueni Evgueni, The source G4IntersectingCone.cc that is in the latest release 10.5 still fails my test. Can you post your log output? I think you are missing something very basic that fails to reproduce my result, but without seeing what you are doing I cannot point it out. -Richard J. Created attachment 538 [details]
Requested log file
Richard, Please find attached a log file produced with modified G4IntersectingCone.cc Evgueni Evgueni, The log file you posted is failing in the tracking for events #1, #2, #3, #4, #5, #7, and #10. It succeeds on #6 for some reason, and on #8. Event #9 has an interaction upstream of the Polycone so it does not test the feature. This demonstrates the bug I am seeking to correct. We should be seeing the incident gamma ray transitioning from Envelope into Shape1 at plane z=-100mm. If it doesn't, that is a serious bug because it completely misses the material in Shape1. My code in G4IntersectingCone.cc-updated is apparently not being correctly linked into your executable when you run this test. Can you add a print statement in G4IntersectingCone::LineHitsCone1 in my modified G4IntersectingCone.cc to make sure that the updated copy is what gets executed in this context? -Richard Jones Richard, Indeed G4IntersectingCone.cc-updated was compiled, but has not been linked. However it does not matter. I've got all necessary information and understood the issue. In case of conical surface there are two extremes: - first, it may become a cylindrical surface - second, it may become a plane Your example demonstrates that current version of G4IntersectCone is not robust enough when the conical surface is almost a plane and should be revised. I'll do it. Evgueni PS. Minor remark: I believe that your example is a bit artificial, not real. Otherwise I would suggest you to get rid of small shift of 1e-6*mm in the definition of zplane[3]. The distance of 1e-6*mm corresponds to few atomic diameters, it has no sense to apply it to definition of real geometrical objects. Created attachment 545 [details]
Fixed G4IntersectingCone.cc
Hi Richard, The source of the problem in the tracking through your geometry was an inaccurate value of the radical caused by loss of precision in calculating the coefficients of the quadratic equation. After applying more accurate procedure of calculation the problem has disappeared, the log file is now identical to the log file produced with G4IntersectingCone.cc-updated. I've attached the modified G4IntersectingCone.cc. A couple of comments: 1. Zero radical. I did not find problems in the current logic of handling radical=0, it covers both "intersection" and "touching". In Geant4 context, "touching" is considered as no intersection. 2. Roundoff error. In computing the relative roundoff error is DBL_EPSILON (in our case EPS), so if radical is near 0, the absolute computation error is order of EPS*b*b. In current code it is used EPS*b, however I left it unchanged. Thank you again for reporting a very interesting algorithmic issue. Evgueni PS. I've made a small presentation relative accurate calculation of radical, you can find it at the following URL: https://indico.cern.ch/event/797876/ |
There are a couple of problems with the handling of rounding errors in the solving the quadratic formula for the intersection between a line and a cone. These produce frequent failures in the tracking of charged particles through my geometry, where I have G4Polycone objects with very steep (but not vertical) faces. In terms of G4IntersectingCone data members, this is the condition where A is of scale 1 and B is of scale 1e-7. In this case, it is very frequent that no intersection is found, even though radical=0 indicates that there is exactly one root. In general there should be no situation under which radical=0 indicates an error -- it should be interpreted to mean that there is one solution which should be returned. Keep in mind that B==0 always gives radical = 0, so it frequently happens when B is much smaller than A. The second problem is with handling the rounding scale, should be order EPS**2, not EPS. Here are my suggested changes. With these, the tracking through my geometry is working again, and the G4 visualization + ray tracing is once again working as expected. $ diff G4IntersectingCone.cc-orig G4IntersectingCone.cc 229c229 < if (radical < -EPS*std::fabs(b)) { return 0; } // No solution --- > if (radical < -EPS*EPS*b*b) { return 0; } // No solution 231c231 < if (radical < EPS*std::fabs(b)) --- > if (radical < EPS*EPS*b*b) 238,239c238,239 < if(B==0.) { return 0; } < if ( std::fabs(x0*ty - y0*tx) < std::fabs(EPS/B) ) --- > //if(B==0.) { return 0; } > //if ( std::fabs(x0*ty - y0*tx) < std::fabs(EPS/(B+1e-99)) ) 333c333 < if (radical < -EPS*std::fabs(b)) { return 0; } // No solution --- > if (radical < -EPS*EPS*b*b) { return 0; } // No solution 335c335 < if (radical < EPS*std::fabs(b)) --- > if (radical < EPS*EPS*b*b) 342c342 < if ( std::fabs(x0*ty - y0*tx) < std::fabs(EPS/B) ) --- > //if ( std::fabs(x0*ty - y0*tx) < std::fabs(EPS/(B+1e-99)) )