| 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 intell 18 // * This code implementation is the intellectual property of the *
19 // * Vanderbilt University Free Electron Laser 19 // * Vanderbilt University Free Electron Laser Center *
20 // * Vanderbilt University, Nashville, TN, USA 20 // * Vanderbilt University, Nashville, TN, USA *
21 // * Development supported by: 21 // * Development supported by: *
22 // * United States MFEL program under grant F 22 // * United States MFEL program under grant FA9550-04-1-0045 *
23 // * and NASA under contract number NNG04CT05P 23 // * and NASA under contract number NNG04CT05P *
24 // * Written by Marcus H. Mendenhall and Rober 24 // * Written by Marcus H. Mendenhall and Robert A. Weller. *
25 // * 25 // * *
26 // * Contributed to the Geant4 Core, January, 26 // * Contributed to the Geant4 Core, January, 2005. *
27 // * 27 // * *
28 // ******************************************* 28 // ********************************************************************
29 // 29 //
30 // Implementation for G4Tet class 30 // Implementation for G4Tet class
31 // 31 //
32 // 03.09.2004 - Marcus Mendenhall, created 32 // 03.09.2004 - Marcus Mendenhall, created
33 // 08.01.2020 - Evgueni Tcherniaev, complete r 33 // 08.01.2020 - Evgueni Tcherniaev, complete revision, speed up
34 // ------------------------------------------- 34 // --------------------------------------------------------------------
35 35
36 #include "G4Tet.hh" 36 #include "G4Tet.hh"
37 37
38 #if !defined(G4GEOM_USE_UTET) 38 #if !defined(G4GEOM_USE_UTET)
39 39
40 #include "G4VoxelLimits.hh" 40 #include "G4VoxelLimits.hh"
41 #include "G4AffineTransform.hh" 41 #include "G4AffineTransform.hh"
42 #include "G4BoundingEnvelope.hh" 42 #include "G4BoundingEnvelope.hh"
43 43
44 #include "G4VPVParameterisation.hh" 44 #include "G4VPVParameterisation.hh"
45 45
46 #include "G4QuickRand.hh" 46 #include "G4QuickRand.hh"
47 47
48 #include "G4VGraphicsScene.hh" 48 #include "G4VGraphicsScene.hh"
49 #include "G4Polyhedron.hh" 49 #include "G4Polyhedron.hh"
50 #include "G4VisExtent.hh" 50 #include "G4VisExtent.hh"
51 51
52 #include "G4AutoLock.hh" 52 #include "G4AutoLock.hh"
53 53
54 namespace 54 namespace
55 { 55 {
56 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE 56 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
57 } 57 }
58 58
59 using namespace CLHEP; 59 using namespace CLHEP;
60 60
61 ////////////////////////////////////////////// 61 ////////////////////////////////////////////////////////////////////////
62 // 62 //
63 // Constructor - create a tetrahedron 63 // Constructor - create a tetrahedron
64 // A Tet has all of its geometrical informatio 64 // A Tet has all of its geometrical information precomputed
65 // 65 //
66 G4Tet::G4Tet(const G4String& pName, 66 G4Tet::G4Tet(const G4String& pName,
67 const G4ThreeVector& p0, 67 const G4ThreeVector& p0,
68 const G4ThreeVector& p1, 68 const G4ThreeVector& p1,
69 const G4ThreeVector& p2, 69 const G4ThreeVector& p2,
70 const G4ThreeVector& p3, G4bool* 70 const G4ThreeVector& p3, G4bool* degeneracyFlag)
71 : G4VSolid(pName) 71 : G4VSolid(pName)
72 { 72 {
73 // Check for degeneracy 73 // Check for degeneracy
74 G4bool degenerate = CheckDegeneracy(p0, p1, 74 G4bool degenerate = CheckDegeneracy(p0, p1, p2, p3);
75 if (degeneracyFlag != nullptr) << 75 if (degeneracyFlag)
76 { 76 {
77 *degeneracyFlag = degenerate; 77 *degeneracyFlag = degenerate;
78 } 78 }
79 else if (degenerate) 79 else if (degenerate)
80 { 80 {
81 std::ostringstream message; 81 std::ostringstream message;
82 message << "Degenerate tetrahedron: " << G 82 message << "Degenerate tetrahedron: " << GetName() << " !\n"
83 << " anchor: " << p0 << "\n" 83 << " anchor: " << p0 << "\n"
84 << " p1 : " << p1 << "\n" 84 << " p1 : " << p1 << "\n"
85 << " p2 : " << p2 << "\n" 85 << " p2 : " << p2 << "\n"
86 << " p3 : " << p3 << "\n" 86 << " p3 : " << p3 << "\n"
87 << " volume: " 87 << " volume: "
88 << std::abs((p1 - p0).cross(p2 - p 88 << std::abs((p1 - p0).cross(p2 - p0).dot(p3 - p0))/6.;
89 G4Exception("G4Tet::G4Tet()", "GeomSolids0 89 G4Exception("G4Tet::G4Tet()", "GeomSolids0002", FatalException, message);
90 } 90 }
91 91
92 // Define surface thickness 92 // Define surface thickness
93 halfTolerance = 0.5 * kCarTolerance; 93 halfTolerance = 0.5 * kCarTolerance;
94 94
95 // Set data members 95 // Set data members
96 Initialize(p0, p1, p2, p3); 96 Initialize(p0, p1, p2, p3);
97 } 97 }
98 98
99 ////////////////////////////////////////////// 99 ////////////////////////////////////////////////////////////////////////
100 // 100 //
101 // Fake default constructor - sets only member 101 // Fake default constructor - sets only member data and allocates memory
102 // for usage restri 102 // for usage restricted to object persistency.
103 // 103 //
104 G4Tet::G4Tet( __void__& a ) 104 G4Tet::G4Tet( __void__& a )
105 : G4VSolid(a) 105 : G4VSolid(a)
106 { 106 {
107 } 107 }
108 108
109 ////////////////////////////////////////////// 109 ////////////////////////////////////////////////////////////////////////
110 // 110 //
111 // Destructor 111 // Destructor
112 // 112 //
113 G4Tet::~G4Tet() 113 G4Tet::~G4Tet()
114 { 114 {
115 delete fpPolyhedron; fpPolyhedron = nullptr; 115 delete fpPolyhedron; fpPolyhedron = nullptr;
116 } 116 }
117 117
118 ////////////////////////////////////////////// 118 ////////////////////////////////////////////////////////////////////////
119 // 119 //
120 // Copy constructor 120 // Copy constructor
121 // 121 //
122 G4Tet::G4Tet(const G4Tet& rhs) 122 G4Tet::G4Tet(const G4Tet& rhs)
123 : G4VSolid(rhs) 123 : G4VSolid(rhs)
124 { 124 {
125 halfTolerance = rhs.halfTolerance; 125 halfTolerance = rhs.halfTolerance;
126 fCubicVolume = rhs.fCubicVolume; 126 fCubicVolume = rhs.fCubicVolume;
127 fSurfaceArea = rhs.fSurfaceArea; 127 fSurfaceArea = rhs.fSurfaceArea;
128 for (G4int i = 0; i < 4; ++i) { fVertex[i] 128 for (G4int i = 0; i < 4; ++i) { fVertex[i] = rhs.fVertex[i]; }
129 for (G4int i = 0; i < 4; ++i) { fNormal[i] 129 for (G4int i = 0; i < 4; ++i) { fNormal[i] = rhs.fNormal[i]; }
130 for (G4int i = 0; i < 4; ++i) { fDist[i] = 130 for (G4int i = 0; i < 4; ++i) { fDist[i] = rhs.fDist[i]; }
131 for (G4int i = 0; i < 4; ++i) { fArea[i] = 131 for (G4int i = 0; i < 4; ++i) { fArea[i] = rhs.fArea[i]; }
132 fBmin = rhs.fBmin; 132 fBmin = rhs.fBmin;
133 fBmax = rhs.fBmax; 133 fBmax = rhs.fBmax;
134 } 134 }
135 135
136 ////////////////////////////////////////////// 136 ////////////////////////////////////////////////////////////////////////
137 // 137 //
138 // Assignment operator 138 // Assignment operator
139 // 139 //
140 G4Tet& G4Tet::operator = (const G4Tet& rhs) 140 G4Tet& G4Tet::operator = (const G4Tet& rhs)
141 { 141 {
142 // Check assignment to self 142 // Check assignment to self
143 // 143 //
144 if (this == &rhs) { return *this; } 144 if (this == &rhs) { return *this; }
145 145
146 // Copy base class data 146 // Copy base class data
147 // 147 //
148 G4VSolid::operator=(rhs); 148 G4VSolid::operator=(rhs);
149 149
150 // Copy data 150 // Copy data
151 // 151 //
152 halfTolerance = rhs.halfTolerance; 152 halfTolerance = rhs.halfTolerance;
153 fCubicVolume = rhs.fCubicVolume; 153 fCubicVolume = rhs.fCubicVolume;
154 fSurfaceArea = rhs.fSurfaceArea; 154 fSurfaceArea = rhs.fSurfaceArea;
155 for (G4int i = 0; i < 4; ++i) { fVertex[i] 155 for (G4int i = 0; i < 4; ++i) { fVertex[i] = rhs.fVertex[i]; }
156 for (G4int i = 0; i < 4; ++i) { fNormal[i] 156 for (G4int i = 0; i < 4; ++i) { fNormal[i] = rhs.fNormal[i]; }
157 for (G4int i = 0; i < 4; ++i) { fDist[i] = 157 for (G4int i = 0; i < 4; ++i) { fDist[i] = rhs.fDist[i]; }
158 for (G4int i = 0; i < 4; ++i) { fArea[i] = 158 for (G4int i = 0; i < 4; ++i) { fArea[i] = rhs.fArea[i]; }
159 fBmin = rhs.fBmin; 159 fBmin = rhs.fBmin;
160 fBmax = rhs.fBmax; 160 fBmax = rhs.fBmax;
161 fRebuildPolyhedron = false; 161 fRebuildPolyhedron = false;
162 delete fpPolyhedron; fpPolyhedron = nullptr 162 delete fpPolyhedron; fpPolyhedron = nullptr;
163 163
164 return *this; 164 return *this;
165 } 165 }
166 166
167 ////////////////////////////////////////////// 167 ////////////////////////////////////////////////////////////////////////
168 // 168 //
169 // Return true if tetrahedron is degenerate 169 // Return true if tetrahedron is degenerate
170 // Tetrahedron is concidered as degenerate in 170 // Tetrahedron is concidered as degenerate in case if its minimal
171 // height is less than degeneracy tolerance 171 // height is less than degeneracy tolerance
172 // 172 //
173 G4bool G4Tet::CheckDegeneracy(const G4ThreeVec 173 G4bool G4Tet::CheckDegeneracy(const G4ThreeVector& p0,
174 const G4ThreeVec 174 const G4ThreeVector& p1,
175 const G4ThreeVec 175 const G4ThreeVector& p2,
176 const G4ThreeVec 176 const G4ThreeVector& p3) const
177 { 177 {
178 G4double hmin = 4. * kCarTolerance; // degen 178 G4double hmin = 4. * kCarTolerance; // degeneracy tolerance
179 179
180 // Calculate volume 180 // Calculate volume
181 G4double vol = std::abs((p1 - p0).cross(p2 - 181 G4double vol = std::abs((p1 - p0).cross(p2 - p0).dot(p3 - p0));
182 182
183 // Calculate face areas squared 183 // Calculate face areas squared
184 G4double ss[4]; 184 G4double ss[4];
185 ss[0] = ((p1 - p0).cross(p2 - p0)).mag2(); 185 ss[0] = ((p1 - p0).cross(p2 - p0)).mag2();
186 ss[1] = ((p2 - p0).cross(p3 - p0)).mag2(); 186 ss[1] = ((p2 - p0).cross(p3 - p0)).mag2();
187 ss[2] = ((p3 - p0).cross(p1 - p0)).mag2(); 187 ss[2] = ((p3 - p0).cross(p1 - p0)).mag2();
188 ss[3] = ((p2 - p1).cross(p3 - p1)).mag2(); 188 ss[3] = ((p2 - p1).cross(p3 - p1)).mag2();
189 189
190 // Find face with max area 190 // Find face with max area
191 G4int k = 0; 191 G4int k = 0;
192 for (G4int i = 1; i < 4; ++i) { if (ss[i] > 192 for (G4int i = 1; i < 4; ++i) { if (ss[i] > ss[k]) k = i; }
193 193
194 // Check: vol^2 / s^2 <= hmin^2 194 // Check: vol^2 / s^2 <= hmin^2
195 return (vol*vol <= ss[k]*hmin*hmin); 195 return (vol*vol <= ss[k]*hmin*hmin);
196 } 196 }
197 197
198 ////////////////////////////////////////////// 198 ////////////////////////////////////////////////////////////////////////
199 // 199 //
200 // Set data members 200 // Set data members
201 // 201 //
202 void G4Tet::Initialize(const G4ThreeVector& p0 202 void G4Tet::Initialize(const G4ThreeVector& p0,
203 const G4ThreeVector& p1 203 const G4ThreeVector& p1,
204 const G4ThreeVector& p2 204 const G4ThreeVector& p2,
205 const G4ThreeVector& p3 205 const G4ThreeVector& p3)
206 { 206 {
207 // Set vertices 207 // Set vertices
208 fVertex[0] = p0; 208 fVertex[0] = p0;
209 fVertex[1] = p1; 209 fVertex[1] = p1;
210 fVertex[2] = p2; 210 fVertex[2] = p2;
211 fVertex[3] = p3; 211 fVertex[3] = p3;
212 212
213 G4ThreeVector norm[4]; 213 G4ThreeVector norm[4];
214 norm[0] = (p2 - p0).cross(p1 - p0); 214 norm[0] = (p2 - p0).cross(p1 - p0);
215 norm[1] = (p3 - p0).cross(p2 - p0); 215 norm[1] = (p3 - p0).cross(p2 - p0);
216 norm[2] = (p1 - p0).cross(p3 - p0); 216 norm[2] = (p1 - p0).cross(p3 - p0);
217 norm[3] = (p2 - p1).cross(p3 - p1); 217 norm[3] = (p2 - p1).cross(p3 - p1);
218 G4double volume = norm[0].dot(p3 - p0); 218 G4double volume = norm[0].dot(p3 - p0);
219 if (volume > 0.) 219 if (volume > 0.)
220 { 220 {
221 for (auto & i : norm) { i = -i; } << 221 for (G4int i = 0; i < 4; ++i) { norm[i] = -norm[i]; }
222 } 222 }
223 223
224 // Set normals to face planes 224 // Set normals to face planes
225 for (G4int i = 0; i < 4; ++i) { fNormal[i] = 225 for (G4int i = 0; i < 4; ++i) { fNormal[i] = norm[i].unit(); }
226 226
227 // Set distances to planes 227 // Set distances to planes
228 for (G4int i = 0; i < 3; ++i) { fDist[i] = f 228 for (G4int i = 0; i < 3; ++i) { fDist[i] = fNormal[i].dot(p0); }
229 fDist[3] = fNormal[3].dot(p1); 229 fDist[3] = fNormal[3].dot(p1);
230 230
231 // Set face areas 231 // Set face areas
232 for (G4int i = 0; i < 4; ++i) { fArea[i] = 0 232 for (G4int i = 0; i < 4; ++i) { fArea[i] = 0.5*norm[i].mag(); }
233 233
234 // Set bounding box 234 // Set bounding box
235 for (G4int i = 0; i < 3; ++i) 235 for (G4int i = 0; i < 3; ++i)
236 { 236 {
237 fBmin[i] = std::min(std::min(std::min(p0[i 237 fBmin[i] = std::min(std::min(std::min(p0[i], p1[i]), p2[i]), p3[i]);
238 fBmax[i] = std::max(std::max(std::max(p0[i 238 fBmax[i] = std::max(std::max(std::max(p0[i], p1[i]), p2[i]), p3[i]);
239 } 239 }
240 240
241 // Set volume and surface area 241 // Set volume and surface area
242 fCubicVolume = std::abs(volume)/6.; 242 fCubicVolume = std::abs(volume)/6.;
243 fSurfaceArea = fArea[0] + fArea[1] + fArea[2 243 fSurfaceArea = fArea[0] + fArea[1] + fArea[2] + fArea[3];
244 } 244 }
245 245
246 ////////////////////////////////////////////// 246 ////////////////////////////////////////////////////////////////////////
247 // 247 //
248 // Set vertices 248 // Set vertices
249 // 249 //
250 void G4Tet::SetVertices(const G4ThreeVector& p 250 void G4Tet::SetVertices(const G4ThreeVector& p0,
251 const G4ThreeVector& p 251 const G4ThreeVector& p1,
252 const G4ThreeVector& p 252 const G4ThreeVector& p2,
253 const G4ThreeVector& p 253 const G4ThreeVector& p3, G4bool* degeneracyFlag)
254 { 254 {
255 // Check for degeneracy 255 // Check for degeneracy
256 G4bool degenerate = CheckDegeneracy(p0, p1, 256 G4bool degenerate = CheckDegeneracy(p0, p1, p2, p3);
257 if (degeneracyFlag != nullptr) << 257 if (degeneracyFlag)
258 { 258 {
259 *degeneracyFlag = degenerate; 259 *degeneracyFlag = degenerate;
260 } 260 }
261 else if (degenerate) 261 else if (degenerate)
262 { 262 {
263 std::ostringstream message; 263 std::ostringstream message;
264 message << "Degenerate tetrahedron is not 264 message << "Degenerate tetrahedron is not permitted: " << GetName() << " !\n"
265 << " anchor: " << p0 << "\n" 265 << " anchor: " << p0 << "\n"
266 << " p1 : " << p1 << "\n" 266 << " p1 : " << p1 << "\n"
267 << " p2 : " << p2 << "\n" 267 << " p2 : " << p2 << "\n"
268 << " p3 : " << p3 << "\n" 268 << " p3 : " << p3 << "\n"
269 << " volume: " 269 << " volume: "
270 << std::abs((p1 - p0).cross(p2 - p 270 << std::abs((p1 - p0).cross(p2 - p0).dot(p3 - p0))/6.;
271 G4Exception("G4Tet::SetVertices()", "GeomS 271 G4Exception("G4Tet::SetVertices()", "GeomSolids0002",
272 FatalException, message); 272 FatalException, message);
273 } 273 }
274 274
275 // Set data members 275 // Set data members
276 Initialize(p0, p1, p2, p3); 276 Initialize(p0, p1, p2, p3);
277 277
278 // Set flag to rebuild polyhedron 278 // Set flag to rebuild polyhedron
279 fRebuildPolyhedron = true; 279 fRebuildPolyhedron = true;
280 } 280 }
281 281
282 ////////////////////////////////////////////// 282 ////////////////////////////////////////////////////////////////////////
283 // 283 //
284 // Return four vertices 284 // Return four vertices
285 // 285 //
286 void G4Tet::GetVertices(G4ThreeVector& p0, 286 void G4Tet::GetVertices(G4ThreeVector& p0,
287 G4ThreeVector& p1, 287 G4ThreeVector& p1,
288 G4ThreeVector& p2, 288 G4ThreeVector& p2,
289 G4ThreeVector& p3) con 289 G4ThreeVector& p3) const
290 { 290 {
291 p0 = fVertex[0]; 291 p0 = fVertex[0];
292 p1 = fVertex[1]; 292 p1 = fVertex[1];
293 p2 = fVertex[2]; 293 p2 = fVertex[2];
294 p3 = fVertex[3]; 294 p3 = fVertex[3];
295 } 295 }
296 296
297 ////////////////////////////////////////////// 297 ////////////////////////////////////////////////////////////////////////
298 // 298 //
299 // Return std::vector of vertices 299 // Return std::vector of vertices
300 // 300 //
301 std::vector<G4ThreeVector> G4Tet::GetVertices( 301 std::vector<G4ThreeVector> G4Tet::GetVertices() const
302 { 302 {
303 std::vector<G4ThreeVector> vertices(4); 303 std::vector<G4ThreeVector> vertices(4);
304 for (G4int i = 0; i < 4; ++i) { vertices[i] 304 for (G4int i = 0; i < 4; ++i) { vertices[i] = fVertex[i]; }
305 return vertices; 305 return vertices;
306 } 306 }
307 307
308 ////////////////////////////////////////////// 308 ////////////////////////////////////////////////////////////////////////
309 // 309 //
310 // Dispatch to parameterisation for replicatio 310 // Dispatch to parameterisation for replication mechanism dimension
311 // computation & modification. 311 // computation & modification.
312 // 312 //
313 void G4Tet::ComputeDimensions(G4VPVParameteris 313 void G4Tet::ComputeDimensions(G4VPVParameterisation* ,
314 const G4int , 314 const G4int ,
315 const G4VPhysica 315 const G4VPhysicalVolume* )
316 { 316 {
317 } 317 }
318 318
319 ////////////////////////////////////////////// 319 ////////////////////////////////////////////////////////////////////////
320 // 320 //
321 // Set bounding box 321 // Set bounding box
322 // 322 //
323 void G4Tet::SetBoundingLimits(const G4ThreeVec 323 void G4Tet::SetBoundingLimits(const G4ThreeVector& pMin,
324 const G4ThreeVec 324 const G4ThreeVector& pMax)
325 { 325 {
326 G4int iout[4] = { 0, 0, 0, 0 }; 326 G4int iout[4] = { 0, 0, 0, 0 };
327 for (G4int i = 0; i < 4; ++i) 327 for (G4int i = 0; i < 4; ++i)
328 { 328 {
329 iout[i] = (G4int)(fVertex[i].x() < pMin.x( << 329 iout[i] = (fVertex[i].x() < pMin.x() ||
330 fVertex[i].y() < pMin.y( << 330 fVertex[i].y() < pMin.y() ||
331 fVertex[i].z() < pMin.z( << 331 fVertex[i].z() < pMin.z() ||
332 fVertex[i].x() > pMax.x( << 332 fVertex[i].x() > pMax.x() ||
333 fVertex[i].y() > pMax.y( << 333 fVertex[i].y() > pMax.y() ||
334 fVertex[i].z() > pMax.z( << 334 fVertex[i].z() > pMax.z());
335 } 335 }
336 if (iout[0] + iout[1] + iout[2] + iout[3] != 336 if (iout[0] + iout[1] + iout[2] + iout[3] != 0)
337 { 337 {
338 std::ostringstream message; 338 std::ostringstream message;
339 message << "Attempt to set bounding box th 339 message << "Attempt to set bounding box that does not encapsulate solid: "
340 << GetName() << " !\n" 340 << GetName() << " !\n"
341 << " Specified bounding box limit 341 << " Specified bounding box limits:\n"
342 << " pmin: " << pMin << "\n" 342 << " pmin: " << pMin << "\n"
343 << " pmax: " << pMax << "\n" 343 << " pmax: " << pMax << "\n"
344 << " Tetrahedron vertices:\n" 344 << " Tetrahedron vertices:\n"
345 << " anchor " << fVertex[0] << << 345 << " anchor " << fVertex[0] << ((iout[0]) ? " is outside\n" : "\n")
346 << " p1 " << fVertex[1] << << 346 << " p1 " << fVertex[1] << ((iout[1]) ? " is outside\n" : "\n")
347 << " p2 " << fVertex[2] << << 347 << " p2 " << fVertex[2] << ((iout[2]) ? " is outside\n" : "\n")
348 << " p3 " << fVertex[3] << << 348 << " p3 " << fVertex[3] << ((iout[3]) ? " is outside" : "");
349 G4Exception("G4Tet::SetBoundingLimits()", 349 G4Exception("G4Tet::SetBoundingLimits()", "GeomSolids0002",
350 FatalException, message); 350 FatalException, message);
351 } 351 }
352 fBmin = pMin; 352 fBmin = pMin;
353 fBmax = pMax; 353 fBmax = pMax;
354 } 354 }
355 355
356 ////////////////////////////////////////////// 356 ////////////////////////////////////////////////////////////////////////
357 // 357 //
358 // Return bounding box 358 // Return bounding box
359 // 359 //
360 void G4Tet::BoundingLimits(G4ThreeVector& pMin 360 void G4Tet::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
361 { 361 {
362 pMin = fBmin; 362 pMin = fBmin;
363 pMax = fBmax; 363 pMax = fBmax;
364 } 364 }
365 365
366 ////////////////////////////////////////////// 366 ////////////////////////////////////////////////////////////////////////
367 // 367 //
368 // Calculate extent under transform and specif 368 // Calculate extent under transform and specified limit
369 // 369 //
370 G4bool G4Tet::CalculateExtent(const EAxis pAxi 370 G4bool G4Tet::CalculateExtent(const EAxis pAxis,
371 const G4VoxelLim 371 const G4VoxelLimits& pVoxelLimit,
372 const G4AffineTr 372 const G4AffineTransform& pTransform,
373 G4double& 373 G4double& pMin, G4double& pMax) const
374 { 374 {
375 G4ThreeVector bmin, bmax; 375 G4ThreeVector bmin, bmax;
376 376
377 // Check bounding box (bbox) 377 // Check bounding box (bbox)
378 // 378 //
379 BoundingLimits(bmin,bmax); 379 BoundingLimits(bmin,bmax);
380 G4BoundingEnvelope bbox(bmin,bmax); 380 G4BoundingEnvelope bbox(bmin,bmax);
381 381
382 // Use simple bounding-box to help in the ca 382 // Use simple bounding-box to help in the case of complex 3D meshes
383 // 383 //
384 return bbox.CalculateExtent(pAxis,pVoxelLimi 384 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
385 385
386 #if 0 386 #if 0
387 // Precise extent computation (disabled by d 387 // Precise extent computation (disabled by default for this shape)
388 // 388 //
389 G4bool exist; 389 G4bool exist;
390 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox 390 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
391 { 391 {
392 return exist = (pMin < pMax) ? true : fals 392 return exist = (pMin < pMax) ? true : false;
393 } 393 }
394 394
395 // Set bounding envelope (benv) and calculat 395 // Set bounding envelope (benv) and calculate extent
396 // 396 //
397 std::vector<G4ThreeVector> vec = GetVertices 397 std::vector<G4ThreeVector> vec = GetVertices();
398 398
399 G4ThreeVectorList anchor(1); 399 G4ThreeVectorList anchor(1);
400 anchor[0].set(vec[0].x(),vec[0].y(),vec[0].z 400 anchor[0].set(vec[0].x(),vec[0].y(),vec[0].z());
401 401
402 G4ThreeVectorList base(3); 402 G4ThreeVectorList base(3);
403 base[0].set(vec[1].x(),vec[1].y(),vec[1].z() 403 base[0].set(vec[1].x(),vec[1].y(),vec[1].z());
404 base[1].set(vec[2].x(),vec[2].y(),vec[2].z() 404 base[1].set(vec[2].x(),vec[2].y(),vec[2].z());
405 base[2].set(vec[3].x(),vec[3].y(),vec[3].z() 405 base[2].set(vec[3].x(),vec[3].y(),vec[3].z());
406 406
407 std::vector<const G4ThreeVectorList *> polyg 407 std::vector<const G4ThreeVectorList *> polygons(2);
408 polygons[0] = &anchor; 408 polygons[0] = &anchor;
409 polygons[1] = &base; 409 polygons[1] = &base;
410 410
411 G4BoundingEnvelope benv(bmin,bmax,polygons); 411 G4BoundingEnvelope benv(bmin,bmax,polygons);
412 return exists = benv.CalculateExtent(pAxis,p 412 return exists = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
413 #endif 413 #endif
414 } 414 }
415 415
416 ////////////////////////////////////////////// 416 ////////////////////////////////////////////////////////////////////////
417 // 417 //
418 // Return whether point inside/outside/on surf 418 // Return whether point inside/outside/on surface
419 // 419 //
420 EInside G4Tet::Inside(const G4ThreeVector& p) 420 EInside G4Tet::Inside(const G4ThreeVector& p) const
421 { 421 {
422 G4double dd[4]; 422 G4double dd[4];
423 for (G4int i = 0; i < 4; ++i) { dd[i] = fNor 423 for (G4int i = 0; i < 4; ++i) { dd[i] = fNormal[i].dot(p) - fDist[i]; }
424 424
425 G4double dist = std::max(std::max(std::max(d 425 G4double dist = std::max(std::max(std::max(dd[0], dd[1]), dd[2]), dd[3]);
426 return (dist <= -halfTolerance) ? 426 return (dist <= -halfTolerance) ?
427 kInside : ((dist <= halfTolerance) ? kSurf 427 kInside : ((dist <= halfTolerance) ? kSurface : kOutside);
428 } 428 }
429 429
430 ////////////////////////////////////////////// 430 ////////////////////////////////////////////////////////////////////////
431 // 431 //
432 // Return unit normal to surface at p 432 // Return unit normal to surface at p
433 // 433 //
434 G4ThreeVector G4Tet::SurfaceNormal( const G4Th 434 G4ThreeVector G4Tet::SurfaceNormal( const G4ThreeVector& p) const
435 { 435 {
436 G4double k[4]; 436 G4double k[4];
437 for (G4int i = 0; i < 4; ++i) 437 for (G4int i = 0; i < 4; ++i)
438 { 438 {
439 k[i] = (G4double)(std::abs(fNormal[i].dot( << 439 k[i] = std::abs(fNormal[i].dot(p) - fDist[i]) <= halfTolerance;
440 } 440 }
441 G4double nsurf = k[0] + k[1] + k[2] + k[3]; 441 G4double nsurf = k[0] + k[1] + k[2] + k[3];
442 G4ThreeVector norm = 442 G4ThreeVector norm =
443 k[0]*fNormal[0] + k[1]*fNormal[1] + k[2]*f 443 k[0]*fNormal[0] + k[1]*fNormal[1] + k[2]*fNormal[2] + k[3]*fNormal[3];
444 444
445 if (nsurf == 1.) return norm; 445 if (nsurf == 1.) return norm;
446 else if (nsurf > 1.) return norm.unit(); // 446 else if (nsurf > 1.) return norm.unit(); // edge or vertex
447 { 447 {
448 #ifdef G4SPECSDEBUG 448 #ifdef G4SPECSDEBUG
449 std::ostringstream message; 449 std::ostringstream message;
450 G4long oldprc = message.precision(16); 450 G4long oldprc = message.precision(16);
451 message << "Point p is not on surface (!?) 451 message << "Point p is not on surface (!?) of solid: "
452 << GetName() << "\n"; 452 << GetName() << "\n";
453 message << "Position:\n"; 453 message << "Position:\n";
454 message << " p.x() = " << p.x()/mm << " 454 message << " p.x() = " << p.x()/mm << " mm\n";
455 message << " p.y() = " << p.y()/mm << " 455 message << " p.y() = " << p.y()/mm << " mm\n";
456 message << " p.z() = " << p.z()/mm << " 456 message << " p.z() = " << p.z()/mm << " mm";
457 G4cout.precision(oldprc); 457 G4cout.precision(oldprc);
458 G4Exception("G4Tet::SurfaceNormal(p)", "Ge 458 G4Exception("G4Tet::SurfaceNormal(p)", "GeomSolids1002",
459 JustWarning, message ); 459 JustWarning, message );
460 DumpInfo(); 460 DumpInfo();
461 #endif 461 #endif
462 return ApproxSurfaceNormal(p); 462 return ApproxSurfaceNormal(p);
463 } 463 }
464 } 464 }
465 465
466 ////////////////////////////////////////////// 466 ////////////////////////////////////////////////////////////////////////
467 // 467 //
468 // Find surface nearest to point and return co 468 // Find surface nearest to point and return corresponding normal
469 // This method normally should not be called 469 // This method normally should not be called
470 // 470 //
471 G4ThreeVector G4Tet::ApproxSurfaceNormal(const 471 G4ThreeVector G4Tet::ApproxSurfaceNormal(const G4ThreeVector& p) const
472 { 472 {
473 G4double dist = -DBL_MAX; 473 G4double dist = -DBL_MAX;
474 G4int iside = 0; 474 G4int iside = 0;
475 for (G4int i = 0; i < 4; ++i) 475 for (G4int i = 0; i < 4; ++i)
476 { 476 {
477 G4double d = fNormal[i].dot(p) - fDist[i]; 477 G4double d = fNormal[i].dot(p) - fDist[i];
478 if (d > dist) { dist = d; iside = i; } 478 if (d > dist) { dist = d; iside = i; }
479 } 479 }
480 return fNormal[iside]; 480 return fNormal[iside];
481 } 481 }
482 482
483 ////////////////////////////////////////////// 483 ////////////////////////////////////////////////////////////////////////
484 // 484 //
485 // Calculate distance to surface from outside, 485 // Calculate distance to surface from outside,
486 // return kInfinity if no intersection 486 // return kInfinity if no intersection
487 // 487 //
488 G4double G4Tet::DistanceToIn(const G4ThreeVect 488 G4double G4Tet::DistanceToIn(const G4ThreeVector& p,
489 const G4ThreeVect 489 const G4ThreeVector& v) const
490 { 490 {
491 G4double tin = -DBL_MAX, tout = DBL_MAX; 491 G4double tin = -DBL_MAX, tout = DBL_MAX;
492 for (G4int i = 0; i < 4; ++i) 492 for (G4int i = 0; i < 4; ++i)
493 { 493 {
494 G4double cosa = fNormal[i].dot(v); 494 G4double cosa = fNormal[i].dot(v);
495 G4double dist = fNormal[i].dot(p) - fDist[ 495 G4double dist = fNormal[i].dot(p) - fDist[i];
496 if (dist >= -halfTolerance) 496 if (dist >= -halfTolerance)
497 { 497 {
498 if (cosa >= 0.) { return kInfinity; } 498 if (cosa >= 0.) { return kInfinity; }
499 tin = std::max(tin, -dist/cosa); 499 tin = std::max(tin, -dist/cosa);
500 } 500 }
501 else if (cosa > 0.) 501 else if (cosa > 0.)
502 { 502 {
503 tout = std::min(tout, -dist/cosa); 503 tout = std::min(tout, -dist/cosa);
504 } 504 }
505 } 505 }
506 506
507 return (tout - tin <= halfTolerance) ? 507 return (tout - tin <= halfTolerance) ?
508 kInfinity : ((tin < halfTolerance) ? 0. : 508 kInfinity : ((tin < halfTolerance) ? 0. : tin);
509 } 509 }
510 510
511 ////////////////////////////////////////////// 511 ////////////////////////////////////////////////////////////////////////
512 // 512 //
513 // Estimate safety distance to surface from ou 513 // Estimate safety distance to surface from outside
514 // 514 //
515 G4double G4Tet::DistanceToIn(const G4ThreeVect 515 G4double G4Tet::DistanceToIn(const G4ThreeVector& p) const
516 { 516 {
517 G4double dd[4]; 517 G4double dd[4];
518 for (G4int i = 0; i < 4; ++i) { dd[i] = fNor 518 for (G4int i = 0; i < 4; ++i) { dd[i] = fNormal[i].dot(p) - fDist[i]; }
519 519
520 G4double dist = std::max(std::max(std::max(d 520 G4double dist = std::max(std::max(std::max(dd[0], dd[1]), dd[2]), dd[3]);
521 return (dist > 0.) ? dist : 0.; 521 return (dist > 0.) ? dist : 0.;
522 } 522 }
523 523
524 ////////////////////////////////////////////// 524 ////////////////////////////////////////////////////////////////////////
525 // 525 //
526 // Calcluate distance to surface from inside 526 // Calcluate distance to surface from inside
527 // 527 //
528 G4double G4Tet::DistanceToOut(const G4ThreeVec 528 G4double G4Tet::DistanceToOut(const G4ThreeVector& p,
529 const G4ThreeVec 529 const G4ThreeVector& v,
530 const G4bool cal 530 const G4bool calcNorm,
531 G4bool* va 531 G4bool* validNorm,
532 G4ThreeVec 532 G4ThreeVector* n) const
533 { 533 {
534 // Calculate distances and cosines 534 // Calculate distances and cosines
535 G4double cosa[4], dist[4]; 535 G4double cosa[4], dist[4];
536 G4int ind[4] = {0}, nside = 0; 536 G4int ind[4] = {0}, nside = 0;
537 for (G4int i = 0; i < 4; ++i) 537 for (G4int i = 0; i < 4; ++i)
538 { 538 {
539 G4double tmp = fNormal[i].dot(v); 539 G4double tmp = fNormal[i].dot(v);
540 cosa[i] = tmp; 540 cosa[i] = tmp;
541 ind[nside] = (G4int)(tmp > 0) * i; << 541 ind[nside] = (tmp > 0) * i;
542 nside += (G4int)(tmp > 0); << 542 nside += (tmp > 0);
543 dist[i] = fNormal[i].dot(p) - fDist[i]; 543 dist[i] = fNormal[i].dot(p) - fDist[i];
544 } 544 }
545 545
546 // Find intersection (in most of cases nside 546 // Find intersection (in most of cases nside == 1)
547 G4double tout = DBL_MAX; 547 G4double tout = DBL_MAX;
548 G4int iside = 0; 548 G4int iside = 0;
549 for (G4int i = 0; i < nside; ++i) 549 for (G4int i = 0; i < nside; ++i)
550 { 550 {
551 G4int k = ind[i]; 551 G4int k = ind[i];
552 // Check: leaving the surface 552 // Check: leaving the surface
553 if (dist[k] >= -halfTolerance) { tout = 0. 553 if (dist[k] >= -halfTolerance) { tout = 0.; iside = k; break; }
554 // Compute distance to intersection 554 // Compute distance to intersection
555 G4double tmp = -dist[k]/cosa[k]; 555 G4double tmp = -dist[k]/cosa[k];
556 if (tmp < tout) { tout = tmp; iside = k; } 556 if (tmp < tout) { tout = tmp; iside = k; }
557 } 557 }
558 558
559 // Set normal, if required, and return dista 559 // Set normal, if required, and return distance to out
560 if (calcNorm) 560 if (calcNorm)
561 { 561 {
562 *validNorm = true; 562 *validNorm = true;
563 *n = fNormal[iside]; 563 *n = fNormal[iside];
564 } 564 }
565 return tout; 565 return tout;
566 } 566 }
567 567
568 ////////////////////////////////////////////// 568 ////////////////////////////////////////////////////////////////////////
569 // 569 //
570 // Calculate safety distance to surface from i 570 // Calculate safety distance to surface from inside
571 // 571 //
572 G4double G4Tet::DistanceToOut(const G4ThreeVec 572 G4double G4Tet::DistanceToOut(const G4ThreeVector& p) const
573 { 573 {
574 G4double dd[4]; 574 G4double dd[4];
575 for (G4int i = 0; i < 4; ++i) { dd[i] = fDis 575 for (G4int i = 0; i < 4; ++i) { dd[i] = fDist[i] - fNormal[i].dot(p); }
576 576
577 G4double dist = std::min(std::min(std::min(d 577 G4double dist = std::min(std::min(std::min(dd[0], dd[1]), dd[2]), dd[3]);
578 return (dist > 0.) ? dist : 0.; 578 return (dist > 0.) ? dist : 0.;
579 } 579 }
580 580
581 ////////////////////////////////////////////// 581 ////////////////////////////////////////////////////////////////////////
582 // 582 //
583 // GetEntityType 583 // GetEntityType
584 // 584 //
585 G4GeometryType G4Tet::GetEntityType() const 585 G4GeometryType G4Tet::GetEntityType() const
586 { 586 {
587 return {"G4Tet"}; << 587 return G4String("G4Tet");
588 } <<
589 <<
590 ////////////////////////////////////////////// <<
591 // <<
592 // IsFaceted <<
593 // <<
594 G4bool G4Tet::IsFaceted() const <<
595 { <<
596 return true; <<
597 } 588 }
598 589
599 ////////////////////////////////////////////// 590 ////////////////////////////////////////////////////////////////////////
600 // 591 //
601 // Make a clone of the object 592 // Make a clone of the object
602 // 593 //
603 G4VSolid* G4Tet::Clone() const 594 G4VSolid* G4Tet::Clone() const
604 { 595 {
605 return new G4Tet(*this); 596 return new G4Tet(*this);
606 } 597 }
607 598
608 ////////////////////////////////////////////// 599 ////////////////////////////////////////////////////////////////////////
609 // 600 //
610 // Stream object contents to an output stream 601 // Stream object contents to an output stream
611 // 602 //
612 std::ostream& G4Tet::StreamInfo(std::ostream& 603 std::ostream& G4Tet::StreamInfo(std::ostream& os) const
613 { 604 {
614 G4long oldprc = os.precision(16); 605 G4long oldprc = os.precision(16);
615 os << "------------------------------------- 606 os << "-----------------------------------------------------------\n"
616 << " *** Dump for solid - " << GetName 607 << " *** Dump for solid - " << GetName() << " ***\n"
617 << " ================================= 608 << " ===================================================\n"
618 << " Solid type: " << GetEntityType() << 609 << " Solid type: " << GetEntityType() << "\n"
619 << " Parameters: \n" 610 << " Parameters: \n"
620 << " anchor: " << fVertex[0]/mm << " m 611 << " anchor: " << fVertex[0]/mm << " mm\n"
621 << " p1 : " << fVertex[1]/mm << " m 612 << " p1 : " << fVertex[1]/mm << " mm\n"
622 << " p2 : " << fVertex[2]/mm << " m 613 << " p2 : " << fVertex[2]/mm << " mm\n"
623 << " p3 : " << fVertex[3]/mm << " m 614 << " p3 : " << fVertex[3]/mm << " mm\n"
624 << "------------------------------------- 615 << "-----------------------------------------------------------\n";
625 os.precision(oldprc); 616 os.precision(oldprc);
626 return os; 617 return os;
627 } 618 }
628 619
629 ////////////////////////////////////////////// 620 ////////////////////////////////////////////////////////////////////////
630 // 621 //
631 // Return random point on the surface 622 // Return random point on the surface
632 // 623 //
633 G4ThreeVector G4Tet::GetPointOnSurface() const 624 G4ThreeVector G4Tet::GetPointOnSurface() const
634 { 625 {
635 constexpr G4int iface[4][3] = { {0,1,2}, {0, 626 constexpr G4int iface[4][3] = { {0,1,2}, {0,2,3}, {0,3,1}, {1,2,3} };
636 627
637 // Select face 628 // Select face
638 G4double select = fSurfaceArea*G4QuickRand() 629 G4double select = fSurfaceArea*G4QuickRand();
639 G4int i = 0; 630 G4int i = 0;
640 i += (G4int)(select > fArea[0]); << 631 for ( ; i < 4; ++i) { if ((select -= fArea[i]) <= 0.) break; }
641 i += (G4int)(select > fArea[0] + fArea[1]); <<
642 i += (G4int)(select > fArea[0] + fArea[1] + <<
643 632
644 // Set selected triangle 633 // Set selected triangle
645 G4ThreeVector p0 = fVertex[iface[i][0]]; 634 G4ThreeVector p0 = fVertex[iface[i][0]];
646 G4ThreeVector e1 = fVertex[iface[i][1]] - p0 635 G4ThreeVector e1 = fVertex[iface[i][1]] - p0;
647 G4ThreeVector e2 = fVertex[iface[i][2]] - p0 636 G4ThreeVector e2 = fVertex[iface[i][2]] - p0;
648 637
649 // Return random point 638 // Return random point
650 G4double r1 = G4QuickRand(); 639 G4double r1 = G4QuickRand();
651 G4double r2 = G4QuickRand(); 640 G4double r2 = G4QuickRand();
652 return (r1 + r2 > 1.) ? 641 return (r1 + r2 > 1.) ?
653 p0 + e1*(1. - r1) + e2*(1. - r2) : p0 + e1 642 p0 + e1*(1. - r1) + e2*(1. - r2) : p0 + e1*r1 + e2*r2;
654 } 643 }
655 644
656 ////////////////////////////////////////////// 645 ////////////////////////////////////////////////////////////////////////
657 // 646 //
658 // Return volume of the tetrahedron 647 // Return volume of the tetrahedron
659 // 648 //
660 G4double G4Tet::GetCubicVolume() 649 G4double G4Tet::GetCubicVolume()
661 { 650 {
662 return fCubicVolume; 651 return fCubicVolume;
663 } 652 }
664 653
665 ////////////////////////////////////////////// 654 ////////////////////////////////////////////////////////////////////////
666 // 655 //
667 // Return surface area of the tetrahedron 656 // Return surface area of the tetrahedron
668 // 657 //
669 G4double G4Tet::GetSurfaceArea() 658 G4double G4Tet::GetSurfaceArea()
670 { 659 {
671 return fSurfaceArea; 660 return fSurfaceArea;
672 } 661 }
673 662
674 ////////////////////////////////////////////// 663 ////////////////////////////////////////////////////////////////////////
675 // 664 //
676 // Methods for visualisation 665 // Methods for visualisation
677 // 666 //
678 void G4Tet::DescribeYourselfTo (G4VGraphicsSce 667 void G4Tet::DescribeYourselfTo (G4VGraphicsScene& scene) const
679 { 668 {
680 scene.AddSolid (*this); 669 scene.AddSolid (*this);
681 } 670 }
682 671
683 ////////////////////////////////////////////// 672 ////////////////////////////////////////////////////////////////////////
684 // 673 //
685 // Return VisExtent 674 // Return VisExtent
686 // 675 //
687 G4VisExtent G4Tet::GetExtent() const 676 G4VisExtent G4Tet::GetExtent() const
688 { 677 {
689 return { fBmin.x(), fBmax.x(), << 678 return G4VisExtent(fBmin.x(), fBmax.x(),
690 fBmin.y(), fBmax.y(), << 679 fBmin.y(), fBmax.y(),
691 fBmin.z(), fBmax.z() }; << 680 fBmin.z(), fBmax.z());
692 } 681 }
693 682
694 ////////////////////////////////////////////// 683 ////////////////////////////////////////////////////////////////////////
695 // 684 //
696 // CreatePolyhedron 685 // CreatePolyhedron
697 // 686 //
698 G4Polyhedron* G4Tet::CreatePolyhedron() const 687 G4Polyhedron* G4Tet::CreatePolyhedron() const
699 { 688 {
700 // Check orientation of vertices 689 // Check orientation of vertices
701 G4ThreeVector v1 = fVertex[1] - fVertex[0]; 690 G4ThreeVector v1 = fVertex[1] - fVertex[0];
702 G4ThreeVector v2 = fVertex[2] - fVertex[0]; 691 G4ThreeVector v2 = fVertex[2] - fVertex[0];
703 G4ThreeVector v3 = fVertex[3] - fVertex[0]; 692 G4ThreeVector v3 = fVertex[3] - fVertex[0];
704 G4bool invert = v1.cross(v2).dot(v3) < 0.; 693 G4bool invert = v1.cross(v2).dot(v3) < 0.;
705 G4int k2 = (invert) ? 3 : 2; 694 G4int k2 = (invert) ? 3 : 2;
706 G4int k3 = (invert) ? 2 : 3; 695 G4int k3 = (invert) ? 2 : 3;
707 696
708 // Set coordinates of vertices 697 // Set coordinates of vertices
709 G4double xyz[4][3]; 698 G4double xyz[4][3];
710 for (G4int i = 0; i < 3; ++i) 699 for (G4int i = 0; i < 3; ++i)
711 { 700 {
712 xyz[0][i] = fVertex[0][i]; 701 xyz[0][i] = fVertex[0][i];
713 xyz[1][i] = fVertex[1][i]; 702 xyz[1][i] = fVertex[1][i];
714 xyz[2][i] = fVertex[k2][i]; 703 xyz[2][i] = fVertex[k2][i];
715 xyz[3][i] = fVertex[k3][i]; 704 xyz[3][i] = fVertex[k3][i];
716 } 705 }
717 706
718 // Create polyhedron 707 // Create polyhedron
719 G4int faces[4][4] = { {1,3,2,0}, {1,4,3,0}, 708 G4int faces[4][4] = { {1,3,2,0}, {1,4,3,0}, {1,2,4,0}, {2,3,4,0} };
720 auto ph = new G4Polyhedron; << 709 G4Polyhedron* ph = new G4Polyhedron;
721 ph->createPolyhedron(4,4,xyz,faces); 710 ph->createPolyhedron(4,4,xyz,faces);
722 711
723 return ph; 712 return ph;
724 } 713 }
725 714
726 ////////////////////////////////////////////// 715 ////////////////////////////////////////////////////////////////////////
727 // 716 //
728 // GetPolyhedron 717 // GetPolyhedron
729 // 718 //
730 G4Polyhedron* G4Tet::GetPolyhedron() const 719 G4Polyhedron* G4Tet::GetPolyhedron() const
731 { 720 {
732 if (fpPolyhedron == nullptr || 721 if (fpPolyhedron == nullptr ||
733 fRebuildPolyhedron || 722 fRebuildPolyhedron ||
734 fpPolyhedron->GetNumberOfRotationStepsAt 723 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
735 fpPolyhedron->GetNumberOfRotationSteps() 724 fpPolyhedron->GetNumberOfRotationSteps())
736 { 725 {
737 G4AutoLock l(&polyhedronMutex); 726 G4AutoLock l(&polyhedronMutex);
738 delete fpPolyhedron; 727 delete fpPolyhedron;
739 fpPolyhedron = CreatePolyhedron(); 728 fpPolyhedron = CreatePolyhedron();
740 fRebuildPolyhedron = false; 729 fRebuildPolyhedron = false;
741 l.unlock(); 730 l.unlock();
742 } 731 }
743 return fpPolyhedron; 732 return fpPolyhedron;
744 } 733 }
745 734
746 #endif 735 #endif
747 736