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