| Geant4 Cross Reference |
>> 1 // This code implementation is the intellectual property of
>> 2 // the GEANT4 collaboration.
1 // 3 //
2 // ******************************************* << 4 // By copying, distributing or modifying the Program (or any work
3 // * License and Disclaimer << 5 // based on the Program) you indicate your acceptance of this statement,
4 // * << 6 // and all its terms.
5 // * The Geant4 software is copyright of th << 7 //
6 // * the Geant4 Collaboration. It is provided << 8 // $Id: G4Para.cc,v 1.4 1999/12/15 14:50:07 gunter Exp $
7 // * conditions of the Geant4 Software License << 9 // GEANT4 tag $Name: geant4-01-01 $
8 // * LICENSE and available at http://cern.ch/ << 10 //
9 // * include a list of copyright holders. << 11 // class G4Para
10 // * <<
11 // * Neither the authors of this software syst <<
12 // * institutes,nor the agencies providing fin <<
13 // * work make any representation or warran <<
14 // * regarding this software system or assum <<
15 // * use. Please see the license in the file <<
16 // * for the full disclaimer and the limitatio <<
17 // * <<
18 // * This code implementation is the result <<
19 // * technical work of the GEANT4 collaboratio <<
20 // * By using, copying, modifying or distri <<
21 // * any work based on the software) you ag <<
22 // * use in resulting scientific publicati <<
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* <<
25 // 12 //
26 // Implementation for G4Para class 13 // Implementation for G4Para class
27 // 14 //
28 // 21.03.95 P.Kent: Modified for `tolerant' ge << 15 // History:
29 // 31.10.96 V.Grichine: Modifications accordin << 16 // 21.03.95 P.Kent Modified for `tolerant' geom
30 // 28.04.05 V.Grichine: new SurfaceNormal acco << 17 // 31.10.96 V.Grichine Modifications according G4Box/Tubs before to commit
31 // 29.05.17 E.Tcherniaev: complete revision, s << 18 // 18.11.99 V. Grichine , kIndefined was added to ESide
32 ////////////////////////////////////////////// <<
33 19
>> 20 #include <math.h>
34 #include "G4Para.hh" 21 #include "G4Para.hh"
35 <<
36 #if !defined(G4GEOM_USE_UPARA) <<
37 <<
38 #include "G4VoxelLimits.hh" 22 #include "G4VoxelLimits.hh"
39 #include "G4AffineTransform.hh" 23 #include "G4AffineTransform.hh"
40 #include "G4BoundingEnvelope.hh" <<
41 #include "Randomize.hh" <<
42 24
43 #include "G4VPVParameterisation.hh" 25 #include "G4VPVParameterisation.hh"
44 26
45 #include "G4VGraphicsScene.hh" 27 #include "G4VGraphicsScene.hh"
>> 28 #include "G4Polyhedron.hh"
>> 29 #include "G4NURBS.hh"
>> 30 #include "G4NURBSbox.hh"
>> 31 #include "G4VisExtent.hh"
46 32
47 using namespace CLHEP; << 33 // Private enum: Not for external use
>> 34
>> 35 enum ESide {kUndefined,kPX,kMX,kPY,kMY,kPZ,kMZ};
48 36
49 ////////////////////////////////////////////// << 37 // used internally for normal routine
50 // <<
51 // Constructor - set & check half widths <<
52 38
53 G4Para::G4Para(const G4String& pName, << 39 enum ENSide {kNZ,kNX,kNY};
54 G4double pDx, G4double pD <<
55 G4double pAlpha, G4double <<
56 : G4CSGSolid(pName), halfCarTolerance(0.5*kC <<
57 { <<
58 SetAllParameters(pDx, pDy, pDz, pAlpha, pThe <<
59 fRebuildPolyhedron = false; // default valu <<
60 } <<
61 40
62 ////////////////////////////////////////////// << 41 /////////////////////////////////////////////////////////////////////
63 // 42 //
64 // Constructor - design of trapezoid based on << 43 // Constructor - check and set half-widths
65 44
66 G4Para::G4Para( const G4String& pName, << 45 void G4Para::SetAllParameters(G4double pDx, G4double pDy, G4double pDz,
67 const G4ThreeVector pt[8] ) << 46 G4double pAlpha, G4double pTheta, G4double pPhi)
68 : G4CSGSolid(pName), halfCarTolerance(0.5*kC <<
69 { 47 {
70 // Find dimensions and trigonometric values << 48 if (pDx>0&&pDy>0&&pDz>0)
71 // <<
72 fDx = (pt[3].x() - pt[2].x())*0.5; <<
73 fDy = (pt[2].y() - pt[1].y())*0.5; <<
74 fDz = pt[7].z(); <<
75 CheckParameters(); // check dimensions <<
76 <<
77 fTalpha = (pt[2].x() + pt[3].x() - pt[1].x() <<
78 fTthetaCphi = (pt[4].x() + fDy*fTalpha + fDx <<
79 fTthetaSphi = (pt[4].y() + fDy)/fDz; <<
80 MakePlanes(); <<
81 <<
82 // Recompute vertices <<
83 // <<
84 G4ThreeVector v[8]; <<
85 G4double DyTalpha = fDy*fTalpha; <<
86 G4double DzTthetaSphi = fDz*fTthetaSphi; <<
87 G4double DzTthetaCphi = fDz*fTthetaCphi; <<
88 v[0].set(-DzTthetaCphi-DyTalpha-fDx, -DzTthe <<
89 v[1].set(-DzTthetaCphi-DyTalpha+fDx, -DzTthe <<
90 v[2].set(-DzTthetaCphi+DyTalpha-fDx, -DzTthe <<
91 v[3].set(-DzTthetaCphi+DyTalpha+fDx, -DzTthe <<
92 v[4].set( DzTthetaCphi-DyTalpha-fDx, DzTthe <<
93 v[5].set( DzTthetaCphi-DyTalpha+fDx, DzTthe <<
94 v[6].set( DzTthetaCphi+DyTalpha-fDx, DzTthe <<
95 v[7].set( DzTthetaCphi+DyTalpha+fDx, DzTthe <<
96 <<
97 // Compare with original vertices <<
98 // <<
99 for (G4int i=0; i<8; ++i) <<
100 { 49 {
101 G4double delx = std::abs(pt[i].x() - v[i]. << 50 fDx=pDx;
102 G4double dely = std::abs(pt[i].y() - v[i]. << 51 fDy=pDy;
103 G4double delz = std::abs(pt[i].z() - v[i]. << 52 fDz=pDz;
104 G4double discrepancy = std::max(std::max(d << 53 fTalpha=tan(pAlpha);
105 if (discrepancy > 0.1*kCarTolerance) << 54 fTthetaCphi=tan(pTheta)*cos(pPhi);
106 { << 55 fTthetaSphi=tan(pTheta)*sin(pPhi);
107 std::ostringstream message; << 56 }
108 G4long oldprc = message.precision(16); << 57 else
109 message << "Invalid vertice coordinates << 58 {
110 << "\nVertix #" << i << ", discr << 59 G4Exception("Error in G4Para::SetAllParameters - Invalid Length Parameters");
111 << "\n original : " << pt[i] <<
112 << "\n recomputed : " << v[i]; <<
113 G4cout.precision(oldprc); <<
114 G4Exception("G4Para::G4Para()", "GeomSol <<
115 FatalException, message); <<
116 <<
117 } <<
118 } 60 }
119 } 61 }
120 62
121 ////////////////////////////////////////////// << 63 ///////////////////////////////////////////////////////////////////////////
122 // <<
123 // Fake default constructor - sets only member <<
124 // for usage restri <<
125 <<
126 G4Para::G4Para( __void__& a ) <<
127 : G4CSGSolid(a), halfCarTolerance(0.5*kCarTo <<
128 { <<
129 SetAllParameters(1., 1., 1., 0., 0., 0.); <<
130 fRebuildPolyhedron = false; // default value <<
131 } <<
132 <<
133 ////////////////////////////////////////////// <<
134 // <<
135 // Destructor <<
136 <<
137 G4Para::~G4Para() = default; <<
138 <<
139 ////////////////////////////////////////////// <<
140 // <<
141 // Copy constructor <<
142 <<
143 G4Para::G4Para(const G4Para& rhs) <<
144 : G4CSGSolid(rhs), halfCarTolerance(rhs.half <<
145 fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz), <<
146 fTthetaCphi(rhs.fTthetaCphi),fTthetaSphi(r <<
147 { <<
148 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs <<
149 } <<
150 <<
151 ////////////////////////////////////////////// <<
152 // 64 //
153 // Assignment operator <<
154 65
155 G4Para& G4Para::operator = (const G4Para& rhs) << 66 G4Para::G4Para(const G4String& pName,G4double pDx, G4double pDy, G4double pDz,
>> 67 G4double pAlpha, G4double pTheta,G4double pPhi) : G4CSGSolid(pName)
156 { 68 {
157 // Check assignment to self << 69 if (pDx>0&&pDy>0&&pDz>0)
158 // << 70 {
159 if (this == &rhs) { return *this; } << 71 SetAllParameters( pDx, pDy, pDz, pAlpha, pTheta, pPhi);
160 << 72 }
161 // Copy base class data << 73 else
162 // << 74 {
163 G4CSGSolid::operator=(rhs); << 75 G4Exception("Error in G4Para::G4Para - Invalid Length Parameters");
164 << 76 }
165 // Copy data <<
166 // <<
167 halfCarTolerance = rhs.halfCarTolerance; <<
168 fDx = rhs.fDx; <<
169 fDy = rhs.fDy; <<
170 fDz = rhs.fDz; <<
171 fTalpha = rhs.fTalpha; <<
172 fTthetaCphi = rhs.fTthetaCphi; <<
173 fTthetaSphi = rhs.fTthetaSphi; <<
174 for (G4int i=0; i<4; ++i) { fPlanes[i] = rh <<
175 <<
176 return *this; <<
177 } 77 }
178 78
179 ////////////////////////////////////////////// << 79 ////////////////////////////////////////////////////////////////////////
180 // 80 //
181 // Set all parameters, as for constructor - se << 81 // Constructor - Design of trapezoid based on 8 G4ThreeVector parameters,
182 << 82 // which are its vertices. Checking of planarity with preparation of
183 void G4Para::SetAllParameters(G4double pDx, G4 << 83 // fPlanes[] and than calculation of other members
184 G4double pAlpha, <<
185 { <<
186 // Reset data of the base class <<
187 fCubicVolume = 0; <<
188 fSurfaceArea = 0; <<
189 fRebuildPolyhedron = true; <<
190 <<
191 // Set parameters <<
192 fDx = pDx; <<
193 fDy = pDy; <<
194 fDz = pDz; <<
195 fTalpha = std::tan(pAlpha); <<
196 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi <<
197 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi <<
198 84
199 CheckParameters(); <<
200 MakePlanes(); <<
201 } <<
202 85
203 ////////////////////////////////////////////// << 86 G4Para::G4Para( const G4String& pName,
204 // << 87 const G4ThreeVector pt[8]) : G4CSGSolid(pName)
205 // Check dimensions <<
206 <<
207 void G4Para::CheckParameters() <<
208 { 88 {
209 if (fDx < 2*kCarTolerance || << 89 if ( pt[0].z()<0 && pt[0].z()==pt[1].z() && pt[0].z()==pt[2].z() && pt[0].z()==pt[3].z()
210 fDy < 2*kCarTolerance || << 90 && pt[4].z()>0 && pt[4].z()==pt[5].z() && pt[4].z()==pt[6].z() && pt[4].z()==pt[7].z()
211 fDz < 2*kCarTolerance) << 91 && (pt[0].z()+pt[4].z())== 0
212 { << 92 && pt[0].y()==pt[1].y() && pt[2].y()==pt[3].y()
213 std::ostringstream message; << 93 && pt[4].y()==pt[5].y() && pt[6].y()==pt[7].y()
214 message << "Invalid (too small or negative << 94 && (pt[0].y()+pt[2].y()+pt[4].y()+pt[6].y())==0 )
215 << GetName() << 95 {
216 << "\n X - " << fDx << 96 fDz = (pt[7]).z() ;
217 << "\n Y - " << fDy << 97
218 << "\n Z - " << fDz; << 98 fDy = ((pt[2]).y()-(pt[1]).y())*0.5 ;
219 G4Exception("G4Para::CheckParameters()", " << 99 fDx = ((pt[1]).x()-(pt[0]).x())*0.5 ;
220 FatalException, message); << 100 fDx = ((pt[3]).x()-(pt[2]).x())*0.5 ;
221 } << 101 fTalpha = ((pt[2]).x()+(pt[3]).x()-(pt[1]).x()-(pt[0]).x())*0.25/fDy ;
222 } << 102
>> 103 // fDy = ((pt[6]).y()-(pt[5]).y())*0.5 ;
>> 104 // fDx = ((pt[5]).x()-(pt[4]).x())*0.5 ;
>> 105 // fDx = ((pt[7]).x()-(pt[6]).x())*0.5 ;
>> 106 // fTalpha = ((pt[6]).x()+(pt[7]).x()-(pt[5]).x()-(pt[4]).x())*0.25/fDy ;
223 107
224 ////////////////////////////////////////////// << 108 fTthetaCphi = ((pt[4]).x()+fDy*fTalpha+fDx)/fDz ;
225 // << 109 fTthetaSphi = ((pt[4]).y()+fDy)/fDz ;
226 // Set side planes <<
227 <<
228 void G4Para::MakePlanes() <<
229 { <<
230 G4ThreeVector vx(1, 0, 0); <<
231 G4ThreeVector vy(fTalpha, 1, 0); <<
232 G4ThreeVector vz(fTthetaCphi, fTthetaSphi, 1 <<
233 <<
234 // Set -Y & +Y planes <<
235 // <<
236 G4ThreeVector ynorm = (vx.cross(vz)).unit(); <<
237 <<
238 fPlanes[0].a = 0.; <<
239 fPlanes[0].b = ynorm.y(); <<
240 fPlanes[0].c = ynorm.z(); <<
241 fPlanes[0].d = fPlanes[0].b*fDy; // point (0 <<
242 <<
243 fPlanes[1].a = 0.; <<
244 fPlanes[1].b = -fPlanes[0].b; <<
245 fPlanes[1].c = -fPlanes[0].c; <<
246 fPlanes[1].d = fPlanes[0].d; <<
247 <<
248 // Set -X & +X planes <<
249 // <<
250 G4ThreeVector xnorm = (vz.cross(vy)).unit(); <<
251 <<
252 fPlanes[2].a = xnorm.x(); <<
253 fPlanes[2].b = xnorm.y(); <<
254 fPlanes[2].c = xnorm.z(); <<
255 fPlanes[2].d = fPlanes[2].a*fDx; // point (f <<
256 <<
257 fPlanes[3].a = -fPlanes[2].a; <<
258 fPlanes[3].b = -fPlanes[2].b; <<
259 fPlanes[3].c = -fPlanes[2].c; <<
260 fPlanes[3].d = fPlanes[2].d; <<
261 } <<
262 110
263 ////////////////////////////////////////////// << 111 }
264 // << 112 else
265 // Get volume << 113 {
>> 114 G4Exception("Error in G4Para::G4Para - Invalid vertice coordinates");
>> 115 }
266 116
267 G4double G4Para::GetCubicVolume() << 117
268 { <<
269 // It is like G4Box, since para transformati <<
270 if (fCubicVolume == 0) <<
271 { <<
272 fCubicVolume = 8*fDx*fDy*fDz; <<
273 } <<
274 return fCubicVolume; <<
275 } 118 }
276 119
277 ////////////////////////////////////////////// 120 //////////////////////////////////////////////////////////////////////////
278 // 121 //
279 // Get surface area <<
280 122
281 G4double G4Para::GetSurfaceArea() << 123 G4Para::~G4Para()
282 { 124 {
283 if(fSurfaceArea == 0) << 125 ;
284 { <<
285 G4ThreeVector vx(fDx, 0, 0); <<
286 G4ThreeVector vy(fDy*fTalpha, fDy, 0); <<
287 G4ThreeVector vz(fDz*fTthetaCphi, fDz*fTth <<
288 <<
289 G4double sxy = fDx*fDy; // (vx.cross(vy)). <<
290 G4double sxz = (vx.cross(vz)).mag(); <<
291 G4double syz = (vy.cross(vz)).mag(); <<
292 <<
293 fSurfaceArea = 8*(sxy+sxz+syz); <<
294 } <<
295 return fSurfaceArea; <<
296 } 126 }
297 127
298 ////////////////////////////////////////////// 128 //////////////////////////////////////////////////////////////////////////
299 // 129 //
300 // Dispatch to parameterisation for replicatio 130 // Dispatch to parameterisation for replication mechanism dimension
301 // computation & modification << 131 // computation & modification.
302 132
303 void G4Para::ComputeDimensions( G4VPVPara << 133 void G4Para::ComputeDimensions(G4VPVParameterisation* p,
304 const G4int n, 134 const G4int n,
305 const G4VPhysi << 135 const G4VPhysicalVolume* pRep)
306 { 136 {
307 p->ComputeDimensions(*this,n,pRep); << 137 p->ComputeDimensions(*this,n,pRep);
308 } 138 }
309 139
310 ////////////////////////////////////////////// <<
311 // <<
312 // Get bounding box <<
313 <<
314 void G4Para::BoundingLimits(G4ThreeVector& pMi <<
315 { <<
316 G4double dz = GetZHalfLength(); <<
317 G4double dx = GetXHalfLength(); <<
318 G4double dy = GetYHalfLength(); <<
319 <<
320 G4double x0 = dz*fTthetaCphi; <<
321 G4double x1 = dy*GetTanAlpha(); <<
322 G4double xmin = <<
323 std::min( <<
324 std::min( <<
325 std::min(-x0-x1-dx,-x0+x1-dx),x0-x1-dx),x0 <<
326 G4double xmax = <<
327 std::max( <<
328 std::max( <<
329 std::max(-x0-x1+dx,-x0+x1+dx),x0-x1+dx),x0 <<
330 <<
331 G4double y0 = dz*fTthetaSphi; <<
332 G4double ymin = std::min(-y0-dy,y0-dy); <<
333 G4double ymax = std::max(-y0+dy,y0+dy); <<
334 <<
335 pMin.set(xmin,ymin,-dz); <<
336 pMax.set(xmax,ymax, dz); <<
337 <<
338 // Check correctness of the bounding box <<
339 // <<
340 if (pMin.x() >= pMax.x() || pMin.y() >= pMax <<
341 { <<
342 std::ostringstream message; <<
343 message << "Bad bounding box (min >= max) <<
344 << GetName() << " !" <<
345 << "\npMin = " << pMin <<
346 << "\npMax = " << pMax; <<
347 G4Exception("G4Para::BoundingLimits()", "G <<
348 JustWarning, message); <<
349 DumpInfo(); <<
350 } <<
351 } <<
352 140
353 ////////////////////////////////////////////// << 141 //////////////////////////////////////////////////////////////
354 // 142 //
355 // Calculate extent under transform and specif 143 // Calculate extent under transform and specified limit
356 144
357 G4bool G4Para::CalculateExtent( const EAxis pA << 145 G4bool G4Para::CalculateExtent(const EAxis pAxis,
358 const G4VoxelL << 146 const G4VoxelLimits& pVoxelLimit,
359 const G4Affine << 147 const G4AffineTransform& pTransform,
360 G4double& << 148 G4double& pMin, G4double& pMax) const
361 { << 149 {
362 G4ThreeVector bmin, bmax; << 150 G4bool flag;
363 G4bool exist; << 151
364 << 152 if (!pTransform.IsRotated())
365 // Check bounding box (bbox) << 153 {
366 // << 154 // Special case handling for unrotated trapezoids
367 BoundingLimits(bmin,bmax); << 155 // Compute z/x/y/ mins and maxs respecting limits, with early returns
368 G4BoundingEnvelope bbox(bmin,bmax); << 156 // if outside limits. Then switch() on pAxis
369 #ifdef G4BBOX_EXTENT << 157 G4int i ;
370 return bbox.CalculateExtent(pAxis,pVoxelLimi << 158 G4double xoffset,xMin,xMax;
371 #endif << 159 G4double yoffset,yMin,yMax;
372 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 160 G4double zoffset,zMin,zMax;
373 { << 161 G4double temp[8] ; // some points for intersection with zMin/zMax
374 return exist = pMin < pMax; << 162
375 } << 163 xoffset=pTransform.NetTranslation().x();
376 << 164 yoffset=pTransform.NetTranslation().y();
377 // Set bounding envelope (benv) and calculat << 165 zoffset=pTransform.NetTranslation().z();
378 // << 166
379 G4double dz = GetZHalfLength(); << 167 G4ThreeVector pt[8]; // vertices after translation
380 G4double dx = GetXHalfLength(); << 168 pt[0]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy*fTalpha-fDx,
381 G4double dy = GetYHalfLength(); << 169 yoffset-fDz*fTthetaSphi-fDy,zoffset-fDz);
382 << 170 pt[1]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy*fTalpha+fDx,
383 G4double x0 = dz*fTthetaCphi; << 171 yoffset-fDz*fTthetaSphi-fDy,zoffset-fDz);
384 G4double x1 = dy*GetTanAlpha(); << 172 pt[2]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy*fTalpha-fDx,
385 G4double y0 = dz*fTthetaSphi; << 173 yoffset-fDz*fTthetaSphi+fDy,zoffset-fDz);
386 << 174 pt[3]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy*fTalpha+fDx,
387 G4ThreeVectorList baseA(4), baseB(4); << 175 yoffset-fDz*fTthetaSphi+fDy,zoffset-fDz);
388 baseA[0].set(-x0-x1-dx,-y0-dy,-dz); << 176 pt[4]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy*fTalpha-fDx,
389 baseA[1].set(-x0-x1+dx,-y0-dy,-dz); << 177 yoffset+fDz*fTthetaSphi-fDy,zoffset+fDz);
390 baseA[2].set(-x0+x1+dx,-y0+dy,-dz); << 178 pt[5]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy*fTalpha+fDx,
391 baseA[3].set(-x0+x1-dx,-y0+dy,-dz); << 179 yoffset+fDz*fTthetaSphi-fDy,zoffset+fDz);
392 << 180 pt[6]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy*fTalpha-fDx,
393 baseB[0].set(+x0-x1-dx, y0-dy, dz); << 181 yoffset+fDz*fTthetaSphi+fDy,zoffset+fDz);
394 baseB[1].set(+x0-x1+dx, y0-dy, dz); << 182 pt[7]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy*fTalpha+fDx,
395 baseB[2].set(+x0+x1+dx, y0+dy, dz); << 183 yoffset+fDz*fTthetaSphi+fDy,zoffset+fDz);
396 baseB[3].set(+x0+x1-dx, y0+dy, dz); << 184 zMin=zoffset-fDz;
397 << 185 zMax=zoffset+fDz;
398 std::vector<const G4ThreeVectorList *> polyg << 186 if (pVoxelLimit.IsZLimited())
399 polygons[0] = &baseA; << 187 {
400 polygons[1] = &baseB; << 188 if (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance
401 << 189 ||zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance)
402 G4BoundingEnvelope benv(bmin,bmax,polygons); << 190 {
403 exist = benv.CalculateExtent(pAxis,pVoxelLim << 191 return false;
404 return exist; << 192 }
405 } << 193 else
406 << 194 {
407 ////////////////////////////////////////////// << 195 if (zMin<pVoxelLimit.GetMinZExtent())
408 // << 196 {
409 // Determine where is point p, inside/on_surfa << 197 zMin=pVoxelLimit.GetMinZExtent();
410 // << 198 }
411 << 199 if (zMax>pVoxelLimit.GetMaxZExtent())
412 EInside G4Para::Inside( const G4ThreeVector& p << 200 {
413 { << 201 zMax=pVoxelLimit.GetMaxZExtent();
414 G4double xx = fPlanes[2].a*p.x()+fPlanes[2]. << 202 }
415 G4double dx = std::abs(xx) + fPlanes[2].d; << 203 }
>> 204 }
>> 205
>> 206 temp[0] = pt[0].y()+(pt[4].y()-pt[0].y())*(zMin-pt[0].z())/(pt[4].z()-pt[0].z()) ;
>> 207 temp[1] = pt[0].y()+(pt[4].y()-pt[0].y())*(zMax-pt[0].z())/(pt[4].z()-pt[0].z()) ;
>> 208 temp[2] = pt[2].y()+(pt[6].y()-pt[2].y())*(zMin-pt[2].z())/(pt[6].z()-pt[2].z()) ;
>> 209 temp[3] = pt[2].y()+(pt[6].y()-pt[2].y())*(zMax-pt[2].z())/(pt[6].z()-pt[2].z()) ;
>> 210 yMax = yoffset - fabs(fDz*fTthetaSphi) - fDy - fDy ;
>> 211 yMin = -yMax ;
>> 212 for(i=0;i<4;i++)
>> 213 {
>> 214 if(temp[i] > yMax) yMax = temp[i] ;
>> 215 if(temp[i] < yMin) yMin = temp[i] ;
>> 216 }
>> 217
>> 218 if (pVoxelLimit.IsYLimited())
>> 219 {
>> 220 if (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance
>> 221 ||yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance)
>> 222 {
>> 223 return false;
>> 224 }
>> 225 else
>> 226 {
>> 227 if (yMin<pVoxelLimit.GetMinYExtent())
>> 228 {
>> 229 yMin=pVoxelLimit.GetMinYExtent();
>> 230 }
>> 231 if (yMax>pVoxelLimit.GetMaxYExtent())
>> 232 {
>> 233 yMax=pVoxelLimit.GetMaxYExtent();
>> 234 }
>> 235 }
>> 236 }
>> 237
>> 238 temp[0] = pt[0].x()+(pt[4].x()-pt[0].x())*(zMin-pt[0].z())/(pt[4].z()-pt[0].z()) ;
>> 239 temp[1] = pt[0].x()+(pt[4].x()-pt[0].x())*(zMax-pt[0].z())/(pt[4].z()-pt[0].z()) ;
>> 240 temp[2] = pt[2].x()+(pt[6].x()-pt[2].x())*(zMin-pt[2].z())/(pt[6].z()-pt[2].z()) ;
>> 241 temp[3] = pt[2].x()+(pt[6].x()-pt[2].x())*(zMax-pt[2].z())/(pt[6].z()-pt[2].z()) ;
>> 242 temp[4] = pt[3].x()+(pt[7].x()-pt[3].x())*(zMin-pt[3].z())/(pt[7].z()-pt[3].z()) ;
>> 243 temp[5] = pt[3].x()+(pt[7].x()-pt[3].x())*(zMax-pt[3].z())/(pt[7].z()-pt[3].z()) ;
>> 244 temp[6] = pt[1].x()+(pt[5].x()-pt[1].x())*(zMin-pt[1].z())/(pt[5].z()-pt[1].z()) ;
>> 245 temp[7] = pt[1].x()+(pt[5].x()-pt[1].x())*(zMax-pt[1].z())/(pt[5].z()-pt[1].z()) ;
>> 246
>> 247 xMax = xoffset - fabs(fDz*fTthetaCphi) - fDx - fDx -fDx - fDx;
>> 248 xMin = -xMax ;
>> 249 for(i=0;i<8;i++)
>> 250 {
>> 251 if(temp[i] > xMax) xMax = temp[i] ;
>> 252 if(temp[i] < xMin) xMin = temp[i] ;
>> 253 }
>> 254 // xMax/Min = f(yMax/Min) ?
>> 255 if (pVoxelLimit.IsXLimited())
>> 256 {
>> 257 if (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance
>> 258 ||xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance)
>> 259 {
>> 260 return false;
>> 261 }
>> 262 else
>> 263 {
>> 264 if (xMin<pVoxelLimit.GetMinXExtent())
>> 265 {
>> 266 xMin=pVoxelLimit.GetMinXExtent();
>> 267 }
>> 268 if (xMax>pVoxelLimit.GetMaxXExtent())
>> 269 {
>> 270 xMax=pVoxelLimit.GetMaxXExtent();
>> 271 }
>> 272 }
>> 273 }
>> 274
>> 275 switch (pAxis)
>> 276 {
>> 277 case kXAxis:
>> 278 pMin=xMin;
>> 279 pMax=xMax;
>> 280 break;
>> 281 case kYAxis:
>> 282 pMin=yMin;
>> 283 pMax=yMax;
>> 284 break;
>> 285 case kZAxis:
>> 286 pMin=zMin;
>> 287 pMax=zMax;
>> 288 break;
>> 289 }
416 290
417 G4double yy = fPlanes[0].b*p.y()+fPlanes[0]. << 291 pMin-=kCarTolerance;
418 G4double dy = std::abs(yy) + fPlanes[0].d; << 292 pMax+=kCarTolerance;
419 G4double dxy = std::max(dx,dy); <<
420 293
421 G4double dz = std::abs(p.z())-fDz; << 294 flag = true;
422 G4double dist = std::max(dxy,dz); << 295 }
>> 296 else
>> 297 {
>> 298 // General rotated case - create and clip mesh to boundaries
423 299
424 if (dist > halfCarTolerance) return kOutside << 300 G4bool existsAfterClip=false;
425 return (dist > -halfCarTolerance) ? kSurface << 301 G4ThreeVectorList *vertices;
426 } <<
427 302
428 ////////////////////////////////////////////// << 303 pMin=+kInfinity;
>> 304 pMax=-kInfinity;
>> 305 // Calculate rotated vertex coordinates
>> 306
>> 307 vertices=CreateRotatedVertices(pTransform);
>> 308 ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
>> 309 ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax);
>> 310 ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
>> 311
>> 312 if (pMin!=kInfinity||pMax!=-kInfinity)
>> 313 {
>> 314 existsAfterClip=true;
>> 315
>> 316 // Add 2*tolerance to avoid precision troubles
>> 317 pMin-=kCarTolerance;
>> 318 pMax+=kCarTolerance;
>> 319
>> 320 }
>> 321 else
>> 322 {
>> 323 // Check for case where completely enveloping clipping volume
>> 324 // If point inside then we are confident that the solid completely
>> 325 // envelopes the clipping volume. Hence set min/max extents according
>> 326 // to clipping volume extents along the specified axis.
>> 327
>> 328 G4ThreeVector clipCentre(
>> 329 (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
>> 330 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
>> 331 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
>> 332
>> 333 if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
>> 334 {
>> 335 existsAfterClip=true;
>> 336 pMin=pVoxelLimit.GetMinExtent(pAxis);
>> 337 pMax=pVoxelLimit.GetMaxExtent(pAxis);
>> 338 }
>> 339 }
>> 340 delete vertices ; // 'new' in the function called
>> 341 flag = existsAfterClip ;
>> 342 }
>> 343 return flag;
>> 344 }
>> 345
>> 346 /////////////////////////////////////////////////////////////////////////////
>> 347 //
>> 348 // Check in p is inside/on surface/outside solid
>> 349
>> 350 EInside G4Para::Inside(const G4ThreeVector& p) const
>> 351 {
>> 352 EInside in=kOutside;
>> 353 G4double xt,yt;
>> 354 // Check
>> 355 if (fabs(p.z())<=fDz-kCarTolerance/2)
>> 356 {
>> 357 yt=fabs(p.y()-fTthetaSphi*p.z());
>> 358 if (yt<=fDy-kCarTolerance/2)
>> 359 {
>> 360 xt=fabs(p.x()-fTthetaCphi*p.z()-fTalpha*yt);
>> 361 if (xt<=fDx-kCarTolerance/2)
>> 362 {
>> 363 in=kInside;
>> 364 }
>> 365 else if (xt<=fDx+kCarTolerance/2)
>> 366 {
>> 367 in=kSurface;
>> 368 }
>> 369 }
>> 370 else if (yt<=fDy+kCarTolerance/2)
>> 371 {
>> 372 xt=fabs(p.x()-fTthetaCphi*p.z()-fTalpha*yt);
>> 373 if (xt<=fDx+kCarTolerance/2)
>> 374 {
>> 375 in=kSurface;
>> 376 }
>> 377 }
>> 378 }
>> 379 else if (fabs(p.z())<=fDz+kCarTolerance/2)
>> 380 {
>> 381 yt=fabs(p.y()-fTthetaSphi*p.z());
>> 382 if (yt<=fDy+kCarTolerance/2)
>> 383 {
>> 384 xt=fabs(p.x()-fTthetaCphi*p.z()-fTalpha*yt);
>> 385 if (xt<=fDx+kCarTolerance/2)
>> 386 {
>> 387 in=kSurface;
>> 388 }
>> 389 }
>> 390 }
>> 391 return in;
>> 392 }
>> 393
>> 394 ///////////////////////////////////////////////////////////////////////////
>> 395 //
>> 396 // Calculate side nearest to p, and return normal
>> 397 // If 2+ sides equidistant, first side's normal returned (arbitrarily)
>> 398
>> 399 G4ThreeVector G4Para::SurfaceNormal( const G4ThreeVector& p) const
>> 400 {
>> 401 ENSide side;
>> 402 G4ThreeVector norm;
>> 403 G4double distx,disty,distz;
>> 404 G4double newpx,newpy,xshift;
>> 405 G4double calpha,salpha; // Sin/Cos(alpha) - needed to recalc G4Parameter
>> 406 G4double tntheta,cosntheta; // tan and cos of normal's theta component
>> 407 G4double ycomp;
>> 408
>> 409 newpx=p.x()-fTthetaCphi*p.z();
>> 410 newpy=p.y()-fTthetaSphi*p.z();
>> 411
>> 412 calpha=1/sqrt(1+fTalpha*fTalpha);
>> 413 if (fTalpha)
>> 414 {
>> 415 salpha=-calpha/fTalpha; // NOTE: actually use MINUS sin(alpha)
>> 416 }
>> 417 else
>> 418 {
>> 419 salpha=0;
>> 420 }
>> 421
>> 422 xshift=newpx*calpha+newpy*salpha;
>> 423
>> 424 distx=fabs(fabs(xshift)-fDx*calpha);
>> 425 disty=fabs(fabs(newpy)-fDy);
>> 426 distz=fabs(fabs(p.z())-fDz);
>> 427
>> 428 if (distx<disty)
>> 429 {
>> 430 if (distx<distz) side=kNX;
>> 431 else side=kNZ;
>> 432 }
>> 433 else
>> 434 {
>> 435 if (disty<distz) side=kNY;
>> 436 else side=kNZ;
>> 437 }
>> 438
>> 439 switch (side)
>> 440 {
>> 441 case kNX:
>> 442 tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
>> 443 if (xshift<0)
>> 444 {
>> 445 cosntheta=-1/sqrt(1+tntheta*tntheta);
>> 446 }
>> 447 else
>> 448 {
>> 449 cosntheta=1/sqrt(1+tntheta*tntheta);
>> 450 }
>> 451 norm=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
>> 452 break;
>> 453 case kNY:
>> 454 if (newpy<0)
>> 455 {
>> 456 ycomp=-1/sqrt(1+fTthetaSphi*fTthetaSphi);
>> 457 }
>> 458 else
>> 459 {
>> 460 ycomp=1/sqrt(1+fTthetaSphi*fTthetaSphi);
>> 461 }
>> 462 norm=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
>> 463 break;
>> 464 case kNZ:
>> 465 // Closest to Z
>> 466 if (p.z()>=0)
>> 467 {
>> 468 norm=G4ThreeVector(0,0,1);
>> 469 }
>> 470 else
>> 471 {
>> 472 norm=G4ThreeVector(0,0,-1);
>> 473 }
>> 474 break;
>> 475 }
>> 476 return norm;
>> 477 }
>> 478
>> 479 //////////////////////////////////////////////////////////////////////////////
>> 480 //
>> 481 // Calculate distance to shape from outside - return kInfinity if no intersection
>> 482 //
>> 483 // ALGORITHM:
>> 484 // For each component, calculate pair of minimum and maximum intersection
>> 485 // values for which the particle is in the extent of the shape
>> 486 // - The smallest (MAX minimum) allowed distance of the pairs is intersect
>> 487 // - Z plane intersectin uses tolerance
>> 488 // - XZ YZ planes use logic & *SLIGHTLY INCORRECT* tolerance
>> 489 // (this saves at least 1 sqrt, 1 multiply and 1 divide... in applicable
>> 490 // cases)
>> 491 // - Note: XZ and YZ planes each divide space into four regions,
>> 492 // characterised by ss1 ss2
>> 493
>> 494 G4double G4Para::DistanceToIn(const G4ThreeVector& p,const G4ThreeVector& v) const
>> 495 {
>> 496 G4double snxt; // snxt = default return value
>> 497 G4double smin,smax;
>> 498 G4double tmin,tmax;
>> 499 G4double yt,vy,xt,vx;
>> 500 G4double max;
>> 501 //
>> 502 // Z Intersection range
>> 503 //
>> 504 if (v.z()>0)
>> 505 {
>> 506 max=fDz-p.z();
>> 507 if (max>kCarTolerance/2)
>> 508 {
>> 509 smax=max/v.z();
>> 510 smin=(-fDz-p.z())/v.z();
>> 511 }
>> 512 else
>> 513 {
>> 514 return snxt=kInfinity;
>> 515 }
>> 516 }
>> 517 else if (v.z()<0)
>> 518 {
>> 519 max=-fDz-p.z();
>> 520 if (max<-kCarTolerance/2)
>> 521 {
>> 522 smax=max/v.z();
>> 523 smin=(fDz-p.z())/v.z();
>> 524 }
>> 525 else
>> 526 {
>> 527 return snxt=kInfinity;
>> 528 }
>> 529 }
>> 530 else
>> 531 {
>> 532 if (fabs(p.z())<=fDz) // Inside
>> 533 {
>> 534 smin=0;
>> 535 smax=kInfinity;
>> 536 }
>> 537 else
>> 538 {
>> 539 return snxt=kInfinity;
>> 540 }
>> 541 }
>> 542
429 // 543 //
430 // Determine side where point is, and return c << 544 // Y G4Parallel planes intersection
431 <<
432 G4ThreeVector G4Para::SurfaceNormal( const G4T <<
433 { <<
434 G4int nsurf = 0; // number of surfaces where <<
435 <<
436 // Check Z faces <<
437 // <<
438 G4double nz = 0; <<
439 G4double dz = std::abs(p.z()) - fDz; <<
440 if (std::abs(dz) <= halfCarTolerance) <<
441 { <<
442 nz = (p.z() < 0) ? -1 : 1; <<
443 ++nsurf; <<
444 } <<
445 <<
446 // Check Y faces <<
447 // <<
448 G4double ny = 0; <<
449 G4double yy = fPlanes[0].b*p.y()+fPlanes[0]. <<
450 if (std::abs(fPlanes[0].d + yy) <= halfCarTo <<
451 { <<
452 ny = fPlanes[0].b; <<
453 nz += fPlanes[0].c; <<
454 ++nsurf; <<
455 } <<
456 else if (std::abs(fPlanes[1].d - yy) <= half <<
457 { <<
458 ny = fPlanes[1].b; <<
459 nz += fPlanes[1].c; <<
460 ++nsurf; <<
461 } <<
462 <<
463 // Check X faces <<
464 // <<
465 G4double nx = 0; <<
466 G4double xx = fPlanes[2].a*p.x()+fPlanes[2]. <<
467 if (std::abs(fPlanes[2].d + xx) <= halfCarTo <<
468 { <<
469 nx = fPlanes[2].a; <<
470 ny += fPlanes[2].b; <<
471 nz += fPlanes[2].c; <<
472 ++nsurf; <<
473 } <<
474 else if (std::abs(fPlanes[3].d - xx) <= half <<
475 { <<
476 nx = fPlanes[3].a; <<
477 ny += fPlanes[3].b; <<
478 nz += fPlanes[3].c; <<
479 ++nsurf; <<
480 } <<
481 <<
482 // Return normal <<
483 // <<
484 if (nsurf == 1) return {nx,ny,nz}; <<
485 else if (nsurf != 0) return G4ThreeVector(nx <<
486 else <<
487 { <<
488 // Point is not on the surface <<
489 // <<
490 #ifdef G4CSGDEBUG <<
491 std::ostringstream message; <<
492 G4int oldprc = message.precision(16); <<
493 message << "Point p is not on surface (!?) <<
494 << GetName() << G4endl; <<
495 message << "Position:\n"; <<
496 message << " p.x() = " << p.x()/mm << " <<
497 message << " p.y() = " << p.y()/mm << " <<
498 message << " p.z() = " << p.z()/mm << " <<
499 G4cout.precision(oldprc) ; <<
500 G4Exception("G4Para::SurfaceNormal(p)", "G <<
501 JustWarning, message ); <<
502 DumpInfo(); <<
503 #endif <<
504 return ApproxSurfaceNormal(p); <<
505 } <<
506 } <<
507 <<
508 ////////////////////////////////////////////// <<
509 // 545 //
510 // Algorithm for SurfaceNormal() following the << 546 yt=p.y()-fTthetaSphi*p.z();
511 // for points not on the surface << 547 vy=v.y()-fTthetaSphi*v.z();
512 <<
513 G4ThreeVector G4Para::ApproxSurfaceNormal( con <<
514 { <<
515 G4double dist = -DBL_MAX; <<
516 G4int iside = 0; <<
517 for (G4int i=0; i<4; ++i) <<
518 { <<
519 G4double d = fPlanes[i].a*p.x() + <<
520 fPlanes[i].b*p.y() + <<
521 fPlanes[i].c*p.z() + fPlanes[ <<
522 if (d > dist) { dist = d; iside = i; } <<
523 } <<
524 548
525 G4double distz = std::abs(p.z()) - fDz; << 549 if (vy>0)
526 if (dist > distz) << 550 {
527 return { fPlanes[iside].a, fPlanes[iside]. << 551 max=fDy-yt;
528 else << 552 if (max>kCarTolerance/2)
529 return { 0, 0, (G4double)((p.z() < 0) ? -1 << 553 {
530 } << 554 tmax=max/vy;
531 << 555 tmin=(-fDy-yt)/vy;
532 ////////////////////////////////////////////// << 556 }
>> 557 else
>> 558 {
>> 559 return snxt=kInfinity;
>> 560 }
>> 561 }
>> 562 else if (vy<0)
>> 563 {
>> 564 max=-fDy-yt;
>> 565 if (max<-kCarTolerance/2)
>> 566 {
>> 567 tmax=max/vy;
>> 568 tmin=(fDy-yt)/vy;
>> 569 }
>> 570 else
>> 571 {
>> 572 return snxt=kInfinity;
>> 573 }
>> 574 }
>> 575 else
>> 576 {
>> 577 if (fabs(yt)<=fDy)
>> 578 {
>> 579 tmin=0;
>> 580 tmax=kInfinity;
>> 581 }
>> 582 else
>> 583 {
>> 584 return snxt=kInfinity;
>> 585 }
>> 586 }
>> 587
>> 588 // Re-Calc valid intersection range
>> 589 if (tmin>smin) smin=tmin;
>> 590 if (tmax<smax) smax=tmax;
>> 591 if (smax<=smin)
>> 592 {
>> 593 return snxt=kInfinity;
>> 594 }
>> 595 else
>> 596 {
533 // 597 //
534 // Calculate distance to shape from outside << 598 // X G4Parallel planes intersection
535 // - return kInfinity if no intersection << 599 //
536 << 600 xt=p.x()-fTthetaCphi*p.z()-fTalpha*yt;
537 G4double G4Para::DistanceToIn(const G4ThreeVec << 601 vx=v.x()-fTthetaCphi*v.z()-fTalpha*vy;
538 const G4ThreeVec << 602 if (vx>0)
539 { << 603 {
540 // Z intersections << 604 max=fDx-xt;
541 // << 605 if (max>kCarTolerance/2)
542 if ((std::abs(p.z()) - fDz) >= -halfCarToler << 606 {
543 return kInfinity; << 607 tmax=max/vx;
544 G4double invz = (-v.z() == 0) ? DBL_MAX : -1 << 608 tmin=(-fDx-xt)/vx;
545 G4double dz = (invz < 0) ? fDz : -fDz; << 609 }
546 G4double tzmin = (p.z() + dz)*invz; << 610 else
547 G4double tzmax = (p.z() - dz)*invz; << 611 {
548 << 612 return snxt=kInfinity;
549 // Y intersections << 613 }
550 // << 614 }
551 G4double tmin0 = tzmin, tmax0 = tzmax; << 615 else if (vx<0)
552 G4double cos0 = fPlanes[0].b*v.y() + fPlanes << 616 {
553 G4double disy = fPlanes[0].b*p.y() + fPlanes << 617 max=-fDx-xt;
554 G4double dis0 = fPlanes[0].d + disy; << 618 if (max<-kCarTolerance/2)
555 if (dis0 >= -halfCarTolerance) << 619 {
556 { << 620 tmax=max/vx;
557 if (cos0 >= 0) return kInfinity; << 621 tmin=(fDx-xt)/vx;
558 G4double tmp = -dis0/cos0; << 622 }
559 if (tmin0 < tmp) tmin0 = tmp; << 623 else
560 } << 624 {
561 else if (cos0 > 0) << 625 return snxt=kInfinity;
562 { << 626 }
563 G4double tmp = -dis0/cos0; << 627 }
564 if (tmax0 > tmp) tmax0 = tmp; << 628 else
565 } << 629 {
566 << 630 if (fabs(xt)<=fDx)
567 G4double tmin1 = tmin0, tmax1 = tmax0; << 631 {
568 G4double cos1 = -cos0; << 632 tmin=0;
569 G4double dis1 = fPlanes[1].d - disy; << 633 tmax=kInfinity;
570 if (dis1 >= -halfCarTolerance) << 634 }
571 { << 635 else
572 if (cos1 >= 0) return kInfinity; << 636 {
573 G4double tmp = -dis1/cos1; << 637 return snxt=kInfinity;
574 if (tmin1 < tmp) tmin1 = tmp; << 638 }
575 } << 639 }
576 else if (cos1 > 0) << 640 if (tmin>smin) smin=tmin;
577 { << 641 if (tmax<smax) smax=tmax;
578 G4double tmp = -dis1/cos1; << 642 }
579 if (tmax1 > tmp) tmax1 = tmp; << 643
580 } << 644 if (smax>0&&smin<smax)
581 << 645 {
582 // X intersections << 646 if (smin>0)
583 // << 647 {
584 G4double tmin2 = tmin1, tmax2 = tmax1; << 648 snxt=smin;
585 G4double cos2 = fPlanes[2].a*v.x() + fPlanes << 649 }
586 G4double disx = fPlanes[2].a*p.x() + fPlanes << 650 else
587 G4double dis2 = fPlanes[2].d + disx; << 651 {
588 if (dis2 >= -halfCarTolerance) << 652 snxt=0;
589 { << 653 }
590 if (cos2 >= 0) return kInfinity; << 654 }
591 G4double tmp = -dis2/cos2; << 655 else
592 if (tmin2 < tmp) tmin2 = tmp; << 656 {
593 } << 657 snxt=kInfinity;
594 else if (cos2 > 0) << 658 }
595 { << 659 return snxt;
596 G4double tmp = -dis2/cos2; <<
597 if (tmax2 > tmp) tmax2 = tmp; <<
598 } <<
599 <<
600 G4double tmin3 = tmin2, tmax3 = tmax2; <<
601 G4double cos3 = -cos2; <<
602 G4double dis3 = fPlanes[3].d - disx; <<
603 if (dis3 >= -halfCarTolerance) <<
604 { <<
605 if (cos3 >= 0) return kInfinity; <<
606 G4double tmp = -dis3/cos3; <<
607 if (tmin3 < tmp) tmin3 = tmp; <<
608 } <<
609 else if (cos3 > 0) <<
610 { <<
611 G4double tmp = -dis3/cos3; <<
612 if (tmax3 > tmp) tmax3 = tmp; <<
613 } <<
614 <<
615 // Find distance <<
616 // <<
617 G4double tmin = tmin3, tmax = tmax3; <<
618 if (tmax <= tmin + halfCarTolerance) return <<
619 return (tmin < halfCarTolerance ) ? 0. : tmi <<
620 } 660 }
621 661
622 ////////////////////////////////////////////// << 662 ////////////////////////////////////////////////////////////////////////////
623 // 663 //
624 // Calculate exact shortest distance to any bo 664 // Calculate exact shortest distance to any boundary from outside
625 // - returns 0 is point inside << 665 // - Returns 0 is point inside
626 666
627 G4double G4Para::DistanceToIn( const G4ThreeVe << 667 G4double G4Para::DistanceToIn(const G4ThreeVector& p) const
628 { 668 {
629 G4double xx = fPlanes[2].a*p.x()+fPlanes[2]. << 669 G4double safe;
630 G4double dx = std::abs(xx) + fPlanes[2].d; << 670 G4double distz1,distz2,disty1,disty2,distx1,distx2;
631 << 671 G4double trany,cosy,tranx,cosx;
632 G4double yy = fPlanes[0].b*p.y()+fPlanes[0]. << 672
633 G4double dy = std::abs(yy) + fPlanes[0].d; << 673 // Z planes
634 G4double dxy = std::max(dx,dy); << 674 distz1=p.z()-fDz;
635 << 675 distz2=-fDz-p.z();
636 G4double dz = std::abs(p.z())-fDz; << 676 if (distz1>distz2)
637 G4double dist = std::max(dxy,dz); << 677 {
>> 678 safe=distz1;
>> 679 }
>> 680 else
>> 681 {
>> 682 safe=distz2;
>> 683 }
>> 684
>> 685 trany=p.y()-fTthetaSphi*p.z(); // Transformed y into `box' system
>> 686 // Transformed x into `box' system
>> 687 cosy=1.0/sqrt(1.0+fTthetaSphi*fTthetaSphi);
>> 688 disty1=(trany-fDy)*cosy;
>> 689 disty2=(-fDy-trany)*cosy;
>> 690
>> 691 if (disty1>safe) safe=disty1;
>> 692 if (disty2>safe) safe=disty2;
638 693
639 return (dist > 0) ? dist : 0.; << 694 tranx=p.x()-fTthetaCphi*p.z()-fTalpha*trany;
>> 695 cosx=1.0/sqrt(1.0+fTalpha*fTalpha+fTthetaCphi*fTthetaCphi);
>> 696 distx1=(tranx-fDx)*cosx;
>> 697 distx2=(-fDx-tranx)*cosx;
>> 698
>> 699 if (distx1>safe) safe=distx1;
>> 700 if (distx2>safe) safe=distx2;
>> 701
>> 702 if (safe<0) safe=0;
>> 703 return safe;
640 } 704 }
641 705
642 ////////////////////////////////////////////// 706 //////////////////////////////////////////////////////////////////////////
643 // 707 //
644 // Calculate distance to surface of shape from << 708 // Calcluate distance to surface of shape from inside
645 // find normal at exit point << 709 // Calculate distance to x/y/z planes - smallest is exiting distance
646 // - when leaving the surface, return 0 <<
647 <<
648 G4double G4Para::DistanceToOut(const G4ThreeVe <<
649 const G4bool ca <<
650 G4bool* v <<
651 { <<
652 // Z intersections <<
653 // <<
654 if ((std::abs(p.z()) - fDz) >= -halfCarToler <<
655 { <<
656 if (calcNorm) <<
657 { <<
658 *validNorm = true; <<
659 n->set(0, 0, (p.z() < 0) ? -1 : 1); <<
660 } <<
661 return 0.; <<
662 } <<
663 G4double vz = v.z(); <<
664 G4double tmax = (vz == 0) ? DBL_MAX : (std:: <<
665 G4int iside = (vz < 0) ? -4 : -2; // little <<
666 <<
667 // Y intersections <<
668 // <<
669 G4double cos0 = fPlanes[0].b*v.y() + fPlanes <<
670 if (cos0 > 0) <<
671 { <<
672 G4double dis0 = fPlanes[0].b*p.y() + fPlan <<
673 if (dis0 >= -halfCarTolerance) <<
674 { <<
675 if (calcNorm) <<
676 { <<
677 *validNorm = true; <<
678 n->set(0, fPlanes[0].b, fPlanes[0].c); <<
679 } <<
680 return 0.; <<
681 } <<
682 G4double tmp = -dis0/cos0; <<
683 if (tmax > tmp) { tmax = tmp; iside = 0; } <<
684 } <<
685 710
686 G4double cos1 = -cos0; << 711 G4double G4Para::DistanceToOut(const G4ThreeVector& p,const G4ThreeVector& v,
687 if (cos1 > 0) << 712 const G4bool calcNorm,
688 { << 713 G4bool *validNorm,G4ThreeVector *n) const
689 G4double dis1 = fPlanes[1].b*p.y() + fPlan << 714 {
690 if (dis1 >= -halfCarTolerance) << 715 ESide side = kUndefined ;
691 { << 716 G4double snxt; // snxt = return value
692 if (calcNorm) << 717 G4double max,tmax;
693 { << 718 G4double yt,vy,xt,vx;
694 *validNorm = true; << 719
695 n->set(0, fPlanes[1].b, fPlanes[1].c); << 720 G4double ycomp,calpha,salpha,tntheta,cosntheta;
696 } << 721 //
697 return 0.; << 722 // Z Intersections
698 } << 723 //
699 G4double tmp = -dis1/cos1; << 724 if (v.z()>0)
700 if (tmax > tmp) { tmax = tmp; iside = 1; } << 725 {
701 } << 726 max=fDz-p.z();
>> 727 if (max>kCarTolerance/2)
>> 728 {
>> 729 snxt=max/v.z();
>> 730 side=kPZ;
>> 731 }
>> 732 else
>> 733 {
>> 734 if (calcNorm)
>> 735 {
>> 736 *validNorm=true;
>> 737 *n=G4ThreeVector(0,0,1);
>> 738 }
>> 739 return snxt=0;
>> 740 }
>> 741 }
>> 742 else if (v.z()<0)
>> 743 {
>> 744 max=-fDz-p.z();
>> 745 if (max<-kCarTolerance/2)
>> 746 {
>> 747 snxt=max/v.z();
>> 748 side=kMZ;
>> 749 }
>> 750 else
>> 751 {
>> 752 if (calcNorm)
>> 753 {
>> 754 *validNorm=true;
>> 755 *n=G4ThreeVector(0,0,-1);
>> 756 }
>> 757 return snxt=0;
>> 758 }
>> 759 }
>> 760 else
>> 761 {
>> 762 snxt=kInfinity;
>> 763 }
702 764
703 // X intersections << 765
704 // << 766 //
705 G4double cos2 = fPlanes[2].a*v.x() + fPlanes << 767 // Y plane intersection
706 if (cos2 > 0) << 768 //
707 { << 769 yt=p.y()-fTthetaSphi*p.z();
708 G4double dis2 = fPlanes[2].a*p.x()+fPlanes << 770 vy=v.y()-fTthetaSphi*v.z();
709 if (dis2 >= -halfCarTolerance) <<
710 { <<
711 if (calcNorm) <<
712 { <<
713 *validNorm = true; <<
714 n->set(fPlanes[2].a, fPlanes[2].b, fP <<
715 } <<
716 return 0.; <<
717 } <<
718 G4double tmp = -dis2/cos2; <<
719 if (tmax > tmp) { tmax = tmp; iside = 2; } <<
720 } <<
721 771
722 G4double cos3 = -cos2; << 772 if (vy>0)
723 if (cos3 > 0) << 773 {
724 { << 774 max=fDy-yt;
725 G4double dis3 = fPlanes[3].a*p.x()+fPlanes << 775 if (max>kCarTolerance/2)
726 if (dis3 >= -halfCarTolerance) << 776 {
727 { << 777 tmax=max/vy;
728 if (calcNorm) << 778 if (tmax<snxt)
729 { << 779 {
730 *validNorm = true; << 780 snxt=tmax;
731 n->set(fPlanes[3].a, fPlanes[3].b, fP << 781 side=kPY;
732 } << 782 }
733 return 0.; << 783 }
734 } << 784 else
735 G4double tmp = -dis3/cos3; << 785 {
736 if (tmax > tmp) { tmax = tmp; iside = 3; } << 786 if (calcNorm)
737 } << 787 {
>> 788 *validNorm=true; // Leaving via plus Y
>> 789 ycomp=1/sqrt(1+fTthetaSphi*fTthetaSphi);
>> 790 *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
>> 791 }
>> 792 return snxt=0;
>> 793 }
>> 794 }
>> 795 else if (vy<0)
>> 796 {
>> 797 max=-fDy-yt;
>> 798 if (max<-kCarTolerance/2)
>> 799 {
>> 800 tmax=max/vy;
>> 801 if (tmax<snxt)
>> 802 {
>> 803 snxt=tmax;
>> 804 side=kMY;
>> 805 }
>> 806 }
>> 807 else
>> 808 {
>> 809 if (calcNorm)
>> 810 {
>> 811 *validNorm=true; // Leaving via minus Y
>> 812 ycomp=-1/sqrt(1+fTthetaSphi*fTthetaSphi);
>> 813 *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
>> 814
>> 815 }
>> 816 return snxt=0;
>> 817 }
>> 818 }
>> 819 //
>> 820 // X plane intersection
>> 821 //
>> 822 xt=p.x()-fTthetaCphi*p.z()-fTalpha*yt;
>> 823 vx=v.x()-fTthetaCphi*v.z()-fTalpha*vy;
>> 824 if (vx>0)
>> 825 {
>> 826 max=fDx-xt;
>> 827 if (max>kCarTolerance/2)
>> 828 {
>> 829 tmax=max/vx;
>> 830 if (tmax<snxt)
>> 831 {
>> 832 snxt=tmax;
>> 833 side=kPX;
>> 834 }
>> 835 }
>> 836 else
>> 837 {
>> 838 if (calcNorm)
>> 839 {
>> 840 *validNorm=true; // Leaving via plus X
>> 841 calpha=1/sqrt(1+fTalpha*fTalpha);
>> 842 if (fTalpha)
>> 843 {
>> 844 salpha=-calpha/fTalpha; // NOTE: actually use MINUS sin(alpha)
>> 845 }
>> 846 else
>> 847 {
>> 848 salpha=0;
>> 849 }
>> 850 tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
>> 851 cosntheta=1/sqrt(1+tntheta*tntheta);
>> 852 *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
>> 853 }
>> 854 return snxt=0;
>> 855 }
>> 856 }
>> 857 else if (vx<0)
>> 858 {
>> 859 max=-fDx-xt;
>> 860 if (max<-kCarTolerance/2)
>> 861 {
>> 862 tmax=max/vx;
>> 863 if (tmax<snxt)
>> 864 {
>> 865 snxt=tmax;
>> 866 side=kMX;
>> 867 }
>> 868 }
>> 869 else
>> 870 {
>> 871 if (calcNorm)
>> 872 {
>> 873 *validNorm=true; // Leaving via minus X
>> 874 calpha=1/sqrt(1+fTalpha*fTalpha);
>> 875 if (fTalpha)
>> 876 {
>> 877 salpha=-calpha/fTalpha; // NOTE: actually use MINUS sin(alpha)
>> 878 }
>> 879 else
>> 880 {
>> 881 salpha=0;
>> 882 }
>> 883 tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
>> 884 cosntheta=-1/sqrt(1+tntheta*tntheta);
>> 885 *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
>> 886 return snxt=0;
>> 887 }
>> 888 }
>> 889 }
738 890
739 // Set normal, if required, and return dista << 891 if (calcNorm)
740 // << 892 {
741 if (calcNorm) << 893 *validNorm=true;
742 { << 894 switch (side)
743 *validNorm = true; << 895 {
744 if (iside < 0) << 896 case kMZ:
745 n->set(0, 0, iside + 3); // (-4+3)=-1, ( << 897 *n=G4ThreeVector(0,0,-1);
746 else << 898 break;
747 n->set(fPlanes[iside].a, fPlanes[iside]. << 899 case kPZ:
748 } << 900 *n=G4ThreeVector(0,0,1);
749 return tmax; << 901 break;
>> 902 case kMY:
>> 903 ycomp=-1/sqrt(1+fTthetaSphi*fTthetaSphi);
>> 904 *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
>> 905 break;
>> 906 case kPY:
>> 907 ycomp=1/sqrt(1+fTthetaSphi*fTthetaSphi);
>> 908 *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
>> 909 break;
>> 910 case kMX:
>> 911 calpha=1/sqrt(1+fTalpha*fTalpha);
>> 912 if (fTalpha)
>> 913 {
>> 914 salpha=-calpha/fTalpha; // NOTE: actually use MINUS sin(alpha)
>> 915 }
>> 916 else
>> 917 {
>> 918 salpha=0;
>> 919 }
>> 920 tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
>> 921 cosntheta=-1/sqrt(1+tntheta*tntheta);
>> 922 *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
>> 923 break;
>> 924 case kPX:
>> 925 calpha=1/sqrt(1+fTalpha*fTalpha);
>> 926 if (fTalpha)
>> 927 {
>> 928 salpha=-calpha/fTalpha; // NOTE: actually use MINUS sin(alpha)
>> 929 }
>> 930 else
>> 931 {
>> 932 salpha=0;
>> 933 }
>> 934 tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
>> 935 cosntheta=1/sqrt(1+tntheta*tntheta);
>> 936 *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
>> 937 break;
>> 938 default:
>> 939 G4Exception("Invalid enum in G4Para::DistanceToOut");
>> 940 break;
>> 941
>> 942 }
>> 943 }
>> 944 return snxt;
750 } 945 }
751 946
752 ////////////////////////////////////////////// << 947 /////////////////////////////////////////////////////////////////////////////
753 // 948 //
754 // Calculate exact shortest distance to any bo 949 // Calculate exact shortest distance to any boundary from inside
755 // - returns 0 is point outside << 950 // - Returns 0 is point outside
756 951
757 G4double G4Para::DistanceToOut( const G4ThreeV << 952 G4double G4Para::DistanceToOut(const G4ThreeVector& p) const
758 { 953 {
759 #ifdef G4CSGDEBUG << 954 G4double safe;
760 if( Inside(p) == kOutside ) << 955 G4double distz1,distz2,disty1,disty2,distx1,distx2;
761 { << 956 G4double trany,cosy,tranx,cosx;
762 std::ostringstream message; << 957
763 G4int oldprc = message.precision(16); << 958 // Z planes
764 message << "Point p is outside (!?) of sol << 959 distz1=fDz-p.z();
765 message << "Position:\n"; << 960 distz2=fDz+p.z();
766 message << " p.x() = " << p.x()/mm << " << 961 if (distz1<distz2)
767 message << " p.y() = " << p.y()/mm << " << 962 {
768 message << " p.z() = " << p.z()/mm << " << 963 safe=distz1;
769 G4cout.precision(oldprc) ; << 964 }
770 G4Exception("G4Para::DistanceToOut(p)", "G << 965 else
771 JustWarning, message ); << 966 {
772 DumpInfo(); << 967 safe=distz2;
773 } << 968 }
774 #endif << 969
775 G4double xx = fPlanes[2].a*p.x()+fPlanes[2]. << 970 trany=p.y()-fTthetaSphi*p.z(); // Transformed y into `box' system
776 G4double dx = std::abs(xx) + fPlanes[2].d; << 971 // Transformed x into `box' system
777 << 972 cosy=1.0/sqrt(1.0+fTthetaSphi*fTthetaSphi);
778 G4double yy = fPlanes[0].b*p.y()+fPlanes[0]. << 973 disty1=(fDy-trany)*cosy;
779 G4double dy = std::abs(yy) + fPlanes[0].d; << 974 disty2=(fDy+trany)*cosy;
780 G4double dxy = std::max(dx,dy); << 975
781 << 976 if (disty1<safe) safe=disty1;
782 G4double dz = std::abs(p.z())-fDz; << 977 if (disty2<safe) safe=disty2;
783 G4double dist = std::max(dxy,dz); <<
784 978
785 return (dist < 0) ? -dist : 0.; << 979 tranx=p.x()-fTthetaCphi*p.z()-fTalpha*trany;
>> 980 cosx=1.0/sqrt(1.0+fTalpha*fTalpha+fTthetaCphi*fTthetaCphi);
>> 981 distx1=(fDx-tranx)*cosx;
>> 982 distx2=(fDx+tranx)*cosx;
>> 983
>> 984 if (distx1<safe) safe=distx1;
>> 985 if (distx2<safe) safe=distx2;
>> 986
>> 987 if (safe<0) safe=0;
>> 988 return safe;
786 } 989 }
787 990
788 ////////////////////////////////////////////// << 991 ////////////////////////////////////////////////////////////////////////////////
789 // 992 //
790 // GetEntityType << 993 // Create a List containing the transformed vertices
791 << 994 // Ordering [0-3] -fDz cross section
792 G4GeometryType G4Para::GetEntityType() const << 995 // [4-7] +fDz cross section such that [0] is below [4],
793 { << 996 // [1] below [5] etc.
794 return {"G4Para"}; << 997 // Note:
>> 998 // Caller has deletion resposibility
>> 999
>> 1000 G4ThreeVectorList*
>> 1001 G4Para::CreateRotatedVertices(const G4AffineTransform& pTransform) const
>> 1002 {
>> 1003 G4ThreeVectorList *vertices;
>> 1004 vertices=new G4ThreeVectorList(8);
>> 1005 if (vertices)
>> 1006 {
>> 1007 G4ThreeVector vertex0(-fDz*fTthetaCphi-fDy*fTalpha-fDx,-fDz*fTthetaSphi-fDy,-fDz);
>> 1008 G4ThreeVector vertex1(-fDz*fTthetaCphi-fDy*fTalpha+fDx,-fDz*fTthetaSphi-fDy,-fDz);
>> 1009 G4ThreeVector vertex2(-fDz*fTthetaCphi+fDy*fTalpha-fDx,-fDz*fTthetaSphi+fDy,-fDz);
>> 1010 G4ThreeVector vertex3(-fDz*fTthetaCphi+fDy*fTalpha+fDx,-fDz*fTthetaSphi+fDy,-fDz);
>> 1011 G4ThreeVector vertex4(+fDz*fTthetaCphi-fDy*fTalpha-fDx,+fDz*fTthetaSphi-fDy,+fDz);
>> 1012 G4ThreeVector vertex5(+fDz*fTthetaCphi-fDy*fTalpha+fDx,+fDz*fTthetaSphi-fDy,+fDz);
>> 1013 G4ThreeVector vertex6(+fDz*fTthetaCphi+fDy*fTalpha-fDx,+fDz*fTthetaSphi+fDy,+fDz);
>> 1014 G4ThreeVector vertex7(+fDz*fTthetaCphi+fDy*fTalpha+fDx,+fDz*fTthetaSphi+fDy,+fDz);
>> 1015
>> 1016 vertices->insert(pTransform.TransformPoint(vertex0));
>> 1017 vertices->insert(pTransform.TransformPoint(vertex1));
>> 1018 vertices->insert(pTransform.TransformPoint(vertex2));
>> 1019 vertices->insert(pTransform.TransformPoint(vertex3));
>> 1020 vertices->insert(pTransform.TransformPoint(vertex4));
>> 1021 vertices->insert(pTransform.TransformPoint(vertex5));
>> 1022 vertices->insert(pTransform.TransformPoint(vertex6));
>> 1023 vertices->insert(pTransform.TransformPoint(vertex7));
>> 1024 }
>> 1025 else
>> 1026 {
>> 1027 G4Exception("G4Para::CreateRotatedVertices Out of memory - Cannot alloc vertices");
>> 1028 }
>> 1029 return vertices;
795 } 1030 }
796 1031
797 ////////////////////////////////////////////// << 1032 ////////////////////////////////////////////////////////////////////////////
798 // 1033 //
799 // IsFaceted << 1034 // Methods for visualisation
800 1035
801 G4bool G4Para::IsFaceted() const << 1036 void G4Para::DescribeYourselfTo (G4VGraphicsScene& scene) const
802 { 1037 {
803 return true; << 1038 scene.AddThis (*this);
>> 1039 // scene.AddThis (*this);
804 } 1040 }
805 1041
806 ////////////////////////////////////////////// << 1042 G4VisExtent G4Para::GetExtent() const
807 // <<
808 // Make a clone of the object <<
809 // <<
810 G4VSolid* G4Para::Clone() const <<
811 { 1043 {
812 return new G4Para(*this); << 1044 G4int i ;
>> 1045 G4double xMin,xMax,yMin,yMax ;
>> 1046 G4ThreeVector pt[8] ;
>> 1047 pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy*fTalpha-fDx,-fDz*fTthetaSphi-fDy,-fDz);
>> 1048 pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy*fTalpha+fDx,-fDz*fTthetaSphi-fDy,-fDz);
>> 1049 pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy*fTalpha-fDx,-fDz*fTthetaSphi+fDy,-fDz);
>> 1050 pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy*fTalpha+fDx,-fDz*fTthetaSphi+fDy,-fDz);
>> 1051 pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy*fTalpha-fDx,+fDz*fTthetaSphi-fDy,+fDz);
>> 1052 pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy*fTalpha+fDx,+fDz*fTthetaSphi-fDy,+fDz);
>> 1053 pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy*fTalpha-fDx,+fDz*fTthetaSphi+fDy,+fDz);
>> 1054 pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy*fTalpha+fDx,+fDz*fTthetaSphi+fDy,+fDz);
>> 1055 xMax = -fabs(fDz*fTthetaCphi)-fDx-fDx-fDx-fDx ;
>> 1056 xMin = -xMax ;
>> 1057 yMax = -fabs(fDz*fTthetaSphi)-fDy-fDy ;
>> 1058 yMin = -yMax ;
>> 1059 for(i=0;i<8;i++)
>> 1060 {
>> 1061 if(pt[i].x() > xMax) xMax = pt[i].x() ;
>> 1062 if(pt[i].x() < xMin) xMin = pt[i].x() ;
>> 1063 if(pt[i].y() > yMax) yMax = pt[i].y() ;
>> 1064 if(pt[i].y() < yMin) yMin = pt[i].y() ;
>> 1065 }
>> 1066 return G4VisExtent (xMin, xMax, yMin, yMax, -fDz, fDz);
813 } 1067 }
814 1068
815 ////////////////////////////////////////////// << 1069 G4Polyhedron* G4Para::CreatePolyhedron () const
816 // <<
817 // Stream object contents to an output stream <<
818 <<
819 std::ostream& G4Para::StreamInfo( std::ostream <<
820 { 1070 {
821 G4double alpha = std::atan(fTalpha); << 1071 G4double phi = atan2(fTthetaSphi, fTthetaCphi);
822 G4double theta = std::atan(std::sqrt(fTtheta << 1072 G4double alpha = atan(fTalpha);
823 fTtheta << 1073 G4double theta = atan(sqrt(fTthetaCphi*fTthetaCphi
824 G4double phi = std::atan2(fTthetaSphi,fTth << 1074 +fTthetaSphi*fTthetaSphi));
825 << 1075
826 G4long oldprc = os.precision(16); << 1076 return new G4PolyhedronPara(fDx, fDy, fDz, alpha, theta, phi);
827 os << "------------------------------------- <<
828 << " *** Dump for solid - " << GetName <<
829 << " ================================= <<
830 << " Solid type: G4Para\n" <<
831 << " Parameters:\n" <<
832 << " half length X: " << fDx/mm << " m <<
833 << " half length Y: " << fDy/mm << " m <<
834 << " half length Z: " << fDz/mm << " m <<
835 << " alpha: " << alpha/degree << "degr <<
836 << " theta: " << theta/degree << "degr <<
837 << " phi: " << phi/degree << "degrees\ <<
838 << "------------------------------------- <<
839 os.precision(oldprc); <<
840 <<
841 return os; <<
842 } 1077 }
843 1078
844 ////////////////////////////////////////////// << 1079 G4NURBS* G4Para::CreateNURBS () const
845 // <<
846 // Return a point randomly and uniformly selec <<
847 <<
848 G4ThreeVector G4Para::GetPointOnSurface() cons <<
849 { 1080 {
850 G4double DyTalpha = fDy*fTalpha; << 1081 // return new G4NURBSbox (fDx, fDy, fDz);
851 G4double DzTthetaSphi = fDz*fTthetaSphi; << 1082 return 0 ;
852 G4double DzTthetaCphi = fDz*fTthetaCphi; <<
853 <<
854 // Set vertices <<
855 // <<
856 G4ThreeVector pt[8]; <<
857 pt[0].set(-DzTthetaCphi-DyTalpha-fDx, -DzTth <<
858 pt[1].set(-DzTthetaCphi-DyTalpha+fDx, -DzTth <<
859 pt[2].set(-DzTthetaCphi+DyTalpha-fDx, -DzTth <<
860 pt[3].set(-DzTthetaCphi+DyTalpha+fDx, -DzTth <<
861 pt[4].set( DzTthetaCphi-DyTalpha-fDx, DzTth <<
862 pt[5].set( DzTthetaCphi-DyTalpha+fDx, DzTth <<
863 pt[6].set( DzTthetaCphi+DyTalpha-fDx, DzTth <<
864 pt[7].set( DzTthetaCphi+DyTalpha+fDx, DzTth <<
865 <<
866 // Set areas (-Z, -Y, +Y, -X, +X, +Z) <<
867 // <<
868 G4ThreeVector vx(fDx, 0, 0); <<
869 G4ThreeVector vy(DyTalpha, fDy, 0); <<
870 G4ThreeVector vz(DzTthetaCphi, DzTthetaSphi, <<
871 <<
872 G4double sxy = fDx*fDy; // (vx.cross(vy)).ma <<
873 G4double sxz = (vx.cross(vz)).mag(); <<
874 G4double syz = (vy.cross(vz)).mag(); <<
875 <<
876 G4double sface[6] = { sxy, syz, syz, sxz, sx <<
877 for (G4int i=1; i<6; ++i) { sface[i] += sfac <<
878 <<
879 // Select face <<
880 // <<
881 G4double select = sface[5]*G4UniformRand(); <<
882 G4int k = 5; <<
883 if (select <= sface[4]) k = 4; <<
884 if (select <= sface[3]) k = 3; <<
885 if (select <= sface[2]) k = 2; <<
886 if (select <= sface[1]) k = 1; <<
887 if (select <= sface[0]) k = 0; <<
888 <<
889 // Generate point <<
890 // <<
891 G4int ip[6][3] = {{0,1,2}, {0,4,1}, {2,3,6}, <<
892 G4double u = G4UniformRand(); <<
893 G4double v = G4UniformRand(); <<
894 return (1.-u-v)*pt[ip[k][0]] + u*pt[ip[k][1] <<
895 } 1083 }
896 1084
897 ////////////////////////////////////////////// <<
898 // 1085 //
899 // Methods for visualisation << 1086 //
>> 1087 /////////////////////////// End of G4Para.cc ///////////////////////////
900 1088
901 void G4Para::DescribeYourselfTo ( G4VGraphicsS <<
902 { <<
903 scene.AddSolid (*this); <<
904 } <<
905 1089
906 G4Polyhedron* G4Para::CreatePolyhedron () cons <<
907 { <<
908 G4double phi = std::atan2(fTthetaSphi, fTthe <<
909 G4double alpha = std::atan(fTalpha); <<
910 G4double theta = std::atan(std::sqrt(fTtheta <<
911 fTtheta <<
912 <<
913 return new G4PolyhedronPara(fDx, fDy, fDz, a <<
914 } <<
915 #endif <<
916 1090