Problem 1809

Summary: G4GenericTrap bug in CalculateExtent
Product: Geant4 Reporter: whit
Component: geometry/solidsAssignee: Evgueni.Tcherniaev
Status: RESOLVED FIXED    
Severity: major CC: Alberto.Ribon, Evgueni.Tcherniaev
Priority: P5    
Version: 10.2   
Hardware: All   
OS: All   
Attachments: Bug fix in G4GenericTrap::CalculateExtent
More detailed warning message for the case where "p is not on surface" in G4GenericTrap::SurfaceNormal(p)

Description whit 2015-12-21 20:20:41 CET
Hello,

The G4GenericTrap has bug in it which I can work around with the adjustment below. CalculateExtent seems to be returning values that are too small.
When I multiply by a factor of two things seem to work just fine.

The runtime error without this fix is:
-------- EEEE ------- G4Exception-START -------- EEEE -------
*** G4Exception : GeomMgt0002
      issued by : G4SmartVoxelHeader::BuildNodes()
PANIC! - Overlapping daughter with mother volume.
         Daughter physical volume collimator2_phys
         is entirely outside mother logical volume collimator_log !!
*** Fatal Exception *** core dump ***
-------- EEEE -------- G4Exception-END --------- EEEE -------

However, this volume is clearly inside of the volume, as visual inspection and  physical placement  overlap check have confirmed.

You will find an example project at https://github.com/whit2333/trap_bug

Cheers,
Whit

PS Also note the FindCLHEP does not work out of the box with the stated version of CLHEP for the latest version 10.2. Perhaps this should be in another bug report.


diff --git a/cmake/Modules/FindCLHEP.cmake b/cmake/Modules/FindCLHEP.cmake
index 8857f9c..ff9d113 100644
--- a/cmake/Modules/FindCLHEP.cmake
+++ b/cmake/Modules/FindCLHEP.cmake
@@ -205,7 +205,7 @@ find_path(CLHEP_INCLUDE_DIR CLHEP/Units/defs.h
 if(CLHEP_INCLUDE_DIR)
     set(CLHEP_VERSION 0)
     file(READ "${CLHEP_INCLUDE_DIR}/CLHEP/Units/defs.h" _CLHEP_DEFS_CONTENTS)
-    string(REGEX REPLACE ".*#define PACKAGE_VERSION \"([0-9.]+).*" "\\1"
+    string(REGEX REPLACE ".*#define CLHEP_UNITS_VERSION \"([0-9.]+).*" "\\1"
         CLHEP_VERSION "${_CLHEP_DEFS_CONTENTS}")
 
     if(NOT CLHEP_FIND_QUIETLY)
diff --git a/source/geometry/solids/specific/src/G4GenericTrap.cc b/source/geometry/solids/specific/src/G4GenericTrap.cc
index 6d0553d..82d4d4f 100644
--- a/source/geometry/solids/specific/src/G4GenericTrap.cc
+++ b/source/geometry/solids/specific/src/G4GenericTrap.cc
@@ -1322,6 +1322,12 @@ G4bool G4GenericTrap::CalculateExtent(const EAxis pAxis,
     }
     pMin-=kCarTolerance;
     pMax+=kCarTolerance;
+    // Artificailly multiply the extent by two because it is erroniously small?
+    //G4cout << " TIMES TWO " << G4endl;
+    pMin*=2.0;
+    pMax*=2.0;
+    //G4cout << " ### pMin " << pMin << G4endl;
+    //G4cout << " ### pMax " << pMax << G4endl;
 
     return true;
   }
@@ -1371,6 +1377,10 @@ G4bool G4GenericTrap::CalculateExtent(const EAxis pAxis,
       }
     }
     delete vertices;
+
+    G4cout << " TIMES TWO " << G4endl;
+    pMin*=2.0;
+    pMax*=2.0;
     return existsAfterClip;
   }
 }
Comment 1 Gabriele Cosmo 2015-12-21 23:00:08 CET
Thanks for the report on G4GenericTrap. We'll verify.
A ticket exists already for CMake and CLHEP; version 2.3.1.0 can be used with release 10.2 until a patch will be issued.
Comment 2 tatiana.nikitina 2016-01-17 21:54:36 CET
Thank you for reporting this problem.
Could you, please, give more details about the G4GenericTrap that you are using?

