GDMLSphere improvement

[lambdam94]

Hello Keith

The test of https://gitlab.cern.ch/VecGeom/VecGeom/-/blob/master/persistency/gdml/gdmls/oneSphere.gdml shows that there are pb with starttheta and deltatheta.

GDML example: oneSphere__l.gdml.txt

toKeith_GDMLSphereCreation.py.txt

Could you please add theses lines to GDMLSphere createGeometry to have the good picture? Thanks, Dam

def createGeometry(self,fp): import math

   # test on rmax
   if not hasattr(fp, 'rmax') :
       return
   if(fp.rmax <= 0.0):
       return

   mul = 1.0
   if hasattr(fp, 'lunit') : mul = GDMLShared.getMult(fp.lunit)

   # main sphere creation
   if hasattr(fp, 'deltaphi') : 
       sphere2 = Part.makeSphere(fp.rmax*mul, 
                           FreeCAD.Vector(0,0,0), 
                           FreeCAD.Vector(0,0,1), 
                           -90.0, 90.0,\
                           getAngleDeg(fp.aunit, fp.deltaphi))
   else: sphere2 = Part.makeSphere(fp.rmax*mul)

   # if starttheta -> cut the upper cone
   if hasattr(fp, 'starttheta') :
       angleTheta_deg = getAngleDeg(fp.aunit, fp.starttheta)
       if(angleTheta_deg>0.0):
           if (angleTheta_deg < 90.0):
               sphere2 = sphere2.cut( Part.makeCone(0.0,
                             fp.rmax*mul/math.cos(getAngleRad(fp.aunit, fp.starttheta)),
                             fp.rmax*mul))
           elif(angleTheta_deg == 90.0):
               Rmax = fp.rmax*mul*2.0 
               sphere2 = sphere2.cut( Part.makeBox(Rmax,Rmax,Rmax,FreeCAD.Vector(-0.5*Rmax,-0.5*Rmax,0.0)))
               del Rmax
           elif(angleTheta_deg <= 180.0):
               sphere2 = sphere2.common( Part.makeCone(0.0,
                             fp.rmax*mul/math.cos((180.0-angleTheta_deg)*math.pi/180.0),
                             fp.rmax*mul,
                             FreeCAD.Vector(0,0,0),
                             FreeCAD.Vector(0,0,-1.0)))
       del angleTheta_deg

   # if deltatheta -> cut the down cone
   if hasattr(fp, 'deltatheta') :
       angleDeltaTheta_deg = getAngleDeg(fp.aunit, fp.deltatheta)
       thetaSum_deg = angleDeltaTheta_deg
       if hasattr(fp, 'starttheta') :
           thetaSum_deg = angleDeltaTheta_deg + getAngleDeg(fp.aunit, fp.starttheta)

       if(thetaSum_deg <180.0):
           if (thetaSum_deg > 90.0):
               sphere2 = sphere2.cut( Part.makeCone(0.0,
                             fp.rmax*mul/math.cos((180.0-thetaSum_deg)*math.pi/180.0),
                             fp.rmax*mul,
                             FreeCAD.Vector(0,0,0),
                             FreeCAD.Vector(0,0,-1.0)))
           elif(thetaSum_deg == 90.0):
               Rmax = fp.rmax*mul*2.0 
               sphere2 = sphere2.cut( Part.makeBox(Rmax,Rmax,Rmax,FreeCAD.Vector(-0.5*Rmax,-0.5*Rmax,-Rmax)))
               del Rmax
           elif(thetaSum_deg >0.0):
               sphere2 = sphere2.common( Part.makeCone(0.0,
                             fp.rmax*mul/math.cos((180.0-angleTheta_deg)*math.pi/180.0),
                             fp.rmax*mul,
                             FreeCAD.Vector(0,0,0),
                             FreeCAD.Vector(0,0,1.0)))
       del thetaSum_deg, angleDeltaTheta_deg

   # if rmin -> cut the rmin sphere
   if hasattr(fp, 'rmin') :
       if(fp.rmin <= 0.0 or fp.rmin>fp.rmax):
           fp.Shape = sphere2 
       else:         
           fp.Shape = sphere2.cut(Part.makeSphere(fp.rmin*mul))

   # if startphi -> rotate startphi angle
   if hasattr(fp, 'startphi') :
       newplace = FreeCAD.Placement(fp.Placement.Base,
                       FreeCAD.Rotation(FreeCAD.Vector(0,0,1), getAngleDeg(fp.aunit, fp.startphi)),
                       FreeCAD.Vector(0,0,0))        
       fp.Placement = newplace                            

   del mul 

[KeithSloan]

Thanks - Have done a variation of the suggested code - If you can buddy check that would be great.

[lambdam94]

Hi Keith

under FreeCAD0.19.AppImage, the import file give this picture (theta is not good)

Under FreeCAD0.18 (official ubuntu distribution) on my old laptop, there is the "angleSectionSolid" pb under Freecad0.18. With the suggested function, we have the "good" picture.

I will take a look on GDMLPolyhedra (perhaps I can give to you a new suggestion how to do without angleSectionSolid !)

Regards, Dam

[KeithSloan]

Can we check the file you are using.

For me oneSphere is a solid sphere

sphere name="world" rmin="0" rmax="100" startphi="0" deltaphi="360" starttheta="0" deltatheta="180" aunit="deg" lunit="cm"

And I just did a pull on the repro to double check.

[lambdam94]

The picture has been done with the GDML example "oneSphere__l.gdml.txt" (file at the beginning of this issue).

The parameters are:

The theta angle ("second solid" angle) is between 25 deg and 115 deg.

Regards Dam

