| 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: G4Trd.cc,v 1.38 2010-10-19 15:42:10 gcosmo Exp $
>> 28 // GEANT4 tag $Name: not supported by cvs2svn $
>> 29 //
>> 30 //
26 // Implementation for G4Trd class 31 // Implementation for G4Trd class
27 // 32 //
28 // 12.01.95 P.Kent: First version << 33 // History:
29 // 28.04.05 V.Grichine: new SurfaceNormal acco << 34 //
30 // 25.05.17 E.Tcherniaev: complete revision, s << 35 // 28.04.05 V.Grichine: new SurfaceNormal according to J. Apostolakis proposal
31 // ------------------------------------------- << 36 // 26.04.05, V.Grichine, new SurfaceNoramal is default
>> 37 // 07.12.04, V.Grichine, SurfaceNoramal with edges/vertices.
>> 38 // 07.05.00, V.Grichine, in d = DistanceToIn(p,v), if d<0.5*kCarTolerance, d=0
>> 39 // ~1996, V.Grichine, 1st implementation based on old code of P.Kent
>> 40 //
>> 41 //////////////////////////////////////////////////////////////////////////////
32 42
33 #include "G4Trd.hh" 43 #include "G4Trd.hh"
34 44
35 #if !defined(G4GEOM_USE_UTRD) << 45 #include "G4VPVParameterisation.hh"
36 <<
37 #include "G4GeomTools.hh" <<
38 <<
39 #include "G4VoxelLimits.hh" 46 #include "G4VoxelLimits.hh"
40 #include "G4AffineTransform.hh" 47 #include "G4AffineTransform.hh"
41 #include "G4BoundingEnvelope.hh" << 48 #include "Randomize.hh"
42 #include "G4QuickRand.hh" <<
43 <<
44 #include "G4VPVParameterisation.hh" <<
45 49
46 #include "G4VGraphicsScene.hh" 50 #include "G4VGraphicsScene.hh"
>> 51 #include "G4Polyhedron.hh"
>> 52 #include "G4NURBS.hh"
>> 53 #include "G4NURBSbox.hh"
47 54
48 using namespace CLHEP; 55 using namespace CLHEP;
49 56
50 ////////////////////////////////////////////// << 57 /////////////////////////////////////////////////////////////////////////
51 // 58 //
52 // Constructor - set & check half widths << 59 // Constructor - check & set half widths
53 60
54 G4Trd::G4Trd(const G4String& pName, << 61 G4Trd::G4Trd( const G4String& pName,
55 G4double pdx1, G4double pdx << 62 G4double pdx1, G4double pdx2,
56 G4double pdy1, G4double pdy << 63 G4double pdy1, G4double pdy2,
57 G4double pdz) << 64 G4double pdz )
58 : G4CSGSolid(pName), halfCarTolerance(0.5*kC << 65 : G4CSGSolid(pName)
59 fDx1(pdx1), fDx2(pdx2), fDy1(pdy1), fDy2(p <<
60 { 66 {
61 CheckParameters(); << 67 CheckAndSetAllParameters (pdx1, pdx2, pdy1, pdy2, pdz);
62 MakePlanes(); <<
63 } 68 }
64 69
65 ////////////////////////////////////////////// << 70 /////////////////////////////////////////////////////////////////////////
>> 71 //
>> 72 // Set and check (coplanarity) of trd parameters
>> 73
>> 74 void G4Trd::CheckAndSetAllParameters ( G4double pdx1, G4double pdx2,
>> 75 G4double pdy1, G4double pdy2,
>> 76 G4double pdz )
>> 77 {
>> 78 if ( pdx1>0&&pdx2>0&&pdy1>0&&pdy2>0&&pdz>0 )
>> 79 {
>> 80 fDx1=pdx1; fDx2=pdx2;
>> 81 fDy1=pdy1; fDy2=pdy2;
>> 82 fDz=pdz;
>> 83 }
>> 84 else
>> 85 {
>> 86 if ( pdx1>=0 && pdx2>=0 && pdy1>=0 && pdy2>=0 && pdz>=0 )
>> 87 {
>> 88 // G4double Minimum_length= (1+per_thousand) * kCarTolerance/2.;
>> 89 // FIX-ME : temporary solution for ZERO or very-small parameters
>> 90 //
>> 91 G4double Minimum_length= kCarTolerance/2.;
>> 92 fDx1=std::max(pdx1,Minimum_length);
>> 93 fDx2=std::max(pdx2,Minimum_length);
>> 94 fDy1=std::max(pdy1,Minimum_length);
>> 95 fDy2=std::max(pdy2,Minimum_length);
>> 96 fDz=std::max(pdz,Minimum_length);
>> 97 }
>> 98 else
>> 99 {
>> 100 std::ostringstream message;
>> 101 message << "Invalid negative dimensions for Solid: " << GetName()
>> 102 << G4endl
>> 103 << " X - " << pdx1 << ", " << pdx2 << G4endl
>> 104 << " Y - " << pdy1 << ", " << pdy2 << G4endl
>> 105 << " Z - " << pdz;
>> 106 G4Exception("G4Trd::CheckAndSetAllParameters()",
>> 107 "GeomSolids0002", FatalException, message);
>> 108 }
>> 109 }
>> 110 fCubicVolume= 0.;
>> 111 fSurfaceArea= 0.;
>> 112 fpPolyhedron = 0;
>> 113 }
>> 114
>> 115 ///////////////////////////////////////////////////////////////////////
66 // 116 //
67 // Fake default constructor - sets only member 117 // Fake default constructor - sets only member data and allocates memory
68 // for usage restri << 118 // for usage restricted to object persistency.
69 // 119 //
70 G4Trd::G4Trd( __void__& a ) 120 G4Trd::G4Trd( __void__& a )
71 : G4CSGSolid(a), halfCarTolerance(0.5*kCarTo << 121 : G4CSGSolid(a), fDx1(0.), fDx2(0.), fDy1(0.), fDy2(0.), fDz(0.)
72 fDx1(1.), fDx2(1.), fDy1(1.), fDy2(1.), fD <<
73 { 122 {
74 MakePlanes(); <<
75 } 123 }
76 124
77 ////////////////////////////////////////////// 125 //////////////////////////////////////////////////////////////////////////
78 // 126 //
79 // Destructor 127 // Destructor
80 128
81 G4Trd::~G4Trd() = default; << 129 G4Trd::~G4Trd()
>> 130 {
>> 131 }
82 132
83 ////////////////////////////////////////////// 133 //////////////////////////////////////////////////////////////////////////
84 // 134 //
85 // Copy constructor 135 // Copy constructor
86 136
87 G4Trd::G4Trd(const G4Trd& rhs) 137 G4Trd::G4Trd(const G4Trd& rhs)
88 : G4CSGSolid(rhs), halfCarTolerance(rhs.half << 138 : G4CSGSolid(rhs), fDx1(rhs.fDx1), fDx2(rhs.fDx2),
89 fDx1(rhs.fDx1), fDx2(rhs.fDx2), << 139 fDy1(rhs.fDy1), fDy2(rhs.fDy2), fDz(rhs.fDz)
90 fDy1(rhs.fDy1), fDy2(rhs.fDy2), fDz(rhs.fD <<
91 fHx(rhs.fHx), fHy(rhs.fHy) <<
92 { 140 {
93 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs <<
94 } 141 }
95 142
96 ////////////////////////////////////////////// 143 //////////////////////////////////////////////////////////////////////////
97 // 144 //
98 // Assignment operator 145 // Assignment operator
99 146
100 G4Trd& G4Trd::operator = (const G4Trd& rhs) << 147 G4Trd& G4Trd::operator = (const G4Trd& rhs)
101 { 148 {
102 // Check assignment to self 149 // Check assignment to self
103 // 150 //
104 if (this == &rhs) { return *this; } 151 if (this == &rhs) { return *this; }
105 152
106 // Copy base class data 153 // Copy base class data
107 // 154 //
108 G4CSGSolid::operator=(rhs); 155 G4CSGSolid::operator=(rhs);
109 156
110 // Copy data 157 // Copy data
111 // 158 //
112 halfCarTolerance = rhs.halfCarTolerance; <<
113 fDx1 = rhs.fDx1; fDx2 = rhs.fDx2; 159 fDx1 = rhs.fDx1; fDx2 = rhs.fDx2;
114 fDy1 = rhs.fDy1; fDy2 = rhs.fDy2; 160 fDy1 = rhs.fDy1; fDy2 = rhs.fDy2;
115 fDz = rhs.fDz; 161 fDz = rhs.fDz;
116 fHx = rhs.fHx; fHy = rhs.fHy; <<
117 for (G4int i=0; i<4; ++i) { fPlanes[i] = rh <<
118 162
119 return *this; 163 return *this;
120 } 164 }
121 165
122 ////////////////////////////////////////////// << 166 ////////////////////////////////////////////////////////////////////////////
123 // 167 //
124 // Set all parameters, as for constructor - se <<
125 <<
126 void G4Trd::SetAllParameters(G4double pdx1, G4 <<
127 G4double pdy1, G4 <<
128 { <<
129 // Reset data of the base class <<
130 fCubicVolume = 0.; <<
131 fSurfaceArea = 0.; <<
132 fRebuildPolyhedron = true; <<
133 <<
134 // Set parameters <<
135 fDx1 = pdx1; fDx2 = pdx2; <<
136 fDy1 = pdy1; fDy2 = pdy2; <<
137 fDz = pdz; <<
138 <<
139 CheckParameters(); <<
140 MakePlanes(); <<
141 } <<
142 <<
143 ////////////////////////////////////////////// <<
144 // 168 //
145 // Check dimensions <<
146 169
147 void G4Trd::CheckParameters() << 170 void G4Trd::SetAllParameters ( G4double pdx1, G4double pdx2, G4double pdy1,
>> 171 G4double pdy2, G4double pdz )
148 { 172 {
149 G4double dmin = 2*kCarTolerance; << 173 CheckAndSetAllParameters (pdx1, pdx2, pdy1, pdy2, pdz);
150 if ((fDx1 < 0 || fDx2 < 0 || fDy1 < 0 || fDy <<
151 (fDx1 < dmin && fDx2 < dmin) || <<
152 (fDy1 < dmin && fDy2 < dmin)) <<
153 { <<
154 std::ostringstream message; <<
155 message << "Invalid (too small or negative <<
156 << GetName() <<
157 << "\n X - " << fDx1 << ", " << f <<
158 << "\n Y - " << fDy1 << ", " << f <<
159 << "\n Z - " << fDz; <<
160 G4Exception("G4Trd::CheckParameters()", "G <<
161 FatalException, message); <<
162 } <<
163 } 174 }
164 175
165 ////////////////////////////////////////////// <<
166 // <<
167 // Set side planes <<
168 <<
169 void G4Trd::MakePlanes() <<
170 { <<
171 G4double dx = fDx1 - fDx2; <<
172 G4double dy = fDy1 - fDy2; <<
173 G4double dz = 2*fDz; <<
174 fHx = std::sqrt(dy*dy + dz*dz); <<
175 fHy = std::sqrt(dx*dx + dz*dz); <<
176 176
177 // Set X planes at -Y & +Y << 177 /////////////////////////////////////////////////////////////////////////
178 // <<
179 fPlanes[0].a = 0.; <<
180 fPlanes[0].b = -dz/fHx; <<
181 fPlanes[0].c = dy/fHx; <<
182 fPlanes[0].d = fPlanes[0].b*fDy1 + fPlanes[0 <<
183 <<
184 fPlanes[1].a = fPlanes[0].a; <<
185 fPlanes[1].b = -fPlanes[0].b; <<
186 fPlanes[1].c = fPlanes[0].c; <<
187 fPlanes[1].d = fPlanes[0].d; <<
188 <<
189 // Set Y planes at -X & +X <<
190 // <<
191 fPlanes[2].a = -dz/fHy; <<
192 fPlanes[2].b = 0.; <<
193 fPlanes[2].c = dx/fHy; <<
194 fPlanes[2].d = fPlanes[2].a*fDx1 + fPlanes[2 <<
195 <<
196 fPlanes[3].a = -fPlanes[2].a; <<
197 fPlanes[3].b = fPlanes[2].b; <<
198 fPlanes[3].c = fPlanes[2].c; <<
199 fPlanes[3].d = fPlanes[2].d; <<
200 } <<
201 <<
202 ////////////////////////////////////////////// <<
203 // <<
204 // Get volume <<
205 <<
206 G4double G4Trd::GetCubicVolume() <<
207 { <<
208 if (fCubicVolume == 0.) <<
209 { <<
210 fCubicVolume = 2*fDz*( (fDx1+fDx2)*(fDy1+f <<
211 (fDx2-fDx1)*(fDy2-f <<
212 } <<
213 return fCubicVolume; <<
214 } <<
215 <<
216 ////////////////////////////////////////////// <<
217 // <<
218 // Get surface area <<
219 <<
220 G4double G4Trd::GetSurfaceArea() <<
221 { <<
222 if (fSurfaceArea == 0.) <<
223 { <<
224 fSurfaceArea = <<
225 4*(fDx1*fDy1 + fDx2*fDy2) + 2*(fDx1+fDx2 <<
226 } <<
227 return fSurfaceArea; <<
228 } <<
229 <<
230 ////////////////////////////////////////////// <<
231 // 178 //
232 // Dispatch to parameterisation for replicatio 179 // Dispatch to parameterisation for replication mechanism dimension
233 // computation & modification << 180 // computation & modification.
234 181
235 void G4Trd::ComputeDimensions( G4VPVPara 182 void G4Trd::ComputeDimensions( G4VPVParameterisation* p,
236 const G4int n, 183 const G4int n,
237 const G4VPhysic 184 const G4VPhysicalVolume* pRep )
238 { 185 {
239 p->ComputeDimensions(*this,n,pRep); 186 p->ComputeDimensions(*this,n,pRep);
240 } 187 }
241 188
242 ////////////////////////////////////////////// <<
243 // <<
244 // Get bounding box <<
245 189
246 void G4Trd::BoundingLimits(G4ThreeVector& pMin << 190 ///////////////////////////////////////////////////////////////////////////
247 { <<
248 G4double dx1 = GetXHalfLength1(); <<
249 G4double dx2 = GetXHalfLength2(); <<
250 G4double dy1 = GetYHalfLength1(); <<
251 G4double dy2 = GetYHalfLength2(); <<
252 G4double dz = GetZHalfLength(); <<
253 <<
254 G4double xmax = std::max(dx1,dx2); <<
255 G4double ymax = std::max(dy1,dy2); <<
256 pMin.set(-xmax,-ymax,-dz); <<
257 pMax.set( xmax, ymax, dz); <<
258 <<
259 // Check correctness of the bounding box <<
260 // <<
261 if (pMin.x() >= pMax.x() || pMin.y() >= pMax <<
262 { <<
263 std::ostringstream message; <<
264 message << "Bad bounding box (min >= max) <<
265 << GetName() << " !" <<
266 << "\npMin = " << pMin <<
267 << "\npMax = " << pMax; <<
268 G4Exception("G4Trd::BoundingLimits()", "Ge <<
269 DumpInfo(); <<
270 } <<
271 } <<
272 <<
273 ////////////////////////////////////////////// <<
274 // 191 //
275 // Calculate extent under transform and specif 192 // Calculate extent under transform and specified limit
276 193
277 G4bool G4Trd::CalculateExtent( const EAxis pAx 194 G4bool G4Trd::CalculateExtent( const EAxis pAxis,
278 const G4VoxelLi 195 const G4VoxelLimits& pVoxelLimit,
279 const G4AffineT 196 const G4AffineTransform& pTransform,
280 G4double& 197 G4double& pMin, G4double& pMax ) const
281 { 198 {
282 G4ThreeVector bmin, bmax; << 199 if (!pTransform.IsRotated())
283 G4bool exist; <<
284 <<
285 // Check bounding box (bbox) <<
286 // <<
287 BoundingLimits(bmin,bmax); <<
288 G4BoundingEnvelope bbox(bmin,bmax); <<
289 #ifdef G4BBOX_EXTENT <<
290 return bbox.CalculateExtent(pAxis,pVoxelLimi <<
291 #endif <<
292 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox <<
293 { 200 {
294 return exist = pMin < pMax; << 201 // Special case handling for unrotated solids
>> 202 // Compute x/y/z mins and maxs respecting limits, with early returns
>> 203 // if outside limits. Then switch() on pAxis
>> 204
>> 205 G4double xoffset,xMin,xMax;
>> 206 G4double yoffset,yMin,yMax;
>> 207 G4double zoffset,zMin,zMax;
>> 208
>> 209 zoffset=pTransform.NetTranslation().z();
>> 210 zMin=zoffset-fDz;
>> 211 zMax=zoffset+fDz;
>> 212 if (pVoxelLimit.IsZLimited())
>> 213 {
>> 214 if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
>> 215 || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
>> 216 {
>> 217 return false;
>> 218 }
>> 219 else
>> 220 {
>> 221 if (zMin<pVoxelLimit.GetMinZExtent())
>> 222 {
>> 223 zMin=pVoxelLimit.GetMinZExtent();
>> 224 }
>> 225 if (zMax>pVoxelLimit.GetMaxZExtent())
>> 226 {
>> 227 zMax=pVoxelLimit.GetMaxZExtent();
>> 228 }
>> 229 }
>> 230 }
>> 231 xoffset=pTransform.NetTranslation().x();
>> 232 if (fDx2 >= fDx1)
>> 233 {
>> 234 xMax = xoffset+(fDx1+fDx2)/2+(zMax-zoffset)*(fDx2-fDx1)/(2*fDz) ;
>> 235 xMin = 2*xoffset - xMax ;
>> 236 }
>> 237 else
>> 238 {
>> 239 xMax = xoffset+(fDx1+fDx2)/2+(zMin-zoffset)*(fDx2-fDx1)/(2*fDz) ;
>> 240 xMin = 2*xoffset - xMax ;
>> 241 }
>> 242 if (pVoxelLimit.IsXLimited())
>> 243 {
>> 244 if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
>> 245 || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
>> 246 {
>> 247 return false;
>> 248 }
>> 249 else
>> 250 {
>> 251 if (xMin<pVoxelLimit.GetMinXExtent())
>> 252 {
>> 253 xMin=pVoxelLimit.GetMinXExtent();
>> 254 }
>> 255 if (xMax>pVoxelLimit.GetMaxXExtent())
>> 256 {
>> 257 xMax=pVoxelLimit.GetMaxXExtent();
>> 258 }
>> 259 }
>> 260 }
>> 261 yoffset= pTransform.NetTranslation().y() ;
>> 262 if(fDy2 >= fDy1)
>> 263 {
>> 264 yMax = yoffset+(fDy2+fDy1)/2+(zMax-zoffset)*(fDy2-fDy1)/(2*fDz) ;
>> 265 yMin = 2*yoffset - yMax ;
>> 266 }
>> 267 else
>> 268 {
>> 269 yMax = yoffset+(fDy2+fDy1)/2+(zMin-zoffset)*(fDy2-fDy1)/(2*fDz) ;
>> 270 yMin = 2*yoffset - yMax ;
>> 271 }
>> 272 if (pVoxelLimit.IsYLimited())
>> 273 {
>> 274 if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
>> 275 || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
>> 276 {
>> 277 return false;
>> 278 }
>> 279 else
>> 280 {
>> 281 if (yMin<pVoxelLimit.GetMinYExtent())
>> 282 {
>> 283 yMin=pVoxelLimit.GetMinYExtent();
>> 284 }
>> 285 if (yMax>pVoxelLimit.GetMaxYExtent())
>> 286 {
>> 287 yMax=pVoxelLimit.GetMaxYExtent();
>> 288 }
>> 289 }
>> 290 }
>> 291
>> 292 switch (pAxis)
>> 293 {
>> 294 case kXAxis:
>> 295 pMin=xMin;
>> 296 pMax=xMax;
>> 297 break;
>> 298 case kYAxis:
>> 299 pMin=yMin;
>> 300 pMax=yMax;
>> 301 break;
>> 302 case kZAxis:
>> 303 pMin=zMin;
>> 304 pMax=zMax;
>> 305 break;
>> 306 default:
>> 307 break;
>> 308 }
>> 309
>> 310 // Add 2*Tolerance to avoid precision troubles ?
>> 311 //
>> 312 pMin-=kCarTolerance;
>> 313 pMax+=kCarTolerance;
>> 314
>> 315 return true;
295 } 316 }
>> 317 else
>> 318 {
>> 319 // General rotated case - create and clip mesh to boundaries
296 320
297 // Set bounding envelope (benv) and calculat << 321 G4bool existsAfterClip=false;
298 // << 322 G4ThreeVectorList *vertices;
299 G4double dx1 = GetXHalfLength1(); << 323
300 G4double dx2 = GetXHalfLength2(); << 324 pMin=+kInfinity;
301 G4double dy1 = GetYHalfLength1(); << 325 pMax=-kInfinity;
302 G4double dy2 = GetYHalfLength2(); << 326
303 G4double dz = GetZHalfLength(); << 327 // Calculate rotated vertex coordinates
304 << 328 //
305 G4ThreeVectorList baseA(4), baseB(4); << 329 vertices=CreateRotatedVertices(pTransform);
306 baseA[0].set(-dx1,-dy1,-dz); << 330 ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
307 baseA[1].set( dx1,-dy1,-dz); << 331 ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax);
308 baseA[2].set( dx1, dy1,-dz); << 332 ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
309 baseA[3].set(-dx1, dy1,-dz); << 333
310 baseB[0].set(-dx2,-dy2, dz); << 334 if (pMin!=kInfinity||pMax!=-kInfinity)
311 baseB[1].set( dx2,-dy2, dz); << 335 {
312 baseB[2].set( dx2, dy2, dz); << 336 existsAfterClip=true;
313 baseB[3].set(-dx2, dy2, dz); << 337
314 << 338 // Add 2*tolerance to avoid precision troubles
315 std::vector<const G4ThreeVectorList *> polyg << 339 //
316 polygons[0] = &baseA; << 340 pMin-=kCarTolerance;
317 polygons[1] = &baseB; << 341 pMax+=kCarTolerance;
318 << 342
319 G4BoundingEnvelope benv(bmin,bmax,polygons); << 343 }
320 exist = benv.CalculateExtent(pAxis,pVoxelLim << 344 else
321 return exist; << 345 {
>> 346 // Check for case where completely enveloping clipping volume
>> 347 // If point inside then we are confident that the solid completely
>> 348 // envelopes the clipping volume. Hence set min/max extents according
>> 349 // to clipping volume extents along the specified axis.
>> 350
>> 351 G4ThreeVector clipCentre(
>> 352 (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
>> 353 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
>> 354 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
>> 355
>> 356 if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
>> 357 {
>> 358 existsAfterClip=true;
>> 359 pMin=pVoxelLimit.GetMinExtent(pAxis);
>> 360 pMax=pVoxelLimit.GetMaxExtent(pAxis);
>> 361 }
>> 362 }
>> 363 delete vertices;
>> 364 return existsAfterClip;
>> 365 }
322 } 366 }
323 367
324 ////////////////////////////////////////////// << 368 ///////////////////////////////////////////////////////////////////
325 // 369 //
326 // Return whether point inside/outside/on surf 370 // Return whether point inside/outside/on surface, using tolerance
327 371
328 EInside G4Trd::Inside( const G4ThreeVector& p 372 EInside G4Trd::Inside( const G4ThreeVector& p ) const
329 { << 373 {
330 G4double dx = fPlanes[3].a*std::abs(p.x())+f << 374 EInside in=kOutside;
331 G4double dy = fPlanes[1].b*std::abs(p.y())+f << 375 G4double x,y,zbase1,zbase2;
332 G4double dxy = std::max(dx,dy); << 376
>> 377 if (std::fabs(p.z())<=fDz-kCarTolerance/2)
>> 378 {
>> 379 zbase1=p.z()+fDz; // Dist from -ve z plane
>> 380 zbase2=fDz-p.z(); // Dist from +ve z plane
333 381
334 G4double dz = std::abs(p.z())-fDz; << 382 // Check whether inside x tolerance
335 G4double dist = std::max(dz,dxy); << 383 //
>> 384 x=0.5*(fDx2*zbase1+fDx1*zbase2)/fDz - kCarTolerance/2;
>> 385 if (std::fabs(p.x())<=x)
>> 386 {
>> 387 y=0.5*((fDy2*zbase1+fDy1*zbase2))/fDz - kCarTolerance/2;
>> 388 if (std::fabs(p.y())<=y)
>> 389 {
>> 390 in=kInside;
>> 391 }
>> 392 else if (std::fabs(p.y())<=y+kCarTolerance)
>> 393 {
>> 394 in=kSurface;
>> 395 }
>> 396 }
>> 397 else if (std::fabs(p.x())<=x+kCarTolerance)
>> 398 {
>> 399 // y = y half width of shape at z of point + tolerant boundary
>> 400 //
>> 401 y=0.5*((fDy2*zbase1+fDy1*zbase2))/fDz + kCarTolerance/2;
>> 402 if (std::fabs(p.y())<=y)
>> 403 {
>> 404 in=kSurface;
>> 405 }
>> 406 }
>> 407 }
>> 408 else if (std::fabs(p.z())<=fDz+kCarTolerance/2)
>> 409 {
>> 410 // Only need to check outer tolerant boundaries
>> 411 //
>> 412 zbase1=p.z()+fDz; // Dist from -ve z plane
>> 413 zbase2=fDz-p.z(); // Dist from +ve z plane
336 414
337 return (dist > halfCarTolerance) ? kOutside << 415 // x = x half width of shape at z of point plus tolerance
338 ((dist > -halfCarTolerance) ? kSurface : k << 416 //
>> 417 x=0.5*(fDx2*zbase1+fDx1*zbase2)/fDz + kCarTolerance/2;
>> 418 if (std::fabs(p.x())<=x)
>> 419 {
>> 420 // y = y half width of shape at z of point
>> 421 //
>> 422 y=0.5*((fDy2*zbase1+fDy1*zbase2))/fDz + kCarTolerance/2;
>> 423 if (std::fabs(p.y())<=y) in=kSurface;
>> 424 }
>> 425 }
>> 426 return in;
339 } 427 }
340 428
341 ////////////////////////////////////////////// 429 //////////////////////////////////////////////////////////////////////////
342 // 430 //
343 // Determine side where point is, and return c << 431 // Calculate side nearest to p, and return normal
>> 432 // If two sides are equidistant, normal of first side (x/y/z)
>> 433 // encountered returned
344 434
345 G4ThreeVector G4Trd::SurfaceNormal( const G4Th 435 G4ThreeVector G4Trd::SurfaceNormal( const G4ThreeVector& p ) const
346 { 436 {
347 G4int nsurf = 0; // number of surfaces where << 437 G4ThreeVector norm, sumnorm(0.,0.,0.);
348 << 438 G4int noSurfaces = 0;
349 // Check Z faces << 439 G4double z = 2.0*fDz, tanx, secx, newpx, widx;
350 // << 440 G4double tany, secy, newpy, widy;
351 G4double nz = 0; << 441 G4double distx, disty, distz, fcos;
352 G4double dz = std::abs(p.z()) - fDz; << 442 G4double delta = 0.5*kCarTolerance;
353 if (std::abs(dz) <= halfCarTolerance) << 443
354 { << 444 tanx = (fDx2 - fDx1)/z;
355 nz = (p.z() < 0) ? -1 : 1; << 445 secx = std::sqrt(1.0+tanx*tanx);
356 ++nsurf; << 446 newpx = std::fabs(p.x())-p.z()*tanx;
>> 447 widx = fDx2 - fDz*tanx;
>> 448
>> 449 tany = (fDy2 - fDy1)/z;
>> 450 secy = std::sqrt(1.0+tany*tany);
>> 451 newpy = std::fabs(p.y())-p.z()*tany;
>> 452 widy = fDy2 - fDz*tany;
>> 453
>> 454 distx = std::fabs(newpx-widx)/secx; // perp. distance to x side
>> 455 disty = std::fabs(newpy-widy)/secy; // to y side
>> 456 distz = std::fabs(std::fabs(p.z())-fDz); // to z side
>> 457
>> 458 fcos = 1.0/secx;
>> 459 G4ThreeVector nX = G4ThreeVector( fcos,0,-tanx*fcos);
>> 460 G4ThreeVector nmX = G4ThreeVector(-fcos,0,-tanx*fcos);
>> 461
>> 462 fcos = 1.0/secy;
>> 463 G4ThreeVector nY = G4ThreeVector(0, fcos,-tany*fcos);
>> 464 G4ThreeVector nmY = G4ThreeVector(0,-fcos,-tany*fcos);
>> 465 G4ThreeVector nZ = G4ThreeVector( 0, 0, 1.0);
>> 466
>> 467 if (distx <= delta)
>> 468 {
>> 469 noSurfaces ++;
>> 470 if ( p.x() >= 0.) sumnorm += nX;
>> 471 else sumnorm += nmX;
>> 472 }
>> 473 if (disty <= delta)
>> 474 {
>> 475 noSurfaces ++;
>> 476 if ( p.y() >= 0.) sumnorm += nY;
>> 477 else sumnorm += nmY;
>> 478 }
>> 479 if (distz <= delta)
>> 480 {
>> 481 noSurfaces ++;
>> 482 if ( p.z() >= 0.) sumnorm += nZ;
>> 483 else sumnorm -= nZ;
357 } 484 }
358 << 485 if ( noSurfaces == 0 )
359 // Check Y faces <<
360 // <<
361 G4double ny = 0; <<
362 G4double dy1 = fPlanes[0].b*p.y(); <<
363 G4double dy2 = fPlanes[0].c*p.z() + fPlanes[ <<
364 if (std::abs(dy2 + dy1) <= halfCarTolerance) <<
365 { 486 {
366 ny += fPlanes[0].b; <<
367 nz += fPlanes[0].c; <<
368 ++nsurf; <<
369 } <<
370 if (std::abs(dy2 - dy1) <= halfCarTolerance) <<
371 { <<
372 ny += fPlanes[1].b; <<
373 nz += fPlanes[1].c; <<
374 ++nsurf; <<
375 } <<
376 <<
377 // Check X faces <<
378 // <<
379 G4double nx = 0; <<
380 G4double dx1 = fPlanes[2].a*p.x(); <<
381 G4double dx2 = fPlanes[2].c*p.z() + fPlanes[ <<
382 if (std::abs(dx2 + dx1) <= halfCarTolerance) <<
383 { <<
384 nx += fPlanes[2].a; <<
385 nz += fPlanes[2].c; <<
386 ++nsurf; <<
387 } <<
388 if (std::abs(dx2 - dx1) <= halfCarTolerance) <<
389 { <<
390 nx += fPlanes[3].a; <<
391 nz += fPlanes[3].c; <<
392 ++nsurf; <<
393 } <<
394 <<
395 // Return normal <<
396 // <<
397 if (nsurf == 1) return {nx,ny,nz}; <<
398 else if (nsurf != 0) return G4ThreeVector(nx <<
399 else <<
400 { <<
401 // Point is not on the surface <<
402 // <<
403 #ifdef G4CSGDEBUG 487 #ifdef G4CSGDEBUG
404 std::ostringstream message; << 488 G4Exception("G4Trd::SurfaceNormal(p)", "GeomSolids1002", JustWarning,
405 G4long oldprc = message.precision(16); << 489 "Point p is not on surface !?" );
406 message << "Point p is not on surface (!?) << 490 #endif
407 << GetName() << G4endl; << 491 norm = ApproxSurfaceNormal(p);
408 message << "Position:\n"; <<
409 message << " p.x() = " << p.x()/mm << " <<
410 message << " p.y() = " << p.y()/mm << " <<
411 message << " p.z() = " << p.z()/mm << " <<
412 G4cout.precision(oldprc) ; <<
413 G4Exception("G4Trd::SurfaceNormal(p)", "Ge <<
414 JustWarning, message ); <<
415 DumpInfo(); <<
416 #endif <<
417 return ApproxSurfaceNormal(p); <<
418 } 492 }
>> 493 else if ( noSurfaces == 1 ) norm = sumnorm;
>> 494 else norm = sumnorm.unit();
>> 495 return norm;
419 } 496 }
420 497
421 ////////////////////////////////////////////// << 498
>> 499 /////////////////////////////////////////////////////////////////////////////
422 // 500 //
423 // Algorithm for SurfaceNormal() following the 501 // Algorithm for SurfaceNormal() following the original specification
424 // for points not on the surface 502 // for points not on the surface
425 503
426 G4ThreeVector G4Trd::ApproxSurfaceNormal( cons 504 G4ThreeVector G4Trd::ApproxSurfaceNormal( const G4ThreeVector& p ) const
427 { 505 {
428 G4double dist = -DBL_MAX; << 506 G4ThreeVector norm;
429 G4int iside = 0; << 507 G4double z,tanx,secx,newpx,widx;
430 for (G4int i=0; i<4; ++i) << 508 G4double tany,secy,newpy,widy;
431 { << 509 G4double distx,disty,distz,fcos;
432 G4double d = fPlanes[i].a*p.x() + << 510
433 fPlanes[i].b*p.y() + << 511 z=2.0*fDz;
434 fPlanes[i].c*p.z() + fPlanes[ << 512
435 if (d > dist) { dist = d; iside = i; } << 513 tanx=(fDx2-fDx1)/z;
436 } << 514 secx=std::sqrt(1.0+tanx*tanx);
>> 515 newpx=std::fabs(p.x())-p.z()*tanx;
>> 516 widx=fDx2-fDz*tanx;
>> 517
>> 518 tany=(fDy2-fDy1)/z;
>> 519 secy=std::sqrt(1.0+tany*tany);
>> 520 newpy=std::fabs(p.y())-p.z()*tany;
>> 521 widy=fDy2-fDz*tany;
>> 522
>> 523 distx=std::fabs(newpx-widx)/secx; // perpendicular distance to x side
>> 524 disty=std::fabs(newpy-widy)/secy; // to y side
>> 525 distz=std::fabs(std::fabs(p.z())-fDz); // to z side
437 526
438 G4double distz = std::abs(p.z()) - fDz; << 527 // find closest side
439 if (dist > distz) << 528 //
440 return { fPlanes[iside].a, fPlanes[iside]. << 529 if (distx<=disty)
>> 530 {
>> 531 if (distx<=distz)
>> 532 {
>> 533 // Closest to X
>> 534 //
>> 535 fcos=1.0/secx;
>> 536 // normal=(+/-std::cos(ang),0,-std::sin(ang))
>> 537 if (p.x()>=0)
>> 538 norm=G4ThreeVector(fcos,0,-tanx*fcos);
>> 539 else
>> 540 norm=G4ThreeVector(-fcos,0,-tanx*fcos);
>> 541 }
>> 542 else
>> 543 {
>> 544 // Closest to Z
>> 545 //
>> 546 if (p.z()>=0)
>> 547 norm=G4ThreeVector(0,0,1);
>> 548 else
>> 549 norm=G4ThreeVector(0,0,-1);
>> 550 }
>> 551 }
441 else 552 else
442 return { 0, 0, (G4double)((p.z() < 0) ? -1 << 553 {
>> 554 if (disty<=distz)
>> 555 {
>> 556 // Closest to Y
>> 557 //
>> 558 fcos=1.0/secy;
>> 559 if (p.y()>=0)
>> 560 norm=G4ThreeVector(0,fcos,-tany*fcos);
>> 561 else
>> 562 norm=G4ThreeVector(0,-fcos,-tany*fcos);
>> 563 }
>> 564 else
>> 565 {
>> 566 // Closest to Z
>> 567 //
>> 568 if (p.z()>=0)
>> 569 norm=G4ThreeVector(0,0,1);
>> 570 else
>> 571 norm=G4ThreeVector(0,0,-1);
>> 572 }
>> 573 }
>> 574 return norm;
443 } 575 }
444 576
445 ////////////////////////////////////////////// << 577 ////////////////////////////////////////////////////////////////////////////
446 // 578 //
447 // Calculate distance to shape from outside 579 // Calculate distance to shape from outside
448 // - return kInfinity if no intersection << 580 // - return kInfinity if no intersection
>> 581 //
>> 582 // ALGORITHM:
>> 583 // For each component, calculate pair of minimum and maximum intersection
>> 584 // values for which the particle is in the extent of the shape
>> 585 // - The smallest (MAX minimum) allowed distance of the pairs is intersect
>> 586 // - Z plane intersectin uses tolerance
>> 587 // - XZ YZ planes use logic & *SLIGHTLY INCORRECT* tolerance
>> 588 // (this saves at least 1 sqrt, 1 multiply and 1 divide... in applicable
>> 589 // cases)
>> 590 // - Note: XZ and YZ planes each divide space into four regions,
>> 591 // characterised by ss1 ss2
>> 592 // NOTE:
>> 593 //
>> 594 // `Inside' safe - meaningful answers given if point is inside the exact
>> 595 // shape.
>> 596
>> 597 G4double G4Trd::DistanceToIn( const G4ThreeVector& p,
>> 598 const G4ThreeVector& v ) const
>> 599 {
>> 600 G4double snxt = kInfinity ; // snxt = default return value
>> 601 G4double smin,smax;
>> 602 G4double s1,s2,tanxz,tanyz,ds1,ds2;
>> 603 G4double ss1,ss2,sn1=0.,sn2=0.,Dist;
449 604
450 G4double G4Trd::DistanceToIn(const G4ThreeVect << 605 if ( v.z() ) // Calculate valid z intersect range
451 const G4ThreeVect << 606 {
452 { << 607 if ( v.z() > 0 ) // Calculate smax: must be +ve or no intersection.
453 // Z intersections << 608 {
454 // << 609 Dist = fDz - p.z() ; // to plane at +dz
455 if ((std::abs(p.z()) - fDz) >= -halfCarToler <<
456 return kInfinity; <<
457 G4double invz = (-v.z() == 0) ? DBL_MAX : -1 <<
458 G4double dz = (invz < 0) ? fDz : -fDz; <<
459 G4double tzmin = (p.z() + dz)*invz; <<
460 G4double tzmax = (p.z() - dz)*invz; <<
461 610
462 // Y intersections << 611 if (Dist >= 0.5*kCarTolerance)
463 // << 612 {
464 G4double tmin0 = tzmin, tmax0 = tzmax; << 613 smax = Dist/v.z() ;
465 G4double ya = fPlanes[0].b*v.y(), yb = fPlan << 614 smin = -(fDz + p.z())/v.z() ;
466 G4double yc = fPlanes[0].b*p.y(), yd = fPlan << 615 }
467 G4double cos0 = yb + ya; << 616 else return snxt ;
468 G4double dis0 = yd + yc; << 617 }
469 if (dis0 >= -halfCarTolerance) << 618 else // v.z <0
>> 619 {
>> 620 Dist=fDz+p.z(); // plane at -dz
>> 621
>> 622 if ( Dist >= 0.5*kCarTolerance )
>> 623 {
>> 624 smax = -Dist/v.z() ;
>> 625 smin = (fDz - p.z())/v.z() ;
>> 626 }
>> 627 else return snxt ;
>> 628 }
>> 629 if (smin < 0 ) smin = 0 ;
>> 630 }
>> 631 else // v.z=0
470 { 632 {
471 if (cos0 >= 0) return kInfinity; << 633 if (std::fabs(p.z()) >= fDz ) return snxt ; // Outside & no intersect
472 G4double tmp = -dis0/cos0; << 634 else
473 if (tmin0 < tmp) tmin0 = tmp; << 635 {
>> 636 smin = 0 ; // Always inside z range
>> 637 smax = kInfinity;
>> 638 }
474 } 639 }
475 else if (cos0 > 0) << 640
>> 641 // Calculate x intersection range
>> 642 //
>> 643 // Calc half width at p.z, and components towards planes
>> 644
>> 645 tanxz = (fDx2 - fDx1)*0.5/fDz ;
>> 646 s1 = 0.5*(fDx1+fDx2) + tanxz*p.z() ; // x half width at p.z
>> 647 ds1 = v.x() - tanxz*v.z() ; // Components of v towards faces at +-x
>> 648 ds2 = v.x() + tanxz*v.z() ;
>> 649 ss1 = s1 - p.x() ; // -delta x to +ve plane
>> 650 // -ve when outside
>> 651 ss2 = -s1 - p.x() ; // -delta x to -ve plane
>> 652 // +ve when outside
>> 653
>> 654 if (ss1 < 0 && ss2 <= 0 )
476 { 655 {
477 G4double tmp = -dis0/cos0; << 656 if (ds1 < 0) // In +ve coord Area
478 if (tmax0 > tmp) tmax0 = tmp; << 657 {
>> 658 sn1 = ss1/ds1 ;
>> 659
>> 660 if ( ds2 < 0 ) sn2 = ss2/ds2 ;
>> 661 else sn2 = kInfinity ;
>> 662 }
>> 663 else return snxt ;
479 } 664 }
>> 665 else if ( ss1 >= 0 && ss2 > 0 )
>> 666 {
>> 667 if ( ds2 > 0 ) // In -ve coord Area
>> 668 {
>> 669 sn1 = ss2/ds2 ;
480 670
481 G4double tmin1 = tmin0, tmax1 = tmax0; << 671 if (ds1 > 0) sn2 = ss1/ds1 ;
482 G4double cos1 = yb - ya; << 672 else sn2 = kInfinity;
483 G4double dis1 = yd - yc; << 673
484 if (dis1 >= -halfCarTolerance) << 674 }
>> 675 else return snxt ;
>> 676 }
>> 677 else if (ss1 >= 0 && ss2 <= 0 )
485 { 678 {
486 if (cos1 >= 0) return kInfinity; << 679 // Inside Area - calculate leaving distance
487 G4double tmp = -dis1/cos1; << 680 // *Don't* use exact distance to side for tolerance
488 if (tmin1 < tmp) tmin1 = tmp; << 681 // = ss1*std::cos(ang xz)
>> 682 // = ss1/std::sqrt(1.0+tanxz*tanxz)
>> 683 sn1 = 0 ;
>> 684
>> 685 if ( ds1 > 0 )
>> 686 {
>> 687 if (ss1 > 0.5*kCarTolerance) sn2 = ss1/ds1 ; // Leave +ve side extent
>> 688 else return snxt ; // Leave immediately by +ve
>> 689 }
>> 690 else sn2 = kInfinity ;
>> 691
>> 692 if ( ds2 < 0 )
>> 693 {
>> 694 if ( ss2 < -0.5*kCarTolerance )
>> 695 {
>> 696 Dist = ss2/ds2 ; // Leave -ve side extent
>> 697 if ( Dist < sn2 ) sn2 = Dist ;
>> 698 }
>> 699 else return snxt ;
>> 700 }
489 } 701 }
490 else if (cos1 > 0) << 702 else if (ss1 < 0 && ss2 > 0 )
491 { 703 {
492 G4double tmp = -dis1/cos1; << 704 // Within +/- plane cross-over areas (not on boundaries ss1||ss2==0)
493 if (tmax1 > tmp) tmax1 = tmp; << 705
>> 706 if ( ds1 >= 0 || ds2 <= 0 )
>> 707 {
>> 708 return snxt ;
>> 709 }
>> 710 else // Will intersect & stay inside
>> 711 {
>> 712 sn1 = ss1/ds1 ;
>> 713 Dist = ss2/ds2 ;
>> 714 if (Dist > sn1 ) sn1 = Dist ;
>> 715 sn2 = kInfinity ;
>> 716 }
494 } 717 }
495 718
496 // X intersections << 719 // Reduce allowed range of distances as appropriate
497 // << 720
498 G4double tmin2 = tmin1, tmax2 = tmax1; << 721 if ( sn1 > smin ) smin = sn1 ;
499 G4double xa = fPlanes[2].a*v.x(), xb = fPlan << 722 if ( sn2 < smax ) smax = sn2 ;
500 G4double xc = fPlanes[2].a*p.x(), xd = fPlan << 723
501 G4double cos2 = xb + xa; << 724 // Check for incompatible ranges (eg z intersects between 50 ->100 and x
502 G4double dis2 = xd + xc; << 725 // only 10-40 -> no intersection)
503 if (dis2 >= -halfCarTolerance) << 726
>> 727 if ( smax < smin ) return snxt ;
>> 728
>> 729 // Calculate valid y intersection range
>> 730 // (repeat of x intersection code)
>> 731
>> 732 tanyz = (fDy2-fDy1)*0.5/fDz ;
>> 733 s2 = 0.5*(fDy1+fDy2) + tanyz*p.z() ; // y half width at p.z
>> 734 ds1 = v.y() - tanyz*v.z() ; // Components of v towards faces at +-y
>> 735 ds2 = v.y() + tanyz*v.z() ;
>> 736 ss1 = s2 - p.y() ; // -delta y to +ve plane
>> 737 ss2 = -s2 - p.y() ; // -delta y to -ve plane
>> 738
>> 739 if ( ss1 < 0 && ss2 <= 0 )
504 { 740 {
505 if (cos2 >= 0) return kInfinity; << 741 if (ds1 < 0 ) // In +ve coord Area
506 G4double tmp = -dis2/cos2; << 742 {
507 if (tmin2 < tmp) tmin2 = tmp; << 743 sn1 = ss1/ds1 ;
>> 744 if ( ds2 < 0 ) sn2 = ss2/ds2 ;
>> 745 else sn2 = kInfinity ;
>> 746 }
>> 747 else return snxt ;
508 } 748 }
509 else if (cos2 > 0) << 749 else if ( ss1 >= 0 && ss2 > 0 )
510 { 750 {
511 G4double tmp = -dis2/cos2; << 751 if ( ds2 > 0 ) // In -ve coord Area
512 if (tmax2 > tmp) tmax2 = tmp; << 752 {
>> 753 sn1 = ss2/ds2 ;
>> 754 if ( ds1 > 0 ) sn2 = ss1/ds1 ;
>> 755 else sn2 = kInfinity ;
>> 756 }
>> 757 else return snxt ;
513 } 758 }
514 << 759 else if (ss1 >= 0 && ss2 <= 0 )
515 G4double tmin3 = tmin2, tmax3 = tmax2; <<
516 G4double cos3 = xb - xa; <<
517 G4double dis3 = xd - xc; <<
518 if (dis3 >= -halfCarTolerance) <<
519 { 760 {
520 if (cos3 >= 0) return kInfinity; << 761 // Inside Area - calculate leaving distance
521 G4double tmp = -dis3/cos3; << 762 // *Don't* use exact distance to side for tolerance
522 if (tmin3 < tmp) tmin3 = tmp; << 763 // = ss1*std::cos(ang yz)
>> 764 // = ss1/std::sqrt(1.0+tanyz*tanyz)
>> 765 sn1 = 0 ;
>> 766
>> 767 if ( ds1 > 0 )
>> 768 {
>> 769 if (ss1 > 0.5*kCarTolerance) sn2 = ss1/ds1 ; // Leave +ve side extent
>> 770 else return snxt ; // Leave immediately by +ve
>> 771 }
>> 772 else sn2 = kInfinity ;
>> 773
>> 774 if ( ds2 < 0 )
>> 775 {
>> 776 if ( ss2 < -0.5*kCarTolerance )
>> 777 {
>> 778 Dist = ss2/ds2 ; // Leave -ve side extent
>> 779 if (Dist < sn2) sn2=Dist;
>> 780 }
>> 781 else return snxt ;
>> 782 }
523 } 783 }
524 else if (cos3 > 0) << 784 else if (ss1 < 0 && ss2 > 0 )
525 { 785 {
526 G4double tmp = -dis3/cos3; << 786 // Within +/- plane cross-over areas (not on boundaries ss1||ss2==0)
527 if (tmax3 > tmp) tmax3 = tmp; << 787
>> 788 if (ds1 >= 0 || ds2 <= 0 )
>> 789 {
>> 790 return snxt ;
>> 791 }
>> 792 else // Will intersect & stay inside
>> 793 {
>> 794 sn1 = ss1/ds1 ;
>> 795 Dist = ss2/ds2 ;
>> 796 if (Dist > sn1 ) sn1 = Dist ;
>> 797 sn2 = kInfinity ;
>> 798 }
528 } 799 }
>> 800
>> 801 // Reduce allowed range of distances as appropriate
529 802
530 // Find distance << 803 if ( sn1 > smin) smin = sn1 ;
531 // << 804 if ( sn2 < smax) smax = sn2 ;
532 G4double tmin = tmin3, tmax = tmax3; << 805
533 if (tmax <= tmin + halfCarTolerance) return << 806 // Check for incompatible ranges (eg x intersects between 50 ->100 and y
534 return (tmin < halfCarTolerance ) ? 0. : tmi << 807 // only 10-40 -> no intersection). Set snxt if ok
>> 808
>> 809 if ( smax > smin ) snxt = smin ;
>> 810 if (snxt < 0.5*kCarTolerance ) snxt = 0.0 ;
>> 811
>> 812 return snxt ;
535 } 813 }
536 814
537 ////////////////////////////////////////////// << 815 /////////////////////////////////////////////////////////////////////////
538 // 816 //
539 // Calculate exact shortest distance to any bo << 817 // Approximate distance to shape
540 // This is the best fast estimation of the sho << 818 // Calculate perpendicular distances to z/x/y surfaces, return largest
541 // - returns 0 if point is inside << 819 // which is the most fast estimation of shortest distance to Trd
>> 820 // - Safe underestimate
>> 821 // - If point within exact shape, return 0
542 822
543 G4double G4Trd::DistanceToIn( const G4ThreeVec 823 G4double G4Trd::DistanceToIn( const G4ThreeVector& p ) const
544 { 824 {
545 G4double dx = fPlanes[3].a*std::abs(p.x())+f << 825 G4double safe=0.0;
546 G4double dy = fPlanes[1].b*std::abs(p.y())+f << 826 G4double tanxz,distx,safx;
547 G4double dxy = std::max(dx,dy); << 827 G4double tanyz,disty,safy;
>> 828 G4double zbase;
>> 829
>> 830 safe=std::fabs(p.z())-fDz;
>> 831 if (safe<0) safe=0; // Also used to ensure x/y distances
>> 832 // POSITIVE
548 833
549 G4double dz = std::abs(p.z())-fDz; << 834 zbase=fDz+p.z();
550 G4double dist = std::max(dz,dxy); <<
551 835
552 return (dist > 0) ? dist : 0.; << 836 // Find distance along x direction to closest x plane
>> 837 //
>> 838 tanxz=(fDx2-fDx1)*0.5/fDz;
>> 839 // widx=fDx1+tanxz*(fDz+p.z()); // x width at p.z
>> 840 // distx=std::fabs(p.x())-widx; // distance to plane
>> 841 distx=std::fabs(p.x())-(fDx1+tanxz*zbase);
>> 842 if (distx>safe)
>> 843 {
>> 844 safx=distx/std::sqrt(1.0+tanxz*tanxz); // vector Dist=Dist*std::cos(ang)
>> 845 if (safx>safe) safe=safx;
>> 846 }
>> 847
>> 848 // Find distance along y direction to slanted wall
>> 849 tanyz=(fDy2-fDy1)*0.5/fDz;
>> 850 // widy=fDy1+tanyz*(fDz+p.z()); // y width at p.z
>> 851 // disty=std::fabs(p.y())-widy; // distance to plane
>> 852 disty=std::fabs(p.y())-(fDy1+tanyz*zbase);
>> 853 if (disty>safe)
>> 854 {
>> 855 safy=disty/std::sqrt(1.0+tanyz*tanyz); // distance along vector
>> 856 if (safy>safe) safe=safy;
>> 857 }
>> 858 return safe;
553 } 859 }
554 860
555 ////////////////////////////////////////////// << 861 ////////////////////////////////////////////////////////////////////////
556 // 862 //
557 // Calculate distance to surface of shape from << 863 // Calcluate distance to surface of shape from inside
558 // find normal at exit point, if required << 864 // Calculate distance to x/y/z planes - smallest is exiting distance
559 // - when leaving the surface, return 0 << 865 // - z planes have std. check for tolerance
560 << 866 // - xz yz planes have check based on distance || to x or y axis
561 G4double G4Trd::DistanceToOut(const G4ThreeVec << 867 // (not corrected for slope of planes)
562 const G4bool cal << 868 // ?BUG? If v.z==0 are there cases when snside not set????
563 G4bool* va << 869
>> 870 G4double G4Trd::DistanceToOut( const G4ThreeVector& p,
>> 871 const G4ThreeVector& v,
>> 872 const G4bool calcNorm,
>> 873 G4bool *validNorm,
>> 874 G4ThreeVector *n ) const
564 { 875 {
565 // Z intersections << 876 ESide side = kUndefined, snside = kUndefined;
566 // << 877 G4double snxt,pdist;
567 if ((std::abs(p.z()) - fDz) >= -halfCarToler << 878 G4double central,ss1,ss2,ds1,ds2,sn=0.,sn2=0.;
>> 879 G4double tanxz=0.,cosxz=0.,tanyz=0.,cosyz=0.;
>> 880
>> 881 if (calcNorm) *validNorm=true; // All normals are valid
>> 882
>> 883 // Calculate z plane intersection
>> 884 if (v.z()>0)
568 { 885 {
569 if (calcNorm) << 886 pdist=fDz-p.z();
>> 887 if (pdist>kCarTolerance/2)
570 { 888 {
571 *validNorm = true; << 889 snxt=pdist/v.z();
572 n->set(0, 0, (p.z() < 0) ? -1 : 1); << 890 side=kPZ;
>> 891 }
>> 892 else
>> 893 {
>> 894 if (calcNorm)
>> 895 {
>> 896 *n=G4ThreeVector(0,0,1);
>> 897 }
>> 898 return snxt=0;
>> 899 }
>> 900 }
>> 901 else if (v.z()<0)
>> 902 {
>> 903 pdist=fDz+p.z();
>> 904 if (pdist>kCarTolerance/2)
>> 905 {
>> 906 snxt=-pdist/v.z();
>> 907 side=kMZ;
>> 908 }
>> 909 else
>> 910 {
>> 911 if (calcNorm)
>> 912 {
>> 913 *n=G4ThreeVector(0,0,-1);
>> 914 }
>> 915 return snxt=0;
573 } 916 }
574 return 0; <<
575 } 917 }
576 G4double vz = v.z(); << 918 else
577 G4double tmax = (vz == 0) ? DBL_MAX : (std:: << 919 {
578 G4int iside = (vz < 0) ? -4 : -2; // little << 920 snxt=kInfinity;
>> 921 }
579 922
580 // Y intersections <<
581 // 923 //
582 G4int i = 0; << 924 // Calculate x intersection
583 for ( ; i<2; ++i) << 925 //
584 { << 926 tanxz=(fDx2-fDx1)*0.5/fDz;
585 G4double cosa = fPlanes[i].b*v.y() + fPlan << 927 central=0.5*(fDx1+fDx2);
586 if (cosa > 0) << 928
>> 929 // +ve plane (1)
>> 930 //
>> 931 ss1=central+tanxz*p.z()-p.x(); // distance || x axis to plane
>> 932 // (+ve if point inside)
>> 933 ds1=v.x()-tanxz*v.z(); // component towards plane at +x
>> 934 // (-ve if +ve -> -ve direction)
>> 935 // -ve plane (2)
>> 936 //
>> 937 ss2=-tanxz*p.z()-p.x()-central; //distance || x axis to plane
>> 938 // (-ve if point inside)
>> 939 ds2=tanxz*v.z()+v.x(); // component towards plane at -x
>> 940
>> 941 if (ss1>0&&ss2<0)
>> 942 {
>> 943 // Normal case - entirely inside region
>> 944 if (ds1<=0&&ds2<0)
>> 945 {
>> 946 if (ss2<-kCarTolerance/2)
>> 947 {
>> 948 sn=ss2/ds2; // Leave by -ve side
>> 949 snside=kMX;
>> 950 }
>> 951 else
>> 952 {
>> 953 sn=0; // Leave immediately by -ve side
>> 954 snside=kMX;
>> 955 }
>> 956 }
>> 957 else if (ds1>0&&ds2>=0)
587 { 958 {
588 G4double dist = fPlanes[i].b*p.y()+fPlan << 959 if (ss1>kCarTolerance/2)
589 if (dist >= -halfCarTolerance) <<
590 { 960 {
591 if (calcNorm) << 961 sn=ss1/ds1; // Leave by +ve side
>> 962 snside=kPX;
>> 963 }
>> 964 else
>> 965 {
>> 966 sn=0; // Leave immediately by +ve side
>> 967 snside=kPX;
>> 968 }
>> 969 }
>> 970 else if (ds1>0&&ds2<0)
>> 971 {
>> 972 if (ss1>kCarTolerance/2)
>> 973 {
>> 974 // sn=ss1/ds1; // Leave by +ve side
>> 975 if (ss2<-kCarTolerance/2)
592 { 976 {
593 *validNorm = true; << 977 sn=ss1/ds1; // Leave by +ve side
594 n->set(0, fPlanes[i].b, fPlanes[i].c << 978 sn2=ss2/ds2;
>> 979 if (sn2<sn)
>> 980 {
>> 981 sn=sn2;
>> 982 snside=kMX;
>> 983 }
>> 984 else
>> 985 {
>> 986 snside=kPX;
>> 987 }
595 } 988 }
596 return 0; << 989 else
>> 990 {
>> 991 sn=0; // Leave immediately by -ve
>> 992 snside=kMX;
>> 993 }
>> 994 }
>> 995 else
>> 996 {
>> 997 sn=0; // Leave immediately by +ve side
>> 998 snside=kPX;
>> 999 }
>> 1000 }
>> 1001 else
>> 1002 {
>> 1003 // Must be || to both
>> 1004 //
>> 1005 sn=kInfinity; // Don't leave by either side
>> 1006 }
>> 1007 }
>> 1008 else if (ss1<=0&&ss2<0)
>> 1009 {
>> 1010 // Outside, in +ve Area
>> 1011
>> 1012 if (ds1>0)
>> 1013 {
>> 1014 sn=0; // Away from shape
>> 1015 // Left by +ve side
>> 1016 snside=kPX;
>> 1017 }
>> 1018 else
>> 1019 {
>> 1020 if (ds2<0)
>> 1021 {
>> 1022 // Ignore +ve plane and use -ve plane intersect
>> 1023 //
>> 1024 sn=ss2/ds2; // Leave by -ve side
>> 1025 snside=kMX;
>> 1026 }
>> 1027 else
>> 1028 {
>> 1029 // Must be || to both -> exit determined by other axes
>> 1030 //
>> 1031 sn=kInfinity; // Don't leave by either side
597 } 1032 }
598 G4double tmp = -dist/cosa; <<
599 if (tmax > tmp) { tmax = tmp; iside = i; <<
600 } 1033 }
601 } 1034 }
>> 1035 else if (ss1>0&&ss2>=0)
>> 1036 {
>> 1037 // Outside, in -ve Area
602 1038
603 // X intersections << 1039 if (ds2<0)
604 // << 1040 {
605 for ( ; i<4; ++i) << 1041 sn=0; // away from shape
>> 1042 // Left by -ve side
>> 1043 snside=kMX;
>> 1044 }
>> 1045 else
>> 1046 {
>> 1047 if (ds1>0)
>> 1048 {
>> 1049 // Ignore +ve plane and use -ve plane intersect
>> 1050 //
>> 1051 sn=ss1/ds1; // Leave by +ve side
>> 1052 snside=kPX;
>> 1053 }
>> 1054 else
>> 1055 {
>> 1056 // Must be || to both -> exit determined by other axes
>> 1057 //
>> 1058 sn=kInfinity; // Don't leave by either side
>> 1059 }
>> 1060 }
>> 1061 }
>> 1062
>> 1063 // Update minimum exit distance
>> 1064
>> 1065 if (sn<snxt)
606 { 1066 {
607 G4double cosa = fPlanes[i].a*v.x()+fPlanes << 1067 snxt=sn;
608 if (cosa > 0) << 1068 side=snside;
>> 1069 }
>> 1070 if (snxt>0)
>> 1071 {
>> 1072 // Calculate y intersection
>> 1073
>> 1074 tanyz=(fDy2-fDy1)*0.5/fDz;
>> 1075 central=0.5*(fDy1+fDy2);
>> 1076
>> 1077 // +ve plane (1)
>> 1078 //
>> 1079 ss1=central+tanyz*p.z()-p.y(); // distance || y axis to plane
>> 1080 // (+ve if point inside)
>> 1081 ds1=v.y()-tanyz*v.z(); // component towards +ve plane
>> 1082 // (-ve if +ve -> -ve direction)
>> 1083 // -ve plane (2)
>> 1084 //
>> 1085 ss2=-tanyz*p.z()-p.y()-central; // distance || y axis to plane
>> 1086 // (-ve if point inside)
>> 1087 ds2=tanyz*v.z()+v.y(); // component towards -ve plane
>> 1088
>> 1089 if (ss1>0&&ss2<0)
>> 1090 {
>> 1091 // Normal case - entirely inside region
>> 1092
>> 1093 if (ds1<=0&&ds2<0)
>> 1094 {
>> 1095 if (ss2<-kCarTolerance/2)
>> 1096 {
>> 1097 sn=ss2/ds2; // Leave by -ve side
>> 1098 snside=kMY;
>> 1099 }
>> 1100 else
>> 1101 {
>> 1102 sn=0; // Leave immediately by -ve side
>> 1103 snside=kMY;
>> 1104 }
>> 1105 }
>> 1106 else if (ds1>0&&ds2>=0)
>> 1107 {
>> 1108 if (ss1>kCarTolerance/2)
>> 1109 {
>> 1110 sn=ss1/ds1; // Leave by +ve side
>> 1111 snside=kPY;
>> 1112 }
>> 1113 else
>> 1114 {
>> 1115 sn=0; // Leave immediately by +ve side
>> 1116 snside=kPY;
>> 1117 }
>> 1118 }
>> 1119 else if (ds1>0&&ds2<0)
>> 1120 {
>> 1121 if (ss1>kCarTolerance/2)
>> 1122 {
>> 1123 // sn=ss1/ds1; // Leave by +ve side
>> 1124 if (ss2<-kCarTolerance/2)
>> 1125 {
>> 1126 sn=ss1/ds1; // Leave by +ve side
>> 1127 sn2=ss2/ds2;
>> 1128 if (sn2<sn)
>> 1129 {
>> 1130 sn=sn2;
>> 1131 snside=kMY;
>> 1132 }
>> 1133 else
>> 1134 {
>> 1135 snside=kPY;
>> 1136 }
>> 1137 }
>> 1138 else
>> 1139 {
>> 1140 sn=0; // Leave immediately by -ve
>> 1141 snside=kMY;
>> 1142 }
>> 1143 }
>> 1144 else
>> 1145 {
>> 1146 sn=0; // Leave immediately by +ve side
>> 1147 snside=kPY;
>> 1148 }
>> 1149 }
>> 1150 else
>> 1151 {
>> 1152 // Must be || to both
>> 1153 //
>> 1154 sn=kInfinity; // Don't leave by either side
>> 1155 }
>> 1156 }
>> 1157 else if (ss1<=0&&ss2<0)
609 { 1158 {
610 G4double dist = fPlanes[i].a*p.x()+fPlan << 1159 // Outside, in +ve Area
611 if (dist >= -halfCarTolerance) << 1160
>> 1161 if (ds1>0)
612 { 1162 {
613 if (calcNorm) << 1163 sn=0; // Away from shape
>> 1164 // Left by +ve side
>> 1165 snside=kPY;
>> 1166 }
>> 1167 else
>> 1168 {
>> 1169 if (ds2<0)
>> 1170 {
>> 1171 // Ignore +ve plane and use -ve plane intersect
>> 1172 //
>> 1173 sn=ss2/ds2; // Leave by -ve side
>> 1174 snside=kMY;
>> 1175 }
>> 1176 else
614 { 1177 {
615 *validNorm = true; << 1178 // Must be || to both -> exit determined by other axes
616 n->set(fPlanes[i].a, fPlanes[i].b, << 1179 //
>> 1180 sn=kInfinity; // Don't leave by either side
617 } 1181 }
618 return 0; <<
619 } 1182 }
620 G4double tmp = -dist/cosa; << 1183 }
621 if (tmax > tmp) { tmax = tmp; iside = i; << 1184 else if (ss1>0&&ss2>=0)
>> 1185 {
>> 1186 // Outside, in -ve Area
>> 1187 if (ds2<0)
>> 1188 {
>> 1189 sn=0; // away from shape
>> 1190 // Left by -ve side
>> 1191 snside=kMY;
>> 1192 }
>> 1193 else
>> 1194 {
>> 1195 if (ds1>0)
>> 1196 {
>> 1197 // Ignore +ve plane and use -ve plane intersect
>> 1198 //
>> 1199 sn=ss1/ds1; // Leave by +ve side
>> 1200 snside=kPY;
>> 1201 }
>> 1202 else
>> 1203 {
>> 1204 // Must be || to both -> exit determined by other axes
>> 1205 //
>> 1206 sn=kInfinity; // Don't leave by either side
>> 1207 }
>> 1208 }
>> 1209 }
>> 1210
>> 1211 // Update minimum exit distance
>> 1212
>> 1213 if (sn<snxt)
>> 1214 {
>> 1215 snxt=sn;
>> 1216 side=snside;
622 } 1217 }
623 } 1218 }
624 1219
625 // Set normal, if required, and return dista <<
626 // <<
627 if (calcNorm) 1220 if (calcNorm)
628 { 1221 {
629 *validNorm = true; << 1222 switch (side)
630 if (iside < 0) << 1223 {
631 n->set(0, 0, iside + 3); // (-4+3)=-1, ( << 1224 case kPX:
632 else << 1225 cosxz=1.0/std::sqrt(1.0+tanxz*tanxz);
633 n->set(fPlanes[iside].a, fPlanes[iside]. << 1226 *n=G4ThreeVector(cosxz,0,-tanxz*cosxz);
>> 1227 break;
>> 1228 case kMX:
>> 1229 cosxz=-1.0/std::sqrt(1.0+tanxz*tanxz);
>> 1230 *n=G4ThreeVector(cosxz,0,tanxz*cosxz);
>> 1231 break;
>> 1232 case kPY:
>> 1233 cosyz=1.0/std::sqrt(1.0+tanyz*tanyz);
>> 1234 *n=G4ThreeVector(0,cosyz,-tanyz*cosyz);
>> 1235 break;
>> 1236 case kMY:
>> 1237 cosyz=-1.0/std::sqrt(1.0+tanyz*tanyz);
>> 1238 *n=G4ThreeVector(0,cosyz,tanyz*cosyz);
>> 1239 break;
>> 1240 case kPZ:
>> 1241 *n=G4ThreeVector(0,0,1);
>> 1242 break;
>> 1243 case kMZ:
>> 1244 *n=G4ThreeVector(0,0,-1);
>> 1245 break;
>> 1246 default:
>> 1247 DumpInfo();
>> 1248 G4Exception("G4Trd::DistanceToOut(p,v,..)",
>> 1249 "GeomSolids1002", JustWarning,
>> 1250 "Undefined side for valid surface normal to solid.");
>> 1251 break;
>> 1252 }
634 } 1253 }
635 return tmax; << 1254 return snxt;
636 } 1255 }
637 1256
638 ////////////////////////////////////////////// << 1257 ///////////////////////////////////////////////////////////////////////////
639 // 1258 //
640 // Calculate exact shortest distance to any bo 1259 // Calculate exact shortest distance to any boundary from inside
641 // - returns 0 if point is outside << 1260 // - Returns 0 is point outside
642 1261
643 G4double G4Trd::DistanceToOut( const G4ThreeVe 1262 G4double G4Trd::DistanceToOut( const G4ThreeVector& p ) const
644 { 1263 {
>> 1264 G4double safe=0.0;
>> 1265 G4double tanxz,xdist,saf1;
>> 1266 G4double tanyz,ydist,saf2;
>> 1267 G4double zbase;
>> 1268
645 #ifdef G4CSGDEBUG 1269 #ifdef G4CSGDEBUG
646 if( Inside(p) == kOutside ) 1270 if( Inside(p) == kOutside )
647 { 1271 {
648 std::ostringstream message; << 1272 G4int oldprc = G4cout.precision(16) ;
649 G4long oldprc = message.precision(16); << 1273 G4cout << G4endl ;
650 message << "Point p is outside (!?) of sol << 1274 DumpInfo();
651 message << "Position:\n"; << 1275 G4cout << "Position:" << G4endl << G4endl ;
652 message << " p.x() = " << p.x()/mm << " << 1276 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
653 message << " p.y() = " << p.y()/mm << " << 1277 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
654 message << " p.z() = " << p.z()/mm << " << 1278 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
655 G4cout.precision(oldprc); << 1279 G4cout.precision(oldprc) ;
656 G4Exception("G4Trd::DistanceToOut(p)", "Ge << 1280 G4Exception("G4Trd::DistanceToOut(p)", "GeomSolids1002", JustWarning,
657 JustWarning, message ); << 1281 "Point p is outside !?" );
658 DumpInfo(); <<
659 } 1282 }
660 #endif 1283 #endif
661 G4double dx = fPlanes[3].a*std::abs(p.x())+f <<
662 G4double dy = fPlanes[1].b*std::abs(p.y())+f <<
663 G4double dxy = std::max(dx,dy); <<
664 1284
665 G4double dz = std::abs(p.z())-fDz; << 1285 safe=fDz-std::fabs(p.z()); // z perpendicular Dist
666 G4double dist = std::max(dz,dxy); << 1286
>> 1287 zbase=fDz+p.z();
667 1288
668 return (dist < 0) ? -dist : 0.; << 1289 // xdist = distance perpendicular to z axis to closest x plane from p
>> 1290 // = (x half width of shape at p.z) - std::fabs(p.x)
>> 1291 //
>> 1292 tanxz=(fDx2-fDx1)*0.5/fDz;
>> 1293 xdist=fDx1+tanxz*zbase-std::fabs(p.x());
>> 1294 saf1=xdist/std::sqrt(1.0+tanxz*tanxz); // x*std::cos(ang_xz) =
>> 1295 // shortest (perpendicular)
>> 1296 // distance to plane
>> 1297 tanyz=(fDy2-fDy1)*0.5/fDz;
>> 1298 ydist=fDy1+tanyz*zbase-std::fabs(p.y());
>> 1299 saf2=ydist/std::sqrt(1.0+tanyz*tanyz);
>> 1300
>> 1301 // Return minimum x/y/z distance
>> 1302 //
>> 1303 if (safe>saf1) safe=saf1;
>> 1304 if (safe>saf2) safe=saf2;
>> 1305
>> 1306 if (safe<0) safe=0;
>> 1307 return safe;
669 } 1308 }
670 1309
671 ////////////////////////////////////////////// << 1310 ////////////////////////////////////////////////////////////////////////////
672 // 1311 //
673 // GetEntityType << 1312 // Create a List containing the transformed vertices
674 << 1313 // Ordering [0-3] -fDz cross section
675 G4GeometryType G4Trd::GetEntityType() const << 1314 // [4-7] +fDz cross section such that [0] is below [4],
676 { << 1315 // [1] below [5] etc.
677 return {"G4Trd"}; << 1316 // Note:
>> 1317 // Caller has deletion resposibility
>> 1318
>> 1319 G4ThreeVectorList*
>> 1320 G4Trd::CreateRotatedVertices( const G4AffineTransform& pTransform ) const
>> 1321 {
>> 1322 G4ThreeVectorList *vertices;
>> 1323 vertices=new G4ThreeVectorList();
>> 1324 if (vertices)
>> 1325 {
>> 1326 vertices->reserve(8);
>> 1327 G4ThreeVector vertex0(-fDx1,-fDy1,-fDz);
>> 1328 G4ThreeVector vertex1(fDx1,-fDy1,-fDz);
>> 1329 G4ThreeVector vertex2(fDx1,fDy1,-fDz);
>> 1330 G4ThreeVector vertex3(-fDx1,fDy1,-fDz);
>> 1331 G4ThreeVector vertex4(-fDx2,-fDy2,fDz);
>> 1332 G4ThreeVector vertex5(fDx2,-fDy2,fDz);
>> 1333 G4ThreeVector vertex6(fDx2,fDy2,fDz);
>> 1334 G4ThreeVector vertex7(-fDx2,fDy2,fDz);
>> 1335
>> 1336 vertices->push_back(pTransform.TransformPoint(vertex0));
>> 1337 vertices->push_back(pTransform.TransformPoint(vertex1));
>> 1338 vertices->push_back(pTransform.TransformPoint(vertex2));
>> 1339 vertices->push_back(pTransform.TransformPoint(vertex3));
>> 1340 vertices->push_back(pTransform.TransformPoint(vertex4));
>> 1341 vertices->push_back(pTransform.TransformPoint(vertex5));
>> 1342 vertices->push_back(pTransform.TransformPoint(vertex6));
>> 1343 vertices->push_back(pTransform.TransformPoint(vertex7));
>> 1344 }
>> 1345 else
>> 1346 {
>> 1347 DumpInfo();
>> 1348 G4Exception("G4Trd::CreateRotatedVertices()",
>> 1349 "GeomSolids0003", FatalException,
>> 1350 "Error in allocation of vertices. Out of memory !");
>> 1351 }
>> 1352 return vertices;
678 } 1353 }
679 1354
680 ////////////////////////////////////////////// 1355 //////////////////////////////////////////////////////////////////////////
681 // 1356 //
682 // IsFaceted << 1357 // GetEntityType
683 1358
684 G4bool G4Trd::IsFaceted() const << 1359 G4GeometryType G4Trd::GetEntityType() const
685 { 1360 {
686 return true; << 1361 return G4String("G4Trd");
687 } 1362 }
688 1363
689 ////////////////////////////////////////////// 1364 //////////////////////////////////////////////////////////////////////////
690 // 1365 //
691 // Make a clone of the object 1366 // Make a clone of the object
692 // 1367 //
693 G4VSolid* G4Trd::Clone() const 1368 G4VSolid* G4Trd::Clone() const
694 { 1369 {
695 return new G4Trd(*this); 1370 return new G4Trd(*this);
696 } 1371 }
697 1372
698 ////////////////////////////////////////////// 1373 //////////////////////////////////////////////////////////////////////////
699 // 1374 //
700 // Stream object contents to an output stream 1375 // Stream object contents to an output stream
701 1376
702 std::ostream& G4Trd::StreamInfo( std::ostream& 1377 std::ostream& G4Trd::StreamInfo( std::ostream& os ) const
703 { 1378 {
704 G4long oldprc = os.precision(16); << 1379 G4int oldprc = os.precision(16);
705 os << "------------------------------------- 1380 os << "-----------------------------------------------------------\n"
706 << " *** Dump for solid - " << GetName 1381 << " *** Dump for solid - " << GetName() << " ***\n"
707 << " ================================= 1382 << " ===================================================\n"
708 << " Solid type: G4Trd\n" 1383 << " Solid type: G4Trd\n"
709 << " Parameters: \n" 1384 << " Parameters: \n"
710 << " half length X, surface -dZ: " << 1385 << " half length X, surface -dZ: " << fDx1/mm << " mm \n"
711 << " half length X, surface +dZ: " << 1386 << " half length X, surface +dZ: " << fDx2/mm << " mm \n"
712 << " half length Y, surface -dZ: " << 1387 << " half length Y, surface -dZ: " << fDy1/mm << " mm \n"
713 << " half length Y, surface +dZ: " << 1388 << " half length Y, surface +dZ: " << fDy2/mm << " mm \n"
714 << " half length Z : " << << 1389 << " half length Z : " << fDz/mm << " mm \n"
715 << "------------------------------------- 1390 << "-----------------------------------------------------------\n";
716 os.precision(oldprc); 1391 os.precision(oldprc);
717 1392
718 return os; 1393 return os;
719 } 1394 }
720 1395
721 ////////////////////////////////////////////// << 1396
>> 1397 ////////////////////////////////////////////////////////////////////////
722 // 1398 //
723 // Return a point randomly and uniformly selec << 1399 // GetPointOnSurface
>> 1400 //
>> 1401 // Return a point (G4ThreeVector) randomly and uniformly
>> 1402 // selected on the solid surface
724 1403
725 G4ThreeVector G4Trd::GetPointOnSurface() const 1404 G4ThreeVector G4Trd::GetPointOnSurface() const
726 { 1405 {
727 // Set areas << 1406 G4double px, py, pz, tgX, tgY, secX, secY, select, sumS, tmp;
728 // << 1407 G4double Sxy1, Sxy2, Sxy, Sxz, Syz;
729 G4double sxz = (fDx1 + fDx2)*fHx; <<
730 G4double syz = (fDy1 + fDy2)*fHy; <<
731 G4double ssurf[6] = { 4.*fDx1*fDy1, sxz, sxz <<
732 ssurf[1] += ssurf[0]; <<
733 ssurf[2] += ssurf[1]; <<
734 ssurf[3] += ssurf[2]; <<
735 ssurf[4] += ssurf[3]; <<
736 ssurf[5] += ssurf[4]; <<
737 1408
738 // Select face << 1409 tgX = 0.5*(fDx2-fDx1)/fDz;
739 // << 1410 secX = std::sqrt(1+tgX*tgX);
740 G4double select = ssurf[5]*G4QuickRand(); << 1411 tgY = 0.5*(fDy2-fDy1)/fDz;
741 G4int k = 5; << 1412 secY = std::sqrt(1+tgY*tgY);
742 k -= (G4int)(select <= ssurf[4]); << 1413
743 k -= (G4int)(select <= ssurf[3]); << 1414 // calculate 0.25 of side surfaces, sumS is 0.25 of total surface
744 k -= (G4int)(select <= ssurf[2]); << 1415
745 k -= (G4int)(select <= ssurf[1]); << 1416 Sxy1 = fDx1*fDy1;
746 k -= (G4int)(select <= ssurf[0]); << 1417 Sxy2 = fDx2*fDy2;
747 << 1418 Sxy = Sxy1 + Sxy2;
748 // Generate point on selected surface << 1419 Sxz = (fDx1 + fDx2)*fDz*secY;
749 // << 1420 Syz = (fDy1 + fDy2)*fDz*secX;
750 G4double u = G4QuickRand(); << 1421 sumS = Sxy + Sxz + Syz;
751 G4double v = G4QuickRand(); << 1422
752 switch(k) << 1423 select = sumS*G4UniformRand();
>> 1424
>> 1425 if( select < Sxy ) // Sxy1 or Sxy2
753 { 1426 {
754 case 0: // base at -Z << 1427 if( select < Sxy1 )
755 { <<
756 return { (2.*u - 1.)*fDx1, (2.*v - 1.)*f <<
757 } <<
758 case 1: // X face at -Y <<
759 { <<
760 if (u + v > 1.) { u = 1. - u; v = 1. - v <<
761 G4ThreeVector p0(-fDx1,-fDy1,-fDz); <<
762 G4ThreeVector p1( fDx2,-fDy2, fDz); <<
763 return (select <= ssurf[0] + fDx1*fHx) ? <<
764 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect <<
765 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect <<
766 } <<
767 case 2: // X face at +Y <<
768 { 1428 {
769 if (u + v > 1.) { u = 1. - u; v = 1. - v << 1429 pz = -fDz;
770 G4ThreeVector p0( fDx1, fDy1,-fDz); << 1430 px = -fDx1 + 2*fDx1*G4UniformRand();
771 G4ThreeVector p1(-fDx2, fDy2, fDz); << 1431 py = -fDy1 + 2*fDy1*G4UniformRand();
772 return (select <= ssurf[1] + fDx1*fHx) ? <<
773 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect <<
774 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect <<
775 } 1432 }
776 case 3: // Y face at -X << 1433 else
777 { 1434 {
778 if (u + v > 1.) { u = 1. - u; v = 1. - v << 1435 pz = fDz;
779 G4ThreeVector p0(-fDx1, fDy1,-fDz); << 1436 px = -fDx2 + 2*fDx2*G4UniformRand();
780 G4ThreeVector p1(-fDx2,-fDy2, fDz); << 1437 py = -fDy2 + 2*fDy2*G4UniformRand();
781 return (select <= ssurf[2] + fDy1*fHy) ? <<
782 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect <<
783 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect <<
784 } <<
785 case 4: // Y face at +X <<
786 { <<
787 if (u + v > 1.) { u = 1. - u; v = 1. - v <<
788 G4ThreeVector p0( fDx1,-fDy1,-fDz); <<
789 G4ThreeVector p1( fDx2, fDy2, fDz); <<
790 return (select <= ssurf[3] + fDy1*fHy) ? <<
791 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect <<
792 (1. - u - v)*p0 + u*p1 + v*G4ThreeVect <<
793 } <<
794 case 5: // base at +Z <<
795 { <<
796 return { (2.*u - 1.)*fDx2, (2.*v - 1.)*f <<
797 } 1438 }
798 } 1439 }
799 return {0., 0., 0.}; << 1440 else if ( ( select - Sxy ) < Sxz ) // Sxz
>> 1441 {
>> 1442 pz = -fDz + 2*fDz*G4UniformRand();
>> 1443 tmp = fDx1 + (pz + fDz)*tgX;
>> 1444 px = -tmp + 2*tmp*G4UniformRand();
>> 1445 tmp = fDy1 + (pz + fDz)*tgY;
>> 1446
>> 1447 if(G4UniformRand() > 0.5) { py = tmp; }
>> 1448 else { py = -tmp; }
>> 1449 }
>> 1450 else // Syz
>> 1451 {
>> 1452 pz = -fDz + 2*fDz*G4UniformRand();
>> 1453 tmp = fDy1 + (pz + fDz)*tgY;
>> 1454 py = -tmp + 2*tmp*G4UniformRand();
>> 1455 tmp = fDx1 + (pz + fDz)*tgX;
>> 1456
>> 1457 if(G4UniformRand() > 0.5) { px = tmp; }
>> 1458 else { px = -tmp; }
>> 1459 }
>> 1460 return G4ThreeVector(px,py,pz);
800 } 1461 }
801 1462
802 ////////////////////////////////////////////// << 1463 ///////////////////////////////////////////////////////////////////////
803 // 1464 //
804 // Methods for visualisation 1465 // Methods for visualisation
805 1466
806 void G4Trd::DescribeYourselfTo ( G4VGraphicsSc 1467 void G4Trd::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
807 { 1468 {
808 scene.AddSolid (*this); 1469 scene.AddSolid (*this);
809 } 1470 }
810 1471
811 G4Polyhedron* G4Trd::CreatePolyhedron () const 1472 G4Polyhedron* G4Trd::CreatePolyhedron () const
812 { 1473 {
813 return new G4PolyhedronTrd2 (fDx1, fDx2, fDy 1474 return new G4PolyhedronTrd2 (fDx1, fDx2, fDy1, fDy2, fDz);
814 } 1475 }
815 1476
816 #endif << 1477 G4NURBS* G4Trd::CreateNURBS () const
>> 1478 {
>> 1479 // return new G4NURBSbox (fDx, fDy, fDz);
>> 1480 return 0;
>> 1481 }
817 1482