ATLAS Offline Software
Loading...
Searching...
No Matches
G4ShiftedCone.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
3*/
4
5//
6// The G4ShiftedCone class copied from standard G4Cons and modified to:
7// 1) have an arbitrary position along Z axis,
8// 2) represent a twopi-cone. Sectors are not supported but the
9// corresponding code kept and just commented out.
10//
11
12//
13// ********************************************************************
14// * License and Disclaimer *
15// * *
16// * The Geant4 software is copyright of the Copyright Holders of *
17// * the Geant4 Collaboration. It is provided under the terms and *
18// * conditions of the Geant4 Software License, included in the file *
19// * LICENSE and available at http://cern.ch/geant4/license . These *
20// * include a list of copyright holders. *
21// * *
22// * Neither the authors of this software system, nor their employing *
23// * institutes,nor the agencies providing financial support for this *
24// * work make any representation or warranty, express or implied, *
25// * regarding this software system or assume any liability for its *
26// * use. Please see the license in the file LICENSE and URL above *
27// * for the full disclaimer and the limitation of liability. *
28// * *
29// * This code implementation is the result of the scientific and *
30// * technical work of the GEANT4 collaboration. *
31// * By using, copying, modifying or distributing the software (or *
32// * any work based on the software) you agree to acknowledge its *
33// * use in resulting scientific publications, and indicate your *
34// * acceptance of all terms of the Geant4 Software license. *
35// ********************************************************************
36//
37//
38//
39//
40// class G4ShiftedCone
41//
42// Implementation for G4ShiftedCone class
43//
44// History:
45// 03.07.2019 A. Sukharev: copied from G4Cons
46// --------------------------------------------------------------------
47
48#include "G4ShiftedCone.h"
49
50#include "G4GeomTools.hh"
51#include "G4VoxelLimits.hh"
52#include "G4AffineTransform.hh"
53#include "G4BoundingEnvelope.hh"
54#include "G4GeometryTolerance.hh"
55
56#include "G4VPVParameterisation.hh"
57
58#include "meshdefs.hh"
59
60#include "Randomize.hh"
61
62#include "G4VGraphicsScene.hh"
63
64using namespace CLHEP;
65
67//
68// Private enum: Not for external use
69namespace {
70 // used by distanceToOut
71 enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kPZ,kMZ};
72
73 // used by normal
74 enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ};
75}
76
78//
79// constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
80// - note if pDPhi>2PI then reset to 2PI
81
82G4ShiftedCone::G4ShiftedCone( const G4String& pName,
83 G4double pZ1, G4double pZ2,
84 G4double pRmin1, G4double pRmax1,
85 G4double pRmin2, G4double pRmax2)
86// G4double pDz,
87// G4double pSPhi, G4double pDPhi)
88 : G4CSGSolid(pName),
89 kRadTolerance (G4GeometryTolerance::GetInstance()->GetRadialTolerance()),
90 kAngTolerance (G4GeometryTolerance::GetInstance()->GetAngularTolerance()),
91 fRmin1(pRmin1), fRmin2(pRmin2),
92 fRmax1(pRmax1), fRmax2(pRmax2),
93 fDz((pZ2 - pZ1) * 0.5), fZshift(pZ1 + fDz),
94 // fSPhi(0.), fDPhi(0.)
95 halfCarTolerance (kCarTolerance*0.5),
96 halfRadTolerance (kRadTolerance*0.5),
97 halfAngTolerance (kAngTolerance*0.5)
98{
99 // Check z-len
100 //
101 if ( fDz < 0 )
102 {
103 std::ostringstream message;
104 message << "Invalid Z half-length for Solid: " << GetName() << G4endl
105 << " hZ = " << fDz;
106 G4Exception("G4ShiftedCone::G4ShiftedCone()", "GeomSolids0002",
107 FatalException, message);
108 }
109
110 // Check radii
111 //
112 if (((pRmin1>=pRmax1) || (pRmin2>=pRmax2) || (pRmin1<0)) && (pRmin2<0))
113 {
114 std::ostringstream message;
115 message << "Invalid values of radii for Solid: " << GetName() << G4endl
116 << " pRmin1 = " << pRmin1 << ", pRmin2 = " << pRmin2
117 << ", pRmax1 = " << pRmax1 << ", pRmax2 = " << pRmax2;
118 G4Exception("G4ShiftedCone::G4ShiftedCone()", "GeomSolids0002",
119 FatalException, message) ;
120 }
121 if( (pRmin1 == 0.0) && (pRmin2 > 0.0) ) { fRmin1 = 1e3*kRadTolerance ; }
122 if( (pRmin2 == 0.0) && (pRmin1 > 0.0) ) { fRmin2 = 1e3*kRadTolerance ; }
123
124 // Check angles
125 //
126// CheckPhiAngles(pSPhi, pDPhi);
127}
128
130//
131// Fake default constructor - sets only member data and allocates memory
132// for usage restricted to object persistency.
133//
134G4ShiftedCone::G4ShiftedCone( __void__& a )
135 : G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.),
136 fRmin1(0.), fRmin2(0.), fRmax1(0.), fRmax2(0.),
137 fDz(0.), fZshift(0.),
138// fSPhi(0.), fDPhi(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.),
139// cosHDPhiIT(0.), sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
140// fPhiFullCone(false),
141 halfCarTolerance(0.), halfRadTolerance(0.), halfAngTolerance(0.)
142{
143}
144
146//
147// Destructor
148
152
154//
155// Copy constructor
156
157G4ShiftedCone::G4ShiftedCone(const G4ShiftedCone& rhs)
158 : G4CSGSolid(rhs), kRadTolerance(rhs.kRadTolerance),
159 kAngTolerance(rhs.kAngTolerance), fRmin1(rhs.fRmin1), fRmin2(rhs.fRmin2),
160 fRmax1(rhs.fRmax1), fRmax2(rhs.fRmax2), fDz(rhs.fDz), fZshift(rhs.fZshift),
161// fSPhi(rhs.fSPhi),
162// fDPhi(rhs.fDPhi), sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
163// cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
164// sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi), sinEPhi(rhs.sinEPhi),
165// cosEPhi(rhs.cosEPhi), fPhiFullCone(rhs.fPhiFullCone),
166 halfCarTolerance(rhs.halfCarTolerance),
167 halfRadTolerance(rhs.halfRadTolerance),
168 halfAngTolerance(rhs.halfAngTolerance)
169{
170}
171
173//
174// Assignment operator
175
176G4ShiftedCone& G4ShiftedCone::operator = (const G4ShiftedCone& rhs)
177{
178 // Check assignment to self
179 //
180 if (this == &rhs) { return *this; }
181
182 // Copy base class data
183 //
184 G4CSGSolid::operator=(rhs);
185
186 // Copy data
187 //
188 kRadTolerance = rhs.kRadTolerance;
189 kAngTolerance = rhs.kAngTolerance;
190 fRmin1 = rhs.fRmin1; fRmin2 = rhs.fRmin2;
191 fRmax1 = rhs.fRmax1; fRmax2 = rhs.fRmax2;
192 fDz = rhs.fDz; fZshift = rhs.fZshift;
193// fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
194// sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi;
195// cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT;
196// sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
197// sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
198// fPhiFullCone = rhs.fPhiFullCone;
199 halfCarTolerance = rhs.halfCarTolerance;
200 halfRadTolerance = rhs.halfRadTolerance;
201 halfAngTolerance = rhs.halfAngTolerance;
202
203 return *this;
204}
205
207//
208// Return whether point inside/outside/on surface
209
210EInside G4ShiftedCone::Inside(const G4ThreeVector& p) const
211{
212 G4double r2, rl, rh, /*pPhi,*/ tolRMin, tolRMax; // rh2, rl2 ;
213 EInside in;
214
215 G4double z = p.z() - fZshift;
216
217 if (std::fabs(z) > fDz + halfCarTolerance ) { return in = kOutside; }
218 else if(std::fabs(z) >= fDz - halfCarTolerance ) { in = kSurface; }
219 else { in = kInside; }
220
221 r2 = p.x()*p.x() + p.y()*p.y() ;
222 rl = 0.5*(fRmin2*(z + fDz) + fRmin1*(fDz - z))/fDz ;
223 rh = 0.5*(fRmax2*(z+fDz)+fRmax1*(fDz-z))/fDz;
224
225 // rh2 = rh*rh;
226
227 tolRMin = rl - halfRadTolerance;
228 if ( tolRMin < 0 ) { tolRMin = 0; }
229 tolRMax = rh + halfRadTolerance;
230
231 if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tolRMax) ) { return in = kOutside; }
232
233 if (rl) { tolRMin = rl + halfRadTolerance; }
234 else { tolRMin = 0.0; }
235 tolRMax = rh - halfRadTolerance;
236
237 if (in == kInside) // else it's kSurface already
238 {
239 if ( (r2 < tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax) ) { in = kSurface; }
240 }
241/*
242 if ( !fPhiFullCone && ((p.x() != 0.0) || (p.y() != 0.0)) )
243 {
244 pPhi = std::atan2(p.y(),p.x()) ;
245
246 if ( pPhi < fSPhi - halfAngTolerance ) { pPhi += twopi; }
247 else if ( pPhi > fSPhi + fDPhi + halfAngTolerance ) { pPhi -= twopi; }
248
249 if ( (pPhi < fSPhi - halfAngTolerance) ||
250 (pPhi > fSPhi + fDPhi + halfAngTolerance) ) { return in = kOutside; }
251
252 else if (in == kInside) // else it's kSurface anyway already
253 {
254 if ( (pPhi < fSPhi + halfAngTolerance) ||
255 (pPhi > fSPhi + fDPhi - halfAngTolerance) ) { in = kSurface; }
256 }
257 }
258 else if ( !fPhiFullCone ) { in = kSurface; }
259*/
260 return in ;
261}
262
264//
265// Dispatch to parameterisation for replication mechanism dimension
266// computation & modification.
267
268void G4ShiftedCone::ComputeDimensions( G4VPVParameterisation*,// p,
269 const G4int ,//n,
270 const G4VPhysicalVolume* )//pRep )
271{
272 std::ostringstream message;
273 message << "ComputeDimensions is not implemented for Solid: " << GetName();
274 G4Exception("G4ShiftedCone::ComputeDimensions()", "GeomSolids0002",
275 FatalException, message) ;
276// p->ComputeDimensions(*this,n,pRep) ;
277}
278
280//
281// Get bounding box
282
283void G4ShiftedCone::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
284{
285// G4double rmin = std::min(GetInnerRadiusMinusZ(),GetInnerRadiusPlusZ());
286 G4double rmax = std::max(GetOuterRadiusMinusZ(),GetOuterRadiusPlusZ());
287
288 // Find bounding box
289 //
290/* if (GetDeltaPhiAngle() < twopi)
291 {
292 G4double dz = GetZHalfLength();
293 G4TwoVector vmin,vmax;
294 G4GeomTools::DiskExtent(rmin,rmax,
295 GetSinStartPhi(),GetCosStartPhi(),
296 GetSinEndPhi(),GetCosEndPhi(),
297 vmin,vmax);
298 pMin.set(vmin.x(),vmin.y(),-dz);
299 pMax.set(vmax.x(),vmax.y(), dz);
300 }
301 else*/
302 {
303 pMin.set(-rmax,-rmax, fZshift - fDz);
304 pMax.set( rmax, rmax, fZshift + fDz);
305 }
306
307 // Check correctness of the bounding box
308 //
309 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
310 {
311 std::ostringstream message;
312 message << "Bad bounding box (min >= max) for solid: "
313 << GetName() << " !"
314 << "\npMin = " << pMin
315 << "\npMax = " << pMax;
316 G4Exception("G4ShiftedCone::BoundingLimits()", "GeomMgt0001",
317 JustWarning, message);
318 DumpInfo();
319 }
320}
321
323//
324// Calculate extent under transform and specified limit
325
326G4bool G4ShiftedCone::CalculateExtent( const EAxis pAxis,
327 const G4VoxelLimits& pVoxelLimit,
328 const G4AffineTransform& pTransform,
329 G4double& pMin,
330 G4double& pMax ) const
331{
332 G4ThreeVector bmin, bmax;
333 G4bool exist;
334
335 // Get bounding box
336 BoundingLimits(bmin,bmax);
337
338 // Check bounding box
339 G4BoundingEnvelope bbox(bmin,bmax);
340#ifdef G4BBOX_EXTENT
341 if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
342#endif
343 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
344 {
345 return exist = (pMin < pMax) ? true : false;
346 }
347
348 // Get parameters of the solid
349 G4double rmin1 = GetInnerRadiusMinusZ();
350 G4double rmax1 = GetOuterRadiusMinusZ();
351 G4double rmin2 = GetInnerRadiusPlusZ();
352 G4double rmax2 = GetOuterRadiusPlusZ();
353 G4double z1 = GetZ1();
354 G4double z2 = GetZ2();
355 G4double dphi = GetDeltaPhiAngle();
356
357 // Find bounding envelope and calculate extent
358 //
359 const G4int NSTEPS = 24; // number of steps for whole circle
360 G4double astep = twopi/NSTEPS; // max angle for one step
361 G4int ksteps = (dphi <= astep) ? 1 : (G4int)((dphi-deg)/astep) + 1;
362 G4double ang = dphi/ksteps;
363
364 G4double sinHalf = std::sin(0.5*ang);
365 G4double cosHalf = std::cos(0.5*ang);
366 G4double sinStep = 2.*sinHalf*cosHalf;
367 G4double cosStep = 1. - 2.*sinHalf*sinHalf;
368 G4double rext1 = rmax1/cosHalf;
369 G4double rext2 = rmax2/cosHalf;
370
371 // bounding envelope for full cone without hole consists of two polygons,
372 // in other cases it is a sequence of quadrilaterals
373 if (rmin1 == 0 && rmin2 == 0 && dphi == twopi)
374 {
375 G4double sinCur = sinHalf;
376 G4double cosCur = cosHalf;
377
378 G4ThreeVectorList baseA(NSTEPS),baseB(NSTEPS);
379 for (G4int k=0; k<NSTEPS; ++k)
380 {
381 baseA[k].set(rext1*cosCur,rext1*sinCur, z1);
382 baseB[k].set(rext2*cosCur,rext2*sinCur, z2);
383
384 G4double sinTmp = sinCur;
385 sinCur = sinCur*cosStep + cosCur*sinStep;
386 cosCur = cosCur*cosStep - sinTmp*sinStep;
387 }
388 std::vector<const G4ThreeVectorList *> polygons(2);
389 polygons[0] = &baseA;
390 polygons[1] = &baseB;
391 G4BoundingEnvelope benv(bmin,bmax,polygons);
392 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
393 }
394 else
395 {
396
397 G4double sinStart = GetSinStartPhi();
398 G4double cosStart = GetCosStartPhi();
399 G4double sinEnd = GetSinEndPhi();
400 G4double cosEnd = GetCosEndPhi();
401 G4double sinCur = sinStart*cosHalf + cosStart*sinHalf;
402 G4double cosCur = cosStart*cosHalf - sinStart*sinHalf;
403
404 // set quadrilaterals
405 G4ThreeVectorList pols[NSTEPS+2];
406 for (G4int k=0; k<ksteps+2; ++k) pols[k].resize(4);
407 pols[0][0].set(rmin2*cosStart,rmin2*sinStart, z2);
408 pols[0][1].set(rmin1*cosStart,rmin1*sinStart, z1);
409 pols[0][2].set(rmax1*cosStart,rmax1*sinStart, z1);
410 pols[0][3].set(rmax2*cosStart,rmax2*sinStart, z2);
411 for (G4int k=1; k<ksteps+1; ++k)
412 {
413 pols[k][0].set(rmin2*cosCur,rmin2*sinCur, z2);
414 pols[k][1].set(rmin1*cosCur,rmin1*sinCur, z1);
415 pols[k][2].set(rext1*cosCur,rext1*sinCur, z1);
416 pols[k][3].set(rext2*cosCur,rext2*sinCur, z2);
417
418 G4double sinTmp = sinCur;
419 sinCur = sinCur*cosStep + cosCur*sinStep;
420 cosCur = cosCur*cosStep - sinTmp*sinStep;
421 }
422 pols[ksteps+1][0].set(rmin2*cosEnd,rmin2*sinEnd, z2);
423 pols[ksteps+1][1].set(rmin1*cosEnd,rmin1*sinEnd, z1);
424 pols[ksteps+1][2].set(rmax1*cosEnd,rmax1*sinEnd, z1);
425 pols[ksteps+1][3].set(rmax2*cosEnd,rmax2*sinEnd, z2);
426
427 // set envelope and calculate extent
428 std::vector<const G4ThreeVectorList *> polygons;
429 polygons.resize(ksteps+2);
430 for (G4int k=0; k<ksteps+2; ++k) polygons[k] = &pols[k];
431 G4BoundingEnvelope benv(bmin,bmax,polygons);
432 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
433 }
434 return exist;
435}
436
438//
439// Return unit normal of surface closest to p
440// - note if point on z axis, ignore phi divided sides
441// - unsafe if point close to z axis a rmin=0 - no explicit checks
442
443G4ThreeVector G4ShiftedCone::SurfaceNormal( const G4ThreeVector& p) const
444{
445 G4int noSurfaces = 0;
446 G4double rho;//, pPhi;
447 G4double distZ, distRMin, distRMax;
448// G4double distSPhi = kInfinity, distEPhi = kInfinity;
449 G4double tanRMin, secRMin, pRMin, widRMin;
450 G4double tanRMax, secRMax, pRMax, widRMax;
451
452 G4ThreeVector norm, sumnorm(0.,0.,0.), nZ = G4ThreeVector(0.,0.,1.);
453 G4ThreeVector nR, nr(0.,0.,0.), nPs, nPe;
454
455 G4double z = p.z() - fZshift;
456
457 distZ = std::fabs(std::fabs(z) - fDz);
458 rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
459
460 tanRMin = (fRmin2 - fRmin1)*0.5/fDz;
461 secRMin = std::sqrt(1 + tanRMin*tanRMin);
462 pRMin = rho - z*tanRMin;
463 widRMin = fRmin2 - fDz*tanRMin;
464 distRMin = std::fabs(pRMin - widRMin)/secRMin;
465
466 tanRMax = (fRmax2 - fRmax1)*0.5/fDz;
467 secRMax = std::sqrt(1+tanRMax*tanRMax);
468 pRMax = rho - z*tanRMax;
469 widRMax = fRmax2 - fDz*tanRMax;
470 distRMax = std::fabs(pRMax - widRMax)/secRMax;
471/*
472 if (!fPhiFullCone) // Protected against (0,0,z)
473 {
474 if ( rho )
475 {
476 pPhi = std::atan2(p.y(),p.x());
477
478 if (pPhi < fSPhi-halfCarTolerance) { pPhi += twopi; }
479 else if (pPhi > fSPhi+fDPhi+halfCarTolerance) { pPhi -= twopi; }
480
481 distSPhi = std::fabs( pPhi - fSPhi );
482 distEPhi = std::fabs( pPhi - fSPhi - fDPhi );
483 }
484 else if( !(fRmin1) || !(fRmin2) )
485 {
486 distSPhi = 0.;
487 distEPhi = 0.;
488 }
489 nPs = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0);
490 nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0);
491 }*/
492 if ( rho > halfCarTolerance )
493 {
494 nR = G4ThreeVector(p.x()/rho/secRMax, p.y()/rho/secRMax, -tanRMax/secRMax);
495 if (fRmin1 || fRmin2)
496 {
497 nr = G4ThreeVector(-p.x()/rho/secRMin,-p.y()/rho/secRMin,tanRMin/secRMin);
498 }
499 }
500
501 if( distRMax <= halfCarTolerance )
502 {
503 noSurfaces ++;
504 sumnorm += nR;
505 }
506 if( (fRmin1 || fRmin2) && (distRMin <= halfCarTolerance) )
507 {
508 noSurfaces ++;
509 sumnorm += nr;
510 }
511/* if( !fPhiFullCone )
512 {
513 if (distSPhi <= halfAngTolerance)
514 {
515 noSurfaces ++;
516 sumnorm += nPs;
517 }
518 if (distEPhi <= halfAngTolerance)
519 {
520 noSurfaces ++;
521 sumnorm += nPe;
522 }
523 }*/
524 if (distZ <= halfCarTolerance)
525 {
526 noSurfaces ++;
527 if ( z >= 0.) { sumnorm += nZ; }
528 else { sumnorm -= nZ; }
529 }
530 if ( noSurfaces == 0 )
531 {
532#ifdef G4CSGDEBUG
533 G4Exception("G4ShiftedCone::SurfaceNormal(p)", "GeomSolids1002",
534 JustWarning, "Point p is not on surface !?" );
535#endif
536 norm = ApproxSurfaceNormal(p);
537 }
538 else if ( noSurfaces == 1 ) { norm = sumnorm; }
539 else { norm = sumnorm.unit(); }
540
541 return norm ;
542}
543
545//
546// Algorithm for SurfaceNormal() following the original specification
547// for points not on the surface
548
549G4ThreeVector G4ShiftedCone::ApproxSurfaceNormal( const G4ThreeVector& p ) const
550{
551 ENorm side ;
552 G4ThreeVector norm ;
553 G4double rho;//, phi ;
554 G4double distZ, distRMin, distRMax;//, distSPhi, distEPhi, distMin ;
555 G4double tanRMin, secRMin, pRMin, widRMin ;
556 G4double tanRMax, secRMax, pRMax, widRMax ;
557
558 G4double z = p.z() - fZshift;
559
560 distZ = std::fabs(std::fabs(z) - fDz) ;
561 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
562
563 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
564 secRMin = std::sqrt(1 + tanRMin*tanRMin) ;
565 pRMin = rho - z*tanRMin ;
566 widRMin = fRmin2 - fDz*tanRMin ;
567 distRMin = std::fabs(pRMin - widRMin)/secRMin ;
568
569 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
570 secRMax = std::sqrt(1+tanRMax*tanRMax) ;
571 pRMax = rho - z*tanRMax ;
572 widRMax = fRmax2 - fDz*tanRMax ;
573 distRMax = std::fabs(pRMax - widRMax)/secRMax ;
574
575 if (distRMin < distRMax) // First minimum
576 {
577 if (distZ < distRMin)
578 {
579// distMin = distZ ;
580 side = kNZ ;
581 }
582 else
583 {
584// distMin = distRMin ;
585 side = kNRMin ;
586 }
587 }
588 else
589 {
590 if (distZ < distRMax)
591 {
592// distMin = distZ ;
593 side = kNZ ;
594 }
595 else
596 {
597// distMin = distRMax ;
598 side = kNRMax ;
599 }
600 }
601/*
602 if ( !fPhiFullCone && rho ) // Protected against (0,0,z)
603 {
604 phi = std::atan2(p.y(),p.x()) ;
605
606 if (phi < 0) { phi += twopi; }
607
608 if (fSPhi < 0) { distSPhi = std::fabs(phi - (fSPhi + twopi))*rho; }
609 else { distSPhi = std::fabs(phi - fSPhi)*rho; }
610
611 distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
612
613 // Find new minimum
614
615 if (distSPhi < distEPhi)
616 {
617 if (distSPhi < distMin) { side = kNSPhi; }
618 }
619 else
620 {
621 if (distEPhi < distMin) { side = kNEPhi; }
622 }
623 }*/
624 switch (side)
625 {
626 case kNRMin: // Inner radius
627 rho *= secRMin ;
628 norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, tanRMin/secRMin) ;
629 break ;
630 case kNRMax: // Outer radius
631 rho *= secRMax ;
632 norm = G4ThreeVector(p.x()/rho, p.y()/rho, -tanRMax/secRMax) ;
633 break ;
634 case kNZ: // +/- dz
635 if (z > 0) { norm = G4ThreeVector(0,0,1); }
636 else { norm = G4ThreeVector(0,0,-1); }
637 break ;
638// case kNSPhi:
639// norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ;
640// break ;
641// case kNEPhi:
642// norm=G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ;
643// break ;
644 default: // Should never reach this case...
645 DumpInfo();
646 G4Exception("G4ShiftedCone::ApproxSurfaceNormal()",
647 "GeomSolids1002", JustWarning,
648 "Undefined side for valid surface normal to solid.");
649 break ;
650 }
651 return norm ;
652}
653
655//
656// Calculate distance to shape from outside, along normalised vector
657// - return kInfinity if no intersection, or intersection distance <= tolerance
658//
659// - Compute the intersection with the z planes
660// - if at valid r, phi, return
661//
662// -> If point is outside cone, compute intersection with rmax1*0.5
663// - if at valid phi,z return
664// - if inside outer cone, handle case when on tolerant outer cone
665// boundary and heading inwards(->0 to in)
666//
667// -> Compute intersection with inner cone, taking largest +ve root
668// - if valid (in z,phi), save intersction
669//
670// -> If phi segmented, compute intersections with phi half planes
671// - return smallest of valid phi intersections and
672// inner radius intersection
673//
674// NOTE:
675// - `if valid' implies tolerant checking of intersection points
676// - z, phi intersection from Tubs
677
678G4double G4ShiftedCone::DistanceToIn( const G4ThreeVector& p,
679 const G4ThreeVector& v ) const
680{
681 G4double snxt = kInfinity ; // snxt = default return value
682 const G4double dRmax = 50*(fRmax1+fRmax2);// 100*(Rmax1+Rmax2)/2.
683
684 G4double tanRMax,secRMax,rMaxAv;//,rMaxOAv ; // Data for cones
685 G4double tanRMin,secRMin,rMinAv;//,rMinOAv ;
686 G4double rout,rin ;
687
688 G4double tolORMin,/*tolORMin2,*/tolIRMin,tolIRMin2 ; // `generous' radii squared
689 G4double /*tolORMax2,*/tolIRMax,tolIRMax2 ;
690 G4double tolODz,tolIDz ;
691
692 G4double /*Dist,*/ sd,xi,yi,zi,ri=0.,risec,rhoi2;//,cosPsi ; // Intersection point vars
693
694 G4double t1,t2,t3,b,c,d ; // Quadratic solver variables
695 G4double nt1,nt2,nt3 ;
696// G4double Comp ;
697
698 G4ThreeVector Normal;
699
700 G4double z = p.z() - fZshift;
701
702 // Cone Precalcs
703
704 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
705 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
706 rMinAv = (fRmin1 + fRmin2)*0.5 ;
707/*
708 if (rMinAv > halfRadTolerance)
709 {
710 rMinOAv = rMinAv - halfRadTolerance ;
711 }
712 else
713 {
714 rMinOAv = 0.0 ;
715 }*/
716 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
717 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
718 rMaxAv = (fRmax1 + fRmax2)*0.5 ;
719// rMaxOAv = rMaxAv + halfRadTolerance ;
720
721 // Intersection with z-surfaces
722
723 tolIDz = fDz - halfCarTolerance ;
724 tolODz = fDz + halfCarTolerance ;
725
726 if (std::fabs(z) >= tolIDz)
727 {
728 if ( z*v.z() < 0 ) // at +Z going in -Z or visa versa
729 {
730 sd = (std::fabs(z) - fDz)/std::fabs(v.z()) ; // Z intersect distance
731
732 if( sd < 0.0 ) { sd = 0.0; } // negative dist -> zero
733
734 xi = p.x() + sd*v.x() ; // Intersection coords
735 yi = p.y() + sd*v.y() ;
736 rhoi2 = xi*xi + yi*yi ;
737
738 // Check validity of intersection
739 // Calculate (outer) tolerant radi^2 at intersecion
740
741 if (v.z() > 0)
742 {
743 tolORMin = fRmin1 - halfRadTolerance*secRMin ;
744 tolIRMin = fRmin1 + halfRadTolerance*secRMin ;
745 tolIRMax = fRmax1 - halfRadTolerance*secRMin ;
746 // tolORMax2 = (fRmax1 + halfRadTolerance*secRMax)*
747 // (fRmax1 + halfRadTolerance*secRMax) ;
748 }
749 else
750 {
751 tolORMin = fRmin2 - halfRadTolerance*secRMin ;
752 tolIRMin = fRmin2 + halfRadTolerance*secRMin ;
753 tolIRMax = fRmax2 - halfRadTolerance*secRMin ;
754 // tolORMax2 = (fRmax2 + halfRadTolerance*secRMax)*
755 // (fRmax2 + halfRadTolerance*secRMax) ;
756 }
757 if ( tolORMin > 0 )
758 {
759 // tolORMin2 = tolORMin*tolORMin ;
760 tolIRMin2 = tolIRMin*tolIRMin ;
761 }
762 else
763 {
764 // tolORMin2 = 0.0 ;
765 tolIRMin2 = 0.0 ;
766 }
767 if ( tolIRMax > 0 ) { tolIRMax2 = tolIRMax*tolIRMax; }
768 else { tolIRMax2 = 0.0; }
769
770 if ( (tolIRMin2 <= rhoi2) && (rhoi2 <= tolIRMax2) )
771 {
772/* if ( !fPhiFullCone && rhoi2 )
773 {
774 // Psi = angle made with central (average) phi of shape
775
776 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
777
778 if (cosPsi >= cosHDPhiIT) { return sd; }
779 }
780 else */
781 {
782 return sd;
783 }
784 }
785 }
786 else // On/outside extent, and heading away -> cannot intersect
787 {
788 return snxt ;
789 }
790 }
791
792// ----> Can not intersect z surfaces
793
794
795// Intersection with outer cone (possible return) and
796// inner cone (must also check phi)
797//
798// Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
799//
800// Intersects with x^2+y^2=(a*z+b)^2
801//
802// where a=tanRMax or tanRMin
803// b=rMaxAv or rMinAv
804//
805// (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a*pz+b))+px^2+py^2-(a*pz+b)^2=0 ;
806// t1 t2 t3
807//
808// \--------u-------/ \-----------v----------/ \---------w--------/
809//
810
811 t1 = 1.0 - v.z()*v.z() ;
812 t2 = p.x()*v.x() + p.y()*v.y() ;
813 t3 = p.x()*p.x() + p.y()*p.y() ;
814 rin = tanRMin*z + rMinAv ;
815 rout = tanRMax*z + rMaxAv ;
816
817 // Outer Cone Intersection
818 // Must be outside/on outer cone for valid intersection
819
820 nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;
821 nt2 = t2 - tanRMax*v.z()*rout ;
822 nt3 = t3 - rout*rout ;
823
824 if (std::fabs(nt1) > kRadTolerance) // Equation quadratic => 2 roots
825 {
826 b = nt2/nt1;
827 c = nt3/nt1;
828 d = b*b-c ;
829 if ( (nt3 > rout*rout*kRadTolerance*kRadTolerance*secRMax*secRMax)
830 || (rout < 0) )
831 {
832 // If outside real cone (should be rho-rout>kRadTolerance*0.5
833 // NOT rho^2 etc) saves a std::sqrt() at expense of accuracy
834
835 if (d >= 0)
836 {
837
838 if ((rout < 0) && (nt3 <= 0))
839 {
840 // Inside `shadow cone' with -ve radius
841 // -> 2nd root could be on real cone
842
843 if (b>0) { sd = c/(-b-std::sqrt(d)); }
844 else { sd = -b + std::sqrt(d); }
845 }
846 else
847 {
848 if ((b <= 0) && (c >= 0)) // both >=0, try smaller root
849 {
850 sd=c/(-b+std::sqrt(d));
851 }
852 else
853 {
854 if ( c <= 0 ) // second >=0
855 {
856 sd = -b + std::sqrt(d) ;
857 if((sd<0) & (sd>-halfRadTolerance)) sd=0;
858 }
859 else // both negative, travel away
860 {
861 return kInfinity ;
862 }
863 }
864 }
865 if ( sd >= 0 ) // If 'forwards'. Check z intersection
866 {
867 if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
868 { // 64 bits systems. Split long distances and recompute
869 G4double fTerm = sd-std::fmod(sd,dRmax);
870 sd = fTerm + DistanceToIn(p+fTerm*v,v);
871 }
872 zi = z + sd*v.z() ;
873
874 if (std::fabs(zi) <= tolODz)
875 {
876 // Z ok. Check phi intersection if reqd
877
878 return sd;
879/* if ( fPhiFullCone ) { return sd; }
880 else
881 {
882 xi = p.x() + sd*v.x() ;
883 yi = p.y() + sd*v.y() ;
884 ri = rMaxAv + zi*tanRMax ;
885 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
886
887 if ( cosPsi >= cosHDPhiIT ) { return sd; }
888 }*/
889 }
890 } // end if (sd>0)
891 }
892 }
893 else
894 {
895 // Inside outer cone
896 // check not inside, and heading through G4ShiftedCone (-> 0 to in)
897
898 if ( ( t3 > (rin + halfRadTolerance*secRMin)*
899 (rin + halfRadTolerance*secRMin) )
900 && (nt2 < 0) && (d >= 0) && (std::fabs(z) <= tolIDz) )
901 {
902 // Inside cones, delta r -ve, inside z extent
903 // Point is on the Surface => check Direction using Normal.dot(v)
904
905 xi = p.x() ;
906 yi = p.y() ;
907 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
908 Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
909/* if ( !fPhiFullCone )
910 {
911 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
912 if ( cosPsi >= cosHDPhiIT )
913 {
914 if ( Normal.dot(v) <= 0 ) { return 0.0; }
915 }
916 }
917 else*/
918 {
919 if ( Normal.dot(v) <= 0 ) { return 0.0; }
920 }
921 }
922 }
923 }
924 else // Single root case
925 {
926 if ( std::fabs(nt2) > kRadTolerance )
927 {
928 sd = -0.5*nt3/nt2 ;
929
930 if ( sd < 0 ) { return kInfinity; } // travel away
931 else // sd >= 0, If 'forwards'. Check z intersection
932 {
933 zi = z + sd*v.z() ;
934
935 if ((std::fabs(zi) <= tolODz) && (nt2 < 0))
936 {
937 // Z ok. Check phi intersection if reqd
938 return sd;
939/* if ( fPhiFullCone ) { return sd; }
940 else
941 {
942 xi = p.x() + sd*v.x() ;
943 yi = p.y() + sd*v.y() ;
944 ri = rMaxAv + zi*tanRMax ;
945 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
946
947 if (cosPsi >= cosHDPhiIT) { return sd; }
948 }*/
949 }
950 }
951 }
952 else // travel || cone surface from its origin
953 {
954 sd = kInfinity ;
955 }
956 }
957
958 // Inner Cone Intersection
959 // o Space is divided into 3 areas:
960 // 1) Radius greater than real inner cone & imaginary cone & outside
961 // tolerance
962 // 2) Radius less than inner or imaginary cone & outside tolarance
963 // 3) Within tolerance of real or imaginary cones
964 // - Extra checks needed for 3's intersections
965 // => lots of duplicated code
966
967 if (rMinAv)
968 {
969 nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
970 nt2 = t2 - tanRMin*v.z()*rin ;
971 nt3 = t3 - rin*rin ;
972
973 if ( nt1 )
974 {
975 if ( nt3 > rin*kRadTolerance*secRMin )
976 {
977 // At radius greater than real & imaginary cones
978 // -> 2nd root, with zi check
979
980 b = nt2/nt1 ;
981 c = nt3/nt1 ;
982 d = b*b-c ;
983 if (d >= 0) // > 0
984 {
985 if(b>0){sd = c/( -b-std::sqrt(d));}
986 else {sd = -b + std::sqrt(d) ;}
987
988 if ( sd >= 0 ) // > 0
989 {
990 if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
991 { // 64 bits systems. Split long distance and recompute
992 G4double fTerm = sd-std::fmod(sd,dRmax);
993 sd = fTerm + DistanceToIn(p+fTerm*v,v);
994 }
995 zi = z + sd*v.z() ;
996
997 if ( std::fabs(zi) <= tolODz )
998 {
999/* if ( !fPhiFullCone )
1000 {
1001 xi = p.x() + sd*v.x() ;
1002 yi = p.y() + sd*v.y() ;
1003 ri = rMinAv + zi*tanRMin ;
1004 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1005
1006 if (cosPsi >= cosHDPhiIT)
1007 {
1008 if ( sd > halfRadTolerance ) { snxt=sd; }
1009 else
1010 {
1011 // Calculate a normal vector in order to check Direction
1012
1013 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1014 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1015 if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1016 }
1017 }
1018 }
1019 else */
1020 {
1021 if ( sd > halfRadTolerance ) { return sd; }
1022 else
1023 {
1024 // Calculate a normal vector in order to check Direction
1025
1026 xi = p.x() + sd*v.x() ;
1027 yi = p.y() + sd*v.y() ;
1028 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1029 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1030 if ( Normal.dot(v) <= 0 ) { return sd; }
1031 }
1032 }
1033 }
1034 }
1035 }
1036 }
1037 else if ( nt3 < -rin*kRadTolerance*secRMin )
1038 {
1039 // Within radius of inner cone (real or imaginary)
1040 // -> Try 2nd root, with checking intersection is with real cone
1041 // -> If check fails, try 1st root, also checking intersection is
1042 // on real cone
1043
1044 b = nt2/nt1 ;
1045 c = nt3/nt1 ;
1046 d = b*b - c ;
1047
1048 if ( d >= 0 ) // > 0
1049 {
1050 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1051 else { sd = -b + std::sqrt(d); }
1052 zi = z + sd*v.z() ;
1053 ri = rMinAv + zi*tanRMin ;
1054
1055 if ( ri > 0 )
1056 {
1057 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd > 0
1058 {
1059 if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1060 { // seen on 64 bits systems. Split and recompute
1061 G4double fTerm = sd-std::fmod(sd,dRmax);
1062 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1063 }
1064/* if ( !fPhiFullCone )
1065 {
1066 xi = p.x() + sd*v.x() ;
1067 yi = p.y() + sd*v.y() ;
1068 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1069
1070 if (cosPsi >= cosHDPhiOT)
1071 {
1072 if ( sd > halfRadTolerance ) { snxt=sd; }
1073 else
1074 {
1075 // Calculate a normal vector in order to check Direction
1076
1077 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1078 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1079 if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1080 }
1081 }
1082 }
1083 else */
1084 {
1085 if( sd > halfRadTolerance ) { return sd; }
1086 else
1087 {
1088 // Calculate a normal vector in order to check Direction
1089
1090 xi = p.x() + sd*v.x() ;
1091 yi = p.y() + sd*v.y() ;
1092 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1093 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1094 if ( Normal.dot(v) <= 0 ) { return sd; }
1095 }
1096 }
1097 }
1098 }
1099 else
1100 {
1101 if (b>0) { sd = -b - std::sqrt(d); }
1102 else { sd = c/(-b+std::sqrt(d)); }
1103 zi = z + sd*v.z() ;
1104 ri = rMinAv + zi*tanRMin ;
1105
1106 if ( (sd >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1107 {
1108 if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1109 { // seen on 64 bits systems. Split and recompute
1110 G4double fTerm = sd-std::fmod(sd,dRmax);
1111 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1112 }
1113/* if ( !fPhiFullCone )
1114 {
1115 xi = p.x() + sd*v.x() ;
1116 yi = p.y() + sd*v.y() ;
1117 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1118
1119 if (cosPsi >= cosHDPhiIT)
1120 {
1121 if ( sd > halfRadTolerance ) { snxt=sd; }
1122 else
1123 {
1124 // Calculate a normal vector in order to check Direction
1125
1126 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1127 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1128 if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1129 }
1130 }
1131 }
1132 else */
1133 {
1134 if ( sd > halfRadTolerance ) { return sd; }
1135 else
1136 {
1137 // Calculate a normal vector in order to check Direction
1138
1139 xi = p.x() + sd*v.x() ;
1140 yi = p.y() + sd*v.y() ;
1141 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1142 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1143 if ( Normal.dot(v) <= 0 ) { return sd; }
1144 }
1145 }
1146 }
1147 }
1148 }
1149 }
1150 else
1151 {
1152 // Within kRadTol*0.5 of inner cone (real OR imaginary)
1153 // ----> Check not travelling through (=>0 to in)
1154 // ----> if not:
1155 // -2nd root with validity check
1156
1157 if ( std::fabs(z) <= tolODz )
1158 {
1159 if ( nt2 > 0 )
1160 {
1161 // Inside inner real cone, heading outwards, inside z range
1162
1163/* if ( !fPhiFullCone )
1164 {
1165 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
1166
1167 if (cosPsi >= cosHDPhiIT) { return 0.0; }
1168 }
1169 else */ { return 0.0; }
1170 }
1171 else
1172 {
1173 // Within z extent, but not travelling through
1174 // -> 2nd root or kInfinity if 1st root on imaginary cone
1175
1176 b = nt2/nt1 ;
1177 c = nt3/nt1 ;
1178 d = b*b - c ;
1179
1180 if ( d >= 0 ) // > 0
1181 {
1182 if (b>0) { sd = -b - std::sqrt(d); }
1183 else { sd = c/(-b+std::sqrt(d)); }
1184 zi = z + sd*v.z() ;
1185 ri = rMinAv + zi*tanRMin ;
1186
1187 if ( ri > 0 ) // 2nd root
1188 {
1189 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1190 else { sd = -b + std::sqrt(d); }
1191
1192 zi = z + sd*v.z() ;
1193
1194 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1195 {
1196 if ( sd>dRmax ) // Avoid rounding errors due to precision issue
1197 { // seen on 64 bits systems. Split and recompute
1198 G4double fTerm = sd-std::fmod(sd,dRmax);
1199 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1200 }
1201/* if ( !fPhiFullCone )
1202 {
1203 xi = p.x() + sd*v.x() ;
1204 yi = p.y() + sd*v.y() ;
1205 ri = rMinAv + zi*tanRMin ;
1206 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1207
1208 if ( cosPsi >= cosHDPhiIT ) { snxt = sd; }
1209 }
1210 else */ { return sd; }
1211 }
1212 }
1213 else { return kInfinity; }
1214 }
1215 }
1216 }
1217 else // 2nd root
1218 {
1219 b = nt2/nt1 ;
1220 c = nt3/nt1 ;
1221 d = b*b - c ;
1222
1223 if ( d > 0 )
1224 {
1225 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1226 else { sd = -b + std::sqrt(d) ; }
1227 zi = z + sd*v.z() ;
1228
1229 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1230 {
1231 if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1232 { // seen on 64 bits systems. Split and recompute
1233 G4double fTerm = sd-std::fmod(sd,dRmax);
1234 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1235 }
1236/* if ( !fPhiFullCone )
1237 {
1238 xi = p.x() + sd*v.x();
1239 yi = p.y() + sd*v.y();
1240 ri = rMinAv + zi*tanRMin ;
1241 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri;
1242
1243 if (cosPsi >= cosHDPhiIT) { snxt = sd; }
1244 }
1245 else */ { return sd; }
1246 }
1247 }
1248 }
1249 }
1250 }
1251 }
1252
1253 // Phi segment intersection
1254 //
1255 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1256 //
1257 // o NOTE: Large duplication of code between sphi & ephi checks
1258 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1259 // intersection check <=0 -> >=0
1260 // -> Should use some form of loop Construct
1261/*
1262 if ( !fPhiFullCone )
1263 {
1264 // First phi surface (starting phi)
1265
1266 Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
1267
1268 if ( Comp < 0 ) // Component in outwards normal dirn
1269 {
1270 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1271
1272 if (Dist < halfCarTolerance)
1273 {
1274 sd = Dist/Comp ;
1275
1276 if ( sd < snxt )
1277 {
1278 if ( sd < 0 ) { sd = 0.0; }
1279
1280 zi = z + sd*v.z() ;
1281
1282 if ( std::fabs(zi) <= tolODz )
1283 {
1284 xi = p.x() + sd*v.x() ;
1285 yi = p.y() + sd*v.y() ;
1286 rhoi2 = xi*xi + yi*yi ;
1287 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
1288 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
1289
1290 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
1291 {
1292 // z and r intersections good - check intersecting with
1293 // correct half-plane
1294
1295 if ((yi*cosCPhi - xi*sinCPhi) <= 0 ) { snxt = sd; }
1296 }
1297 }
1298 }
1299 }
1300 }
1301
1302 // Second phi surface (Ending phi)
1303
1304 Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
1305
1306 if ( Comp < 0 ) // Component in outwards normal dirn
1307 {
1308 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1309 if (Dist < halfCarTolerance)
1310 {
1311 sd = Dist/Comp ;
1312
1313 if ( sd < snxt )
1314 {
1315 if ( sd < 0 ) { sd = 0.0; }
1316
1317 zi = z + sd*v.z() ;
1318
1319 if (std::fabs(zi) <= tolODz)
1320 {
1321 xi = p.x() + sd*v.x() ;
1322 yi = p.y() + sd*v.y() ;
1323 rhoi2 = xi*xi + yi*yi ;
1324 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
1325 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
1326
1327 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
1328 {
1329 // z and r intersections good - check intersecting with
1330 // correct half-plane
1331
1332 if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 ) { snxt = sd; }
1333 }
1334 }
1335 }
1336 }
1337 }
1338 }*/
1339 if (snxt < halfCarTolerance) { snxt = 0.; }
1340
1341 return snxt ;
1342}
1343
1345//
1346// Calculate distance (<= actual) to closest surface of shape from outside
1347// - Calculate distance to z, radial planes
1348// - Only to phi planes if outside phi extent
1349// - Return 0 if point inside
1350
1351G4double G4ShiftedCone::DistanceToIn(const G4ThreeVector& p) const
1352{
1353 G4double safe=0.0, rho, safeR1, safeR2, safeZ;//, safePhi, cosPsi ;
1354 G4double tanRMin, secRMin, pRMin ;
1355 G4double tanRMax, secRMax, pRMax ;
1356
1357 G4double z = p.z() - fZshift;
1358
1359 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1360 safeZ = std::fabs(z) - fDz ;
1361
1362 if ( fRmin1 || fRmin2 )
1363 {
1364 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
1365 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1366 pRMin = tanRMin*z + (fRmin1 + fRmin2)*0.5 ;
1367 safeR1 = (pRMin - rho)/secRMin ;
1368
1369 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1370 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1371 pRMax = tanRMax*z + (fRmax1 + fRmax2)*0.5 ;
1372 safeR2 = (rho - pRMax)/secRMax ;
1373
1374 if ( safeR1 > safeR2) { safe = safeR1; }
1375 else { safe = safeR2; }
1376 }
1377 else
1378 {
1379 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1380 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1381 pRMax = tanRMax*z + (fRmax1 + fRmax2)*0.5 ;
1382 safe = (rho - pRMax)/secRMax ;
1383 }
1384 if ( safeZ > safe ) { safe = safeZ; }
1385
1386/* if ( !fPhiFullCone && rho )
1387 {
1388 // Psi=angle from central phi to point
1389
1390 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
1391
1392 if ( cosPsi < std::cos(fDPhi*0.5) ) // Point lies outside phi range
1393 {
1394 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0.0 )
1395 {
1396 safePhi = std::fabs(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi));
1397 }
1398 else
1399 {
1400 safePhi = std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
1401 }
1402 if ( safePhi > safe ) { safe = safePhi; }
1403 }
1404 }*/
1405 if ( safe < 0.0 ) { safe = 0.0; }
1406
1407 return safe ;
1408}
1409
1411//
1412// Calculate distance to surface of shape from 'inside', allowing for tolerance
1413// - Only Calc rmax intersection if no valid rmin intersection
1414
1415G4double G4ShiftedCone::DistanceToOut( const G4ThreeVector& p,
1416 const G4ThreeVector& v,
1417 const G4bool calcNorm,
1418 G4bool *validNorm,
1419 G4ThreeVector *n) const
1420{
1421 ESide side = kNull, sider = kNull;//, sidephi = kNull;
1422
1423 G4double snxt,srd,/*sphi,*/pdist ;
1424
1425 G4double tanRMax, secRMax, rMaxAv ; // Data for outer cone
1426 G4double tanRMin, secRMin, rMinAv ; // Data for inner cone
1427
1428 G4double t1, t2, t3, rout, rin, nt1, nt2, nt3 ;
1429 G4double b, c, d, sr2, sr3 ;
1430
1431 // Vars for intersection within tolerance
1432
1433 ESide sidetol = kNull ;
1434 G4double slentol = kInfinity ;
1435
1436 // Vars for phi intersection:
1437
1438// G4double pDistS, compS, pDistE, compE, sphi2, vphi ;
1439 G4double zi, ri, deltaRoi2, xi, yi, risec ;
1440
1441 // Z plane intersection
1442
1443 G4double z = p.z() - fZshift;
1444
1445 if ( v.z() > 0.0 )
1446 {
1447 pdist = fDz - z ;
1448
1449 if (pdist > halfCarTolerance)
1450 {
1451 snxt = pdist/v.z() ;
1452 side = kPZ ;
1453 }
1454 else
1455 {
1456 if (calcNorm)
1457 {
1458 *n = G4ThreeVector(0,0,1) ;
1459 *validNorm = true ;
1460 }
1461 return snxt = 0.0;
1462 }
1463 }
1464 else if ( v.z() < 0.0 )
1465 {
1466 pdist = fDz + z ;
1467
1468 if ( pdist > halfCarTolerance)
1469 {
1470 snxt = -pdist/v.z() ;
1471 side = kMZ ;
1472 }
1473 else
1474 {
1475 if ( calcNorm )
1476 {
1477 *n = G4ThreeVector(0,0,-1) ;
1478 *validNorm = true ;
1479 }
1480 return snxt = 0.0 ;
1481 }
1482 }
1483 else // Travel perpendicular to z axis
1484 {
1485 snxt = kInfinity ;
1486 side = kNull ;
1487 }
1488
1489 // Radial Intersections
1490 //
1491 // Intersection with outer cone (possible return) and
1492 // inner cone (must also check phi)
1493 //
1494 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1495 //
1496 // Intersects with x^2+y^2=(a*z+b)^2
1497 //
1498 // where a=tanRMax or tanRMin
1499 // b=rMaxAv or rMinAv
1500 //
1501 // (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a*pz+b))+px^2+py^2-(a*pz+b)^2=0 ;
1502 // t1 t2 t3
1503 //
1504 // \--------u-------/ \-----------v----------/ \---------w--------/
1505
1506 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1507 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1508 rMaxAv = (fRmax1 + fRmax2)*0.5 ;
1509
1510
1511 t1 = 1.0 - v.z()*v.z() ; // since v normalised
1512 t2 = p.x()*v.x() + p.y()*v.y() ;
1513 t3 = p.x()*p.x() + p.y()*p.y() ;
1514 rout = tanRMax*z + rMaxAv ;
1515
1516 nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;
1517 nt2 = t2 - tanRMax*v.z()*rout ;
1518 nt3 = t3 - rout*rout ;
1519
1520 if (v.z() > 0.0)
1521 {
1522 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
1523 - fRmax2*(fRmax2 + kRadTolerance*secRMax);
1524 }
1525 else if ( v.z() < 0.0 )
1526 {
1527 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
1528 - fRmax1*(fRmax1 + kRadTolerance*secRMax);
1529 }
1530 else
1531 {
1532 deltaRoi2 = 1.0;
1533 }
1534
1535 if ( nt1 && (deltaRoi2 > 0.0) )
1536 {
1537 // Equation quadratic => 2 roots : second root must be leaving
1538
1539 b = nt2/nt1 ;
1540 c = nt3/nt1 ;
1541 d = b*b - c ;
1542
1543 if ( d >= 0 )
1544 {
1545 // Check if on outer cone & heading outwards
1546 // NOTE: Should use rho-rout>-kRadTolerance*0.5
1547
1548 if (nt3 > -halfRadTolerance && nt2 >= 0 )
1549 {
1550 if (calcNorm)
1551 {
1552 risec = std::sqrt(t3)*secRMax ;
1553 *validNorm = true ;
1554 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
1555 }
1556 return snxt=0 ;
1557 }
1558 else
1559 {
1560 sider = kRMax ;
1561 if (b>0) { srd = -b - std::sqrt(d); }
1562 else { srd = c/(-b+std::sqrt(d)) ; }
1563
1564 zi = z + srd*v.z() ;
1565 ri = tanRMax*zi + rMaxAv ;
1566
1567 if ((ri >= 0) && (-halfRadTolerance <= srd) && (srd <= halfRadTolerance))
1568 {
1569 // An intersection within the tolerance
1570 // we will Store it in case it is good -
1571 //
1572 slentol = srd ;
1573 sidetol = kRMax ;
1574 }
1575 if ( (ri < 0) || (srd < halfRadTolerance) )
1576 {
1577 // Safety: if both roots -ve ensure that srd cannot `win'
1578 // distance to out
1579
1580 if (b>0) { sr2 = c/(-b-std::sqrt(d)); }
1581 else { sr2 = -b + std::sqrt(d); }
1582 zi = z + sr2*v.z() ;
1583 ri = tanRMax*zi + rMaxAv ;
1584
1585 if ((ri >= 0) && (sr2 > halfRadTolerance))
1586 {
1587 srd = sr2;
1588 }
1589 else
1590 {
1591 srd = kInfinity ;
1592
1593 if( (-halfRadTolerance <= sr2) && ( sr2 <= halfRadTolerance) )
1594 {
1595 // An intersection within the tolerance.
1596 // Storing it in case it is good.
1597
1598 slentol = sr2 ;
1599 sidetol = kRMax ;
1600 }
1601 }
1602 }
1603 }
1604 }
1605 else
1606 {
1607 // No intersection with outer cone & not parallel
1608 // -> already outside, no intersection
1609
1610 if ( calcNorm )
1611 {
1612 risec = std::sqrt(t3)*secRMax;
1613 *validNorm = true;
1614 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
1615 }
1616 return snxt = 0.0 ;
1617 }
1618 }
1619 else if ( nt2 && (deltaRoi2 > 0.0) )
1620 {
1621 // Linear case (only one intersection) => point outside outer cone
1622
1623 if ( calcNorm )
1624 {
1625 risec = std::sqrt(t3)*secRMax;
1626 *validNorm = true;
1627 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
1628 }
1629 return snxt = 0.0 ;
1630 }
1631 else
1632 {
1633 // No intersection -> parallel to outer cone
1634 // => Z or inner cone intersection
1635
1636 srd = kInfinity ;
1637 }
1638
1639 // Check possible intersection within tolerance
1640
1641 if ( slentol <= halfCarTolerance )
1642 {
1643 // An intersection within the tolerance was found.
1644 // We must accept it only if the momentum points outwards.
1645 //
1646 // G4ThreeVector ptTol ; // The point of the intersection
1647 // ptTol= p + slentol*v ;
1648 // ri=tanRMax*zi+rMaxAv ;
1649 //
1650 // Calculate a normal vector, as below
1651
1652 xi = p.x() + slentol*v.x();
1653 yi = p.y() + slentol*v.y();
1654 risec = std::sqrt(xi*xi + yi*yi)*secRMax;
1655 G4ThreeVector Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax);
1656
1657 if ( Normal.dot(v) > 0 ) // We will leave the Cone immediatelly
1658 {
1659 if ( calcNorm )
1660 {
1661 *n = Normal.unit() ;
1662 *validNorm = true ;
1663 }
1664 return snxt = 0.0 ;
1665 }
1666 else // On the surface, but not heading out so we ignore this intersection
1667 { // (as it is within tolerance).
1668 slentol = kInfinity ;
1669 }
1670 }
1671
1672 // Inner Cone intersection
1673
1674 if ( fRmin1 || fRmin2 )
1675 {
1676 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
1677 nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
1678
1679 if ( nt1 )
1680 {
1681 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1682 rMinAv = (fRmin1 + fRmin2)*0.5 ;
1683 rin = tanRMin*z + rMinAv ;
1684 nt2 = t2 - tanRMin*v.z()*rin ;
1685 nt3 = t3 - rin*rin ;
1686
1687 // Equation quadratic => 2 roots : first root must be leaving
1688
1689 b = nt2/nt1 ;
1690 c = nt3/nt1 ;
1691 d = b*b - c ;
1692
1693 if ( d >= 0.0 )
1694 {
1695 // NOTE: should be rho-rin<kRadTolerance*0.5,
1696 // but using squared versions for efficiency
1697
1698 if (nt3 < kRadTolerance*(rin + kRadTolerance*0.25))
1699 {
1700 if ( nt2 < 0.0 )
1701 {
1702 if (calcNorm) { *validNorm = false; }
1703 return snxt = 0.0;
1704 }
1705 }
1706 else
1707 {
1708 if (b>0) { sr2 = -b - std::sqrt(d); }
1709 else { sr2 = c/(-b+std::sqrt(d)); }
1710 zi = z + sr2*v.z() ;
1711 ri = tanRMin*zi + rMinAv ;
1712
1713 if( (ri>=0.0)&&(-halfRadTolerance<=sr2)&&(sr2<=halfRadTolerance) )
1714 {
1715 // An intersection within the tolerance
1716 // storing it in case it is good.
1717
1718 slentol = sr2 ;
1719 sidetol = kRMax ;
1720 }
1721 if( (ri<0) || (sr2 < halfRadTolerance) )
1722 {
1723 if (b>0) { sr3 = c/(-b-std::sqrt(d)); }
1724 else { sr3 = -b + std::sqrt(d) ; }
1725
1726 // Safety: if both roots -ve ensure that srd cannot `win'
1727 // distancetoout
1728
1729 if ( sr3 > halfRadTolerance )
1730 {
1731 if( sr3 < srd )
1732 {
1733 zi = z + sr3*v.z() ;
1734 ri = tanRMin*zi + rMinAv ;
1735
1736 if ( ri >= 0.0 )
1737 {
1738 srd=sr3 ;
1739 sider=kRMin ;
1740 }
1741 }
1742 }
1743 else if ( sr3 > -halfRadTolerance )
1744 {
1745 // Intersection in tolerance. Store to check if it's good
1746
1747 slentol = sr3 ;
1748 sidetol = kRMin ;
1749 }
1750 }
1751 else if ( (sr2 < srd) && (sr2 > halfCarTolerance) )
1752 {
1753 srd = sr2 ;
1754 sider = kRMin ;
1755 }
1756 else if (sr2 > -halfCarTolerance)
1757 {
1758 // Intersection in tolerance. Store to check if it's good
1759
1760 slentol = sr2 ;
1761 sidetol = kRMin ;
1762 }
1763 if( slentol <= halfCarTolerance )
1764 {
1765 // An intersection within the tolerance was found.
1766 // We must accept it only if the momentum points outwards.
1767
1768 G4ThreeVector Normal ;
1769
1770 // Calculate a normal vector, as below
1771
1772 xi = p.x() + slentol*v.x() ;
1773 yi = p.y() + slentol*v.y() ;
1774 if( sidetol==kRMax )
1775 {
1776 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
1777 Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
1778 }
1779 else
1780 {
1781 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1782 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1783 }
1784 if( Normal.dot(v) > 0 )
1785 {
1786 // We will leave the cone immediately
1787
1788 if( calcNorm )
1789 {
1790 *n = Normal.unit() ;
1791 *validNorm = true ;
1792 }
1793 return snxt = 0.0 ;
1794 }
1795 else
1796 {
1797 // On the surface, but not heading out so we ignore this
1798 // intersection (as it is within tolerance).
1799
1800 slentol = kInfinity ;
1801 }
1802 }
1803 }
1804 }
1805 }
1806 }
1807
1808 // Linear case => point outside inner cone ---> outer cone intersect
1809 //
1810 // Phi Intersection
1811/*
1812 if ( !fPhiFullCone )
1813 {
1814 // add angle calculation with correction
1815 // of the difference in domain of atan2 and Sphi
1816
1817 vphi = std::atan2(v.y(),v.x()) ;
1818
1819 if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; }
1820 else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; }
1821
1822 if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
1823 {
1824 // pDist -ve when inside
1825
1826 pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
1827 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1828
1829 // Comp -ve when in direction of outwards normal
1830
1831 compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
1832 compE = sinEPhi*v.x() - cosEPhi*v.y() ;
1833
1834 sidephi = kNull ;
1835
1836 if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
1837 && (pDistE <= halfCarTolerance) ) )
1838 || ( (fDPhi > pi) && !((pDistS > halfCarTolerance)
1839 && (pDistE > halfCarTolerance) ) ) )
1840 {
1841 // Inside both phi *full* planes
1842 if ( compS < 0 )
1843 {
1844 sphi = pDistS/compS ;
1845 if (sphi >= -halfCarTolerance)
1846 {
1847 xi = p.x() + sphi*v.x() ;
1848 yi = p.y() + sphi*v.y() ;
1849
1850 // Check intersecting with correct half-plane
1851 // (if not -> no intersect)
1852 //
1853 if ( (std::fabs(xi)<=kCarTolerance)
1854 && (std::fabs(yi)<=kCarTolerance) )
1855 {
1856 sidephi= kSPhi;
1857 if ( ( fSPhi-halfAngTolerance <= vphi )
1858 && ( fSPhi+fDPhi+halfAngTolerance >=vphi ) )
1859 {
1860 sphi = kInfinity;
1861 }
1862 }
1863 else
1864 if ( (yi*cosCPhi-xi*sinCPhi)>=0 )
1865 {
1866 sphi = kInfinity ;
1867 }
1868 else
1869 {
1870 sidephi = kSPhi ;
1871 if ( pDistS > -halfCarTolerance )
1872 {
1873 sphi = 0.0 ; // Leave by sphi immediately
1874 }
1875 }
1876 }
1877 else
1878 {
1879 sphi = kInfinity ;
1880 }
1881 }
1882 else
1883 {
1884 sphi = kInfinity ;
1885 }
1886
1887 if ( compE < 0 )
1888 {
1889 sphi2 = pDistE/compE ;
1890
1891 // Only check further if < starting phi intersection
1892 //
1893 if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
1894 {
1895 xi = p.x() + sphi2*v.x() ;
1896 yi = p.y() + sphi2*v.y() ;
1897
1898 // Check intersecting with correct half-plane
1899
1900 if ( (std::fabs(xi)<=kCarTolerance)
1901 && (std::fabs(yi)<=kCarTolerance) )
1902 {
1903 // Leaving via ending phi
1904
1905 if(!( (fSPhi-halfAngTolerance <= vphi)
1906 && (fSPhi+fDPhi+halfAngTolerance >= vphi) ) )
1907 {
1908 sidephi = kEPhi ;
1909 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; }
1910 else { sphi = 0.0; }
1911 }
1912 }
1913 else // Check intersecting with correct half-plane
1914 if ( yi*cosCPhi-xi*sinCPhi >= 0 )
1915 {
1916 // Leaving via ending phi
1917
1918 sidephi = kEPhi ;
1919 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; }
1920 else { sphi = 0.0; }
1921 }
1922 }
1923 }
1924 }
1925 else
1926 {
1927 sphi = kInfinity ;
1928 }
1929 }
1930 else
1931 {
1932 // On z axis + travel not || to z axis -> if phi of vector direction
1933 // within phi of shape, Step limited by rmax, else Step =0
1934
1935 if ( (fSPhi-halfAngTolerance <= vphi)
1936 && (vphi <= fSPhi+fDPhi+halfAngTolerance) )
1937 {
1938 sphi = kInfinity ;
1939 }
1940 else
1941 {
1942 sidephi = kSPhi ; // arbitrary
1943 sphi = 0.0 ;
1944 }
1945 }
1946 if ( sphi < snxt ) // Order intersecttions
1947 {
1948 snxt=sphi ;
1949 side=sidephi ;
1950 }
1951 }
1952*/
1953 if ( srd < snxt ) // Order intersections
1954 {
1955 snxt = srd ;
1956 side = sider ;
1957 }
1958 if (calcNorm)
1959 {
1960 switch(side)
1961 { // Note: returned vector not normalised
1962 case kRMax: // (divide by frmax for unit vector)
1963 xi = p.x() + snxt*v.x() ;
1964 yi = p.y() + snxt*v.y() ;
1965 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
1966 *n = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
1967 *validNorm = true ;
1968 break ;
1969 case kRMin:
1970 *validNorm = false ; // Rmin is inconvex
1971 break ;
1972/* case kSPhi:
1973 if ( fDPhi <= pi )
1974 {
1975 *n = G4ThreeVector(sinSPhi, -cosSPhi, 0);
1976 *validNorm = true ;
1977 }
1978 else
1979 {
1980 *validNorm = false ;
1981 }
1982 break ;
1983 case kEPhi:
1984 if ( fDPhi <= pi )
1985 {
1986 *n = G4ThreeVector(-sinEPhi, cosEPhi, 0);
1987 *validNorm = true ;
1988 }
1989 else
1990 {
1991 *validNorm = false ;
1992 }
1993 break ;*/
1994 case kPZ:
1995 *n = G4ThreeVector(0,0,1) ;
1996 *validNorm = true ;
1997 break ;
1998 case kMZ:
1999 *n = G4ThreeVector(0,0,-1) ;
2000 *validNorm = true ;
2001 break ;
2002 default:
2003 G4cout << G4endl ;
2004 DumpInfo();
2005 std::ostringstream message;
2006 G4int oldprc = message.precision(16) ;
2007 message << "Undefined side for valid surface normal to solid."
2008 << G4endl
2009 << "Position:" << G4endl << G4endl
2010 << "p.x() = " << p.x()/mm << " mm" << G4endl
2011 << "p.y() = " << p.y()/mm << " mm" << G4endl
2012 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
2013 << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
2014 << " mm" << G4endl << G4endl ;
2015 if( p.x() != 0. || p.y() != 0.)
2016 {
2017 message << "point phi = " << std::atan2(p.y(),p.x())/degree
2018 << " degree" << G4endl << G4endl ;
2019 }
2020 message << "Direction:" << G4endl << G4endl
2021 << "v.x() = " << v.x() << G4endl
2022 << "v.y() = " << v.y() << G4endl
2023 << "v.z() = " << v.z() << G4endl<< G4endl
2024 << "Proposed distance :" << G4endl<< G4endl
2025 << "snxt = " << snxt/mm << " mm" << G4endl ;
2026 message.precision(oldprc) ;
2027 G4Exception("G4ShiftedCone::DistanceToOut(p,v,..)","GeomSolids1002",
2028 JustWarning, message) ;
2029 break ;
2030 }
2031 }
2032 if (snxt < halfCarTolerance) { snxt = 0.; }
2033
2034 return snxt ;
2035}
2036
2038//
2039// Calculate distance (<=actual) to closest surface of shape from inside
2040
2041G4double G4ShiftedCone::DistanceToOut(const G4ThreeVector& p) const
2042{
2043 G4double safe=0.0, rho, safeR1, safeR2, safeZ;//, safePhi;
2044 G4double tanRMin, secRMin, pRMin;
2045 G4double tanRMax, secRMax, pRMax;
2046
2047#ifdef G4CSGDEBUG
2048 if( Inside(p) == kOutside )
2049 {
2050 G4int oldprc=G4cout.precision(16) ;
2051 G4cout << G4endl ;
2052 DumpInfo();
2053 G4cout << "Position:" << G4endl << G4endl ;
2054 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
2055 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
2056 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
2057 G4cout << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
2058 << " mm" << G4endl << G4endl ;
2059 if( (p.x() != 0.) || (p.x() != 0.) )
2060 {
2061 G4cout << "point phi = " << std::atan2(p.y(),p.x())/degree
2062 << " degree" << G4endl << G4endl ;
2063 }
2064 G4cout.precision(oldprc) ;
2065 G4Exception("G4ShiftedCone::DistanceToOut(p)", "GeomSolids1002",
2066 JustWarning, "Point p is outside !?" );
2067 }
2068#endif
2069
2070 G4double z = p.z() - fZshift;
2071
2072 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
2073 safeZ = fDz - std::fabs(z) ;
2074
2075 if (fRmin1 || fRmin2)
2076 {
2077 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
2078 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
2079 pRMin = tanRMin*z + (fRmin1 + fRmin2)*0.5 ;
2080 safeR1 = (rho - pRMin)/secRMin ;
2081 }
2082 else
2083 {
2084 safeR1 = kInfinity ;
2085 }
2086
2087 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
2088 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
2089 pRMax = tanRMax*z + (fRmax1+fRmax2)*0.5 ;
2090 safeR2 = (pRMax - rho)/secRMax ;
2091
2092 if (safeR1 < safeR2) { safe = safeR1; }
2093 else { safe = safeR2; }
2094 if (safeZ < safe) { safe = safeZ ; }
2095
2096 // Check if phi divided, Calc distances closest phi plane
2097/*
2098 if (!fPhiFullCone)
2099 {
2100 // Above/below central phi of G4ShiftedCone?
2101
2102 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
2103 {
2104 safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
2105 }
2106 else
2107 {
2108 safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
2109 }
2110 if (safePhi < safe) { safe = safePhi; }
2111 }*/
2112 if ( safe < 0 ) { safe = 0; }
2113
2114 return safe ;
2115}
2116
2118//
2119// GetEntityType
2120
2121G4GeometryType G4ShiftedCone::GetEntityType() const
2122{
2123 return G4String("G4ShiftedCone");
2124}
2125
2127//
2128// Make a clone of the object
2129//
2130G4VSolid* G4ShiftedCone::Clone() const
2131{
2132 return new G4ShiftedCone(*this);
2133}
2134
2136//
2137// Stream object contents to an output stream
2138
2139std::ostream& G4ShiftedCone::StreamInfo(std::ostream& os) const
2140{
2141 G4int oldprc = os.precision(16);
2142 os << "-----------------------------------------------------------\n"
2143 << " *** Dump for solid - " << GetName() << " ***\n"
2144 << " ===================================================\n"
2145 << " Solid type: G4ShiftedCone\n"
2146 << " Parameters: \n"
2147 << " inside -fDz radius: " << fRmin1/mm << " mm \n"
2148 << " outside -fDz radius: " << fRmax1/mm << " mm \n"
2149 << " inside +fDz radius: " << fRmin2/mm << " mm \n"
2150 << " outside +fDz radius: " << fRmax2/mm << " mm \n"
2151 << " Z1 : " << GetZ1()/mm << " mm \n"
2152 << " Z2 : " << GetZ2()/mm << " mm \n"
2153// << " starting angle of segment: " << fSPhi/degree << " degrees \n"
2154// << " delta angle of segment : " << fDPhi/degree << " degrees \n"
2155 << "-----------------------------------------------------------\n";
2156 os.precision(oldprc);
2157
2158 return os;
2159}
2160
2161
2162
2164//
2165// GetPointOnSurface
2166
2167G4ThreeVector G4ShiftedCone::GetPointOnSurface() const
2168{
2169 // declare working variables
2170 //
2171 G4double rone = (fRmax1-fRmax2)/(2.*fDz);
2172 G4double rtwo = (fRmin1-fRmin2)/(2.*fDz);
2173 G4double qone = (fRmax1 == fRmax2) ? 0. : fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2);
2174 G4double qtwo = (fRmin1 == fRmin2) ? 0. : fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2);
2175
2176 G4double slin = std::hypot(fRmin1-fRmin2, 2.*fDz);
2177 G4double slout = std::hypot(fRmax1-fRmax2, 2.*fDz);
2178 G4double Aone = 0.5*GetDeltaPhiAngle()*(fRmax2 + fRmax1)*slout; // outer surface
2179 G4double Atwo = 0.5*GetDeltaPhiAngle()*(fRmin2 + fRmin1)*slin; // inner surface
2180 G4double Athree = 0.5*GetDeltaPhiAngle()*(fRmax1*fRmax1-fRmin1*fRmin1); // base at -Dz
2181 G4double Afour = 0.5*GetDeltaPhiAngle()*(fRmax2*fRmax2-fRmin2*fRmin2); // base at +Dz
2182 G4double Afive = fDz*(fRmax1-fRmin1+fRmax2-fRmin2); // phi section
2183
2184 G4double phi = G4RandFlat::shoot(GetStartPhiAngle(),GetStartPhiAngle() + GetDeltaPhiAngle());
2185 G4double cosu = std::cos(phi);
2186 G4double sinu = std::sin(phi);
2187 G4double rRand1 = GetRadiusInRing(fRmin1, fRmax1);
2188 G4double rRand2 = GetRadiusInRing(fRmin2, fRmax2);
2189
2190 G4bool fPhiFullCone = true;
2191 if ( (GetStartPhiAngle() == 0.) && fPhiFullCone ) { Afive = 0.; }
2192 G4double chose = G4RandFlat::shoot(0.,Aone+Atwo+Athree+Afour+2.*Afive);
2193
2194 if( (chose >= 0.) && (chose < Aone) ) // outer surface
2195 {
2196 if(fRmax1 != fRmax2)
2197 {
2198 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2199 return G4ThreeVector (rone*cosu*(qone-zRand),
2200 rone*sinu*(qone-zRand), zRand + fZshift);
2201 }
2202 else
2203 {
2204 return G4ThreeVector(fRmax1*cosu, fRmax2*sinu,
2205 G4RandFlat::shoot(-1.*fDz,fDz) + fZshift);
2206 }
2207 }
2208 else if( (chose >= Aone) && (chose < Aone + Atwo) ) // inner surface
2209 {
2210 if(fRmin1 != fRmin2)
2211 {
2212 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2213 return G4ThreeVector (rtwo*cosu*(qtwo-zRand),
2214 rtwo*sinu*(qtwo-zRand), zRand + fZshift);
2215 }
2216 else
2217 {
2218 return G4ThreeVector(fRmin1*cosu, fRmin2*sinu,
2219 G4RandFlat::shoot(-1.*fDz,fDz) + fZshift);
2220 }
2221 }
2222 else if( (chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree) ) // base at -Dz
2223 {
2224 return G4ThreeVector (rRand1*cosu, rRand1*sinu, -1*fDz + fZshift);
2225 }
2226 else if( (chose >= Aone + Atwo + Athree)
2227 && (chose < Aone + Atwo + Athree + Afour) ) // base at +Dz
2228 {
2229 return G4ThreeVector (rRand2*cosu,rRand2*sinu,fDz + fZshift);
2230 }
2231 else if( (chose >= Aone + Atwo + Athree + Afour) // SPhi section
2232 && (chose < Aone + Atwo + Athree + Afour + Afive) )
2233 {
2234 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2235 rRand1 = G4RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
2236 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
2237 return G4ThreeVector (rRand1*GetCosStartPhi(),
2238 rRand1*GetSinStartPhi(), zRand + fZshift);
2239 }
2240 else // SPhi+DPhi section
2241 {
2242 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2243 rRand1 = G4RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
2244 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
2245 return G4ThreeVector (rRand1*GetCosEndPhi(),
2246 rRand1*GetSinEndPhi(), zRand + fZshift);
2247 }
2248}
2249
2251//
2252// Methods for visualisation
2253
2254void G4ShiftedCone::DescribeYourselfTo (G4VGraphicsScene& scene) const
2255{
2256 scene.AddSolid (*this);
2257}
2258
2259G4Polyhedron* G4ShiftedCone::CreatePolyhedron () const
2260{
2261 G4double rmin[2] = { GetRmin1(), GetRmin2() };
2262 G4double rmax[2] = { GetRmax1(), GetRmax2() };
2263 G4double z[2] = { GetZ1(), GetZ2() };
2264 return new G4PolyhedronPcon(
2265 GetStartPhiAngle(), GetDeltaPhiAngle(), 2, z, rmin, rmax
2266 );
2267}
phi
Scalar phi() const
phi method
Definition AmgMatrixBasePlugin.h:67
a
static Double_t a
Definition LArPhysWaveHECTool.cxx:38
deg
#define deg
Definition SbPolyhedron.cxx:17
z
#define z
twopi
constexpr double twopi
Definition VertexPointEstimator.cxx:16
G4ShiftedCone::Inside
EInside Inside(const G4ThreeVector &p) const
Definition G4ShiftedCone.cxx:210
G4ShiftedCone::DistanceToOut
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
Definition G4ShiftedCone.cxx:1415
G4ShiftedCone::halfRadTolerance
G4double halfRadTolerance
Definition G4ShiftedCone.h:240
G4ShiftedCone::DistanceToIn
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition G4ShiftedCone.cxx:678
G4ShiftedCone::GetRmin2
G4double GetRmin2() const
G4ShiftedCone::GetOuterRadiusMinusZ
G4double GetOuterRadiusMinusZ() const
G4ShiftedCone::DescribeYourselfTo
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition G4ShiftedCone.cxx:2254
G4ShiftedCone::GetDeltaPhiAngle
G4double GetDeltaPhiAngle() const
G4ShiftedCone::GetInnerRadiusMinusZ
G4double GetInnerRadiusMinusZ() const
G4ShiftedCone::operator=
G4ShiftedCone & operator=(const G4ShiftedCone &rhs)
Definition G4ShiftedCone.cxx:176
G4ShiftedCone::GetRmax2
G4double GetRmax2() const
G4ShiftedCone::CreatePolyhedron
G4Polyhedron * CreatePolyhedron() const
Definition G4ShiftedCone.cxx:2259
G4ShiftedCone::GetRmax1
G4double GetRmax1() const
G4ShiftedCone::GetSinEndPhi
G4double GetSinEndPhi() const
G4ShiftedCone::ComputeDimensions
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition G4ShiftedCone.cxx:268
G4ShiftedCone::GetEntityType
G4GeometryType GetEntityType() const
Definition G4ShiftedCone.cxx:2121
G4ShiftedCone::halfAngTolerance
G4double halfAngTolerance
Definition G4ShiftedCone.h:240
G4ShiftedCone::GetStartPhiAngle
G4double GetStartPhiAngle() const
G4ShiftedCone::ApproxSurfaceNormal
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
Definition G4ShiftedCone.cxx:549
G4ShiftedCone::GetPointOnSurface
G4ThreeVector GetPointOnSurface() const
Definition G4ShiftedCone.cxx:2167
G4ShiftedCone::GetInnerRadiusPlusZ
G4double GetInnerRadiusPlusZ() const
G4ShiftedCone::G4ShiftedCone
G4ShiftedCone(const G4String &pName, G4double pZ1, G4double pZ2, G4double pRmin1, G4double pRmax1, G4double pRmin2, G4double pRmax2)
Definition G4ShiftedCone.cxx:82
G4ShiftedCone::GetZ1
G4double GetZ1() const
G4ShiftedCone::halfCarTolerance
G4double halfCarTolerance
Definition G4ShiftedCone.h:240
G4ShiftedCone::kRadTolerance
G4double kRadTolerance
Definition G4ShiftedCone.h:222
G4ShiftedCone::GetZ2
G4double GetZ2() const
G4ShiftedCone::BoundingLimits
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
Definition G4ShiftedCone.cxx:283
G4ShiftedCone::SurfaceNormal
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition G4ShiftedCone.cxx:443
G4ShiftedCone::CalculateExtent
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
Definition G4ShiftedCone.cxx:326
G4ShiftedCone::Clone
G4VSolid * Clone() const
Definition G4ShiftedCone.cxx:2130
G4ShiftedCone::GetOuterRadiusPlusZ
G4double GetOuterRadiusPlusZ() const
G4ShiftedCone::GetCosEndPhi
G4double GetCosEndPhi() const
G4ShiftedCone::kAngTolerance
G4double kAngTolerance
Definition G4ShiftedCone.h:222
G4ShiftedCone::GetCosStartPhi
G4double GetCosStartPhi() const
G4ShiftedCone::GetSinStartPhi
G4double GetSinStartPhi() const
G4ShiftedCone::GetRmin1
G4double GetRmin1() const
G4ShiftedCone::StreamInfo
std::ostream & StreamInfo(std::ostream &os) const
Definition G4ShiftedCone.cxx:2139
Generated on for ATLAS Offline Software by doxygen 1.14.0