Thank you in advance,
Tatiana.
Comment 3 whit 2016-01-18 18:07:42 CET
(In reply to comment #2)
> Thank you for reporting this problem.
> Could you, please, give more details about the G4GenericTrap that you are
> using?
> 
> Thank you in advance,
> Tatiana.

I have provided a complete example above. You will fine the definition of the G4GenericTrap here: https://github.com/whit2333/trap_bug/blob/master/src/B1DetectorConstruction.cc#L294
Comment 4 tatiana.nikitina 2016-01-22 11:47:29 CET
Hello,

Thank you for the details. I could extract G4GenericTrap from your example and I'm testing it.

But I cannot run your example, it is not compiling on Linux. What machine do you use for run your test? 

Thank you in advance.
Comment 5 whit 2016-01-22 17:21:37 CET
(In reply to comment #4)
> Hello,
> 
> Thank you for the details. I could extract G4GenericTrap from your example and
> I'm testing it.
> 
> But I cannot run your example, it is not compiling on Linux. What machine do
> you use for run your test? 
> 
> Thank you in advance.

I am using Linux. What is your specific problem compiling?
Comment 6 tatiana.nikitina 2016-01-23 16:38:35 CET
Thank you.
The problem can be in differences in used C++.
The compiler is complying about wrong use of 'array' in  B1ParallelWorldConstruction.cc.

In any case the test with extracted G4GenericTrap shows an error in CalculateExtent(), CalculateExtent() is return centered values of extent, but in your case is not centered. I'm working on the fix.
Comment 7 whit 2016-01-25 22:15:26 CET
(In reply to comment #6)
> Thank you.
> The problem can be in differences in used C++.
> The compiler is complying about wrong use of 'array' in 
> B1ParallelWorldConstruction.cc.
> 
> In any case the test with extracted G4GenericTrap shows an error in
> CalculateExtent(), CalculateExtent() is return centered values of extent, but
> in your case is not centered. I'm working on the fix.

I also occasionally get the following warning when running my simulation:

-------- WWWW ------- G4Exception-START -------- WWWW -------
*** G4Exception : GeomSolids1002
      issued by : G4GenericTrap::SurfaceNormal(p)
Point p is not on surface !?
*** This is just a warning message. ***
-------- WWWW -------- G4Exception-END --------- WWWW -------

This might help you debug. It does not happen a lot.
Comment 8 Evgueni.Tcherniaev 2016-02-02 17:33:08 CET
Created attachment 380 [details]
Bug fix in G4GenericTrap::CalculateExtent
Comment 9 Evgueni.Tcherniaev 2016-02-02 17:37:08 CET
Hello,

Please find attached an bug fix in the G4GenericTrap::CalculateExtent. Let us know if you still have problem in this solid.

Best regards,
Evgueni
Comment 10 whit 2016-02-07 19:03:24 CET
(In reply to comment #9)
> Hello,
> 
> Please find attached an bug fix in the G4GenericTrap::CalculateExtent. Let us
> know if you still have problem in this solid.
> 
> Best regards,
> Evgueni

Thank you. The false overlaps have gone away, but I still get the following warning:

-------- WWWW ------- G4Exception-START -------- WWWW -------
*** G4Exception : GeomSolids1002
      issued by : G4GenericTrap::SurfaceNormal(p)
Point p is not on surface !?
*** This is just a warning message. ***
-------- WWWW -------- G4Exception-END --------- WWWW -------

This does not happen frequently but it is clearly a bug since the surface normals should be easily calculated because they are just planes.
Comment 11 Evgueni.Tcherniaev 2016-02-09 17:16:44 CET
Created attachment 382 [details]
More detailed warning message for the case where "p is not on surface" in G4GenericTrap::SurfaceNormal(p)
Comment 12 Evgueni.Tcherniaev 2016-02-09 17:49:45 CET
> Thank you. The false overlaps have gone away, but I still get the following
> warning:
> 
> -------- WWWW ------- G4Exception-START -------- WWWW -------
> *** G4Exception : GeomSolids1002
>       issued by : G4GenericTrap::SurfaceNormal(p)
> Point p is not on surface !?
> *** This is just a warning message. ***
> -------- WWWW -------- G4Exception-END --------- WWWW -------
> 
> This does not happen frequently but it is clearly a bug since the surface
> normals should be easily calculated because they are just planes.

I attach the code for G4GenericTrap.cc that produces more details output for
the case where SurfaceNormal() is called for a point situated outside the surface
of the G4GenericTrap. Please send us back the warning message that we would
be able to reproduce the problem.

By the way, the G4GenericTrap in your case is twisted. To create a solid with plane surfaces
some vertices should be shifted, for example as follows:

 double trap_width = 195.259*mm;
 std::vector< G4TwoVector> trap_points;
 trap_points.push_back(G4TwoVector(-103.734*mm, 179.672 *mm));
 trap_points.push_back(G4TwoVector(-1625.95*mm, 3172.94 *mm));
 trap_points.push_back(G4TwoVector(1625.95 *mm, 3172.94 *mm));
 trap_points.push_back(G4TwoVector(103.734 *mm, 179.672 *mm));
 trap_points.push_back(G4TwoVector(-103.734*mm, -3.29559*mm)); // was -2.42984*mm
 trap_points.push_back(G4TwoVector(-1683.74*mm, 3103.61 *mm));
 trap_points.push_back(G4TwoVector(1683.74 *mm, 3103.61 *mm));
 trap_points.push_back(G4TwoVector(103.734 *mm, -3.29559*mm)); // was -2.42984*mm
 
Let us know whether the problem also exists for such untwisted G4GenericTrap().

Thank you
Comment 13 Evgueni.Tcherniaev 2016-02-26 17:00:25 CET
The problem with incorrect calculation in G4GenericTrap::CalculateExtent() has been fixed. Corrected code will be available in the next Geant4 release.
Comment 14 whit 2016-02-26 17:55:41 CET
(In reply to comment #13)
> The problem with incorrect calculation in G4GenericTrap::CalculateExtent() has
> been fixed. Corrected code will be available in the next Geant4 release.

Did the surface normal problem also get fixed or should it be filed as another bug?
Comment 15 Evgueni.Tcherniaev 2016-02-26 18:52:25 CET
No, it was not fixed. We cannot reproduce the problem. In the second attachment there is a special G4GenericTrap.cc which should produce detailed output when the problem happens. Without this it is difficult to figure out what is the reason of the warning.

Thanks
Evgueni

(In reply to comment #14)
> (In reply to comment #13)
> > The problem with incorrect calculation in G4GenericTrap::CalculateExtent() has
> > been fixed. Corrected code will be available in the next Geant4 release.
> 
> Did the surface normal problem also get fixed or should it be filed as another
> bug?