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