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 11.2.1)


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