Getting wrong result for volume of G4MultiUnion

Munther_Hindi · July 22, 2024, 1:39pm

Geant4 Version: geant4-v11.2.1
Operating System: Linux Fedora 40
Compiler/Version: 14.1.1
CMake Version: 3.28.2

The following gdml file listing (sorry I could not figure out how to upload the file) has a G4MultiUnion of 10 cones. Each cone as a base radius of 1 mm and height of 2 mm, so a volume of 2.094 mm^3. The geant application shows the G4MultiUnion volume of 15.875 mm^3 instead of the expected 10x2.094=20.94 mm^3. The application reports the correct number of nodes (10) and correct volume of each cone (2.094 mm^3). The Multiunion solid volume was reported wrong from both the application and the command /vis/ATree with a verbosity of 14.

Can any one confirm this please?

<?xml version='1.0' encoding='utf-8'?>
<gdml xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="http://service-spi.web.cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd">
  <define>
    <constant name="HALFPI" value="pi/2."/>
    <constant name="PI" value="1.*pi"/>
    <constant name="TWOPI" value="2.*pi"/>
    <position name="P-MultiUnion-PathArray_01" unit="mm" x="10.0" y="0.0" z="-6.0"/>
    <rotation name="identity" x="0" y="0" z="0"/>
    <position name="P-MultiUnion-PathArray_12" unit="mm" x="-1.74" y="9.84" z="-3.77"/>
    <position name="P-MultiUnion-PathArray_23" unit="mm" x="-9.39" y="-3.4" z="-1.55"/>
    <position name="P-MultiUnion-PathArray_34" unit="mm" x="5.00" y="-8.66" z="0.66"/>
    <position name="P-MultiUnion-PathArray_45" unit="mm" x="7.66" y="6.42" z="2.88"/>
    <position name="P-MultiUnion-PathArray_56" unit="mm" x="-7.66" y="6.42" z="5.11"/>
    <position name="P-MultiUnion-PathArray_67" unit="mm" x="-5.00" y="-8.66" z="7.33"/>
    <position name="P-MultiUnion-PathArray_78" unit="mm" x="9.39" y="-3.42" z="9.55"/>
    <position name="P-MultiUnion-PathArray_89" unit="mm" x="1.74" y="9.84" z="11.77"/>
    <position name="P-MultiUnion-PathArray_910" unit="mm" x="-10.0" y="0.0" z="14.00"/>
    <position name="center" x="0" y="0" z="0" unit="mm"/>
  </define>
  <materials/>
  <solids>
    <box name="WorldBox" x="33.80" y="33.80" z="57.2" lunit="mm"/>
    <cone name="Cone" rmax1="0.0" rmax2="1.0" deltaphi="360.0" aunit="deg" z="2.0" lunit="mm"/>
    <multiUnion name="MultiUnion-PathArray">
      <multiUnionNode name="MultiUnion-PathArray_0">
        <solid ref="Cone"/>
        <positionref ref="P-MultiUnion-PathArray_01"/>
        <rotationref ref="identity"/>
      </multiUnionNode>
      <multiUnionNode name="MultiUnion-PathArray_1">
        <solid ref="Cone"/>
        <positionref ref="P-MultiUnion-PathArray_12"/>
        <rotationref ref="identity"/>
      </multiUnionNode>
      <multiUnionNode name="MultiUnion-PathArray_2">
        <solid ref="Cone"/>
        <positionref ref="P-MultiUnion-PathArray_23"/>
        <rotationref ref="identity"/>
      </multiUnionNode>
      <multiUnionNode name="MultiUnion-PathArray_3">
        <solid ref="Cone"/>
        <positionref ref="P-MultiUnion-PathArray_34"/>
        <rotationref ref="identity"/>
      </multiUnionNode>
      <multiUnionNode name="MultiUnion-PathArray_4">
        <solid ref="Cone"/>
        <positionref ref="P-MultiUnion-PathArray_45"/>
        <rotationref ref="identity"/>
      </multiUnionNode>
      <multiUnionNode name="MultiUnion-PathArray_5">
        <solid ref="Cone"/>
        <positionref ref="P-MultiUnion-PathArray_56"/>
        <rotationref ref="identity"/>
      </multiUnionNode>
      <multiUnionNode name="MultiUnion-PathArray_6">
        <solid ref="Cone"/>
        <positionref ref="P-MultiUnion-PathArray_67"/>
        <rotationref ref="identity"/>
      </multiUnionNode>
      <multiUnionNode name="MultiUnion-PathArray_7">
        <solid ref="Cone"/>
        <positionref ref="P-MultiUnion-PathArray_78"/>
        <rotationref ref="identity"/>
      </multiUnionNode>
      <multiUnionNode name="MultiUnion-PathArray_8">
        <solid ref="Cone"/>
        <positionref ref="P-MultiUnion-PathArray_89"/>
        <rotationref ref="identity"/>
      </multiUnionNode>
      <multiUnionNode name="MultiUnion-PathArray_9">
        <solid ref="Cone"/>
        <positionref ref="P-MultiUnion-PathArray_910"/>
        <rotationref ref="identity"/>
      </multiUnionNode>
    </multiUnion>
  </solids>
  <structure>
    <volume name="Part">
      <materialref ref="G4_STAINLESS-STEEL"/>
      <solidref ref="MultiUnion-PathArray"/>
      <auxiliary auxtype="Color" auxvalue="#cccccc00"/>
    </volume>
    <volume name="worldVOL">
      <materialref ref="G4_AIR"/>
      <solidref ref="WorldBox"/>
      <physvol name="PV-Part">
        <volumeref ref="Part"/>
        <positionref ref="center"/>
        <rotationref ref="identity"/>
      </physvol>
    </volume>
  </structure>
  <setup name="Default" version="1.0">
    <world ref="worldVOL"/>
  </setup>
