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