| 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 // $Id: G4Paraboloid.cc 104316 2017-05-24 13:04:23Z gcosmo $
>> 27 //
>> 28 // class G4Paraboloid
>> 29 //
26 // Implementation for G4Paraboloid class 30 // Implementation for G4Paraboloid class
27 // 31 //
28 // Author : Lukas Lindroos (CERN), July 2007 32 // Author : Lukas Lindroos (CERN), July 2007
29 // Revised: Tatiana Nikitina (CERN) 33 // Revised: Tatiana Nikitina (CERN)
30 // ------------------------------------------- 34 // --------------------------------------------------------------------
31 35
32 #include "globals.hh" 36 #include "globals.hh"
33 37
34 #include "G4Paraboloid.hh" 38 #include "G4Paraboloid.hh"
35 39
36 #if !(defined(G4GEOM_USE_UPARABOLOID) && defin 40 #if !(defined(G4GEOM_USE_UPARABOLOID) && defined(G4GEOM_USE_SYS_USOLIDS))
37 41
38 #include "G4VoxelLimits.hh" 42 #include "G4VoxelLimits.hh"
39 #include "G4AffineTransform.hh" 43 #include "G4AffineTransform.hh"
40 #include "G4BoundingEnvelope.hh" 44 #include "G4BoundingEnvelope.hh"
41 45
42 #include "meshdefs.hh" 46 #include "meshdefs.hh"
43 47
44 #include "Randomize.hh" 48 #include "Randomize.hh"
45 49
46 #include "G4VGraphicsScene.hh" 50 #include "G4VGraphicsScene.hh"
47 #include "G4VisExtent.hh" 51 #include "G4VisExtent.hh"
48 52
49 #include "G4AutoLock.hh" 53 #include "G4AutoLock.hh"
50 54
51 namespace 55 namespace
52 { 56 {
53 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE 57 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
54 } 58 }
55 59
56 using namespace CLHEP; 60 using namespace CLHEP;
57 61
58 ////////////////////////////////////////////// 62 ///////////////////////////////////////////////////////////////////////////////
59 // 63 //
60 // constructor - check parameters 64 // constructor - check parameters
61 // << 65
62 G4Paraboloid::G4Paraboloid(const G4String& pNa 66 G4Paraboloid::G4Paraboloid(const G4String& pName,
63 G4double pDz, 67 G4double pDz,
64 G4double pR1, 68 G4double pR1,
65 G4double pR2) 69 G4double pR2)
66 : G4VSolid(pName) << 70 : G4VSolid(pName), fRebuildPolyhedron(false), fpPolyhedron(0),
>> 71 fSurfaceArea(0.), fCubicVolume(0.)
>> 72
67 { 73 {
68 if( (pDz <= 0.) || (pR2 <= pR1) || (pR1 < 0. 74 if( (pDz <= 0.) || (pR2 <= pR1) || (pR1 < 0.) )
69 { 75 {
70 std::ostringstream message; 76 std::ostringstream message;
71 message << "Invalid dimensions. Negative I 77 message << "Invalid dimensions. Negative Input Values or R1>=R2 - "
72 << GetName(); 78 << GetName();
73 G4Exception("G4Paraboloid::G4Paraboloid()" 79 G4Exception("G4Paraboloid::G4Paraboloid()", "GeomSolids0002",
74 FatalErrorInArgument, message, 80 FatalErrorInArgument, message,
75 "Z half-length must be larger 81 "Z half-length must be larger than zero or R1>=R2.");
76 } 82 }
77 83
78 r1 = pR1; 84 r1 = pR1;
79 r2 = pR2; 85 r2 = pR2;
80 dz = pDz; 86 dz = pDz;
81 87
82 // r1^2 = k1 * (-dz) + k2 88 // r1^2 = k1 * (-dz) + k2
83 // r2^2 = k1 * ( dz) + k2 89 // r2^2 = k1 * ( dz) + k2
84 // => r1^2 + r2^2 = k2 + k2 => k2 = (r2^2 + 90 // => r1^2 + r2^2 = k2 + k2 => k2 = (r2^2 + r1^2) / 2
85 // and r2^2 - r1^2 = k1 * dz - k1 * (-dz) => 91 // and r2^2 - r1^2 = k1 * dz - k1 * (-dz) => k1 = (r2^2 - r1^2) / 2 / dz
86 92
87 k1 = (r2 * r2 - r1 * r1) / 2 / dz; 93 k1 = (r2 * r2 - r1 * r1) / 2 / dz;
88 k2 = (r2 * r2 + r1 * r1) / 2; 94 k2 = (r2 * r2 + r1 * r1) / 2;
89 } 95 }
90 96
91 ////////////////////////////////////////////// 97 ///////////////////////////////////////////////////////////////////////////////
92 // 98 //
93 // Fake default constructor - sets only member 99 // Fake default constructor - sets only member data and allocates memory
94 // for usage restri 100 // for usage restricted to object persistency.
95 // 101 //
96 G4Paraboloid::G4Paraboloid( __void__& a ) 102 G4Paraboloid::G4Paraboloid( __void__& a )
97 : G4VSolid(a), dz(0.), r1(0.), r2(0.), k1(0. << 103 : G4VSolid(a), fRebuildPolyhedron(false), fpPolyhedron(0),
>> 104 fSurfaceArea(0.), fCubicVolume(0.),
>> 105 dz(0.), r1(0.), r2(0.), k1(0.), k2(0.)
98 { 106 {
99 } 107 }
100 108
101 ////////////////////////////////////////////// 109 ///////////////////////////////////////////////////////////////////////////////
102 // 110 //
103 // Destructor 111 // Destructor
104 // << 112
105 G4Paraboloid::~G4Paraboloid() 113 G4Paraboloid::~G4Paraboloid()
106 { 114 {
107 delete fpPolyhedron; fpPolyhedron = nullptr; << 115 delete fpPolyhedron; fpPolyhedron = 0;
108 } 116 }
109 117
110 ////////////////////////////////////////////// 118 ///////////////////////////////////////////////////////////////////////////////
111 // 119 //
112 // Copy constructor 120 // Copy constructor
113 // << 121
114 G4Paraboloid::G4Paraboloid(const G4Paraboloid& 122 G4Paraboloid::G4Paraboloid(const G4Paraboloid& rhs)
115 : G4VSolid(rhs), << 123 : G4VSolid(rhs), fRebuildPolyhedron(false), fpPolyhedron(0),
116 fSurfaceArea(rhs.fSurfaceArea), fCubicVolu 124 fSurfaceArea(rhs.fSurfaceArea), fCubicVolume(rhs.fCubicVolume),
117 dz(rhs.dz), r1(rhs.r1), r2(rhs.r2), k1(rhs 125 dz(rhs.dz), r1(rhs.r1), r2(rhs.r2), k1(rhs.k1), k2(rhs.k2)
118 { 126 {
119 } 127 }
120 128
>> 129
121 ////////////////////////////////////////////// 130 ///////////////////////////////////////////////////////////////////////////////
122 // 131 //
123 // Assignment operator 132 // Assignment operator
124 // << 133
125 G4Paraboloid& G4Paraboloid::operator = (const 134 G4Paraboloid& G4Paraboloid::operator = (const G4Paraboloid& rhs)
126 { 135 {
127 // Check assignment to self 136 // Check assignment to self
128 // 137 //
129 if (this == &rhs) { return *this; } 138 if (this == &rhs) { return *this; }
130 139
131 // Copy base class data 140 // Copy base class data
132 // 141 //
133 G4VSolid::operator=(rhs); 142 G4VSolid::operator=(rhs);
134 143
135 // Copy data 144 // Copy data
136 // 145 //
137 fSurfaceArea = rhs.fSurfaceArea; fCubicVolu 146 fSurfaceArea = rhs.fSurfaceArea; fCubicVolume = rhs.fCubicVolume;
138 dz = rhs.dz; r1 = rhs.r1; r2 = rhs.r2; k1 = 147 dz = rhs.dz; r1 = rhs.r1; r2 = rhs.r2; k1 = rhs.k1; k2 = rhs.k2;
139 fRebuildPolyhedron = false; 148 fRebuildPolyhedron = false;
140 delete fpPolyhedron; fpPolyhedron = nullptr << 149 delete fpPolyhedron; fpPolyhedron = 0;
141 150
142 return *this; 151 return *this;
143 } 152 }
144 153
>> 154 /////////////////////////////////////////////////////////////////////////
>> 155 //
>> 156 // Dispatch to parameterisation for replication mechanism dimension
>> 157 // computation & modification.
>> 158
>> 159 //void ComputeDimensions( G4VPVParamerisation p,
>> 160 // const G4Int n,
>> 161 // const G4VPhysicalVolume* pRep )
>> 162 //{
>> 163 // p->ComputeDimensions(*this,n,pRep) ;
>> 164 //}
>> 165
145 ////////////////////////////////////////////// 166 ///////////////////////////////////////////////////////////////////////////////
146 // 167 //
147 // Get bounding box 168 // Get bounding box
148 // << 169
149 void G4Paraboloid::BoundingLimits(G4ThreeVecto 170 void G4Paraboloid::BoundingLimits(G4ThreeVector& pMin,
150 G4ThreeVecto 171 G4ThreeVector& pMax) const
151 { 172 {
152 pMin.set(-r2,-r2,-dz); 173 pMin.set(-r2,-r2,-dz);
153 pMax.set( r2, r2, dz); 174 pMax.set( r2, r2, dz);
154 175
155 // Check correctness of the bounding box 176 // Check correctness of the bounding box
156 // 177 //
157 if (pMin.x() >= pMax.x() || pMin.y() >= pMax 178 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
158 { 179 {
159 std::ostringstream message; 180 std::ostringstream message;
160 message << "Bad bounding box (min >= max) 181 message << "Bad bounding box (min >= max) for solid: "
161 << GetName() << " !" 182 << GetName() << " !"
162 << "\npMin = " << pMin 183 << "\npMin = " << pMin
163 << "\npMax = " << pMax; 184 << "\npMax = " << pMax;
164 G4Exception("G4Paraboloid::BoundingLimits( 185 G4Exception("G4Paraboloid::BoundingLimits()", "GeomMgt0001",
165 JustWarning, message); 186 JustWarning, message);
166 DumpInfo(); 187 DumpInfo();
167 } 188 }
168 } 189 }
169 190
170 ////////////////////////////////////////////// 191 ///////////////////////////////////////////////////////////////////////////////
171 // 192 //
172 // Calculate extent under transform and specif 193 // Calculate extent under transform and specified limit
173 // << 194
174 G4bool 195 G4bool
175 G4Paraboloid::CalculateExtent(const EAxis pAxi 196 G4Paraboloid::CalculateExtent(const EAxis pAxis,
176 const G4VoxelLim << 197 const G4VoxelLimits& pVoxelLimit,
177 const G4AffineTr << 198 const G4AffineTransform& pTransform,
178 G4double& << 199 G4double& pMin, G4double& pMax) const
179 { 200 {
180 G4ThreeVector bmin, bmax; 201 G4ThreeVector bmin, bmax;
181 202
182 // Get bounding box 203 // Get bounding box
183 BoundingLimits(bmin,bmax); 204 BoundingLimits(bmin,bmax);
184 205
185 // Find extent 206 // Find extent
186 G4BoundingEnvelope bbox(bmin,bmax); 207 G4BoundingEnvelope bbox(bmin,bmax);
187 return bbox.CalculateExtent(pAxis,pVoxelLimi 208 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
188 } 209 }
189 210
190 ////////////////////////////////////////////// 211 ///////////////////////////////////////////////////////////////////////////////
191 // 212 //
192 // Return whether point inside/outside/on surf 213 // Return whether point inside/outside/on surface
193 // << 214
194 EInside G4Paraboloid::Inside(const G4ThreeVect 215 EInside G4Paraboloid::Inside(const G4ThreeVector& p) const
195 { 216 {
196 // First check is the point is above or bel 217 // First check is the point is above or below the solid.
197 // 218 //
198 if(std::fabs(p.z()) > dz + 0.5 * kCarToleran 219 if(std::fabs(p.z()) > dz + 0.5 * kCarTolerance) { return kOutside; }
199 220
200 G4double rho2 = p.perp2(), 221 G4double rho2 = p.perp2(),
201 rhoSurfTimesTol2 = (k1 * p.z() + k 222 rhoSurfTimesTol2 = (k1 * p.z() + k2) * sqr(kCarTolerance),
202 A = rho2 - ((k1 *p.z() + k2) + 0.25 223 A = rho2 - ((k1 *p.z() + k2) + 0.25 * kCarTolerance * kCarTolerance);
203 224
204 if(A < 0 && sqr(A) > rhoSurfTimesTol2) 225 if(A < 0 && sqr(A) > rhoSurfTimesTol2)
205 { 226 {
206 // Actually checking rho < radius of parab 227 // Actually checking rho < radius of paraboloid at z = p.z().
207 // We're either inside or in lower/upper c 228 // We're either inside or in lower/upper cutoff area.
208 229
209 if(std::fabs(p.z()) > dz - 0.5 * kCarToler 230 if(std::fabs(p.z()) > dz - 0.5 * kCarTolerance)
210 { 231 {
211 // We're in the upper/lower cutoff area, 232 // We're in the upper/lower cutoff area, sides have a paraboloid shape
212 // maybe further checks should be made t 233 // maybe further checks should be made to make these nicer
213 234
214 return kSurface; 235 return kSurface;
215 } 236 }
216 else 237 else
217 { 238 {
218 return kInside; 239 return kInside;
219 } 240 }
220 } 241 }
221 else if(A <= 0 || sqr(A) < rhoSurfTimesTol2) 242 else if(A <= 0 || sqr(A) < rhoSurfTimesTol2)
222 { 243 {
223 // We're in the parabolic surface. 244 // We're in the parabolic surface.
224 245
225 return kSurface; 246 return kSurface;
226 } 247 }
227 else 248 else
228 { 249 {
229 return kOutside; 250 return kOutside;
230 } 251 }
231 } 252 }
232 253
233 ////////////////////////////////////////////// 254 ///////////////////////////////////////////////////////////////////////////////
234 // 255 //
235 // SurfaceNormal << 256
236 // <<
237 G4ThreeVector G4Paraboloid::SurfaceNormal( con 257 G4ThreeVector G4Paraboloid::SurfaceNormal( const G4ThreeVector& p) const
238 { 258 {
239 G4ThreeVector n(0, 0, 0); 259 G4ThreeVector n(0, 0, 0);
240 if(std::fabs(p.z()) > dz + 0.5 * kCarToleran 260 if(std::fabs(p.z()) > dz + 0.5 * kCarTolerance)
241 { 261 {
242 // If above or below just return normal ve 262 // If above or below just return normal vector for the cutoff plane.
243 263
244 n = G4ThreeVector(0, 0, p.z()/std::fabs(p. 264 n = G4ThreeVector(0, 0, p.z()/std::fabs(p.z()));
245 } 265 }
246 else if(std::fabs(p.z()) > dz - 0.5 * kCarTo 266 else if(std::fabs(p.z()) > dz - 0.5 * kCarTolerance)
247 { 267 {
248 // This means we're somewhere in the plane 268 // This means we're somewhere in the plane z = dz or z = -dz.
249 // (As far as the program is concerned any 269 // (As far as the program is concerned anyway.
250 270
251 if(p.z() < 0) // Are we in upper or lower 271 if(p.z() < 0) // Are we in upper or lower plane?
252 { 272 {
253 if(p.perp2() > sqr(r1 + 0.5 * kCarTolera 273 if(p.perp2() > sqr(r1 + 0.5 * kCarTolerance))
254 { 274 {
255 n = G4ThreeVector(p.x(), p.y(), -k1 / 275 n = G4ThreeVector(p.x(), p.y(), -k1 / 2).unit();
256 } 276 }
257 else if(r1 < 0.5 * kCarTolerance 277 else if(r1 < 0.5 * kCarTolerance
258 || p.perp2() > sqr(r1 - 0.5 * kCarT 278 || p.perp2() > sqr(r1 - 0.5 * kCarTolerance))
259 { 279 {
260 n = G4ThreeVector(p.x(), p.y(), 0.).un 280 n = G4ThreeVector(p.x(), p.y(), 0.).unit()
261 + G4ThreeVector(0., 0., -1.).unit(); 281 + G4ThreeVector(0., 0., -1.).unit();
262 n = n.unit(); 282 n = n.unit();
263 } 283 }
264 else 284 else
265 { 285 {
266 n = G4ThreeVector(0., 0., -1.); 286 n = G4ThreeVector(0., 0., -1.);
267 } 287 }
268 } 288 }
269 else 289 else
270 { 290 {
271 if(p.perp2() > sqr(r2 + 0.5 * kCarTolera 291 if(p.perp2() > sqr(r2 + 0.5 * kCarTolerance))
272 { 292 {
273 n = G4ThreeVector(p.x(), p.y(), 0.).un 293 n = G4ThreeVector(p.x(), p.y(), 0.).unit();
274 } 294 }
275 else if(r2 < 0.5 * kCarTolerance 295 else if(r2 < 0.5 * kCarTolerance
276 || p.perp2() > sqr(r2 - 0.5 * kCarT 296 || p.perp2() > sqr(r2 - 0.5 * kCarTolerance))
277 { 297 {
278 n = G4ThreeVector(p.x(), p.y(), 0.).un 298 n = G4ThreeVector(p.x(), p.y(), 0.).unit()
279 + G4ThreeVector(0., 0., 1.).unit(); 299 + G4ThreeVector(0., 0., 1.).unit();
280 n = n.unit(); 300 n = n.unit();
281 } 301 }
282 else 302 else
283 { 303 {
284 n = G4ThreeVector(0., 0., 1.); 304 n = G4ThreeVector(0., 0., 1.);
285 } 305 }
286 } 306 }
287 } 307 }
288 else 308 else
289 { 309 {
290 G4double rho2 = p.perp2(); 310 G4double rho2 = p.perp2();
291 G4double rhoSurfTimesTol2 = (k1 * p.z() + 311 G4double rhoSurfTimesTol2 = (k1 * p.z() + k2) * sqr(kCarTolerance);
292 G4double A = rho2 - ((k1 *p.z() + k2) 312 G4double A = rho2 - ((k1 *p.z() + k2)
293 + 0.25 * kCarTolerance * kCarTo 313 + 0.25 * kCarTolerance * kCarTolerance);
294 314
295 if(A < 0 && sqr(A) > rhoSurfTimesTol2) 315 if(A < 0 && sqr(A) > rhoSurfTimesTol2)
296 { 316 {
297 // Actually checking rho < radius of par 317 // Actually checking rho < radius of paraboloid at z = p.z().
298 // We're inside. 318 // We're inside.
299 319
300 if(p.mag2() != 0) { n = p.unit(); } 320 if(p.mag2() != 0) { n = p.unit(); }
301 } 321 }
302 else if(A <= 0 || sqr(A) < rhoSurfTimesTol 322 else if(A <= 0 || sqr(A) < rhoSurfTimesTol2)
303 { 323 {
304 // We're in the parabolic surface. 324 // We're in the parabolic surface.
305 325
306 n = G4ThreeVector(p.x(), p.y(), - k1 / 2 326 n = G4ThreeVector(p.x(), p.y(), - k1 / 2).unit();
307 } 327 }
308 else 328 else
309 { 329 {
310 n = G4ThreeVector(p.x(), p.y(), - k1 / 2 330 n = G4ThreeVector(p.x(), p.y(), - k1 / 2).unit();
311 } 331 }
312 } 332 }
313 333
314 if(n.mag2() == 0) 334 if(n.mag2() == 0)
315 { 335 {
316 std::ostringstream message; 336 std::ostringstream message;
317 message << "No normal defined for this poi 337 message << "No normal defined for this point p." << G4endl
318 << " p = " << 1 / mm * p 338 << " p = " << 1 / mm * p << " mm";
319 G4Exception("G4Paraboloid::SurfaceNormal(p 339 G4Exception("G4Paraboloid::SurfaceNormal(p)", "GeomSolids1002",
320 JustWarning, message); 340 JustWarning, message);
321 } 341 }
322 return n; 342 return n;
323 } 343 }
324 344
325 ////////////////////////////////////////////// 345 ///////////////////////////////////////////////////////////////////////////////
326 // 346 //
327 // Calculate distance to shape from outside, a 347 // Calculate distance to shape from outside, along normalised vector
328 // - return kInfinity if no intersection 348 // - return kInfinity if no intersection
329 // 349 //
330 // << 350
331 G4double G4Paraboloid::DistanceToIn( const G4T 351 G4double G4Paraboloid::DistanceToIn( const G4ThreeVector& p,
332 const G4T << 352 const G4ThreeVector& v ) const
333 { 353 {
334 G4double rho2 = p.perp2(), paraRho2 = std::f 354 G4double rho2 = p.perp2(), paraRho2 = std::fabs(k1 * p.z() + k2);
335 G4double tol2 = kCarTolerance*kCarTolerance; 355 G4double tol2 = kCarTolerance*kCarTolerance;
336 G4double tolh = 0.5*kCarTolerance; 356 G4double tolh = 0.5*kCarTolerance;
337 357
338 if((r2 != 0.0) && p.z() > - tolh + dz) << 358 if(r2 && p.z() > - tolh + dz)
339 { 359 {
340 // If the points is above check for inters 360 // If the points is above check for intersection with upper edge.
341 361
342 if(v.z() < 0) 362 if(v.z() < 0)
343 { 363 {
344 G4double intersection = (dz - p.z()) / v 364 G4double intersection = (dz - p.z()) / v.z(); // With plane z = dz.
345 if(sqr(p.x() + v.x()*intersection) 365 if(sqr(p.x() + v.x()*intersection)
346 + sqr(p.y() + v.y()*intersection) < sqr 366 + sqr(p.y() + v.y()*intersection) < sqr(r2 + 0.5 * kCarTolerance))
347 { 367 {
348 if(p.z() < tolh + dz) 368 if(p.z() < tolh + dz)
349 { return 0; } 369 { return 0; }
350 else 370 else
351 { return intersection; } 371 { return intersection; }
352 } 372 }
353 } 373 }
354 else // Direction away, no possibility of 374 else // Direction away, no possibility of intersection
355 { 375 {
356 return kInfinity; 376 return kInfinity;
357 } 377 }
358 } 378 }
359 else if((r1 != 0.0) && p.z() < tolh - dz) << 379 else if(r1 && p.z() < tolh - dz)
360 { 380 {
361 // If the points is belove check for inter 381 // If the points is belove check for intersection with lower edge.
362 382
363 if(v.z() > 0) 383 if(v.z() > 0)
364 { 384 {
365 G4double intersection = (-dz - p.z()) / 385 G4double intersection = (-dz - p.z()) / v.z(); // With plane z = -dz.
366 if(sqr(p.x() + v.x()*intersection) 386 if(sqr(p.x() + v.x()*intersection)
367 + sqr(p.y() + v.y()*intersection) < sqr 387 + sqr(p.y() + v.y()*intersection) < sqr(r1 + 0.5 * kCarTolerance))
368 { 388 {
369 if(p.z() > -tolh - dz) 389 if(p.z() > -tolh - dz)
370 { 390 {
371 return 0; 391 return 0;
372 } 392 }
373 else 393 else
374 { 394 {
375 return intersection; 395 return intersection;
376 } 396 }
377 } 397 }
378 } 398 }
379 else // Direction away, no possibility of 399 else // Direction away, no possibility of intersection
380 { 400 {
381 return kInfinity; 401 return kInfinity;
382 } 402 }
383 } 403 }
384 404
385 G4double A = k1 / 2 * v.z() - p.x() * v.x() 405 G4double A = k1 / 2 * v.z() - p.x() * v.x() - p.y() * v.y(),
386 vRho2 = v.perp2(), intersection, 406 vRho2 = v.perp2(), intersection,
387 B = (k1 * p.z() + k2 - rho2) * vRho 407 B = (k1 * p.z() + k2 - rho2) * vRho2;
388 408
389 if ( ( (rho2 > paraRho2) && (sqr(rho2-paraRh 409 if ( ( (rho2 > paraRho2) && (sqr(rho2-paraRho2-0.25*tol2) > tol2*paraRho2) )
390 || (p.z() < - dz+kCarTolerance) 410 || (p.z() < - dz+kCarTolerance)
391 || (p.z() > dz-kCarTolerance) ) // Make su 411 || (p.z() > dz-kCarTolerance) ) // Make sure it's safely outside.
392 { 412 {
393 // Is there a problem with squaring rho tw 413 // Is there a problem with squaring rho twice?
394 414
395 if(vRho2<tol2) // Needs to be treated sepe 415 if(vRho2<tol2) // Needs to be treated seperately.
396 { 416 {
397 intersection = ((rho2 - k2)/k1 - p.z())/ 417 intersection = ((rho2 - k2)/k1 - p.z())/v.z();
398 if(intersection < 0) { return kInfinity; 418 if(intersection < 0) { return kInfinity; }
399 else if(std::fabs(p.z() + v.z() * inters 419 else if(std::fabs(p.z() + v.z() * intersection) <= dz)
400 { 420 {
401 return intersection; 421 return intersection;
402 } 422 }
403 else 423 else
404 { 424 {
405 return kInfinity; 425 return kInfinity;
406 } 426 }
407 } 427 }
408 else if(A*A + B < 0) // No real intersecti 428 else if(A*A + B < 0) // No real intersections.
409 { 429 {
410 return kInfinity; 430 return kInfinity;
411 } 431 }
412 else 432 else
413 { 433 {
414 intersection = (A - std::sqrt(B + sqr(A) 434 intersection = (A - std::sqrt(B + sqr(A))) / vRho2;
415 if(intersection < 0) 435 if(intersection < 0)
416 { 436 {
417 return kInfinity; 437 return kInfinity;
418 } 438 }
419 else if(std::fabs(p.z() + intersection * 439 else if(std::fabs(p.z() + intersection * v.z()) < dz + tolh)
420 { 440 {
421 return intersection; 441 return intersection;
422 } 442 }
423 else 443 else
424 { 444 {
425 return kInfinity; 445 return kInfinity;
426 } 446 }
427 } 447 }
428 } 448 }
429 else if(sqr(rho2 - paraRho2 - .25 * tol2) <= 449 else if(sqr(rho2 - paraRho2 - .25 * tol2) <= tol2 * paraRho2)
430 { 450 {
431 // If this is true we're somewhere in the 451 // If this is true we're somewhere in the border.
432 452
433 G4ThreeVector normal(p.x(), p.y(), -k1/2); 453 G4ThreeVector normal(p.x(), p.y(), -k1/2);
434 if(normal.dot(v) <= 0) 454 if(normal.dot(v) <= 0)
435 { return 0; } 455 { return 0; }
436 else 456 else
437 { return kInfinity; } 457 { return kInfinity; }
438 } 458 }
439 else 459 else
440 { 460 {
441 std::ostringstream message; 461 std::ostringstream message;
442 if(Inside(p) == kInside) 462 if(Inside(p) == kInside)
443 { 463 {
444 message << "Point p is inside! - " << Ge 464 message << "Point p is inside! - " << GetName() << G4endl;
445 } 465 }
446 else 466 else
447 { 467 {
448 message << "Likely a problem in this fun 468 message << "Likely a problem in this function, for solid: " << GetName()
449 << G4endl; 469 << G4endl;
450 } 470 }
451 message << " p = " << p * (1/mm) 471 message << " p = " << p * (1/mm) << " mm" << G4endl
452 << " v = " << v * (1/mm) 472 << " v = " << v * (1/mm) << " mm";
453 G4Exception("G4Paraboloid::DistanceToIn(p, 473 G4Exception("G4Paraboloid::DistanceToIn(p,v)", "GeomSolids1002",
454 JustWarning, message); 474 JustWarning, message);
455 return 0; 475 return 0;
456 } 476 }
457 } 477 }
458 478
459 ////////////////////////////////////////////// 479 ///////////////////////////////////////////////////////////////////////////////
460 // 480 //
461 // Calculate distance (<= actual) to closest s 481 // Calculate distance (<= actual) to closest surface of shape from outside
462 // - Return zero if point inside << 482 // - Return 0 if point inside
463 // << 483
464 G4double G4Paraboloid::DistanceToIn(const G4Th 484 G4double G4Paraboloid::DistanceToIn(const G4ThreeVector& p) const
465 { 485 {
466 G4double safz = -dz+std::fabs(p.z()); 486 G4double safz = -dz+std::fabs(p.z());
467 if(safz<0.) { safz=0.; } << 487 if(safz<0) { safz=0; }
468 G4double safr = kInfinity; 488 G4double safr = kInfinity;
469 489
470 G4double rho = p.x()*p.x()+p.y()*p.y(); 490 G4double rho = p.x()*p.x()+p.y()*p.y();
471 G4double paraRho = (p.z()-k2)/k1; 491 G4double paraRho = (p.z()-k2)/k1;
472 G4double sqrho = std::sqrt(rho); 492 G4double sqrho = std::sqrt(rho);
473 493
474 if(paraRho<0.) << 494 if(paraRho<0)
475 { 495 {
476 safr=sqrho-r2; 496 safr=sqrho-r2;
477 if(safr>safz) { safz=safr; } 497 if(safr>safz) { safz=safr; }
478 return safz; 498 return safz;
479 } 499 }
480 500
481 G4double sqprho = std::sqrt(paraRho); 501 G4double sqprho = std::sqrt(paraRho);
482 G4double dRho = sqrho-sqprho; 502 G4double dRho = sqrho-sqprho;
483 if(dRho<0.) { return safz; } << 503 if(dRho<0) { return safz; }
484 504
485 G4double talf = -2.*k1*sqprho; 505 G4double talf = -2.*k1*sqprho;
486 G4double tmp = 1+talf*talf; 506 G4double tmp = 1+talf*talf;
487 if(tmp<0.) { return safz; } 507 if(tmp<0.) { return safz; }
488 508
489 G4double salf = talf/std::sqrt(tmp); 509 G4double salf = talf/std::sqrt(tmp);
490 safr = std::fabs(dRho*salf); 510 safr = std::fabs(dRho*salf);
491 if(safr>safz) { safz=safr; } 511 if(safr>safz) { safz=safr; }
492 512
493 return safz; 513 return safz;
494 } 514 }
495 515
496 ////////////////////////////////////////////// 516 ///////////////////////////////////////////////////////////////////////////////
497 // 517 //
498 // Calculate distance to surface of shape from 518 // Calculate distance to surface of shape from 'inside'
499 // << 519
500 G4double G4Paraboloid::DistanceToOut(const G4T 520 G4double G4Paraboloid::DistanceToOut(const G4ThreeVector& p,
501 const G4Th 521 const G4ThreeVector& v,
502 const G4bo 522 const G4bool calcNorm,
503 G4bo << 523 G4bool *validNorm,
504 G4Th << 524 G4ThreeVector *n ) const
505 { 525 {
506 G4double rho2 = p.perp2(), paraRho2 = std::f 526 G4double rho2 = p.perp2(), paraRho2 = std::fabs(k1 * p.z() + k2);
507 G4double vRho2 = v.perp2(), intersection; 527 G4double vRho2 = v.perp2(), intersection;
508 G4double tol2 = kCarTolerance*kCarTolerance; 528 G4double tol2 = kCarTolerance*kCarTolerance;
509 G4double tolh = 0.5*kCarTolerance; 529 G4double tolh = 0.5*kCarTolerance;
510 530
511 if(calcNorm) { *validNorm = false; } 531 if(calcNorm) { *validNorm = false; }
512 532
513 // We have that the particle p follows the l 533 // We have that the particle p follows the line x = p + s * v
514 // meaning x = p.x() + s * v.x(), y = p.y() 534 // meaning x = p.x() + s * v.x(), y = p.y() + s * v.y() and
515 // z = p.z() + s * v.z() 535 // z = p.z() + s * v.z()
516 // The equation for all points on the surfac 536 // The equation for all points on the surface (surface expanded for
517 // to include all z) x^2 + y^2 = k1 * z + k2 537 // to include all z) x^2 + y^2 = k1 * z + k2 => .. =>
518 // => s = (A +- std::sqrt(A^2 + B)) / vRho2 538 // => s = (A +- std::sqrt(A^2 + B)) / vRho2
519 // where: 539 // where:
520 // 540 //
521 G4double A = k1 / 2 * v.z() - p.x() * v.x() 541 G4double A = k1 / 2 * v.z() - p.x() * v.x() - p.y() * v.y();
522 // 542 //
523 // and: 543 // and:
524 // 544 //
525 G4double B = (-rho2 + paraRho2) * vRho2; 545 G4double B = (-rho2 + paraRho2) * vRho2;
526 546
527 if ( rho2 < paraRho2 && sqr(rho2 - paraRho2 547 if ( rho2 < paraRho2 && sqr(rho2 - paraRho2 - 0.25 * tol2) > tol2 * paraRho2
528 && std::fabs(p.z()) < dz - kCarTolerance) 548 && std::fabs(p.z()) < dz - kCarTolerance)
529 { 549 {
530 // Make sure it's safely inside. 550 // Make sure it's safely inside.
531 551
532 if(v.z() > 0) 552 if(v.z() > 0)
533 { 553 {
534 // It's heading upwards, check where it 554 // It's heading upwards, check where it colides with the plane z = dz.
535 // When it does, is that in the surface 555 // When it does, is that in the surface of the paraboloid.
536 // z = p.z() + variable * v.z() for all 556 // z = p.z() + variable * v.z() for all points the particle can go.
537 // => variable = (z - p.z()) / v.z() so 557 // => variable = (z - p.z()) / v.z() so intersection must be:
538 558
539 intersection = (dz - p.z()) / v.z(); 559 intersection = (dz - p.z()) / v.z();
540 G4ThreeVector ip = p + intersection * v; 560 G4ThreeVector ip = p + intersection * v; // Point of intersection.
541 561
542 if(ip.perp2() < sqr(r2 + kCarTolerance)) 562 if(ip.perp2() < sqr(r2 + kCarTolerance))
543 { 563 {
544 if(calcNorm) 564 if(calcNorm)
545 { 565 {
546 *n = G4ThreeVector(0, 0, 1); 566 *n = G4ThreeVector(0, 0, 1);
547 if(r2 < tolh || ip.perp2() > sqr(r2 567 if(r2 < tolh || ip.perp2() > sqr(r2 - tolh))
548 { 568 {
549 *n += G4ThreeVector(ip.x(), ip.y() 569 *n += G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
550 *n = n->unit(); 570 *n = n->unit();
551 } 571 }
552 *validNorm = true; 572 *validNorm = true;
553 } 573 }
554 return intersection; 574 return intersection;
555 } 575 }
556 } 576 }
557 else if(v.z() < 0) 577 else if(v.z() < 0)
558 { 578 {
559 // It's heading downwards, check were it 579 // It's heading downwards, check were it colides with the plane z = -dz.
560 // When it does, is that in the surface 580 // When it does, is that in the surface of the paraboloid.
561 // z = p.z() + variable * v.z() for all 581 // z = p.z() + variable * v.z() for all points the particle can go.
562 // => variable = (z - p.z()) / v.z() so 582 // => variable = (z - p.z()) / v.z() so intersection must be:
563 583
564 intersection = (-dz - p.z()) / v.z(); 584 intersection = (-dz - p.z()) / v.z();
565 G4ThreeVector ip = p + intersection * v; 585 G4ThreeVector ip = p + intersection * v; // Point of intersection.
566 586
567 if(ip.perp2() < sqr(r1 + tolh)) 587 if(ip.perp2() < sqr(r1 + tolh))
568 { 588 {
569 if(calcNorm) 589 if(calcNorm)
570 { 590 {
571 *n = G4ThreeVector(0, 0, -1); 591 *n = G4ThreeVector(0, 0, -1);
572 if(r1 < tolh || ip.perp2() > sqr(r1 592 if(r1 < tolh || ip.perp2() > sqr(r1 - tolh))
573 { 593 {
574 *n += G4ThreeVector(ip.x(), ip.y() 594 *n += G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
575 *n = n->unit(); 595 *n = n->unit();
576 } 596 }
577 *validNorm = true; 597 *validNorm = true;
578 } 598 }
579 return intersection; 599 return intersection;
580 } 600 }
581 } 601 }
582 602
583 // Now check for collisions with paraboloi 603 // Now check for collisions with paraboloid surface.
584 604
585 if(vRho2 == 0) // Needs to be treated sepe 605 if(vRho2 == 0) // Needs to be treated seperately.
586 { 606 {
587 intersection = ((rho2 - k2)/k1 - p.z())/ 607 intersection = ((rho2 - k2)/k1 - p.z())/v.z();
588 if(calcNorm) 608 if(calcNorm)
589 { 609 {
590 G4ThreeVector intersectionP = p + v * 610 G4ThreeVector intersectionP = p + v * intersection;
591 *n = G4ThreeVector(intersectionP.x(), 611 *n = G4ThreeVector(intersectionP.x(), intersectionP.y(), -k1/2);
592 *n = n->unit(); 612 *n = n->unit();
593 613
594 *validNorm = true; 614 *validNorm = true;
595 } 615 }
596 return intersection; 616 return intersection;
597 } 617 }
598 else if( ((A <= 0) && (B >= sqr(A) * (sqr( 618 else if( ((A <= 0) && (B >= sqr(A) * (sqr(vRho2) - 1))) || (A >= 0))
599 { 619 {
600 // intersection = (A + std::sqrt(B + sqr 620 // intersection = (A + std::sqrt(B + sqr(A))) / vRho2;
601 // The above calculation has a precision 621 // The above calculation has a precision problem:
602 // known problem of solving quadratic eq 622 // known problem of solving quadratic equation with small A
603 623
604 A = A/vRho2; 624 A = A/vRho2;
605 B = (k1 * p.z() + k2 - rho2)/vRho2; 625 B = (k1 * p.z() + k2 - rho2)/vRho2;
606 intersection = B/(-A + std::sqrt(B + sqr 626 intersection = B/(-A + std::sqrt(B + sqr(A)));
607 if(calcNorm) 627 if(calcNorm)
608 { 628 {
609 G4ThreeVector intersectionP = p + v * 629 G4ThreeVector intersectionP = p + v * intersection;
610 *n = G4ThreeVector(intersectionP.x(), 630 *n = G4ThreeVector(intersectionP.x(), intersectionP.y(), -k1/2);
611 *n = n->unit(); 631 *n = n->unit();
612 *validNorm = true; 632 *validNorm = true;
613 } 633 }
614 return intersection; 634 return intersection;
615 } 635 }
616 std::ostringstream message; 636 std::ostringstream message;
617 message << "There is no intersection betwe 637 message << "There is no intersection between given line and solid!"
618 << G4endl 638 << G4endl
619 << " p = " << p << G4endl 639 << " p = " << p << G4endl
620 << " v = " << v; 640 << " v = " << v;
621 G4Exception("G4Paraboloid::DistanceToOut(p 641 G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "GeomSolids1002",
622 JustWarning, message); 642 JustWarning, message);
623 643
624 return kInfinity; 644 return kInfinity;
625 } 645 }
626 else if ( (rho2 < paraRho2 + kCarTolerance 646 else if ( (rho2 < paraRho2 + kCarTolerance
627 || sqr(rho2 - paraRho2 - 0.25 * tol2) 647 || sqr(rho2 - paraRho2 - 0.25 * tol2) < tol2 * paraRho2 )
628 && std::fabs(p.z()) < dz + tolh) 648 && std::fabs(p.z()) < dz + tolh)
629 { 649 {
630 // If this is true we're somewhere in the 650 // If this is true we're somewhere in the border.
631 651
632 G4ThreeVector normal = G4ThreeVector (p.x( 652 G4ThreeVector normal = G4ThreeVector (p.x(), p.y(), -k1/2);
633 653
634 if(std::fabs(p.z()) > dz - tolh) 654 if(std::fabs(p.z()) > dz - tolh)
635 { 655 {
636 // We're in the lower or upper edge 656 // We're in the lower or upper edge
637 // 657 //
638 if( ((v.z() > 0) && (p.z() > 0)) || ((v. 658 if( ((v.z() > 0) && (p.z() > 0)) || ((v.z() < 0) && (p.z() < 0)) )
639 { // If we're heading out of 659 { // If we're heading out of the object that is treated here
640 if(calcNorm) 660 if(calcNorm)
641 { 661 {
642 *validNorm = true; 662 *validNorm = true;
643 if(p.z() > 0) 663 if(p.z() > 0)
644 { *n = G4ThreeVector(0, 0, 1); } 664 { *n = G4ThreeVector(0, 0, 1); }
645 else 665 else
646 { *n = G4ThreeVector(0, 0, -1); } 666 { *n = G4ThreeVector(0, 0, -1); }
647 } 667 }
648 return 0; 668 return 0;
649 } 669 }
650 670
651 if(v.z() == 0) 671 if(v.z() == 0)
652 { 672 {
653 // Case where we're moving inside the 673 // Case where we're moving inside the surface needs to be
654 // treated separately. 674 // treated separately.
655 // Distance until it goes out through 675 // Distance until it goes out through a side is returned.
656 676
657 G4double r = (p.z() > 0)? r2 : r1; 677 G4double r = (p.z() > 0)? r2 : r1;
658 G4double pDotV = p.dot(v); 678 G4double pDotV = p.dot(v);
659 A = vRho2 * ( sqr(r) - sqr(p.x()) - sq 679 A = vRho2 * ( sqr(r) - sqr(p.x()) - sqr(p.y()));
660 intersection = (-pDotV + std::sqrt(A + 680 intersection = (-pDotV + std::sqrt(A + sqr(pDotV))) / vRho2;
661 681
662 if(calcNorm) 682 if(calcNorm)
663 { 683 {
664 *validNorm = true; 684 *validNorm = true;
665 685
666 *n = (G4ThreeVector(0, 0, p.z()/std: 686 *n = (G4ThreeVector(0, 0, p.z()/std::fabs(p.z()))
667 + G4ThreeVector(p.x() + v.x() * 687 + G4ThreeVector(p.x() + v.x() * intersection, p.y() + v.y()
668 * intersection, -k1/2).unit()).u 688 * intersection, -k1/2).unit()).unit();
669 } 689 }
670 return intersection; 690 return intersection;
671 } 691 }
672 } 692 }
673 // 693 //
674 // Problem in the Logic :: Following condi 694 // Problem in the Logic :: Following condition for point on upper surface
675 // and Vz<0 will 695 // and Vz<0 will return 0 (Problem #1015), but
676 // it has to retur 696 // it has to return intersection with parabolic
677 // surface or with 697 // surface or with lower plane surface (z = -dz)
678 // The logic has to be :: If not found in 698 // The logic has to be :: If not found intersection until now,
679 // do not exit but continue to search for 699 // do not exit but continue to search for possible intersection.
680 // Only for point situated on both borders 700 // Only for point situated on both borders (Z and parabolic)
681 // this condition has to be taken into acc 701 // this condition has to be taken into account and done later
682 // 702 //
683 // 703 //
684 // else if(normal.dot(v) >= 0) 704 // else if(normal.dot(v) >= 0)
685 // { 705 // {
686 // if(calcNorm) 706 // if(calcNorm)
687 // { 707 // {
688 // *validNorm = true; 708 // *validNorm = true;
689 // *n = normal.unit(); 709 // *n = normal.unit();
690 // } 710 // }
691 // return 0; 711 // return 0;
692 // } 712 // }
693 713
694 if(v.z() > 0) 714 if(v.z() > 0)
695 { 715 {
696 // Check for collision with upper edge. 716 // Check for collision with upper edge.
697 717
698 intersection = (dz - p.z()) / v.z(); 718 intersection = (dz - p.z()) / v.z();
699 G4ThreeVector ip = p + intersection * v; 719 G4ThreeVector ip = p + intersection * v;
700 720
701 if(ip.perp2() < sqr(r2 - tolh)) 721 if(ip.perp2() < sqr(r2 - tolh))
702 { 722 {
703 if(calcNorm) 723 if(calcNorm)
704 { 724 {
705 *validNorm = true; 725 *validNorm = true;
706 *n = G4ThreeVector(0, 0, 1); 726 *n = G4ThreeVector(0, 0, 1);
707 } 727 }
708 return intersection; 728 return intersection;
709 } 729 }
710 else if(ip.perp2() < sqr(r2 + tolh)) 730 else if(ip.perp2() < sqr(r2 + tolh))
711 { 731 {
712 if(calcNorm) 732 if(calcNorm)
713 { 733 {
714 *validNorm = true; 734 *validNorm = true;
715 *n = G4ThreeVector(0, 0, 1) 735 *n = G4ThreeVector(0, 0, 1)
716 + G4ThreeVector(ip.x(), ip.y(), - 736 + G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
717 *n = n->unit(); 737 *n = n->unit();
718 } 738 }
719 return intersection; 739 return intersection;
720 } 740 }
721 } 741 }
722 if( v.z() < 0) 742 if( v.z() < 0)
723 { 743 {
724 // Check for collision with lower edge. 744 // Check for collision with lower edge.
725 745
726 intersection = (-dz - p.z()) / v.z(); 746 intersection = (-dz - p.z()) / v.z();
727 G4ThreeVector ip = p + intersection * v; 747 G4ThreeVector ip = p + intersection * v;
728 748
729 if(ip.perp2() < sqr(r1 - tolh)) 749 if(ip.perp2() < sqr(r1 - tolh))
730 { 750 {
731 if(calcNorm) 751 if(calcNorm)
732 { 752 {
733 *validNorm = true; 753 *validNorm = true;
734 *n = G4ThreeVector(0, 0, -1); 754 *n = G4ThreeVector(0, 0, -1);
735 } 755 }
736 return intersection; 756 return intersection;
737 } 757 }
738 else if(ip.perp2() < sqr(r1 + tolh)) 758 else if(ip.perp2() < sqr(r1 + tolh))
739 { 759 {
740 if(calcNorm) 760 if(calcNorm)
741 { 761 {
742 *validNorm = true; 762 *validNorm = true;
743 *n = G4ThreeVector(0, 0, -1) 763 *n = G4ThreeVector(0, 0, -1)
744 + G4ThreeVector(ip.x(), ip.y(), - 764 + G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
745 *n = n->unit(); 765 *n = n->unit();
746 } 766 }
747 return intersection; 767 return intersection;
748 } 768 }
749 } 769 }
750 770
751 // Note: comparison with zero below would 771 // Note: comparison with zero below would not be correct !
752 // 772 //
753 if(std::fabs(vRho2) > tol2) // precision e 773 if(std::fabs(vRho2) > tol2) // precision error in the calculation of
754 { // intersectio 774 { // intersection = (A+std::sqrt(B+sqr(A)))/vRho2
755 A = A/vRho2; 775 A = A/vRho2;
756 B = (k1 * p.z() + k2 - rho2); 776 B = (k1 * p.z() + k2 - rho2);
757 if(std::fabs(B)>kCarTolerance) 777 if(std::fabs(B)>kCarTolerance)
758 { 778 {
759 B = (B)/vRho2; 779 B = (B)/vRho2;
760 intersection = B/(-A + std::sqrt(B + s 780 intersection = B/(-A + std::sqrt(B + sqr(A)));
761 } 781 }
762 else // Point is On 782 else // Point is On both borders: Z and parabolic
763 { // solution de 783 { // solution depends on normal.dot(v) sign
764 if(normal.dot(v) >= 0) 784 if(normal.dot(v) >= 0)
765 { 785 {
766 if(calcNorm) 786 if(calcNorm)
767 { 787 {
768 *validNorm = true; 788 *validNorm = true;
769 *n = normal.unit(); 789 *n = normal.unit();
770 } 790 }
771 return 0; 791 return 0;
772 } 792 }
773 intersection = 2.*A; 793 intersection = 2.*A;
774 } 794 }
775 } 795 }
776 else 796 else
777 { 797 {
778 intersection = ((rho2 - k2) / k1 - p.z() 798 intersection = ((rho2 - k2) / k1 - p.z()) / v.z();
779 } 799 }
780 800
781 if(calcNorm) 801 if(calcNorm)
782 { 802 {
783 *validNorm = true; 803 *validNorm = true;
784 *n = G4ThreeVector(p.x() + intersection 804 *n = G4ThreeVector(p.x() + intersection * v.x(), p.y()
785 + intersection * v.y(), - k1 / 2); 805 + intersection * v.y(), - k1 / 2);
786 *n = n->unit(); 806 *n = n->unit();
787 } 807 }
788 return intersection; 808 return intersection;
789 } 809 }
790 else 810 else
791 { 811 {
792 #ifdef G4SPECSDEBUG 812 #ifdef G4SPECSDEBUG
793 if(kOutside == Inside(p)) 813 if(kOutside == Inside(p))
794 { 814 {
795 G4Exception("G4Paraboloid::DistanceToOut 815 G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "GeomSolids1002",
796 JustWarning, "Point p is out 816 JustWarning, "Point p is outside!");
797 } 817 }
798 else 818 else
799 G4Exception("G4Paraboloid::DistanceToOut 819 G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "GeomSolids1002",
800 JustWarning, "There's an err 820 JustWarning, "There's an error in this functions code.");
801 #endif 821 #endif
802 return kInfinity; 822 return kInfinity;
803 } 823 }
804 return 0; 824 return 0;
805 } 825 }
806 826
807 ////////////////////////////////////////////// 827 ///////////////////////////////////////////////////////////////////////////////
808 // 828 //
809 // Calculate distance (<=actual) to closest su 829 // Calculate distance (<=actual) to closest surface of shape from inside
810 // << 830
811 G4double G4Paraboloid::DistanceToOut(const G4T 831 G4double G4Paraboloid::DistanceToOut(const G4ThreeVector& p) const
812 { 832 {
813 G4double safe=0.0,rho,safeR,safeZ ; 833 G4double safe=0.0,rho,safeR,safeZ ;
814 G4double tanRMax,secRMax,pRMax ; 834 G4double tanRMax,secRMax,pRMax ;
815 835
816 #ifdef G4SPECSDEBUG 836 #ifdef G4SPECSDEBUG
817 if( Inside(p) == kOutside ) 837 if( Inside(p) == kOutside )
818 { 838 {
819 G4cout << G4endl ; 839 G4cout << G4endl ;
820 DumpInfo(); 840 DumpInfo();
821 std::ostringstream message; 841 std::ostringstream message;
822 G4long oldprc = message.precision(16); << 842 G4int oldprc = message.precision(16);
823 message << "Point p is outside !?" << G4e 843 message << "Point p is outside !?" << G4endl
824 << "Position:" << G4endl 844 << "Position:" << G4endl
825 << " p.x() = " << p.x()/mm << 845 << " p.x() = " << p.x()/mm << " mm" << G4endl
826 << " p.y() = " << p.y()/mm << 846 << " p.y() = " << p.y()/mm << " mm" << G4endl
827 << " p.z() = " << p.z()/mm << 847 << " p.z() = " << p.z()/mm << " mm";
828 message.precision(oldprc) ; 848 message.precision(oldprc) ;
829 G4Exception("G4Paraboloid::DistanceToOut( 849 G4Exception("G4Paraboloid::DistanceToOut(p)", "GeomSolids1002",
830 JustWarning, message); 850 JustWarning, message);
831 } 851 }
832 #endif 852 #endif
833 853
834 rho = p.perp(); 854 rho = p.perp();
835 safeZ = dz - std::fabs(p.z()) ; 855 safeZ = dz - std::fabs(p.z()) ;
836 856
837 tanRMax = (r2 - r1)*0.5/dz ; 857 tanRMax = (r2 - r1)*0.5/dz ;
838 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ; 858 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
839 pRMax = tanRMax*p.z() + (r1+r2)*0.5 ; 859 pRMax = tanRMax*p.z() + (r1+r2)*0.5 ;
840 safeR = (pRMax - rho)/secRMax ; 860 safeR = (pRMax - rho)/secRMax ;
841 861
842 if (safeZ < safeR) { safe = safeZ; } 862 if (safeZ < safeR) { safe = safeZ; }
843 else { safe = safeR; } 863 else { safe = safeR; }
844 if ( safe < 0.5 * kCarTolerance ) { safe = 0 864 if ( safe < 0.5 * kCarTolerance ) { safe = 0; }
845 return safe ; 865 return safe ;
846 } 866 }
847 867
848 ////////////////////////////////////////////// 868 //////////////////////////////////////////////////////////////////////////
849 // 869 //
850 // G4EntityType 870 // G4EntityType
851 // << 871
852 G4GeometryType G4Paraboloid::GetEntityType() c 872 G4GeometryType G4Paraboloid::GetEntityType() const
853 { 873 {
854 return {"G4Paraboloid"}; << 874 return G4String("G4Paraboloid");
855 } 875 }
856 876
857 ////////////////////////////////////////////// 877 //////////////////////////////////////////////////////////////////////////
858 // 878 //
859 // Make a clone of the object 879 // Make a clone of the object
860 // << 880
861 G4VSolid* G4Paraboloid::Clone() const 881 G4VSolid* G4Paraboloid::Clone() const
862 { 882 {
863 return new G4Paraboloid(*this); 883 return new G4Paraboloid(*this);
864 } 884 }
865 885
866 ////////////////////////////////////////////// 886 //////////////////////////////////////////////////////////////////////////
867 // 887 //
868 // Stream object contents to an output stream 888 // Stream object contents to an output stream
869 // << 889
870 std::ostream& G4Paraboloid::StreamInfo( std::o 890 std::ostream& G4Paraboloid::StreamInfo( std::ostream& os ) const
871 { 891 {
872 G4long oldprc = os.precision(16); << 892 G4int oldprc = os.precision(16);
873 os << "------------------------------------- 893 os << "-----------------------------------------------------------\n"
874 << " *** Dump for solid - " << GetName 894 << " *** Dump for solid - " << GetName() << " ***\n"
875 << " ================================= 895 << " ===================================================\n"
876 << " Solid type: G4Paraboloid\n" 896 << " Solid type: G4Paraboloid\n"
877 << " Parameters: \n" 897 << " Parameters: \n"
878 << " z half-axis: " << dz/mm << " mm 898 << " z half-axis: " << dz/mm << " mm \n"
879 << " radius at -dz: " << r1/mm << " mm 899 << " radius at -dz: " << r1/mm << " mm \n"
880 << " radius at dz: " << r2/mm << " mm 900 << " radius at dz: " << r2/mm << " mm \n"
881 << "------------------------------------- 901 << "-----------------------------------------------------------\n";
882 os.precision(oldprc); 902 os.precision(oldprc);
883 903
884 return os; 904 return os;
885 } 905 }
886 906
887 ////////////////////////////////////////////// 907 ////////////////////////////////////////////////////////////////////
888 // 908 //
889 // GetPointOnSurface 909 // GetPointOnSurface
890 // << 910
891 G4ThreeVector G4Paraboloid::GetPointOnSurface( 911 G4ThreeVector G4Paraboloid::GetPointOnSurface() const
892 { 912 {
893 G4double A = (fSurfaceArea == 0)? CalculateS 913 G4double A = (fSurfaceArea == 0)? CalculateSurfaceArea(): fSurfaceArea;
894 G4double z = G4RandFlat::shoot(0.,1.); 914 G4double z = G4RandFlat::shoot(0.,1.);
895 G4double phi = G4RandFlat::shoot(0., twopi); 915 G4double phi = G4RandFlat::shoot(0., twopi);
896 if(pi*(sqr(r1) + sqr(r2))/A >= z) 916 if(pi*(sqr(r1) + sqr(r2))/A >= z)
897 { 917 {
898 G4double rho; 918 G4double rho;
899 if(pi * sqr(r1) / A > z) 919 if(pi * sqr(r1) / A > z)
900 { 920 {
901 rho = r1 * std::sqrt(G4RandFlat::shoot(0 921 rho = r1 * std::sqrt(G4RandFlat::shoot(0., 1.));
902 return { rho * std::cos(phi), rho * std: << 922 return G4ThreeVector(rho * std::cos(phi), rho * std::sin(phi), -dz);
903 } 923 }
904 else 924 else
905 { 925 {
906 rho = r2 * std::sqrt(G4RandFlat::shoot(0 926 rho = r2 * std::sqrt(G4RandFlat::shoot(0., 1));
907 return { rho * std::cos(phi), rho * std: << 927 return G4ThreeVector(rho * std::cos(phi), rho * std::sin(phi), dz);
908 } 928 }
909 } 929 }
910 else 930 else
911 { 931 {
912 z = G4RandFlat::shoot(0., 1.)*2*dz - dz; 932 z = G4RandFlat::shoot(0., 1.)*2*dz - dz;
913 return { std::sqrt(z*k1 + k2)*std::cos(phi << 933 return G4ThreeVector(std::sqrt(z*k1 + k2)*std::cos(phi),
914 std::sqrt(z*k1 + k2)*std::sin(phi << 934 std::sqrt(z*k1 + k2)*std::sin(phi), z);
915 } 935 }
916 } 936 }
917 937
918 ////////////////////////////////////////////// 938 /////////////////////////////////////////////////////////////////////////////
919 // 939 //
920 // Methods for visualisation 940 // Methods for visualisation
921 // << 941
922 void G4Paraboloid::DescribeYourselfTo (G4VGrap 942 void G4Paraboloid::DescribeYourselfTo (G4VGraphicsScene& scene) const
923 { 943 {
924 scene.AddSolid(*this); 944 scene.AddSolid(*this);
925 } 945 }
926 946
927 G4Polyhedron* G4Paraboloid::CreatePolyhedron ( 947 G4Polyhedron* G4Paraboloid::CreatePolyhedron () const
928 { 948 {
929 return new G4PolyhedronParaboloid(r1, r2, dz 949 return new G4PolyhedronParaboloid(r1, r2, dz, 0., twopi);
930 } 950 }
931 951
>> 952
932 G4Polyhedron* G4Paraboloid::GetPolyhedron () c 953 G4Polyhedron* G4Paraboloid::GetPolyhedron () const
933 { 954 {
934 if (fpPolyhedron == nullptr || << 955 if (!fpPolyhedron ||
935 fRebuildPolyhedron || 956 fRebuildPolyhedron ||
936 fpPolyhedron->GetNumberOfRotationStepsAt 957 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
937 fpPolyhedron->GetNumberOfRotationSteps() 958 fpPolyhedron->GetNumberOfRotationSteps())
938 { 959 {
939 G4AutoLock l(&polyhedronMutex); 960 G4AutoLock l(&polyhedronMutex);
940 delete fpPolyhedron; 961 delete fpPolyhedron;
941 fpPolyhedron = CreatePolyhedron(); 962 fpPolyhedron = CreatePolyhedron();
942 fRebuildPolyhedron = false; 963 fRebuildPolyhedron = false;
943 l.unlock(); 964 l.unlock();
944 } 965 }
945 return fpPolyhedron; 966 return fpPolyhedron;
946 } 967 }
947 968
948 #endif // !defined(G4GEOM_USE_UPARABOLOID) || 969 #endif // !defined(G4GEOM_USE_UPARABOLOID) || !defined(G4GEOM_USE_SYS_USOLIDS)
949 970