[lambdam94]

In freecad makeSphere, the angle parameters cut the sphere by a plan cf. https://wiki.freecadweb.org/Part_Sphere It is not easy to do a "solid" angle sphere (gdml definition) with this function. A cut or a common with a sphere can do that "solide" angle.

[KeithSloan]

Okay - these files are not valid as far as the Geant4 version of gdmlview. Have had to add materials before it works, but now I get the same visual as you. Thanks for the heads up on Part_Sphere. My tests were with units of 90 so did not show up my lack of understanding.

[lambdam94]

Hi Keith our messages probably crossed! I propose you a mix of our codes which works on FreeCad0.18

Regards Dam

` #new code for GDMLSphere

def createGeometry(self,fp): import math

   # test on rmax
   if not hasattr(fp, 'rmax') :
       return
   if(fp.rmax <= 0.0):
       return

   mul = 1.0
   if hasattr(fp, 'lunit') : mul = GDMLShared.getMult(fp.lunit)

   spos = FreeCAD.Vector(0,0,0)
   sdir = FreeCAD.Vector(0,0,1)

   # main sphere creation
   if hasattr(fp, 'deltaphi') : 
       sphere2 = Part.makeSphere(fp.rmax*mul, spos, sdir,
                       -90.0, 90.0,\
                       getAngleDeg(fp.aunit, fp.deltaphi))

       if hasattr(fp, 'startphi') :
           (sphere2).rotate(spos, sdir,getAngleDeg(fp.aunit, fp.startphi))

   else: sphere2 = Part.makeSphere(fp.rmax*mul)

   # if starttheta -> cut the upper cone
   if hasattr(fp, 'starttheta') :
       angleTheta_deg = getAngleDeg(fp.aunit, fp.starttheta)
       if(angleTheta_deg>0.0):
           if (angleTheta_deg < 90.0):
               sphere2 = sphere2.cut( Part.makeCone(0.0,
                         fp.rmax*mul/math.cos(getAngleRad(fp.aunit, fp.starttheta)),
                         fp.rmax*mul))
           elif(angleTheta_deg == 90.0):
               Rmax = fp.rmax*mul*2.0 
               sphere2 = sphere2.cut( Part.makeBox(Rmax,Rmax,Rmax,FreeCAD.Vector(-0.5*Rmax,-0.5*Rmax,0.0)))
               del Rmax
           elif(angleTheta_deg <= 180.0):
               sphere2 = sphere2.common( Part.makeCone(0.0,
                         fp.rmax*mul/math.cos((180.0-angleTheta_deg)*math.pi/180.0),
                         fp.rmax*mul,
                         spos,
                         FreeCAD.Vector(0,0,-1.0)))
       del angleTheta_deg

   # if deltatheta -> cut the down cone
   if hasattr(fp, 'deltatheta') :
       angleDeltaTheta_deg = getAngleDeg(fp.aunit, fp.deltatheta)
       thetaSum_deg = angleDeltaTheta_deg
       if hasattr(fp, 'starttheta') :
           thetaSum_deg = angleDeltaTheta_deg + getAngleDeg(fp.aunit, fp.starttheta)

       if(thetaSum_deg <180.0):
           if (thetaSum_deg > 90.0):
               sphere2 = sphere2.cut( Part.makeCone(0.0,
                         fp.rmax*mul/math.cos((180.0-thetaSum_deg)*math.pi/180.0),
                         fp.rmax*mul,
                         spos,
                         FreeCAD.Vector(0,0,-1.0)))
           elif(thetaSum_deg == 90.0):
               Rmax = fp.rmax*mul*2.0 
               sphere2 = sphere2.cut( Part.makeBox(Rmax,Rmax,Rmax,FreeCAD.Vector(-0.5*Rmax,-0.5*Rmax,-Rmax)))
               del Rmax
       elif(thetaSum_deg >0.0):
           sphere2 = sphere2.common( Part.makeCone(0.0,
                         fp.rmax*mul/math.cos((180.0-angleTheta_deg)*math.pi/180.0),
                         fp.rmax*mul,
                         spos, sdir))

       del thetaSum_deg, angleDeltaTheta_deg

   # if rmin -> cut the rmin sphere
   if hasattr(fp, 'rmin') :
       if(fp.rmin <= 0.0 or fp.rmin>fp.rmax):
           fp.Shape = sphere2 
       else:         
           fp.Shape = sphere2.cut(Part.makeSphere(fp.rmin*mul))
   else: fp.Shape = sphere2 

   # delete old var -> free memory
   del mul, spos, sdir

`

[KeithSloan]

Why are you using expressions like math.cos((180.0-angleTheta_deg)*math.pi/180.0)? Python math wants radians so why not use angleDeltaTheta_rad = getAngleRad(fp.aunit, fp.deltatheta) so math.cos(pi-angleDeltaTheta_rad)

[lambdam94]

I am agree with you. We can call also getAngleRad. But I think that we should have only one angle definition chosen in rad or chosen in deg. You are right, rad is more useful for cos.

Thus " if(thetaSum_deg <180.0): " should replace by " if(thetaSum_rad <2.0*math.pi): " etc.

[KeithSloan]

Please Try latest commit

[lambdam94]

That seems good! Thanks. For starttheta=90deg, there is one pb ("1.0/0.0"). For deltatheta<90deg, thetadeltaRad should be replace by deltathetaRad.

[KeithSloan]

Not sure I am clear where these corrections are. Please could you give line numbers. Thanks

[lambdam94]

hello Keith, There was probably a pb during the commit process because there are two § about starttheta. I add some comments on the "Changes to Tube & Sphere" version. Regards Dam