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