Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/CSG/src/G4Para.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /geometry/solids/CSG/src/G4Para.cc (Version 11.3.0) and /geometry/solids/CSG/src/G4Para.cc (Version 10.5)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 //
                                                   >>  27 //
                                                   >>  28 // class G4Para
                                                   >>  29 //
 26 // Implementation for G4Para class                 30 // Implementation for G4Para class
 27 //                                                 31 //
 28 // 21.03.95 P.Kent: Modified for `tolerant' ge <<  32 // History:
 29 // 31.10.96 V.Grichine: Modifications accordin <<  33 //
 30 // 28.04.05 V.Grichine: new SurfaceNormal acco << 
 31 // 29.05.17 E.Tcherniaev: complete revision, s     34 // 29.05.17 E.Tcherniaev: complete revision, speed-up
                                                   >>  35 // 23.09.16 E.Tcherniaev: use G4BoundingEnvelope for CalculateExtent(),
                                                   >>  36 //                      removed CreateRotatedVertices()
                                                   >>  37 // 23.10.05 V.Grichine: bug fixed in DistanceToOut(p,v,...) for the v.x()<0 case
                                                   >>  38 // 28.04.05 V.Grichine: new SurfaceNormal according to J. Apostolakis proposal
                                                   >>  39 // 30.11.04 V.Grichine: modifications in SurfaceNormal for edges/vertices and
                                                   >>  40 //                      in constructor with vertices
                                                   >>  41 // 14.02.02 V.Grichine: bug fixed in Inside according to proposal of D.Wright
                                                   >>  42 // 18.11.99 V.Grichine: kUndef was added to ESide
                                                   >>  43 // 31.10.96 V.Grichine: Modifications according G4Box/Tubs before to commit
                                                   >>  44 // 21.03.95 P.Kent: Modified for `tolerant' geom
                                                   >>  45 //
 32 //////////////////////////////////////////////     46 ////////////////////////////////////////////////////////////////////////////
 33                                                    47 
 34 #include "G4Para.hh"                               48 #include "G4Para.hh"
 35                                                    49 
 36 #if !defined(G4GEOM_USE_UPARA)                     50 #if !defined(G4GEOM_USE_UPARA)
 37                                                    51 
 38 #include "G4VoxelLimits.hh"                        52 #include "G4VoxelLimits.hh"
 39 #include "G4AffineTransform.hh"                    53 #include "G4AffineTransform.hh"
 40 #include "G4BoundingEnvelope.hh"                   54 #include "G4BoundingEnvelope.hh"
 41 #include "Randomize.hh"                            55 #include "Randomize.hh"
 42                                                    56 
 43 #include "G4VPVParameterisation.hh"                57 #include "G4VPVParameterisation.hh"
 44                                                    58 
 45 #include "G4VGraphicsScene.hh"                     59 #include "G4VGraphicsScene.hh"
 46                                                    60 
 47 using namespace CLHEP;                             61 using namespace CLHEP;
 48                                                    62 
 49 //////////////////////////////////////////////     63 //////////////////////////////////////////////////////////////////////////
 50 //                                                 64 //
 51 //  Constructor - set & check half widths          65 //  Constructor - set & check half widths
 52                                                    66 
 53 G4Para::G4Para(const G4String& pName,              67 G4Para::G4Para(const G4String& pName,
 54                      G4double pDx, G4double pD     68                      G4double pDx, G4double pDy, G4double pDz,
 55                      G4double pAlpha, G4double     69                      G4double pAlpha, G4double pTheta, G4double pPhi)
 56   : G4CSGSolid(pName), halfCarTolerance(0.5*kC     70   : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance)
 57 {                                                  71 {
 58   SetAllParameters(pDx, pDy, pDz, pAlpha, pThe     72   SetAllParameters(pDx, pDy, pDz, pAlpha, pTheta, pPhi);
 59   fRebuildPolyhedron = false;  // default valu     73   fRebuildPolyhedron = false;  // default value for G4CSGSolid
 60 }                                                  74 }
 61                                                    75 
 62 //////////////////////////////////////////////     76 //////////////////////////////////////////////////////////////////////////
 63 //                                                 77 //
 64 // Constructor - design of trapezoid based on      78 // Constructor - design of trapezoid based on 8 vertices
 65                                                    79 
 66 G4Para::G4Para( const G4String& pName,             80 G4Para::G4Para( const G4String& pName,
 67                 const G4ThreeVector pt[8] )        81                 const G4ThreeVector pt[8] )
 68   : G4CSGSolid(pName), halfCarTolerance(0.5*kC     82   : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance)
 69 {                                                  83 {
 70   // Find dimensions and trigonometric values      84   // Find dimensions and trigonometric values
 71   //                                               85   // 
 72   fDx = (pt[3].x() - pt[2].x())*0.5;               86   fDx = (pt[3].x() - pt[2].x())*0.5;
 73   fDy = (pt[2].y() - pt[1].y())*0.5;               87   fDy = (pt[2].y() - pt[1].y())*0.5;
 74   fDz = pt[7].z();                                 88   fDz = pt[7].z();
 75   CheckParameters(); // check dimensions           89   CheckParameters(); // check dimensions
 76                                                    90 
 77   fTalpha = (pt[2].x() + pt[3].x() - pt[1].x()     91   fTalpha = (pt[2].x() + pt[3].x() - pt[1].x() - pt[0].x())*0.25/fDy;
 78   fTthetaCphi = (pt[4].x() + fDy*fTalpha + fDx     92   fTthetaCphi = (pt[4].x() + fDy*fTalpha + fDx)/fDz;
 79   fTthetaSphi = (pt[4].y() + fDy)/fDz;             93   fTthetaSphi = (pt[4].y() + fDy)/fDz;
 80   MakePlanes();                                    94   MakePlanes();
 81                                                    95 
 82   // Recompute vertices                            96   // Recompute vertices
 83   //                                               97   //
 84   G4ThreeVector v[8];                              98   G4ThreeVector v[8];
 85   G4double DyTalpha = fDy*fTalpha;                 99   G4double DyTalpha = fDy*fTalpha;
 86   G4double DzTthetaSphi = fDz*fTthetaSphi;        100   G4double DzTthetaSphi = fDz*fTthetaSphi;
 87   G4double DzTthetaCphi = fDz*fTthetaCphi;        101   G4double DzTthetaCphi = fDz*fTthetaCphi;
 88   v[0].set(-DzTthetaCphi-DyTalpha-fDx, -DzTthe    102   v[0].set(-DzTthetaCphi-DyTalpha-fDx, -DzTthetaSphi-fDy, -fDz);
 89   v[1].set(-DzTthetaCphi-DyTalpha+fDx, -DzTthe    103   v[1].set(-DzTthetaCphi-DyTalpha+fDx, -DzTthetaSphi-fDy, -fDz);
 90   v[2].set(-DzTthetaCphi+DyTalpha-fDx, -DzTthe    104   v[2].set(-DzTthetaCphi+DyTalpha-fDx, -DzTthetaSphi+fDy, -fDz);
 91   v[3].set(-DzTthetaCphi+DyTalpha+fDx, -DzTthe    105   v[3].set(-DzTthetaCphi+DyTalpha+fDx, -DzTthetaSphi+fDy, -fDz);
 92   v[4].set( DzTthetaCphi-DyTalpha-fDx,  DzTthe    106   v[4].set( DzTthetaCphi-DyTalpha-fDx,  DzTthetaSphi-fDy,  fDz);
 93   v[5].set( DzTthetaCphi-DyTalpha+fDx,  DzTthe    107   v[5].set( DzTthetaCphi-DyTalpha+fDx,  DzTthetaSphi-fDy,  fDz);
 94   v[6].set( DzTthetaCphi+DyTalpha-fDx,  DzTthe    108   v[6].set( DzTthetaCphi+DyTalpha-fDx,  DzTthetaSphi+fDy,  fDz);
 95   v[7].set( DzTthetaCphi+DyTalpha+fDx,  DzTthe    109   v[7].set( DzTthetaCphi+DyTalpha+fDx,  DzTthetaSphi+fDy,  fDz);
 96                                                   110 
 97   // Compare with original vertices               111   // Compare with original vertices
 98   //                                              112   //
 99   for (G4int i=0; i<8; ++i)                       113   for (G4int i=0; i<8; ++i)
100   {                                               114   {
101     G4double delx = std::abs(pt[i].x() - v[i].    115     G4double delx = std::abs(pt[i].x() - v[i].x());
102     G4double dely = std::abs(pt[i].y() - v[i].    116     G4double dely = std::abs(pt[i].y() - v[i].y());
103     G4double delz = std::abs(pt[i].z() - v[i].    117     G4double delz = std::abs(pt[i].z() - v[i].z());
104     G4double discrepancy = std::max(std::max(d    118     G4double discrepancy = std::max(std::max(delx,dely),delz);
105     if (discrepancy > 0.1*kCarTolerance)          119     if (discrepancy > 0.1*kCarTolerance)
106     {                                             120     {
107       std::ostringstream message;                 121       std::ostringstream message;
108       G4long oldprc = message.precision(16);   << 122       G4int oldprc = message.precision(16);
109       message << "Invalid vertice coordinates     123       message << "Invalid vertice coordinates for Solid: " << GetName()
110               << "\nVertix #" << i << ", discr    124               << "\nVertix #" << i << ", discrepancy = " << discrepancy
111               << "\n  original   : " << pt[i]     125               << "\n  original   : " << pt[i]
112               << "\n  recomputed : " << v[i];     126               << "\n  recomputed : " << v[i];
113       G4cout.precision(oldprc);                   127       G4cout.precision(oldprc);
114       G4Exception("G4Para::G4Para()", "GeomSol    128       G4Exception("G4Para::G4Para()", "GeomSolids0002",
115                   FatalException, message);       129                   FatalException, message);
116                                                   130 
117     }                                             131     }
118   }                                               132   }
119 }                                                 133 }
120                                                   134 
121 //////////////////////////////////////////////    135 //////////////////////////////////////////////////////////////////////////
122 //                                                136 //
123 // Fake default constructor - sets only member    137 // Fake default constructor - sets only member data and allocates memory
124 //                            for usage restri    138 //                            for usage restricted to object persistency
125                                                   139 
126 G4Para::G4Para( __void__& a )                     140 G4Para::G4Para( __void__& a )
127   : G4CSGSolid(a), halfCarTolerance(0.5*kCarTo    141   : G4CSGSolid(a), halfCarTolerance(0.5*kCarTolerance)
128 {                                                 142 {
129   SetAllParameters(1., 1., 1., 0., 0., 0.);       143   SetAllParameters(1., 1., 1., 0., 0., 0.);
130   fRebuildPolyhedron = false; // default value    144   fRebuildPolyhedron = false; // default value for G4CSGSolid
131 }                                                 145 }
132                                                   146 
133 //////////////////////////////////////////////    147 //////////////////////////////////////////////////////////////////////////
134 //                                                148 //
135 // Destructor                                     149 // Destructor
136                                                   150 
137 G4Para::~G4Para() = default;                   << 151 G4Para::~G4Para()
                                                   >> 152 {
                                                   >> 153 }
138                                                   154 
139 //////////////////////////////////////////////    155 //////////////////////////////////////////////////////////////////////////
140 //                                                156 //
141 // Copy constructor                               157 // Copy constructor
142                                                   158 
143 G4Para::G4Para(const G4Para& rhs)                 159 G4Para::G4Para(const G4Para& rhs)
144   : G4CSGSolid(rhs), halfCarTolerance(rhs.half    160   : G4CSGSolid(rhs), halfCarTolerance(rhs.halfCarTolerance),
145     fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz),     161     fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz), fTalpha(rhs.fTalpha),
146     fTthetaCphi(rhs.fTthetaCphi),fTthetaSphi(r    162     fTthetaCphi(rhs.fTthetaCphi),fTthetaSphi(rhs.fTthetaSphi)
147 {                                                 163 {
148   for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs    164   for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; }
149 }                                                 165 }
150                                                   166 
151 //////////////////////////////////////////////    167 //////////////////////////////////////////////////////////////////////////
152 //                                                168 //
153 // Assignment operator                            169 // Assignment operator
154                                                   170 
155 G4Para& G4Para::operator = (const G4Para& rhs)    171 G4Para& G4Para::operator = (const G4Para& rhs)
156 {                                                 172 {
157    // Check assignment to self                    173    // Check assignment to self
158    //                                             174    //
159    if (this == &rhs)  { return *this; }           175    if (this == &rhs)  { return *this; }
160                                                   176 
161    // Copy base class data                        177    // Copy base class data
162    //                                             178    //
163    G4CSGSolid::operator=(rhs);                    179    G4CSGSolid::operator=(rhs);
164                                                   180 
165    // Copy data                                   181    // Copy data
166    //                                             182    //
167    halfCarTolerance = rhs.halfCarTolerance;       183    halfCarTolerance = rhs.halfCarTolerance;
168    fDx = rhs.fDx;                                 184    fDx = rhs.fDx;
169    fDy = rhs.fDy;                                 185    fDy = rhs.fDy;
170    fDz = rhs.fDz;                                 186    fDz = rhs.fDz;
171    fTalpha = rhs.fTalpha;                         187    fTalpha = rhs.fTalpha;
172    fTthetaCphi = rhs.fTthetaCphi;                 188    fTthetaCphi = rhs.fTthetaCphi;
173    fTthetaSphi = rhs.fTthetaSphi;                 189    fTthetaSphi = rhs.fTthetaSphi;
174    for (G4int i=0; i<4; ++i) { fPlanes[i] = rh    190    for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; }
175                                                   191 
176    return *this;                                  192    return *this;
177 }                                                 193 }
178                                                   194 
179 //////////////////////////////////////////////    195 //////////////////////////////////////////////////////////////////////////
180 //                                                196 //
181 // Set all parameters, as for constructor - se    197 // Set all parameters, as for constructor - set and check half-widths
182                                                   198 
183 void G4Para::SetAllParameters(G4double pDx, G4    199 void G4Para::SetAllParameters(G4double pDx, G4double pDy, G4double pDz,
184                               G4double pAlpha,    200                               G4double pAlpha, G4double pTheta, G4double pPhi)
185 {                                                 201 {
186   // Reset data of the base class                 202   // Reset data of the base class
187   fCubicVolume = 0;                               203   fCubicVolume = 0;
188   fSurfaceArea = 0;                               204   fSurfaceArea = 0;
189   fRebuildPolyhedron = true;                      205   fRebuildPolyhedron = true;
190                                                   206 
191   // Set parameters                               207   // Set parameters
192   fDx = pDx;                                      208   fDx = pDx;
193   fDy = pDy;                                      209   fDy = pDy;
194   fDz = pDz;                                      210   fDz = pDz;
195   fTalpha = std::tan(pAlpha);                     211   fTalpha = std::tan(pAlpha);
196   fTthetaCphi = std::tan(pTheta)*std::cos(pPhi    212   fTthetaCphi = std::tan(pTheta)*std::cos(pPhi);
197   fTthetaSphi = std::tan(pTheta)*std::sin(pPhi    213   fTthetaSphi = std::tan(pTheta)*std::sin(pPhi);
198                                                   214 
199   CheckParameters();                              215   CheckParameters();
200   MakePlanes();                                   216   MakePlanes();
201 }                                                 217 }
202                                                   218 
203 //////////////////////////////////////////////    219 //////////////////////////////////////////////////////////////////////////
204 //                                                220 //
205 // Check dimensions                               221 // Check dimensions
206                                                   222 
207 void G4Para::CheckParameters()                    223 void G4Para::CheckParameters()
208 {                                                 224 {
209   if (fDx < 2*kCarTolerance ||                    225   if (fDx < 2*kCarTolerance ||
210       fDy < 2*kCarTolerance ||                    226       fDy < 2*kCarTolerance ||
211       fDz < 2*kCarTolerance)                      227       fDz < 2*kCarTolerance)
212   {                                               228   {
213     std::ostringstream message;                   229     std::ostringstream message;
214     message << "Invalid (too small or negative    230     message << "Invalid (too small or negative) dimensions for Solid: "
215             << GetName()                          231             << GetName()
216             << "\n  X - " << fDx                  232             << "\n  X - " << fDx
217             << "\n  Y - " << fDy                  233             << "\n  Y - " << fDy
218             << "\n  Z - " << fDz;                 234             << "\n  Z - " << fDz;
219     G4Exception("G4Para::CheckParameters()", "    235     G4Exception("G4Para::CheckParameters()", "GeomSolids0002",
220                 FatalException, message);         236                 FatalException, message);
221   }                                               237   }
222 }                                                 238 }
223                                                   239 
224 //////////////////////////////////////////////    240 //////////////////////////////////////////////////////////////////////////
225 //                                                241 //
226 // Set side planes                                242 // Set side planes
227                                                   243 
228 void G4Para::MakePlanes()                         244 void G4Para::MakePlanes()
229 {                                                 245 {
230   G4ThreeVector vx(1, 0, 0);                      246   G4ThreeVector vx(1, 0, 0);
231   G4ThreeVector vy(fTalpha, 1, 0);                247   G4ThreeVector vy(fTalpha, 1, 0);
232   G4ThreeVector vz(fTthetaCphi, fTthetaSphi, 1    248   G4ThreeVector vz(fTthetaCphi, fTthetaSphi, 1);
233                                                   249 
234   // Set -Y & +Y planes                           250   // Set -Y & +Y planes
235   //                                              251   //
236   G4ThreeVector ynorm = (vx.cross(vz)).unit();    252   G4ThreeVector ynorm = (vx.cross(vz)).unit();
237                                                   253 
238   fPlanes[0].a = 0.;                              254   fPlanes[0].a = 0.;
239   fPlanes[0].b = ynorm.y();                       255   fPlanes[0].b = ynorm.y();
240   fPlanes[0].c = ynorm.z();                       256   fPlanes[0].c = ynorm.z();
241   fPlanes[0].d = fPlanes[0].b*fDy; // point (0    257   fPlanes[0].d = fPlanes[0].b*fDy; // point (0,fDy,0) is on plane
242                                                   258 
243   fPlanes[1].a =  0.;                             259   fPlanes[1].a =  0.;
244   fPlanes[1].b = -fPlanes[0].b;                   260   fPlanes[1].b = -fPlanes[0].b;
245   fPlanes[1].c = -fPlanes[0].c;                   261   fPlanes[1].c = -fPlanes[0].c;
246   fPlanes[1].d =  fPlanes[0].d;                   262   fPlanes[1].d =  fPlanes[0].d;
247                                                   263 
248   // Set -X & +X planes                           264   // Set -X & +X planes
249   //                                              265   //
250   G4ThreeVector xnorm = (vz.cross(vy)).unit();    266   G4ThreeVector xnorm = (vz.cross(vy)).unit();
251                                                   267 
252   fPlanes[2].a = xnorm.x();                       268   fPlanes[2].a = xnorm.x();
253   fPlanes[2].b = xnorm.y();                       269   fPlanes[2].b = xnorm.y();
254   fPlanes[2].c = xnorm.z();                       270   fPlanes[2].c = xnorm.z();
255   fPlanes[2].d = fPlanes[2].a*fDx; // point (f    271   fPlanes[2].d = fPlanes[2].a*fDx; // point (fDx,0,0) is on plane
256                                                   272 
257   fPlanes[3].a = -fPlanes[2].a;                   273   fPlanes[3].a = -fPlanes[2].a;
258   fPlanes[3].b = -fPlanes[2].b;                   274   fPlanes[3].b = -fPlanes[2].b;
259   fPlanes[3].c = -fPlanes[2].c;                   275   fPlanes[3].c = -fPlanes[2].c;
260   fPlanes[3].d =  fPlanes[2].d;                   276   fPlanes[3].d =  fPlanes[2].d;
261 }                                                 277 }
262                                                   278 
263 //////////////////////////////////////////////    279 //////////////////////////////////////////////////////////////////////////
264 //                                                280 //
265 // Get volume                                     281 // Get volume
266                                                   282 
267 G4double G4Para::GetCubicVolume()                 283 G4double G4Para::GetCubicVolume()
268 {                                                 284 {
269   // It is like G4Box, since para transformati    285   // It is like G4Box, since para transformations keep the volume to be const
270   if (fCubicVolume == 0)                          286   if (fCubicVolume == 0)
271   {                                               287   {
272     fCubicVolume = 8*fDx*fDy*fDz;                 288     fCubicVolume = 8*fDx*fDy*fDz;
273   }                                               289   }
274   return fCubicVolume;                            290   return fCubicVolume;
275 }                                                 291 }
276                                                   292 
277 //////////////////////////////////////////////    293 //////////////////////////////////////////////////////////////////////////
278 //                                                294 //
279 // Get surface area                               295 // Get surface area
280                                                   296 
281 G4double G4Para::GetSurfaceArea()                 297 G4double G4Para::GetSurfaceArea()
282 {                                                 298 {
283   if(fSurfaceArea == 0)                           299   if(fSurfaceArea == 0)
284   {                                               300   {
285     G4ThreeVector vx(fDx, 0, 0);                  301     G4ThreeVector vx(fDx, 0, 0);
286     G4ThreeVector vy(fDy*fTalpha, fDy, 0);        302     G4ThreeVector vy(fDy*fTalpha, fDy, 0);
287     G4ThreeVector vz(fDz*fTthetaCphi, fDz*fTth    303     G4ThreeVector vz(fDz*fTthetaCphi, fDz*fTthetaSphi, fDz);
288                                                   304 
289     G4double sxy = fDx*fDy; // (vx.cross(vy)).    305     G4double sxy = fDx*fDy; // (vx.cross(vy)).mag();
290     G4double sxz = (vx.cross(vz)).mag();          306     G4double sxz = (vx.cross(vz)).mag();
291     G4double syz = (vy.cross(vz)).mag();          307     G4double syz = (vy.cross(vz)).mag();
292                                                   308 
293     fSurfaceArea = 8*(sxy+sxz+syz);               309     fSurfaceArea = 8*(sxy+sxz+syz);
294   }                                               310   }
295   return fSurfaceArea;                            311   return fSurfaceArea;
296 }                                                 312 }
297                                                   313 
298 //////////////////////////////////////////////    314 //////////////////////////////////////////////////////////////////////////
299 //                                                315 //
300 // Dispatch to parameterisation for replicatio    316 // Dispatch to parameterisation for replication mechanism dimension
301 // computation & modification                     317 // computation & modification
302                                                   318 
303 void G4Para::ComputeDimensions(      G4VPVPara    319 void G4Para::ComputeDimensions(      G4VPVParameterisation* p,
304                                 const G4int n,    320                                 const G4int n,
305                                 const G4VPhysi    321                                 const G4VPhysicalVolume* pRep )
306 {                                                 322 {
307   p->ComputeDimensions(*this,n,pRep);             323   p->ComputeDimensions(*this,n,pRep);
308 }                                                 324 }
309                                                   325 
310 //////////////////////////////////////////////    326 //////////////////////////////////////////////////////////////////////////
311 //                                                327 //
312 // Get bounding box                               328 // Get bounding box
313                                                   329 
314 void G4Para::BoundingLimits(G4ThreeVector& pMi    330 void G4Para::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
315 {                                                 331 {
316   G4double dz = GetZHalfLength();                 332   G4double dz = GetZHalfLength();
317   G4double dx = GetXHalfLength();                 333   G4double dx = GetXHalfLength();
318   G4double dy = GetYHalfLength();                 334   G4double dy = GetYHalfLength();
319                                                   335 
320   G4double x0 = dz*fTthetaCphi;                   336   G4double x0 = dz*fTthetaCphi;
321   G4double x1 = dy*GetTanAlpha();                 337   G4double x1 = dy*GetTanAlpha();
322   G4double xmin =                                 338   G4double xmin =
323     std::min(                                     339     std::min(
324     std::min(                                     340     std::min(
325     std::min(-x0-x1-dx,-x0+x1-dx),x0-x1-dx),x0    341     std::min(-x0-x1-dx,-x0+x1-dx),x0-x1-dx),x0+x1-dx);
326   G4double xmax =                                 342   G4double xmax =
327     std::max(                                     343     std::max(
328     std::max(                                     344     std::max(
329     std::max(-x0-x1+dx,-x0+x1+dx),x0-x1+dx),x0    345     std::max(-x0-x1+dx,-x0+x1+dx),x0-x1+dx),x0+x1+dx);
330                                                   346 
331   G4double y0 = dz*fTthetaSphi;                   347   G4double y0 = dz*fTthetaSphi;
332   G4double ymin = std::min(-y0-dy,y0-dy);         348   G4double ymin = std::min(-y0-dy,y0-dy);
333   G4double ymax = std::max(-y0+dy,y0+dy);         349   G4double ymax = std::max(-y0+dy,y0+dy);
334                                                   350 
335   pMin.set(xmin,ymin,-dz);                        351   pMin.set(xmin,ymin,-dz);
336   pMax.set(xmax,ymax, dz);                        352   pMax.set(xmax,ymax, dz);
337                                                   353 
338   // Check correctness of the bounding box        354   // Check correctness of the bounding box
339   //                                              355   //
340   if (pMin.x() >= pMax.x() || pMin.y() >= pMax    356   if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
341   {                                               357   {
342     std::ostringstream message;                   358     std::ostringstream message;
343     message << "Bad bounding box (min >= max)     359     message << "Bad bounding box (min >= max) for solid: "
344             << GetName() << " !"                  360             << GetName() << " !"
345             << "\npMin = " << pMin                361             << "\npMin = " << pMin
346             << "\npMax = " << pMax;               362             << "\npMax = " << pMax;
347     G4Exception("G4Para::BoundingLimits()", "G    363     G4Exception("G4Para::BoundingLimits()", "GeomMgt0001",
348                 JustWarning, message);            364                 JustWarning, message);
349     DumpInfo();                                   365     DumpInfo();
350   }                                               366   }
351 }                                                 367 }
352                                                   368 
353 //////////////////////////////////////////////    369 //////////////////////////////////////////////////////////////////////////
354 //                                                370 //
355 // Calculate extent under transform and specif    371 // Calculate extent under transform and specified limit
356                                                   372 
357 G4bool G4Para::CalculateExtent( const EAxis pA    373 G4bool G4Para::CalculateExtent( const EAxis pAxis,
358                                 const G4VoxelL    374                                 const G4VoxelLimits& pVoxelLimit,
359                                 const G4Affine    375                                 const G4AffineTransform& pTransform,
360                                      G4double&    376                                      G4double& pMin, G4double& pMax ) const
361 {                                                 377 {
362   G4ThreeVector bmin, bmax;                       378   G4ThreeVector bmin, bmax;
363   G4bool exist;                                   379   G4bool exist;
364                                                   380 
365   // Check bounding box (bbox)                    381   // Check bounding box (bbox)
366   //                                              382   //
367   BoundingLimits(bmin,bmax);                      383   BoundingLimits(bmin,bmax);
368   G4BoundingEnvelope bbox(bmin,bmax);             384   G4BoundingEnvelope bbox(bmin,bmax);
369 #ifdef G4BBOX_EXTENT                              385 #ifdef G4BBOX_EXTENT
370   return bbox.CalculateExtent(pAxis,pVoxelLimi << 386   if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
371 #endif                                            387 #endif
372   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox    388   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
373   {                                               389   {
374     return exist = pMin < pMax;                << 390     return exist = (pMin < pMax) ? true : false;
375   }                                               391   }
376                                                   392 
377   // Set bounding envelope (benv) and calculat    393   // Set bounding envelope (benv) and calculate extent
378   //                                              394   //
379   G4double dz = GetZHalfLength();                 395   G4double dz = GetZHalfLength();
380   G4double dx = GetXHalfLength();                 396   G4double dx = GetXHalfLength();
381   G4double dy = GetYHalfLength();                 397   G4double dy = GetYHalfLength();
382                                                   398 
383   G4double x0 = dz*fTthetaCphi;                   399   G4double x0 = dz*fTthetaCphi;
384   G4double x1 = dy*GetTanAlpha();                 400   G4double x1 = dy*GetTanAlpha();
385   G4double y0 = dz*fTthetaSphi;                   401   G4double y0 = dz*fTthetaSphi;
386                                                   402 
387   G4ThreeVectorList baseA(4), baseB(4);           403   G4ThreeVectorList baseA(4), baseB(4);
388   baseA[0].set(-x0-x1-dx,-y0-dy,-dz);             404   baseA[0].set(-x0-x1-dx,-y0-dy,-dz);
389   baseA[1].set(-x0-x1+dx,-y0-dy,-dz);             405   baseA[1].set(-x0-x1+dx,-y0-dy,-dz);
390   baseA[2].set(-x0+x1+dx,-y0+dy,-dz);             406   baseA[2].set(-x0+x1+dx,-y0+dy,-dz);
391   baseA[3].set(-x0+x1-dx,-y0+dy,-dz);             407   baseA[3].set(-x0+x1-dx,-y0+dy,-dz);
392                                                   408 
393   baseB[0].set(+x0-x1-dx, y0-dy, dz);             409   baseB[0].set(+x0-x1-dx, y0-dy, dz);
394   baseB[1].set(+x0-x1+dx, y0-dy, dz);             410   baseB[1].set(+x0-x1+dx, y0-dy, dz);
395   baseB[2].set(+x0+x1+dx, y0+dy, dz);             411   baseB[2].set(+x0+x1+dx, y0+dy, dz);
396   baseB[3].set(+x0+x1-dx, y0+dy, dz);             412   baseB[3].set(+x0+x1-dx, y0+dy, dz);
397                                                   413 
398   std::vector<const G4ThreeVectorList *> polyg    414   std::vector<const G4ThreeVectorList *> polygons(2);
399   polygons[0] = &baseA;                           415   polygons[0] = &baseA;
400   polygons[1] = &baseB;                           416   polygons[1] = &baseB;
401                                                   417 
402   G4BoundingEnvelope benv(bmin,bmax,polygons);    418   G4BoundingEnvelope benv(bmin,bmax,polygons);
403   exist = benv.CalculateExtent(pAxis,pVoxelLim    419   exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
404   return exist;                                   420   return exist;
405 }                                                 421 }
406                                                   422 
407 //////////////////////////////////////////////    423 //////////////////////////////////////////////////////////////////////////
408 //                                                424 //
409 // Determine where is point p, inside/on_surfa    425 // Determine where is point p, inside/on_surface/outside
410 //                                                426 //
411                                                   427 
412 EInside G4Para::Inside( const G4ThreeVector& p    428 EInside G4Para::Inside( const G4ThreeVector& p ) const
413 {                                                 429 {
414   G4double xx = fPlanes[2].a*p.x()+fPlanes[2].    430   G4double xx = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z();
415   G4double dx = std::abs(xx) + fPlanes[2].d;      431   G4double dx = std::abs(xx) + fPlanes[2].d;
416                                                   432 
417   G4double yy = fPlanes[0].b*p.y()+fPlanes[0].    433   G4double yy = fPlanes[0].b*p.y()+fPlanes[0].c*p.z();
418   G4double dy = std::abs(yy) + fPlanes[0].d;      434   G4double dy = std::abs(yy) + fPlanes[0].d;
419   G4double dxy = std::max(dx,dy);                 435   G4double dxy = std::max(dx,dy);
420                                                   436 
421   G4double dz = std::abs(p.z())-fDz;              437   G4double dz = std::abs(p.z())-fDz;
422   G4double dist = std::max(dxy,dz);               438   G4double dist = std::max(dxy,dz);
423                                                   439 
424   if (dist > halfCarTolerance) return kOutside    440   if (dist > halfCarTolerance) return kOutside;
425   return (dist > -halfCarTolerance) ? kSurface    441   return (dist > -halfCarTolerance) ? kSurface : kInside;
426 }                                                 442 }
427                                                   443 
428 //////////////////////////////////////////////    444 //////////////////////////////////////////////////////////////////////////
429 //                                                445 //
430 // Determine side where point is, and return c    446 // Determine side where point is, and return corresponding normal
431                                                   447 
432 G4ThreeVector G4Para::SurfaceNormal( const G4T    448 G4ThreeVector G4Para::SurfaceNormal( const G4ThreeVector& p ) const
433 {                                                 449 {
434   G4int nsurf = 0; // number of surfaces where    450   G4int nsurf = 0; // number of surfaces where p is placed
435                                                   451 
436   // Check Z faces                                452   // Check Z faces
437   //                                              453   //
438   G4double nz = 0;                                454   G4double nz = 0;
439   G4double dz = std::abs(p.z()) - fDz;            455   G4double dz = std::abs(p.z()) - fDz;
440   if (std::abs(dz) <= halfCarTolerance)           456   if (std::abs(dz) <= halfCarTolerance)
441   {                                               457   {
442     nz = (p.z() < 0) ? -1 : 1;                    458     nz = (p.z() < 0) ? -1 : 1;
443     ++nsurf;                                      459     ++nsurf;
444   }                                               460   }
445                                                   461 
446   // Check Y faces                                462   // Check Y faces
447   //                                              463   //
448   G4double ny = 0;                                464   G4double ny = 0;
449   G4double yy = fPlanes[0].b*p.y()+fPlanes[0].    465   G4double yy = fPlanes[0].b*p.y()+fPlanes[0].c*p.z();
450   if (std::abs(fPlanes[0].d + yy) <= halfCarTo    466   if (std::abs(fPlanes[0].d + yy) <= halfCarTolerance)
451   {                                               467   {
452     ny  = fPlanes[0].b;                           468     ny  = fPlanes[0].b;
453     nz += fPlanes[0].c;                           469     nz += fPlanes[0].c;
454     ++nsurf;                                      470     ++nsurf;
455   }                                               471   }
456   else if (std::abs(fPlanes[1].d - yy) <= half    472   else if (std::abs(fPlanes[1].d - yy) <= halfCarTolerance)
457   {                                               473   {
458     ny  = fPlanes[1].b;                           474     ny  = fPlanes[1].b;
459     nz += fPlanes[1].c;                           475     nz += fPlanes[1].c;
460     ++nsurf;                                      476     ++nsurf;
461   }                                               477   }
462                                                   478 
463   // Check X faces                                479   // Check X faces
464   //                                              480   //
465   G4double nx = 0;                                481   G4double nx = 0;
466   G4double xx = fPlanes[2].a*p.x()+fPlanes[2].    482   G4double xx = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z();
467   if (std::abs(fPlanes[2].d + xx) <= halfCarTo    483   if (std::abs(fPlanes[2].d + xx) <= halfCarTolerance)
468   {                                               484   {
469     nx  = fPlanes[2].a;                           485     nx  = fPlanes[2].a;
470     ny += fPlanes[2].b;                           486     ny += fPlanes[2].b;
471     nz += fPlanes[2].c;                           487     nz += fPlanes[2].c;
472     ++nsurf;                                      488     ++nsurf;
473   }                                               489   }
474   else if (std::abs(fPlanes[3].d - xx) <= half    490   else if (std::abs(fPlanes[3].d - xx) <= halfCarTolerance)
475   {                                               491   {
476     nx  = fPlanes[3].a;                           492     nx  = fPlanes[3].a;
477     ny += fPlanes[3].b;                           493     ny += fPlanes[3].b;
478     nz += fPlanes[3].c;                           494     nz += fPlanes[3].c;
479     ++nsurf;                                      495     ++nsurf;
480   }                                               496   }
481                                                   497 
482   // Return normal                                498   // Return normal
483   //                                              499   //
484   if (nsurf == 1)      return {nx,ny,nz};      << 500   if (nsurf == 1)      return G4ThreeVector(nx,ny,nz);
485   else if (nsurf != 0) return G4ThreeVector(nx    501   else if (nsurf != 0) return G4ThreeVector(nx,ny,nz).unit(); // edge or corner
486   else                                            502   else
487   {                                               503   {
488     // Point is not on the surface                504     // Point is not on the surface
489     //                                            505     //
490 #ifdef G4CSGDEBUG                                 506 #ifdef G4CSGDEBUG
491     std::ostringstream message;                   507     std::ostringstream message;
492     G4int oldprc = message.precision(16);         508     G4int oldprc = message.precision(16);
493     message << "Point p is not on surface (!?)    509     message << "Point p is not on surface (!?) of solid: "
494             << GetName() << G4endl;               510             << GetName() << G4endl;
495     message << "Position:\n";                     511     message << "Position:\n";
496     message << "   p.x() = " << p.x()/mm << "     512     message << "   p.x() = " << p.x()/mm << " mm\n";
497     message << "   p.y() = " << p.y()/mm << "     513     message << "   p.y() = " << p.y()/mm << " mm\n";
498     message << "   p.z() = " << p.z()/mm << "     514     message << "   p.z() = " << p.z()/mm << " mm";
499     G4cout.precision(oldprc) ;                    515     G4cout.precision(oldprc) ;
500     G4Exception("G4Para::SurfaceNormal(p)", "G    516     G4Exception("G4Para::SurfaceNormal(p)", "GeomSolids1002",
501                 JustWarning, message );           517                 JustWarning, message );
502     DumpInfo();                                   518     DumpInfo();
503 #endif                                            519 #endif
504     return ApproxSurfaceNormal(p);                520     return ApproxSurfaceNormal(p);
505   }                                               521   }
506 }                                                 522 }
507                                                   523 
508 //////////////////////////////////////////////    524 //////////////////////////////////////////////////////////////////////////
509 //                                                525 //
510 // Algorithm for SurfaceNormal() following the    526 // Algorithm for SurfaceNormal() following the original specification
511 // for points not on the surface                  527 // for points not on the surface
512                                                   528 
513 G4ThreeVector G4Para::ApproxSurfaceNormal( con    529 G4ThreeVector G4Para::ApproxSurfaceNormal( const G4ThreeVector& p ) const
514 {                                                 530 {
515   G4double dist = -DBL_MAX;                       531   G4double dist = -DBL_MAX;
516   G4int iside = 0;                                532   G4int iside = 0;
517   for (G4int i=0; i<4; ++i)                       533   for (G4int i=0; i<4; ++i)
518   {                                               534   {
519     G4double d = fPlanes[i].a*p.x() +             535     G4double d = fPlanes[i].a*p.x() +
520                  fPlanes[i].b*p.y() +             536                  fPlanes[i].b*p.y() +
521                  fPlanes[i].c*p.z() + fPlanes[    537                  fPlanes[i].c*p.z() + fPlanes[i].d;
522     if (d > dist) { dist = d; iside = i; }        538     if (d > dist) { dist = d; iside = i; }
523   }                                               539   }
524                                                   540 
525   G4double distz = std::abs(p.z()) - fDz;         541   G4double distz = std::abs(p.z()) - fDz;
526   if (dist > distz)                               542   if (dist > distz)
527     return { fPlanes[iside].a, fPlanes[iside]. << 543     return G4ThreeVector(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c);
528   else                                            544   else
529     return { 0, 0, (G4double)((p.z() < 0) ? -1 << 545     return G4ThreeVector(0, 0, (p.z() < 0) ? -1 : 1);
530 }                                                 546 }
531                                                   547 
532 //////////////////////////////////////////////    548 //////////////////////////////////////////////////////////////////////////
533 //                                                549 //
534 // Calculate distance to shape from outside       550 // Calculate distance to shape from outside
535 //  - return kInfinity if no intersection         551 //  - return kInfinity if no intersection
536                                                   552 
537 G4double G4Para::DistanceToIn(const G4ThreeVec    553 G4double G4Para::DistanceToIn(const G4ThreeVector& p,
538                               const G4ThreeVec    554                               const G4ThreeVector& v ) const
539 {                                                 555 {
540   // Z intersections                              556   // Z intersections
541   //                                              557   //
542   if ((std::abs(p.z()) - fDz) >= -halfCarToler    558   if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() >= 0)
543     return kInfinity;                             559     return kInfinity;
544   G4double invz = (-v.z() == 0) ? DBL_MAX : -1    560   G4double invz = (-v.z() == 0) ? DBL_MAX : -1./v.z();
545   G4double dz = (invz < 0) ? fDz : -fDz;          561   G4double dz = (invz < 0) ? fDz : -fDz; 
546   G4double tzmin = (p.z() + dz)*invz;             562   G4double tzmin = (p.z() + dz)*invz;
547   G4double tzmax = (p.z() - dz)*invz;             563   G4double tzmax = (p.z() - dz)*invz;
548                                                   564 
549   // Y intersections                              565   // Y intersections
550   //                                              566   //
551   G4double tmin0 = tzmin, tmax0 = tzmax;          567   G4double tmin0 = tzmin, tmax0 = tzmax;
552   G4double cos0 = fPlanes[0].b*v.y() + fPlanes    568   G4double cos0 = fPlanes[0].b*v.y() + fPlanes[0].c*v.z();
553   G4double disy = fPlanes[0].b*p.y() + fPlanes    569   G4double disy = fPlanes[0].b*p.y() + fPlanes[0].c*p.z();
554   G4double dis0 = fPlanes[0].d + disy;            570   G4double dis0 = fPlanes[0].d + disy;
555   if (dis0 >= -halfCarTolerance)                  571   if (dis0 >= -halfCarTolerance)
556   {                                               572   {
557     if (cos0 >= 0) return kInfinity;              573     if (cos0 >= 0) return kInfinity;
558     G4double tmp  = -dis0/cos0;                   574     G4double tmp  = -dis0/cos0;
559     if (tmin0 < tmp) tmin0 = tmp;                 575     if (tmin0 < tmp) tmin0 = tmp;
560   }                                               576   }
561   else if (cos0 > 0)                              577   else if (cos0 > 0)
562   {                                               578   {
563     G4double tmp  = -dis0/cos0;                   579     G4double tmp  = -dis0/cos0;
564     if (tmax0 > tmp) tmax0 = tmp;                 580     if (tmax0 > tmp) tmax0 = tmp;
565   }                                               581   }
566                                                   582 
567   G4double tmin1 = tmin0, tmax1 = tmax0;          583   G4double tmin1 = tmin0, tmax1 = tmax0;
568   G4double cos1 = -cos0;                          584   G4double cos1 = -cos0;
569   G4double dis1 = fPlanes[1].d - disy;            585   G4double dis1 = fPlanes[1].d - disy;
570   if (dis1 >= -halfCarTolerance)                  586   if (dis1 >= -halfCarTolerance)
571   {                                               587   {
572     if (cos1 >= 0) return kInfinity;              588     if (cos1 >= 0) return kInfinity;
573     G4double tmp  = -dis1/cos1;                   589     G4double tmp  = -dis1/cos1;
574     if (tmin1 < tmp) tmin1 = tmp;                 590     if (tmin1 < tmp) tmin1 = tmp;
575   }                                               591   }
576   else if (cos1 > 0)                              592   else if (cos1 > 0)
577   {                                               593   {
578     G4double tmp  = -dis1/cos1;                   594     G4double tmp  = -dis1/cos1;
579     if (tmax1 > tmp) tmax1 = tmp;                 595     if (tmax1 > tmp) tmax1 = tmp;
580   }                                               596   }
581                                                   597 
582   // X intersections                              598   // X intersections
583   //                                              599   //
584   G4double tmin2 = tmin1, tmax2 = tmax1;          600   G4double tmin2 = tmin1, tmax2 = tmax1;
585   G4double cos2 = fPlanes[2].a*v.x() + fPlanes    601   G4double cos2 = fPlanes[2].a*v.x() + fPlanes[2].b*v.y() + fPlanes[2].c*v.z();
586   G4double disx = fPlanes[2].a*p.x() + fPlanes    602   G4double disx = fPlanes[2].a*p.x() + fPlanes[2].b*p.y() + fPlanes[2].c*p.z();
587   G4double dis2 = fPlanes[2].d + disx;            603   G4double dis2 = fPlanes[2].d + disx;
588   if (dis2 >= -halfCarTolerance)                  604   if (dis2 >= -halfCarTolerance)
589   {                                               605   {
590     if (cos2 >= 0) return kInfinity;              606     if (cos2 >= 0) return kInfinity;
591     G4double tmp  = -dis2/cos2;                   607     G4double tmp  = -dis2/cos2;
592     if (tmin2 < tmp) tmin2 = tmp;                 608     if (tmin2 < tmp) tmin2 = tmp;
593   }                                               609   }
594   else if (cos2 > 0)                              610   else if (cos2 > 0)
595   {                                               611   {
596     G4double tmp  = -dis2/cos2;                   612     G4double tmp  = -dis2/cos2;
597     if (tmax2 > tmp) tmax2 = tmp;                 613     if (tmax2 > tmp) tmax2 = tmp;
598   }                                               614   }
599                                                   615 
600   G4double tmin3 = tmin2, tmax3 = tmax2;          616   G4double tmin3 = tmin2, tmax3 = tmax2;
601   G4double cos3 = -cos2;                          617   G4double cos3 = -cos2;
602   G4double dis3 = fPlanes[3].d - disx;            618   G4double dis3 = fPlanes[3].d - disx;
603   if (dis3 >= -halfCarTolerance)                  619   if (dis3 >= -halfCarTolerance)
604   {                                               620   {
605     if (cos3 >= 0) return kInfinity;              621     if (cos3 >= 0) return kInfinity;
606     G4double tmp  = -dis3/cos3;                   622     G4double tmp  = -dis3/cos3;
607     if (tmin3 < tmp) tmin3 = tmp;                 623     if (tmin3 < tmp) tmin3 = tmp;
608   }                                               624   }
609   else if (cos3 > 0)                              625   else if (cos3 > 0)
610   {                                               626   {
611     G4double tmp  = -dis3/cos3;                   627     G4double tmp  = -dis3/cos3;
612     if (tmax3 > tmp) tmax3 = tmp;                 628     if (tmax3 > tmp) tmax3 = tmp;
613   }                                               629   }
614                                                   630 
615   // Find distance                                631   // Find distance
616   //                                              632   //
617   G4double tmin = tmin3, tmax = tmax3;            633   G4double tmin = tmin3, tmax = tmax3;
618   if (tmax <= tmin + halfCarTolerance) return     634   if (tmax <= tmin + halfCarTolerance) return kInfinity; // touch or no hit
619   return (tmin < halfCarTolerance ) ? 0. : tmi    635   return (tmin < halfCarTolerance ) ? 0. : tmin;
620 }                                                 636 }
621                                                   637 
622 //////////////////////////////////////////////    638 //////////////////////////////////////////////////////////////////////////
623 //                                                639 //
624 // Calculate exact shortest distance to any bo    640 // Calculate exact shortest distance to any boundary from outside
625 // - returns 0 is point inside                    641 // - returns 0 is point inside
626                                                   642 
627 G4double G4Para::DistanceToIn( const G4ThreeVe    643 G4double G4Para::DistanceToIn( const G4ThreeVector& p ) const
628 {                                                 644 {
629   G4double xx = fPlanes[2].a*p.x()+fPlanes[2].    645   G4double xx = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z();
630   G4double dx = std::abs(xx) + fPlanes[2].d;      646   G4double dx = std::abs(xx) + fPlanes[2].d;
631                                                   647 
632   G4double yy = fPlanes[0].b*p.y()+fPlanes[0].    648   G4double yy = fPlanes[0].b*p.y()+fPlanes[0].c*p.z();
633   G4double dy = std::abs(yy) + fPlanes[0].d;      649   G4double dy = std::abs(yy) + fPlanes[0].d;
634   G4double dxy = std::max(dx,dy);                 650   G4double dxy = std::max(dx,dy);
635                                                   651 
636   G4double dz = std::abs(p.z())-fDz;              652   G4double dz = std::abs(p.z())-fDz;
637   G4double dist = std::max(dxy,dz);               653   G4double dist = std::max(dxy,dz);
638                                                   654 
639   return (dist > 0) ? dist : 0.;                  655   return (dist > 0) ? dist : 0.;
640 }                                                 656 }
641                                                   657 
642 //////////////////////////////////////////////    658 //////////////////////////////////////////////////////////////////////////
643 //                                                659 //
644 // Calculate distance to surface of shape from    660 // Calculate distance to surface of shape from inside and, if required,
645 // find normal at exit point                      661 // find normal at exit point
646 // - when leaving the surface, return 0           662 // - when leaving the surface, return 0
647                                                   663 
648 G4double G4Para::DistanceToOut(const G4ThreeVe    664 G4double G4Para::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
649                                const G4bool ca    665                                const G4bool calcNorm,
650                                      G4bool* v << 666                                      G4bool *validNorm, G4ThreeVector *n) const
651 {                                                 667 {
652   // Z intersections                              668   // Z intersections
653   //                                              669   //
654   if ((std::abs(p.z()) - fDz) >= -halfCarToler    670   if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() > 0)
655   {                                               671   {
656     if (calcNorm)                                 672     if (calcNorm)
657     {                                             673     {
658       *validNorm = true;                          674       *validNorm = true;
659       n->set(0, 0, (p.z() < 0) ? -1 : 1);         675       n->set(0, 0, (p.z() < 0) ? -1 : 1);
660     }                                             676     }
661     return 0.;                                    677     return 0.;
662   }                                               678   }
663   G4double vz = v.z();                            679   G4double vz = v.z();
664   G4double tmax = (vz == 0) ? DBL_MAX : (std::    680   G4double tmax = (vz == 0) ? DBL_MAX : (std::copysign(fDz,vz) - p.z())/vz;
665   G4int iside = (vz < 0) ? -4 : -2; // little     681   G4int iside = (vz < 0) ? -4 : -2; // little trick: (-4+3)=-1, (-2+3)=+1
666                                                   682 
667   // Y intersections                              683   // Y intersections
668   //                                              684   //
669   G4double cos0 = fPlanes[0].b*v.y() + fPlanes    685   G4double cos0 = fPlanes[0].b*v.y() + fPlanes[0].c*v.z();
670   if (cos0 > 0)                                   686   if (cos0 > 0)
671   {                                               687   {
672     G4double dis0 = fPlanes[0].b*p.y() + fPlan    688     G4double dis0 = fPlanes[0].b*p.y() + fPlanes[0].c*p.z() + fPlanes[0].d;
673     if (dis0 >= -halfCarTolerance)                689     if (dis0 >= -halfCarTolerance)
674     {                                             690     {
675       if (calcNorm)                               691       if (calcNorm)
676       {                                           692       {
677         *validNorm = true;                        693         *validNorm = true;
678         n->set(0, fPlanes[0].b, fPlanes[0].c);    694         n->set(0, fPlanes[0].b, fPlanes[0].c);
679       }                                           695       }
680       return 0.;                                  696       return 0.;
681     }                                             697     }
682     G4double tmp = -dis0/cos0;                    698     G4double tmp = -dis0/cos0;
683     if (tmax > tmp) { tmax = tmp; iside = 0; }    699     if (tmax > tmp) { tmax = tmp; iside = 0; }
684   }                                               700   }
685                                                   701 
686   G4double cos1 = -cos0;                          702   G4double cos1 = -cos0;
687   if (cos1 > 0)                                   703   if (cos1 > 0)
688   {                                               704   {
689     G4double dis1 = fPlanes[1].b*p.y() + fPlan    705     G4double dis1 = fPlanes[1].b*p.y() + fPlanes[1].c*p.z() + fPlanes[1].d;
690     if (dis1 >= -halfCarTolerance)                706     if (dis1 >= -halfCarTolerance)
691     {                                             707     {
692       if (calcNorm)                               708       if (calcNorm)
693       {                                           709       {
694         *validNorm = true;                        710         *validNorm = true;
695         n->set(0, fPlanes[1].b, fPlanes[1].c);    711         n->set(0, fPlanes[1].b, fPlanes[1].c);
696       }                                           712       }
697       return 0.;                                  713       return 0.;
698     }                                             714     }
699     G4double tmp = -dis1/cos1;                    715     G4double tmp = -dis1/cos1;
700     if (tmax > tmp) { tmax = tmp; iside = 1; }    716     if (tmax > tmp) { tmax = tmp; iside = 1; }
701   }                                               717   }
702                                                   718 
703   // X intersections                              719   // X intersections
704   //                                              720   //
705   G4double cos2 = fPlanes[2].a*v.x() + fPlanes    721   G4double cos2 = fPlanes[2].a*v.x() + fPlanes[2].b*v.y() + fPlanes[2].c*v.z();
706   if (cos2 > 0)                                   722   if (cos2 > 0)
707   {                                               723   {
708     G4double dis2 = fPlanes[2].a*p.x()+fPlanes    724     G4double dis2 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z()+fPlanes[2].d;
709     if (dis2 >= -halfCarTolerance)                725     if (dis2 >= -halfCarTolerance)
710     {                                             726     {
711       if (calcNorm)                               727       if (calcNorm)
712       {                                           728       {
713          *validNorm = true;                       729          *validNorm = true;
714          n->set(fPlanes[2].a, fPlanes[2].b, fP    730          n->set(fPlanes[2].a, fPlanes[2].b, fPlanes[2].c);
715       }                                           731       }
716       return 0.;                                  732       return 0.;
717     }                                             733     }
718     G4double tmp = -dis2/cos2;                    734     G4double tmp = -dis2/cos2;
719     if (tmax > tmp) { tmax = tmp; iside = 2; }    735     if (tmax > tmp) { tmax = tmp; iside = 2; }
720   }                                               736   }
721                                                   737 
722   G4double cos3 = -cos2;                          738   G4double cos3 = -cos2;
723   if (cos3 > 0)                                   739   if (cos3 > 0)
724   {                                               740   {
725     G4double dis3 = fPlanes[3].a*p.x()+fPlanes    741     G4double dis3 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y()+fPlanes[3].c*p.z()+fPlanes[3].d;
726     if (dis3 >= -halfCarTolerance)                742     if (dis3 >= -halfCarTolerance)
727     {                                             743     {
728       if (calcNorm)                               744       if (calcNorm)
729       {                                           745       {
730          *validNorm = true;                       746          *validNorm = true;
731          n->set(fPlanes[3].a, fPlanes[3].b, fP    747          n->set(fPlanes[3].a, fPlanes[3].b, fPlanes[3].c);
732       }                                           748       }
733       return 0.;                                  749       return 0.;
734     }                                             750     }
735     G4double tmp = -dis3/cos3;                    751     G4double tmp = -dis3/cos3;
736     if (tmax > tmp) { tmax = tmp; iside = 3; }    752     if (tmax > tmp) { tmax = tmp; iside = 3; }
737   }                                               753   }
738                                                   754 
739   // Set normal, if required, and return dista    755   // Set normal, if required, and return distance
740   //                                              756   //
741   if (calcNorm)                                   757   if (calcNorm) 
742   {                                               758   {
743     *validNorm = true;                            759     *validNorm = true;
744     if (iside < 0)                                760     if (iside < 0)
745       n->set(0, 0, iside + 3); // (-4+3)=-1, (    761       n->set(0, 0, iside + 3); // (-4+3)=-1, (-2+3)=+1
746     else                                          762     else
747       n->set(fPlanes[iside].a, fPlanes[iside].    763       n->set(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c);
748   }                                               764   }
749   return tmax;                                    765   return tmax;
750 }                                                 766 }
751                                                   767 
752 //////////////////////////////////////////////    768 //////////////////////////////////////////////////////////////////////////
753 //                                                769 //
754 // Calculate exact shortest distance to any bo    770 // Calculate exact shortest distance to any boundary from inside
755 // - returns 0 is point outside                   771 // - returns 0 is point outside
756                                                   772 
757 G4double G4Para::DistanceToOut( const G4ThreeV    773 G4double G4Para::DistanceToOut( const G4ThreeVector& p ) const
758 {                                                 774 {
759 #ifdef G4CSGDEBUG                                 775 #ifdef G4CSGDEBUG
760   if( Inside(p) == kOutside )                     776   if( Inside(p) == kOutside )
761   {                                               777   {
762     std::ostringstream message;                   778     std::ostringstream message;
763     G4int oldprc = message.precision(16);         779     G4int oldprc = message.precision(16);
764     message << "Point p is outside (!?) of sol    780     message << "Point p is outside (!?) of solid: " << GetName() << G4endl;
765     message << "Position:\n";                     781     message << "Position:\n";
766     message << "   p.x() = " << p.x()/mm << "     782     message << "   p.x() = " << p.x()/mm << " mm\n";
767     message << "   p.y() = " << p.y()/mm << "     783     message << "   p.y() = " << p.y()/mm << " mm\n";
768     message << "   p.z() = " << p.z()/mm << "     784     message << "   p.z() = " << p.z()/mm << " mm";
769     G4cout.precision(oldprc) ;                    785     G4cout.precision(oldprc) ;
770     G4Exception("G4Para::DistanceToOut(p)", "G    786     G4Exception("G4Para::DistanceToOut(p)", "GeomSolids1002",
771                 JustWarning, message );           787                 JustWarning, message );
772     DumpInfo();                                   788     DumpInfo();
773     }                                             789     }
774 #endif                                            790 #endif
775   G4double xx = fPlanes[2].a*p.x()+fPlanes[2].    791   G4double xx = fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z();
776   G4double dx = std::abs(xx) + fPlanes[2].d;      792   G4double dx = std::abs(xx) + fPlanes[2].d;
777                                                   793 
778   G4double yy = fPlanes[0].b*p.y()+fPlanes[0].    794   G4double yy = fPlanes[0].b*p.y()+fPlanes[0].c*p.z();
779   G4double dy = std::abs(yy) + fPlanes[0].d;      795   G4double dy = std::abs(yy) + fPlanes[0].d;
780   G4double dxy = std::max(dx,dy);                 796   G4double dxy = std::max(dx,dy);
781                                                   797 
782   G4double dz = std::abs(p.z())-fDz;              798   G4double dz = std::abs(p.z())-fDz;
783   G4double dist = std::max(dxy,dz);               799   G4double dist = std::max(dxy,dz);
784                                                   800 
785   return (dist < 0) ? -dist : 0.;                 801   return (dist < 0) ? -dist : 0.;
786 }                                                 802 }
787                                                   803 
788 //////////////////////////////////////////////    804 //////////////////////////////////////////////////////////////////////////
789 //                                                805 //
790 // GetEntityType                                  806 // GetEntityType
791                                                   807 
792 G4GeometryType G4Para::GetEntityType() const      808 G4GeometryType G4Para::GetEntityType() const
793 {                                                 809 {
794   return {"G4Para"};                           << 810   return G4String("G4Para");
795 }                                              << 
796                                                << 
797 ////////////////////////////////////////////// << 
798 //                                             << 
799 // IsFaceted                                   << 
800                                                << 
801 G4bool G4Para::IsFaceted() const               << 
802 {                                              << 
803   return true;                                 << 
804 }                                                 811 }
805                                                   812 
806 //////////////////////////////////////////////    813 //////////////////////////////////////////////////////////////////////////
807 //                                                814 //
808 // Make a clone of the object                     815 // Make a clone of the object
809 //                                                816 //
810 G4VSolid* G4Para::Clone() const                   817 G4VSolid* G4Para::Clone() const
811 {                                                 818 {
812   return new G4Para(*this);                       819   return new G4Para(*this);
813 }                                                 820 }
814                                                   821 
815 //////////////////////////////////////////////    822 //////////////////////////////////////////////////////////////////////////
816 //                                                823 //
817 // Stream object contents to an output stream     824 // Stream object contents to an output stream
818                                                   825 
819 std::ostream& G4Para::StreamInfo( std::ostream    826 std::ostream& G4Para::StreamInfo( std::ostream& os ) const
820 {                                                 827 {
821   G4double alpha = std::atan(fTalpha);            828   G4double alpha = std::atan(fTalpha);
822   G4double theta = std::atan(std::sqrt(fTtheta    829   G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi +
823                                        fTtheta    830                                        fTthetaSphi*fTthetaSphi));
824   G4double phi   = std::atan2(fTthetaSphi,fTth    831   G4double phi   = std::atan2(fTthetaSphi,fTthetaCphi);
                                                   >> 832   G4String signDegree = "\u00B0";
825                                                   833 
826   G4long oldprc = os.precision(16);            << 834   G4int oldprc = os.precision(16);
827   os << "-------------------------------------    835   os << "-----------------------------------------------------------\n"
828      << "    *** Dump for solid - " << GetName    836      << "    *** Dump for solid - " << GetName() << " ***\n"
829      << "    =================================    837      << "    ===================================================\n"
830      << " Solid type: G4Para\n"                   838      << " Solid type: G4Para\n"
831      << " Parameters:\n"                          839      << " Parameters:\n"
832      << "    half length X: " << fDx/mm << " m    840      << "    half length X: " << fDx/mm << " mm\n"
833      << "    half length Y: " << fDy/mm << " m    841      << "    half length Y: " << fDy/mm << " mm\n"
834      << "    half length Z: " << fDz/mm << " m    842      << "    half length Z: " << fDz/mm << " mm\n"
835      << "    alpha: " << alpha/degree << "degr << 843      << "    alpha: " << alpha/degree << signDegree << "\n"
836      << "    theta: " << theta/degree << "degr << 844      << "    theta: " << theta/degree << signDegree << "\n"
837      << "    phi: " << phi/degree << "degrees\ << 845      << "    phi: " << phi/degree << signDegree << "\n"
838      << "-------------------------------------    846      << "-----------------------------------------------------------\n";
839   os.precision(oldprc);                           847   os.precision(oldprc);
840                                                   848 
841   return os;                                      849   return os;
842 }                                                 850 }
843                                                   851 
844 //////////////////////////////////////////////    852 //////////////////////////////////////////////////////////////////////////
845 //                                                853 //
846 // Return a point randomly and uniformly selec    854 // Return a point randomly and uniformly selected on the solid surface
847                                                   855 
848 G4ThreeVector G4Para::GetPointOnSurface() cons    856 G4ThreeVector G4Para::GetPointOnSurface() const
849 {                                                 857 {
850   G4double DyTalpha = fDy*fTalpha;                858   G4double DyTalpha = fDy*fTalpha;
851   G4double DzTthetaSphi = fDz*fTthetaSphi;        859   G4double DzTthetaSphi = fDz*fTthetaSphi;
852   G4double DzTthetaCphi = fDz*fTthetaCphi;        860   G4double DzTthetaCphi = fDz*fTthetaCphi;
853                                                   861 
854   // Set vertices                                 862   // Set vertices
855   //                                              863   //
856   G4ThreeVector pt[8];                            864   G4ThreeVector pt[8];
857   pt[0].set(-DzTthetaCphi-DyTalpha-fDx, -DzTth    865   pt[0].set(-DzTthetaCphi-DyTalpha-fDx, -DzTthetaSphi-fDy, -fDz);
858   pt[1].set(-DzTthetaCphi-DyTalpha+fDx, -DzTth    866   pt[1].set(-DzTthetaCphi-DyTalpha+fDx, -DzTthetaSphi-fDy, -fDz);
859   pt[2].set(-DzTthetaCphi+DyTalpha-fDx, -DzTth    867   pt[2].set(-DzTthetaCphi+DyTalpha-fDx, -DzTthetaSphi+fDy, -fDz);
860   pt[3].set(-DzTthetaCphi+DyTalpha+fDx, -DzTth    868   pt[3].set(-DzTthetaCphi+DyTalpha+fDx, -DzTthetaSphi+fDy, -fDz);
861   pt[4].set( DzTthetaCphi-DyTalpha-fDx,  DzTth    869   pt[4].set( DzTthetaCphi-DyTalpha-fDx,  DzTthetaSphi-fDy,  fDz);
862   pt[5].set( DzTthetaCphi-DyTalpha+fDx,  DzTth    870   pt[5].set( DzTthetaCphi-DyTalpha+fDx,  DzTthetaSphi-fDy,  fDz);
863   pt[6].set( DzTthetaCphi+DyTalpha-fDx,  DzTth    871   pt[6].set( DzTthetaCphi+DyTalpha-fDx,  DzTthetaSphi+fDy,  fDz);
864   pt[7].set( DzTthetaCphi+DyTalpha+fDx,  DzTth    872   pt[7].set( DzTthetaCphi+DyTalpha+fDx,  DzTthetaSphi+fDy,  fDz);
865                                                   873 
866   // Set areas (-Z, -Y, +Y, -X, +X, +Z)           874   // Set areas (-Z, -Y, +Y, -X, +X, +Z)
867   //                                              875   //
868   G4ThreeVector vx(fDx, 0, 0);                    876   G4ThreeVector vx(fDx, 0, 0);
869   G4ThreeVector vy(DyTalpha, fDy, 0);             877   G4ThreeVector vy(DyTalpha, fDy, 0);
870   G4ThreeVector vz(DzTthetaCphi, DzTthetaSphi,    878   G4ThreeVector vz(DzTthetaCphi, DzTthetaSphi, fDz);
871                                                   879 
872   G4double sxy = fDx*fDy; // (vx.cross(vy)).ma    880   G4double sxy = fDx*fDy; // (vx.cross(vy)).mag();
873   G4double sxz = (vx.cross(vz)).mag();            881   G4double sxz = (vx.cross(vz)).mag();
874   G4double syz = (vy.cross(vz)).mag();            882   G4double syz = (vy.cross(vz)).mag();
875                                                   883   
876   G4double sface[6] = { sxy, syz, syz, sxz, sx    884   G4double sface[6] = { sxy, syz, syz, sxz, sxz, sxy };
877   for (G4int i=1; i<6; ++i) { sface[i] += sfac    885   for (G4int i=1; i<6; ++i) { sface[i] += sface[i-1]; }
878                                                   886 
879   // Select face                                  887   // Select face
880   //                                              888   //
881   G4double select = sface[5]*G4UniformRand();     889   G4double select = sface[5]*G4UniformRand();
882   G4int k = 5;                                    890   G4int k = 5;
883   if (select <= sface[4]) k = 4;                  891   if (select <= sface[4]) k = 4;
884   if (select <= sface[3]) k = 3;                  892   if (select <= sface[3]) k = 3;
885   if (select <= sface[2]) k = 2;                  893   if (select <= sface[2]) k = 2;
886   if (select <= sface[1]) k = 1;                  894   if (select <= sface[1]) k = 1;
887   if (select <= sface[0]) k = 0;                  895   if (select <= sface[0]) k = 0;
888                                                   896 
889   // Generate point                               897   // Generate point
890   //                                              898   //
891   G4int ip[6][3] = {{0,1,2}, {0,4,1}, {2,3,6},    899   G4int ip[6][3] = {{0,1,2}, {0,4,1}, {2,3,6}, {0,2,4}, {1,5,3}, {4,6,5}};
892   G4double u = G4UniformRand();                   900   G4double u = G4UniformRand();
893   G4double v = G4UniformRand();                   901   G4double v = G4UniformRand();
894   return (1.-u-v)*pt[ip[k][0]] + u*pt[ip[k][1]    902   return (1.-u-v)*pt[ip[k][0]] + u*pt[ip[k][1]] + v*pt[ip[k][2]];
895 }                                                 903 }
896                                                   904 
897 //////////////////////////////////////////////    905 //////////////////////////////////////////////////////////////////////////
898 //                                                906 //
899 // Methods for visualisation                      907 // Methods for visualisation
900                                                   908 
901 void G4Para::DescribeYourselfTo ( G4VGraphicsS    909 void G4Para::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
902 {                                                 910 {
903   scene.AddSolid (*this);                         911   scene.AddSolid (*this);
904 }                                                 912 }
905                                                   913 
906 G4Polyhedron* G4Para::CreatePolyhedron () cons    914 G4Polyhedron* G4Para::CreatePolyhedron () const
907 {                                                 915 {
908   G4double phi = std::atan2(fTthetaSphi, fTthe    916   G4double phi = std::atan2(fTthetaSphi, fTthetaCphi);
909   G4double alpha = std::atan(fTalpha);            917   G4double alpha = std::atan(fTalpha);
910   G4double theta = std::atan(std::sqrt(fTtheta    918   G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi +
911                                        fTtheta    919                                        fTthetaSphi*fTthetaSphi));
912                                                   920     
913   return new G4PolyhedronPara(fDx, fDy, fDz, a    921   return new G4PolyhedronPara(fDx, fDy, fDz, alpha, theta, phi);
914 }                                                 922 }
915 #endif                                            923 #endif
916                                                   924