</gdml>

evc · July 22, 2024, 6:57pm

Volume of boolean solids is estimated using the Monte-Carlo method. In case of G4MultiUnion the number of random points used for the estimation is only 10000, that gives rather big estimation error.

In my case estimated volume is equal to 25.7972 mm^3. When I increased the number of random points to 1 million, the returned value became 20.8858 mm^3.

BTW, what was the reason to use G4MultiUnion for the cones which do not intersect each other?

Munther_Hindi · July 22, 2024, 8:55pm

Thanks Evgueni for the quick feedback. I suspected that might be the cause, but I did not expect the answer to be so far off.

The reason for using the G4MultiUnion is a little involved. The GDML file is produced by a GDML exporter for the FreeCAD open source CAD application. The exporter was written largley by Keith Sloane with a lot of help from me. FreeCAD allows the easy construction of arrays of objects and arrays of arrays of objects (as many depths as one wishes, but in practice one usually uses orthogonal arrays of polar arrays or vise versa). For better or worse, I decided to implement the exported array as G4MultiUnion, since it invloves only exporting one solid and one volume, with each node placement derived from the array parameters. For sure this is not the only way to implement array exports, but that’s what we have now. That is the reason for seeing MultiUnion of non-overlapping solids in this gdml.

For completeness, I should mention what led me to discover the mismatch: I am in the process of developing tools to debug the export of new and/or more complicated structures from the FreeCAD side to gdml. We would like to rely on more than just a visual comparison of the geomtery loaded into geant4 with that we see in FreeCAD. In the past I compared the volumes calculated in geant to those calculated in FreeCAD (which led me to discover and report a bug in the calculation of one of the volumes in geant4; this has been corrected in release 11+), but that is now not sufficent, since our volumes may be correct but their positioning and rotation in the exported gdml might be wrong and in large geometries might escape visual inspection. To that end I am calculating the center of mass and moments of intertia of the geometries in both FreeCAD and geant and comparing the two. That’s when I noticed that the array volume in FreeCAD was not the same as its exported G4MultiUnion in geant. By the way, for boolean solids I am also using a MonteCarlo calculation to calculate the center of mass in geant. For now I am using 10^6 simulations for each boolean, which seems to be fast enough.

Given your feedback, I think I will reply on a more accurate way to caclulate the volumes of booleans on the geant side. The geometry engine in FreeCAD uses openCascade which seems to give “exact” results for the volume and center of mass of booleans. I thought geant also uses openCascade, so am not sure why one has to rely on MC for the volume calculation.

Thanks gain for your feedback. Much appreciated.

Munther

Munther_Hindi · July 22, 2024, 9:37pm

I think I understand why the big difference in volumes: I assume the Monte Carlo generates the points in bounding volume of the G4MultiUnion and since in this case the bounding volume is much larger than the volume of the nodes, only a small fraction of the MC points land in in the nodes and hence the poor statistics and large difference.

evc · July 23, 2024, 7:24am

Yes, the accuracy of volume calculation depends not only on the number of points, but also on the size of the bounding box. In any case, 10000 points is too few to expect good accuracy and should be corrected. The default number of points for other boolean solids is 1 million.

I would like to note that GetCubicVolume() is not the only function that returns the volume of a solid. The volume of any solid can be estimated with G4VSolid::GetCubicVolume(), which uses 1 million points, and with G4VSolid::EstimateCubicVolume(npoints, epsilon), which provides a possibility to specify any number of points.

Munther_Hindi · July 23, 2024, 3:04pm

Thanks Evgueni. Good to know about EstimateCubicVolume. I probably can come up with an algorithm that gives accurate volumes for G4MultiUnions depending on how close or far apart the individual nodes are.

Thanks again.