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