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