ATLAS Offline Software
Loading...
Searching...
No Matches
G4ShiftedCone.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2026 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 //coverity[DEADCODE]
645 default: // Should never reach this case...
646 DumpInfo();
647 G4Exception("G4ShiftedCone::ApproxSurfaceNormal()",
648 "GeomSolids1002", JustWarning,
649 "Undefined side for valid surface normal to solid.");
650 break ;
651 }
652 return norm ;
653}
654
656//
657// Calculate distance to shape from outside, along normalised vector
658// - return kInfinity if no intersection, or intersection distance <= tolerance
659//
660// - Compute the intersection with the z planes
661// - if at valid r, phi, return
662//
663// -> If point is outside cone, compute intersection with rmax1*0.5
664// - if at valid phi,z return
665// - if inside outer cone, handle case when on tolerant outer cone
666// boundary and heading inwards(->0 to in)
667//
668// -> Compute intersection with inner cone, taking largest +ve root
669// - if valid (in z,phi), save intersction
670//
671// -> If phi segmented, compute intersections with phi half planes
672// - return smallest of valid phi intersections and
673// inner radius intersection
674//
675// NOTE:
676// - `if valid' implies tolerant checking of intersection points
677// - z, phi intersection from Tubs
678
679G4double G4ShiftedCone::DistanceToIn( const G4ThreeVector& p,
680 const G4ThreeVector& v ) const
681{
682 G4double snxt = kInfinity ; // snxt = default return value
683 const G4double dRmax = 50*(fRmax1+fRmax2);// 100*(Rmax1+Rmax2)/2.
684
685 G4double tanRMax,secRMax,rMaxAv;//,rMaxOAv ; // Data for cones
686 G4double tanRMin,secRMin,rMinAv;//,rMinOAv ;
687 G4double rout,rin ;
688
689 G4double tolORMin,/*tolORMin2,*/tolIRMin,tolIRMin2 ; // `generous' radii squared
690 G4double /*tolORMax2,*/tolIRMax,tolIRMax2 ;
691 G4double tolODz,tolIDz ;
692
693 G4double /*Dist,*/ sd,xi,yi,zi,ri=0.,risec,rhoi2;//,cosPsi ; // Intersection point vars
694
695 G4double t1,t2,t3,b,c,d ; // Quadratic solver variables
696 G4double nt1,nt2,nt3 ;
697// G4double Comp ;
698
699 G4ThreeVector Normal;
700
701 G4double z = p.z() - fZshift;
702
703 // Cone Precalcs
704
705 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
706 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
707 rMinAv = (fRmin1 + fRmin2)*0.5 ;
708/*
709 if (rMinAv > halfRadTolerance)
710 {
711 rMinOAv = rMinAv - halfRadTolerance ;
712 }
713 else
714 {
715 rMinOAv = 0.0 ;
716 }*/
717 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
718 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
719 rMaxAv = (fRmax1 + fRmax2)*0.5 ;
720// rMaxOAv = rMaxAv + halfRadTolerance ;
721
722 // Intersection with z-surfaces
723
724 tolIDz = fDz - halfCarTolerance ;
725 tolODz = fDz + halfCarTolerance ;
726
727 if (std::fabs(z) >= tolIDz)
728 {
729 if ( z*v.z() < 0 ) // at +Z going in -Z or visa versa
730 {
731 sd = (std::fabs(z) - fDz)/std::fabs(v.z()) ; // Z intersect distance
732
733 if( sd < 0.0 ) { sd = 0.0; } // negative dist -> zero
734
735 xi = p.x() + sd*v.x() ; // Intersection coords
736 yi = p.y() + sd*v.y() ;
737 rhoi2 = xi*xi + yi*yi ;
738
739 // Check validity of intersection
740 // Calculate (outer) tolerant radi^2 at intersecion
741
742 if (v.z() > 0)
743 {
744 tolORMin = fRmin1 - halfRadTolerance*secRMin ;
745 tolIRMin = fRmin1 + halfRadTolerance*secRMin ;
746 tolIRMax = fRmax1 - halfRadTolerance*secRMin ;
747 // tolORMax2 = (fRmax1 + halfRadTolerance*secRMax)*
748 // (fRmax1 + halfRadTolerance*secRMax) ;
749 }
750 else
751 {
752 tolORMin = fRmin2 - halfRadTolerance*secRMin ;
753 tolIRMin = fRmin2 + halfRadTolerance*secRMin ;
754 tolIRMax = fRmax2 - halfRadTolerance*secRMin ;
755 // tolORMax2 = (fRmax2 + halfRadTolerance*secRMax)*
756 // (fRmax2 + halfRadTolerance*secRMax) ;
757 }
758 if ( tolORMin > 0 )
759 {
760 // tolORMin2 = tolORMin*tolORMin ;
761 tolIRMin2 = tolIRMin*tolIRMin ;
762 }
763 else
764 {
765 // tolORMin2 = 0.0 ;
766 tolIRMin2 = 0.0 ;
767 }
768 if ( tolIRMax > 0 ) { tolIRMax2 = tolIRMax*tolIRMax; }
769 else { tolIRMax2 = 0.0; }
770
771 if ( (tolIRMin2 <= rhoi2) && (rhoi2 <= tolIRMax2) )
772 {
773/* if ( !fPhiFullCone && rhoi2 )
774 {
775 // Psi = angle made with central (average) phi of shape
776
777 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
778
779 if (cosPsi >= cosHDPhiIT) { return sd; }
780 }
781 else */
782 {
783 return sd;
784 }
785 }
786 }
787 else // On/outside extent, and heading away -> cannot intersect
788 {
789 return snxt ;
790 }
791 }
792
793// ----> Can not intersect z surfaces
794
795
796// Intersection with outer cone (possible return) and
797// inner cone (must also check phi)
798//
799// Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
800//
801// Intersects with x^2+y^2=(a*z+b)^2
802//
803// where a=tanRMax or tanRMin
804// b=rMaxAv or rMinAv
805//
806// (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 ;
807// t1 t2 t3
808//
809// \--------u-------/ \-----------v----------/ \---------w--------/
810//
811
812 t1 = 1.0 - v.z()*v.z() ;
813 t2 = p.x()*v.x() + p.y()*v.y() ;
814 t3 = p.x()*p.x() + p.y()*p.y() ;
815 rin = tanRMin*z + rMinAv ;
816 rout = tanRMax*z + rMaxAv ;
817
818 // Outer Cone Intersection
819 // Must be outside/on outer cone for valid intersection
820
821 nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;
822 nt2 = t2 - tanRMax*v.z()*rout ;
823 nt3 = t3 - rout*rout ;
824
825 if (std::fabs(nt1) > kRadTolerance) // Equation quadratic => 2 roots
826 {
827 b = nt2/nt1;
828 c = nt3/nt1;
829 d = b*b-c ;
830 if ( (nt3 > rout*rout*kRadTolerance*kRadTolerance*secRMax*secRMax)
831 || (rout < 0) )
832 {
833 // If outside real cone (should be rho-rout>kRadTolerance*0.5
834 // NOT rho^2 etc) saves a std::sqrt() at expense of accuracy
835
836 if (d >= 0)
837 {
838
839 if ((rout < 0) && (nt3 <= 0))
840 {
841 // Inside `shadow cone' with -ve radius
842 // -> 2nd root could be on real cone
843
844 if (b>0) { sd = c/(-b-std::sqrt(d)); }
845 else { sd = -b + std::sqrt(d); }
846 }
847 else
848 {
849 if ((b <= 0) && (c >= 0)) // both >=0, try smaller root
850 {
851 sd=c/(-b+std::sqrt(d));
852 }
853 else
854 {
855 if ( c <= 0 ) // second >=0
856 {
857 sd = -b + std::sqrt(d) ;
858 if((sd<0) & (sd>-halfRadTolerance)) sd=0;
859 }
860 else // both negative, travel away
861 {
862 return kInfinity ;
863 }
864 }
865 }
866 if ( sd >= 0 ) // If 'forwards'. Check z intersection
867 {
868 if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
869 { // 64 bits systems. Split long distances and recompute
870 G4double fTerm = sd-std::fmod(sd,dRmax);
871 sd = fTerm + DistanceToIn(p+fTerm*v,v);
872 }
873 zi = z + sd*v.z() ;
874
875 if (std::fabs(zi) <= tolODz)
876 {
877 // Z ok. Check phi intersection if reqd
878
879 return sd;
880/* if ( fPhiFullCone ) { return sd; }
881 else
882 {
883 xi = p.x() + sd*v.x() ;
884 yi = p.y() + sd*v.y() ;
885 ri = rMaxAv + zi*tanRMax ;
886 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
887
888 if ( cosPsi >= cosHDPhiIT ) { return sd; }
889 }*/
890 }
891 } // end if (sd>0)
892 }
893 }
894 else
895 {
896 // Inside outer cone
897 // check not inside, and heading through G4ShiftedCone (-> 0 to in)
898
899 if ( ( t3 > (rin + halfRadTolerance*secRMin)*
900 (rin + halfRadTolerance*secRMin) )
901 && (nt2 < 0) && (d >= 0) && (std::fabs(z) <= tolIDz) )
902 {
903 // Inside cones, delta r -ve, inside z extent
904 // Point is on the Surface => check Direction using Normal.dot(v)
905
906 xi = p.x() ;
907 yi = p.y() ;
908 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
909 Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
910/* if ( !fPhiFullCone )
911 {
912 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
913 if ( cosPsi >= cosHDPhiIT )
914 {
915 if ( Normal.dot(v) <= 0 ) { return 0.0; }
916 }
917 }
918 else*/
919 {
920 if ( Normal.dot(v) <= 0 ) { return 0.0; }
921 }
922 }
923 }
924 }
925 else // Single root case
926 {
927 if ( std::fabs(nt2) > kRadTolerance )
928 {
929 sd = -0.5*nt3/nt2 ;
930
931 if ( sd < 0 ) { return kInfinity; } // travel away
932 else // sd >= 0, If 'forwards'. Check z intersection
933 {
934 zi = z + sd*v.z() ;
935
936 if ((std::fabs(zi) <= tolODz) && (nt2 < 0))
937 {
938 // Z ok. Check phi intersection if reqd
939 return sd;
940/* if ( fPhiFullCone ) { return sd; }
941 else
942 {
943 xi = p.x() + sd*v.x() ;
944 yi = p.y() + sd*v.y() ;
945 ri = rMaxAv + zi*tanRMax ;
946 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
947
948 if (cosPsi >= cosHDPhiIT) { return sd; }
949 }*/
950 }
951 }
952 }
953 }
954
955 // Inner Cone Intersection
956 // o Space is divided into 3 areas:
957 // 1) Radius greater than real inner cone & imaginary cone & outside
958 // tolerance
959 // 2) Radius less than inner or imaginary cone & outside tolarance
960 // 3) Within tolerance of real or imaginary cones
961 // - Extra checks needed for 3's intersections
962 // => lots of duplicated code
963
964 if (rMinAv)
965 {
966 nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
967 nt2 = t2 - tanRMin*v.z()*rin ;
968 nt3 = t3 - rin*rin ;
969
970 if ( nt1 )
971 {
972 if ( nt3 > rin*kRadTolerance*secRMin )
973 {
974 // At radius greater than real & imaginary cones
975 // -> 2nd root, with zi check
976
977 b = nt2/nt1 ;
978 c = nt3/nt1 ;
979 d = b*b-c ;
980 if (d >= 0) // > 0
981 {
982 if(b>0){sd = c/( -b-std::sqrt(d));}
983 else {sd = -b + std::sqrt(d) ;}
984
985 if ( sd >= 0 ) // > 0
986 {
987 if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
988 { // 64 bits systems. Split long distance and recompute
989 G4double fTerm = sd-std::fmod(sd,dRmax);
990 sd = fTerm + DistanceToIn(p+fTerm*v,v);
991 }
992 zi = z + sd*v.z() ;
993
994 if ( std::fabs(zi) <= tolODz )
995 {
996/* if ( !fPhiFullCone )
997 {
998 xi = p.x() + sd*v.x() ;
999 yi = p.y() + sd*v.y() ;
1000 ri = rMinAv + zi*tanRMin ;
1001 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1002
1003 if (cosPsi >= cosHDPhiIT)
1004 {
1005 if ( sd > halfRadTolerance ) { snxt=sd; }
1006 else
1007 {
1008 // Calculate a normal vector in order to check Direction
1009
1010 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1011 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1012 if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1013 }
1014 }
1015 }
1016 else */
1017 {
1018 if ( sd > halfRadTolerance ) { return sd; }
1019 else
1020 {
1021 // Calculate a normal vector in order to check Direction
1022
1023 xi = p.x() + sd*v.x() ;
1024 yi = p.y() + sd*v.y() ;
1025 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1026 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1027 if ( Normal.dot(v) <= 0 ) { return sd; }
1028 }
1029 }
1030 }
1031 }
1032 }
1033 }
1034 else if ( nt3 < -rin*kRadTolerance*secRMin )
1035 {
1036 // Within radius of inner cone (real or imaginary)
1037 // -> Try 2nd root, with checking intersection is with real cone
1038 // -> If check fails, try 1st root, also checking intersection is
1039 // on real cone
1040
1041 b = nt2/nt1 ;
1042 c = nt3/nt1 ;
1043 d = b*b - c ;
1044
1045 if ( d >= 0 ) // > 0
1046 {
1047 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1048 else { sd = -b + std::sqrt(d); }
1049 zi = z + sd*v.z() ;
1050 ri = rMinAv + zi*tanRMin ;
1051
1052 if ( ri > 0 )
1053 {
1054 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd > 0
1055 {
1056 if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1057 { // seen on 64 bits systems. Split and recompute
1058 G4double fTerm = sd-std::fmod(sd,dRmax);
1059 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1060 }
1061/* if ( !fPhiFullCone )
1062 {
1063 xi = p.x() + sd*v.x() ;
1064 yi = p.y() + sd*v.y() ;
1065 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1066
1067 if (cosPsi >= cosHDPhiOT)
1068 {
1069 if ( sd > halfRadTolerance ) { snxt=sd; }
1070 else
1071 {
1072 // Calculate a normal vector in order to check Direction
1073
1074 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1075 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1076 if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1077 }
1078 }
1079 }
1080 else */
1081 {
1082 if( sd > halfRadTolerance ) { return sd; }
1083 else
1084 {
1085 // Calculate a normal vector in order to check Direction
1086
1087 xi = p.x() + sd*v.x() ;
1088 yi = p.y() + sd*v.y() ;
1089 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1090 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1091 if ( Normal.dot(v) <= 0 ) { return sd; }
1092 }
1093 }
1094 }
1095 }
1096 else
1097 {
1098 if (b>0) { sd = -b - std::sqrt(d); }
1099 else { sd = c/(-b+std::sqrt(d)); }
1100 zi = z + sd*v.z() ;
1101 ri = rMinAv + zi*tanRMin ;
1102
1103 if ( (sd >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1104 {
1105 if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1106 { // seen on 64 bits systems. Split and recompute
1107 G4double fTerm = sd-std::fmod(sd,dRmax);
1108 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1109 }
1110/* if ( !fPhiFullCone )
1111 {
1112 xi = p.x() + sd*v.x() ;
1113 yi = p.y() + sd*v.y() ;
1114 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1115
1116 if (cosPsi >= cosHDPhiIT)
1117 {
1118 if ( sd > halfRadTolerance ) { snxt=sd; }
1119 else
1120 {
1121 // Calculate a normal vector in order to check Direction
1122
1123 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1124 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1125 if ( Normal.dot(v) <= 0 ) { snxt = sd; }
1126 }
1127 }
1128 }
1129 else */
1130 {
1131 if ( sd > halfRadTolerance ) { return sd; }
1132 else
1133 {
1134 // Calculate a normal vector in order to check Direction
1135
1136 xi = p.x() + sd*v.x() ;
1137 yi = p.y() + sd*v.y() ;
1138 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1139 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1140 if ( Normal.dot(v) <= 0 ) { return sd; }
1141 }
1142 }
1143 }
1144 }
1145 }
1146 }
1147 else
1148 {
1149 // Within kRadTol*0.5 of inner cone (real OR imaginary)
1150 // ----> Check not travelling through (=>0 to in)
1151 // ----> if not:
1152 // -2nd root with validity check
1153
1154 if ( std::fabs(z) <= tolODz )
1155 {
1156 if ( nt2 > 0 )
1157 {
1158 // Inside inner real cone, heading outwards, inside z range
1159
1160/* if ( !fPhiFullCone )
1161 {
1162 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
1163
1164 if (cosPsi >= cosHDPhiIT) { return 0.0; }
1165 }
1166 else */ { return 0.0; }
1167 }
1168 else
1169 {
1170 // Within z extent, but not travelling through
1171 // -> 2nd root or kInfinity if 1st root on imaginary cone
1172
1173 b = nt2/nt1 ;
1174 c = nt3/nt1 ;
1175 d = b*b - c ;
1176
1177 if ( d >= 0 ) // > 0
1178 {
1179 if (b>0) { sd = -b - std::sqrt(d); }
1180 else { sd = c/(-b+std::sqrt(d)); }
1181 zi = z + sd*v.z() ;
1182 ri = rMinAv + zi*tanRMin ;
1183
1184 if ( ri > 0 ) // 2nd root
1185 {
1186 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1187 else { sd = -b + std::sqrt(d); }
1188
1189 zi = z + sd*v.z() ;
1190
1191 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1192 {
1193 if ( sd>dRmax ) // Avoid rounding errors due to precision issue
1194 { // seen on 64 bits systems. Split and recompute
1195 G4double fTerm = sd-std::fmod(sd,dRmax);
1196 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1197 }
1198/* if ( !fPhiFullCone )
1199 {
1200 xi = p.x() + sd*v.x() ;
1201 yi = p.y() + sd*v.y() ;
1202 ri = rMinAv + zi*tanRMin ;
1203 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1204
1205 if ( cosPsi >= cosHDPhiIT ) { snxt = sd; }
1206 }
1207 else */ { return sd; }
1208 }
1209 }
1210 else { return kInfinity; }
1211 }
1212 }
1213 }
1214 else // 2nd root
1215 {
1216 b = nt2/nt1 ;
1217 c = nt3/nt1 ;
1218 d = b*b - c ;
1219
1220 if ( d > 0 )
1221 {
1222 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1223 else { sd = -b + std::sqrt(d) ; }
1224 zi = z + sd*v.z() ;
1225
1226 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1227 {
1228 if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1229 { // seen on 64 bits systems. Split and recompute
1230 G4double fTerm = sd-std::fmod(sd,dRmax);
1231 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1232 }
1233/* if ( !fPhiFullCone )
1234 {
1235 xi = p.x() + sd*v.x();
1236 yi = p.y() + sd*v.y();
1237 ri = rMinAv + zi*tanRMin ;
1238 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri;
1239
1240 if (cosPsi >= cosHDPhiIT) { snxt = sd; }
1241 }
1242 else */ { return sd; }
1243 }
1244 }
1245 }
1246 }
1247 }
1248 }
1249
1250 // Phi segment intersection
1251 //
1252 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1253 //
1254 // o NOTE: Large duplication of code between sphi & ephi checks
1255 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1256 // intersection check <=0 -> >=0
1257 // -> Should use some form of loop Construct
1258/*
1259 if ( !fPhiFullCone )
1260 {
1261 // First phi surface (starting phi)
1262
1263 Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
1264
1265 if ( Comp < 0 ) // Component in outwards normal dirn
1266 {
1267 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1268
1269 if (Dist < halfCarTolerance)
1270 {
1271 sd = Dist/Comp ;
1272
1273 if ( sd < snxt )
1274 {
1275 if ( sd < 0 ) { sd = 0.0; }
1276
1277 zi = z + sd*v.z() ;
1278
1279 if ( std::fabs(zi) <= tolODz )
1280 {
1281 xi = p.x() + sd*v.x() ;
1282 yi = p.y() + sd*v.y() ;
1283 rhoi2 = xi*xi + yi*yi ;
1284 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
1285 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
1286
1287 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
1288 {
1289 // z and r intersections good - check intersecting with
1290 // correct half-plane
1291
1292 if ((yi*cosCPhi - xi*sinCPhi) <= 0 ) { snxt = sd; }
1293 }
1294 }
1295 }
1296 }
1297 }
1298
1299 // Second phi surface (Ending phi)
1300
1301 Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
1302
1303 if ( Comp < 0 ) // Component in outwards normal dirn
1304 {
1305 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1306 if (Dist < halfCarTolerance)
1307 {
1308 sd = Dist/Comp ;
1309
1310 if ( sd < snxt )
1311 {
1312 if ( sd < 0 ) { sd = 0.0; }
1313
1314 zi = z + sd*v.z() ;
1315
1316 if (std::fabs(zi) <= tolODz)
1317 {
1318 xi = p.x() + sd*v.x() ;
1319 yi = p.y() + sd*v.y() ;
1320 rhoi2 = xi*xi + yi*yi ;
1321 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
1322 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
1323
1324 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
1325 {
1326 // z and r intersections good - check intersecting with
1327 // correct half-plane
1328
1329 if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 ) { snxt = sd; }
1330 }
1331 }
1332 }
1333 }
1334 }
1335 }*/
1336 if (snxt < halfCarTolerance) { snxt = 0.; }
1337
1338 return snxt ;
1339}
1340
1342//
1343// Calculate distance (<= actual) to closest surface of shape from outside
1344// - Calculate distance to z, radial planes
1345// - Only to phi planes if outside phi extent
1346// - Return 0 if point inside
1347
1348G4double G4ShiftedCone::DistanceToIn(const G4ThreeVector& p) const
1349{
1350 G4double safe=0.0, rho, safeR1, safeR2, safeZ;//, safePhi, cosPsi ;
1351 G4double tanRMin, secRMin, pRMin ;
1352 G4double tanRMax, secRMax, pRMax ;
1353
1354 G4double z = p.z() - fZshift;
1355
1356 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1357 safeZ = std::fabs(z) - fDz ;
1358
1359 if ( fRmin1 || fRmin2 )
1360 {
1361 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
1362 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1363 pRMin = tanRMin*z + (fRmin1 + fRmin2)*0.5 ;
1364 safeR1 = (pRMin - rho)/secRMin ;
1365
1366 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1367 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1368 pRMax = tanRMax*z + (fRmax1 + fRmax2)*0.5 ;
1369 safeR2 = (rho - pRMax)/secRMax ;
1370
1371 if ( safeR1 > safeR2) { safe = safeR1; }
1372 else { safe = safeR2; }
1373 }
1374 else
1375 {
1376 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1377 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1378 pRMax = tanRMax*z + (fRmax1 + fRmax2)*0.5 ;
1379 safe = (rho - pRMax)/secRMax ;
1380 }
1381 if ( safeZ > safe ) { safe = safeZ; }
1382
1383/* if ( !fPhiFullCone && rho )
1384 {
1385 // Psi=angle from central phi to point
1386
1387 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
1388
1389 if ( cosPsi < std::cos(fDPhi*0.5) ) // Point lies outside phi range
1390 {
1391 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0.0 )
1392 {
1393 safePhi = std::fabs(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi));
1394 }
1395 else
1396 {
1397 safePhi = std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
1398 }
1399 if ( safePhi > safe ) { safe = safePhi; }
1400 }
1401 }*/
1402 if ( safe < 0.0 ) { safe = 0.0; }
1403
1404 return safe ;
1405}
1406
1408//
1409// Calculate distance to surface of shape from 'inside', allowing for tolerance
1410// - Only Calc rmax intersection if no valid rmin intersection
1411
1412G4double G4ShiftedCone::DistanceToOut( const G4ThreeVector& p,
1413 const G4ThreeVector& v,
1414 const G4bool calcNorm,
1415 G4bool *validNorm,
1416 G4ThreeVector *n) const
1417{
1418 ESide side = kNull, sider = kNull;//, sidephi = kNull;
1419
1420 G4double snxt,srd,/*sphi,*/pdist ;
1421
1422 G4double tanRMax, secRMax, rMaxAv ; // Data for outer cone
1423 G4double tanRMin, secRMin, rMinAv ; // Data for inner cone
1424
1425 G4double t1, t2, t3, rout, rin, nt1, nt2, nt3 ;
1426 G4double b, c, d, sr2, sr3 ;
1427
1428 // Vars for intersection within tolerance
1429
1430 ESide sidetol = kNull ;
1431 G4double slentol = kInfinity ;
1432
1433 // Vars for phi intersection:
1434
1435// G4double pDistS, compS, pDistE, compE, sphi2, vphi ;
1436 G4double zi, ri, deltaRoi2, xi, yi, risec ;
1437
1438 // Z plane intersection
1439
1440 G4double z = p.z() - fZshift;
1441
1442 if ( v.z() > 0.0 )
1443 {
1444 pdist = fDz - z ;
1445
1446 if (pdist > halfCarTolerance)
1447 {
1448 snxt = pdist/v.z() ;
1449 side = kPZ ;
1450 }
1451 else
1452 {
1453 if (calcNorm)
1454 {
1455 *n = G4ThreeVector(0,0,1) ;
1456 *validNorm = true ;
1457 }
1458 return snxt = 0.0;
1459 }
1460 }
1461 else if ( v.z() < 0.0 )
1462 {
1463 pdist = fDz + z ;
1464
1465 if ( pdist > halfCarTolerance)
1466 {
1467 snxt = -pdist/v.z() ;
1468 side = kMZ ;
1469 }
1470 else
1471 {
1472 if ( calcNorm )
1473 {
1474 *n = G4ThreeVector(0,0,-1) ;
1475 *validNorm = true ;
1476 }
1477 return snxt = 0.0 ;
1478 }
1479 }
1480 else // Travel perpendicular to z axis
1481 {
1482 snxt = kInfinity ;
1483 side = kNull ;
1484 }
1485
1486 // Radial Intersections
1487 //
1488 // Intersection with outer cone (possible return) and
1489 // inner cone (must also check phi)
1490 //
1491 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1492 //
1493 // Intersects with x^2+y^2=(a*z+b)^2
1494 //
1495 // where a=tanRMax or tanRMin
1496 // b=rMaxAv or rMinAv
1497 //
1498 // (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 ;
1499 // t1 t2 t3
1500 //
1501 // \--------u-------/ \-----------v----------/ \---------w--------/
1502
1503 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1504 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1505 rMaxAv = (fRmax1 + fRmax2)*0.5 ;
1506
1507
1508 t1 = 1.0 - v.z()*v.z() ; // since v normalised
1509 t2 = p.x()*v.x() + p.y()*v.y() ;
1510 t3 = p.x()*p.x() + p.y()*p.y() ;
1511 rout = tanRMax*z + rMaxAv ;
1512
1513 nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;
1514 nt2 = t2 - tanRMax*v.z()*rout ;
1515 nt3 = t3 - rout*rout ;
1516
1517 if (v.z() > 0.0)
1518 {
1519 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
1520 - fRmax2*(fRmax2 + kRadTolerance*secRMax);
1521 }
1522 else if ( v.z() < 0.0 )
1523 {
1524 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
1525 - fRmax1*(fRmax1 + kRadTolerance*secRMax);
1526 }
1527 else
1528 {
1529 deltaRoi2 = 1.0;
1530 }
1531
1532 if ( nt1 && (deltaRoi2 > 0.0) )
1533 {
1534 // Equation quadratic => 2 roots : second root must be leaving
1535
1536 b = nt2/nt1 ;
1537 c = nt3/nt1 ;
1538 d = b*b - c ;
1539
1540 if ( d >= 0 )
1541 {
1542 // Check if on outer cone & heading outwards
1543 // NOTE: Should use rho-rout>-kRadTolerance*0.5
1544
1545 if (nt3 > -halfRadTolerance && nt2 >= 0 )
1546 {
1547 if (calcNorm)
1548 {
1549 risec = std::sqrt(t3)*secRMax ;
1550 *validNorm = true ;
1551 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
1552 }
1553 return snxt=0 ;
1554 }
1555 else
1556 {
1557 sider = kRMax ;
1558 if (b>0) { srd = -b - std::sqrt(d); }
1559 else { srd = c/(-b+std::sqrt(d)) ; }
1560
1561 zi = z + srd*v.z() ;
1562 ri = tanRMax*zi + rMaxAv ;
1563
1564 if ((ri >= 0) && (-halfRadTolerance <= srd) && (srd <= halfRadTolerance))
1565 {
1566 // An intersection within the tolerance
1567 // we will Store it in case it is good -
1568 //
1569 slentol = srd ;
1570 sidetol = kRMax ;
1571 }
1572 if ( (ri < 0) || (srd < halfRadTolerance) )
1573 {
1574 // Safety: if both roots -ve ensure that srd cannot `win'
1575 // distance to out
1576
1577 if (b>0) { sr2 = c/(-b-std::sqrt(d)); }
1578 else { sr2 = -b + std::sqrt(d); }
1579 zi = z + sr2*v.z() ;
1580 ri = tanRMax*zi + rMaxAv ;
1581
1582 if ((ri >= 0) && (sr2 > halfRadTolerance))
1583 {
1584 srd = sr2;
1585 }
1586 else
1587 {
1588 srd = kInfinity ;
1589
1590 if( (-halfRadTolerance <= sr2) && ( sr2 <= halfRadTolerance) )
1591 {
1592 // An intersection within the tolerance.
1593 // Storing it in case it is good.
1594
1595 slentol = sr2 ;
1596 sidetol = kRMax ;
1597 }
1598 }
1599 }
1600 }
1601 }
1602 else
1603 {
1604 // No intersection with outer cone & not parallel
1605 // -> already outside, no intersection
1606
1607 if ( calcNorm )
1608 {
1609 risec = std::sqrt(t3)*secRMax;
1610 *validNorm = true;
1611 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
1612 }
1613 return snxt = 0.0 ;
1614 }
1615 }
1616 else if ( nt2 && (deltaRoi2 > 0.0) )
1617 {
1618 // Linear case (only one intersection) => point outside outer cone
1619
1620 if ( calcNorm )
1621 {
1622 risec = std::sqrt(t3)*secRMax;
1623 *validNorm = true;
1624 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
1625 }
1626 return snxt = 0.0 ;
1627 }
1628 else
1629 {
1630 // No intersection -> parallel to outer cone
1631 // => Z or inner cone intersection
1632
1633 srd = kInfinity ;
1634 }
1635
1636 // Check possible intersection within tolerance
1637
1638 if ( slentol <= halfCarTolerance )
1639 {
1640 // An intersection within the tolerance was found.
1641 // We must accept it only if the momentum points outwards.
1642 //
1643 // G4ThreeVector ptTol ; // The point of the intersection
1644 // ptTol= p + slentol*v ;
1645 // ri=tanRMax*zi+rMaxAv ;
1646 //
1647 // Calculate a normal vector, as below
1648
1649 xi = p.x() + slentol*v.x();
1650 yi = p.y() + slentol*v.y();
1651 risec = std::sqrt(xi*xi + yi*yi)*secRMax;
1652 G4ThreeVector Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax);
1653
1654 if ( Normal.dot(v) > 0 ) // We will leave the Cone immediatelly
1655 {
1656 if ( calcNorm )
1657 {
1658 *n = Normal.unit() ;
1659 *validNorm = true ;
1660 }
1661 return snxt = 0.0 ;
1662 }
1663 else // On the surface, but not heading out so we ignore this intersection
1664 { // (as it is within tolerance).
1665 slentol = kInfinity ;
1666 }
1667 }
1668
1669 // Inner Cone intersection
1670
1671 if ( fRmin1 || fRmin2 )
1672 {
1673 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
1674 nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
1675
1676 if ( nt1 )
1677 {
1678 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1679 rMinAv = (fRmin1 + fRmin2)*0.5 ;
1680 rin = tanRMin*z + rMinAv ;
1681 nt2 = t2 - tanRMin*v.z()*rin ;
1682 nt3 = t3 - rin*rin ;
1683
1684 // Equation quadratic => 2 roots : first root must be leaving
1685
1686 b = nt2/nt1 ;
1687 c = nt3/nt1 ;
1688 d = b*b - c ;
1689
1690 if ( d >= 0.0 )
1691 {
1692 // NOTE: should be rho-rin<kRadTolerance*0.5,
1693 // but using squared versions for efficiency
1694
1695 if (nt3 < kRadTolerance*(rin + kRadTolerance*0.25))
1696 {
1697 if ( nt2 < 0.0 )
1698 {
1699 if (calcNorm) { *validNorm = false; }
1700 return snxt = 0.0;
1701 }
1702 }
1703 else
1704 {
1705 if (b>0) { sr2 = -b - std::sqrt(d); }
1706 else { sr2 = c/(-b+std::sqrt(d)); }
1707 zi = z + sr2*v.z() ;
1708 ri = tanRMin*zi + rMinAv ;
1709
1710 if( (ri>=0.0)&&(-halfRadTolerance<=sr2)&&(sr2<=halfRadTolerance) )
1711 {
1712 // An intersection within the tolerance
1713 // storing it in case it is good.
1714
1715 slentol = sr2 ;
1716 sidetol = kRMax ;
1717 }
1718 if( (ri<0) || (sr2 < halfRadTolerance) )
1719 {
1720 if (b>0) { sr3 = c/(-b-std::sqrt(d)); }
1721 else { sr3 = -b + std::sqrt(d) ; }
1722
1723 // Safety: if both roots -ve ensure that srd cannot `win'
1724 // distancetoout
1725
1726 if ( sr3 > halfRadTolerance )
1727 {
1728 if( sr3 < srd )
1729 {
1730 zi = z + sr3*v.z() ;
1731 ri = tanRMin*zi + rMinAv ;
1732
1733 if ( ri >= 0.0 )
1734 {
1735 srd=sr3 ;
1736 sider=kRMin ;
1737 }
1738 }
1739 }
1740 else if ( sr3 > -halfRadTolerance )
1741 {
1742 // Intersection in tolerance. Store to check if it's good
1743
1744 slentol = sr3 ;
1745 sidetol = kRMin ;
1746 }
1747 }
1748 else if ( (sr2 < srd) && (sr2 > halfCarTolerance) )
1749 {
1750 srd = sr2 ;
1751 sider = kRMin ;
1752 }
1753 else if (sr2 > -halfCarTolerance)
1754 {
1755 // Intersection in tolerance. Store to check if it's good
1756
1757 slentol = sr2 ;
1758 sidetol = kRMin ;
1759 }
1760 if( slentol <= halfCarTolerance )
1761 {
1762 // An intersection within the tolerance was found.
1763 // We must accept it only if the momentum points outwards.
1764
1765 G4ThreeVector Normal ;
1766
1767 // Calculate a normal vector, as below
1768
1769 xi = p.x() + slentol*v.x() ;
1770 yi = p.y() + slentol*v.y() ;
1771 if( sidetol==kRMax )
1772 {
1773 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
1774 Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
1775 }
1776 else
1777 {
1778 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1779 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1780 }
1781 if( Normal.dot(v) > 0 )
1782 {
1783 // We will leave the cone immediately
1784
1785 if( calcNorm )
1786 {
1787 *n = Normal.unit() ;
1788 *validNorm = true ;
1789 }
1790 return snxt = 0.0 ;
1791 }
1792 else
1793 {
1794 // On the surface, but not heading out so we ignore this
1795 // intersection (as it is within tolerance).
1796
1797 slentol = kInfinity ;
1798 }
1799 }
1800 }
1801 }
1802 }
1803 }
1804
1805 // Linear case => point outside inner cone ---> outer cone intersect
1806 //
1807 // Phi Intersection
1808/*
1809 if ( !fPhiFullCone )
1810 {
1811 // add angle calculation with correction
1812 // of the difference in domain of atan2 and Sphi
1813
1814 vphi = std::atan2(v.y(),v.x()) ;
1815
1816 if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; }
1817 else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; }
1818
1819 if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
1820 {
1821 // pDist -ve when inside
1822
1823 pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
1824 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1825
1826 // Comp -ve when in direction of outwards normal
1827
1828 compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
1829 compE = sinEPhi*v.x() - cosEPhi*v.y() ;
1830
1831 sidephi = kNull ;
1832
1833 if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
1834 && (pDistE <= halfCarTolerance) ) )
1835 || ( (fDPhi > pi) && !((pDistS > halfCarTolerance)
1836 && (pDistE > halfCarTolerance) ) ) )
1837 {
1838 // Inside both phi *full* planes
1839 if ( compS < 0 )
1840 {
1841 sphi = pDistS/compS ;
1842 if (sphi >= -halfCarTolerance)
1843 {
1844 xi = p.x() + sphi*v.x() ;
1845 yi = p.y() + sphi*v.y() ;
1846
1847 // Check intersecting with correct half-plane
1848 // (if not -> no intersect)
1849 //
1850 if ( (std::fabs(xi)<=kCarTolerance)
1851 && (std::fabs(yi)<=kCarTolerance) )
1852 {
1853 sidephi= kSPhi;
1854 if ( ( fSPhi-halfAngTolerance <= vphi )
1855 && ( fSPhi+fDPhi+halfAngTolerance >=vphi ) )
1856 {
1857 sphi = kInfinity;
1858 }
1859 }
1860 else
1861 if ( (yi*cosCPhi-xi*sinCPhi)>=0 )
1862 {
1863 sphi = kInfinity ;
1864 }
1865 else
1866 {
1867 sidephi = kSPhi ;
1868 if ( pDistS > -halfCarTolerance )
1869 {
1870 sphi = 0.0 ; // Leave by sphi immediately
1871 }
1872 }
1873 }
1874 else
1875 {
1876 sphi = kInfinity ;
1877 }
1878 }
1879 else
1880 {
1881 sphi = kInfinity ;
1882 }
1883
1884 if ( compE < 0 )
1885 {
1886 sphi2 = pDistE/compE ;
1887
1888 // Only check further if < starting phi intersection
1889 //
1890 if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
1891 {
1892 xi = p.x() + sphi2*v.x() ;
1893 yi = p.y() + sphi2*v.y() ;
1894
1895 // Check intersecting with correct half-plane
1896
1897 if ( (std::fabs(xi)<=kCarTolerance)
1898 && (std::fabs(yi)<=kCarTolerance) )
1899 {
1900 // Leaving via ending phi
1901
1902 if(!( (fSPhi-halfAngTolerance <= vphi)
1903 && (fSPhi+fDPhi+halfAngTolerance >= vphi) ) )
1904 {
1905 sidephi = kEPhi ;
1906 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; }
1907 else { sphi = 0.0; }
1908 }
1909 }
1910 else // Check intersecting with correct half-plane
1911 if ( yi*cosCPhi-xi*sinCPhi >= 0 )
1912 {
1913 // Leaving via ending phi
1914
1915 sidephi = kEPhi ;
1916 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; }
1917 else { sphi = 0.0; }
1918 }
1919 }
1920 }
1921 }
1922 else
1923 {
1924 sphi = kInfinity ;
1925 }
1926 }
1927 else
1928 {
1929 // On z axis + travel not || to z axis -> if phi of vector direction
1930 // within phi of shape, Step limited by rmax, else Step =0
1931
1932 if ( (fSPhi-halfAngTolerance <= vphi)
1933 && (vphi <= fSPhi+fDPhi+halfAngTolerance) )
1934 {
1935 sphi = kInfinity ;
1936 }
1937 else
1938 {
1939 sidephi = kSPhi ; // arbitrary
1940 sphi = 0.0 ;
1941 }
1942 }
1943 if ( sphi < snxt ) // Order intersecttions
1944 {
1945 snxt=sphi ;
1946 side=sidephi ;
1947 }
1948 }
1949*/
1950 if ( srd < snxt ) // Order intersections
1951 {
1952 snxt = srd ;
1953 side = sider ;
1954 }
1955 if (calcNorm)
1956 {
1957 switch(side)
1958 { // Note: returned vector not normalised
1959 case kRMax: // (divide by frmax for unit vector)
1960 xi = p.x() + snxt*v.x() ;
1961 yi = p.y() + snxt*v.y() ;
1962 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
1963 *n = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
1964 *validNorm = true ;
1965 break ;
1966 case kRMin:
1967 *validNorm = false ; // Rmin is inconvex
1968 break ;
1969/* case kSPhi:
1970 if ( fDPhi <= pi )
1971 {
1972 *n = G4ThreeVector(sinSPhi, -cosSPhi, 0);
1973 *validNorm = true ;
1974 }
1975 else
1976 {
1977 *validNorm = false ;
1978 }
1979 break ;
1980 case kEPhi:
1981 if ( fDPhi <= pi )
1982 {
1983 *n = G4ThreeVector(-sinEPhi, cosEPhi, 0);
1984 *validNorm = true ;
1985 }
1986 else
1987 {
1988 *validNorm = false ;
1989 }
1990 break ;*/
1991 case kPZ:
1992 *n = G4ThreeVector(0,0,1) ;
1993 *validNorm = true ;
1994 break ;
1995 case kMZ:
1996 *n = G4ThreeVector(0,0,-1) ;
1997 *validNorm = true ;
1998 break ;
1999 default:
2000 G4cout << G4endl ;
2001 DumpInfo();
2002 std::ostringstream message;
2003 G4int oldprc = message.precision(16) ;
2004 message << "Undefined side for valid surface normal to solid."
2005 << G4endl
2006 << "Position:" << G4endl << G4endl
2007 << "p.x() = " << p.x()/mm << " mm" << G4endl
2008 << "p.y() = " << p.y()/mm << " mm" << G4endl
2009 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
2010 << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
2011 << " mm" << G4endl << G4endl ;
2012 if( p.x() != 0. || p.y() != 0.)
2013 {
2014 message << "point phi = " << std::atan2(p.y(),p.x())/degree
2015 << " degree" << G4endl << G4endl ;
2016 }
2017 message << "Direction:" << G4endl << G4endl
2018 << "v.x() = " << v.x() << G4endl
2019 << "v.y() = " << v.y() << G4endl
2020 << "v.z() = " << v.z() << G4endl<< G4endl
2021 << "Proposed distance :" << G4endl<< G4endl
2022 << "snxt = " << snxt/mm << " mm" << G4endl ;
2023 message.precision(oldprc) ;
2024 G4Exception("G4ShiftedCone::DistanceToOut(p,v,..)","GeomSolids1002",
2025 JustWarning, message) ;
2026 break ;
2027 }
2028 }
2029 if (snxt < halfCarTolerance) { snxt = 0.; }
2030
2031 return snxt ;
2032}
2033
2035//
2036// Calculate distance (<=actual) to closest surface of shape from inside
2037
2038G4double G4ShiftedCone::DistanceToOut(const G4ThreeVector& p) const
2039{
2040 G4double safe=0.0, rho, safeR1, safeR2, safeZ;//, safePhi;
2041 G4double tanRMin, secRMin, pRMin;
2042 G4double tanRMax, secRMax, pRMax;
2043
2044#ifdef G4CSGDEBUG
2045 if( Inside(p) == kOutside )
2046 {
2047 G4int oldprc=G4cout.precision(16) ;
2048 G4cout << G4endl ;
2049 DumpInfo();
2050 G4cout << "Position:" << G4endl << G4endl ;
2051 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
2052 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
2053 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
2054 G4cout << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
2055 << " mm" << G4endl << G4endl ;
2056 if( (p.x() != 0.) || (p.x() != 0.) )
2057 {
2058 G4cout << "point phi = " << std::atan2(p.y(),p.x())/degree
2059 << " degree" << G4endl << G4endl ;
2060 }
2061 G4cout.precision(oldprc) ;
2062 G4Exception("G4ShiftedCone::DistanceToOut(p)", "GeomSolids1002",
2063 JustWarning, "Point p is outside !?" );
2064 }
2065#endif
2066
2067 G4double z = p.z() - fZshift;
2068
2069 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
2070 safeZ = fDz - std::fabs(z) ;
2071
2072 if (fRmin1 || fRmin2)
2073 {
2074 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
2075 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
2076 pRMin = tanRMin*z + (fRmin1 + fRmin2)*0.5 ;
2077 safeR1 = (rho - pRMin)/secRMin ;
2078 }
2079 else
2080 {
2081 safeR1 = kInfinity ;
2082 }
2083
2084 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
2085 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
2086 pRMax = tanRMax*z + (fRmax1+fRmax2)*0.5 ;
2087 safeR2 = (pRMax - rho)/secRMax ;
2088
2089 if (safeR1 < safeR2) { safe = safeR1; }
2090 else { safe = safeR2; }
2091 if (safeZ < safe) { safe = safeZ ; }
2092
2093 // Check if phi divided, Calc distances closest phi plane
2094/*
2095 if (!fPhiFullCone)
2096 {
2097 // Above/below central phi of G4ShiftedCone?
2098
2099 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
2100 {
2101 safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
2102 }
2103 else
2104 {
2105 safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
2106 }
2107 if (safePhi < safe) { safe = safePhi; }
2108 }*/
2109 if ( safe < 0 ) { safe = 0; }
2110
2111 return safe ;
2112}
2113
2115//
2116// GetEntityType
2117
2118G4GeometryType G4ShiftedCone::GetEntityType() const
2119{
2120 return G4String("G4ShiftedCone");
2121}
2122
2124//
2125// Make a clone of the object
2126//
2127G4VSolid* G4ShiftedCone::Clone() const
2128{
2129 return new G4ShiftedCone(*this);
2130}
2131
2133//
2134// Stream object contents to an output stream
2135
2136std::ostream& G4ShiftedCone::StreamInfo(std::ostream& os) const
2137{
2138 G4int oldprc = os.precision(16);
2139 os << "-----------------------------------------------------------\n"
2140 << " *** Dump for solid - " << GetName() << " ***\n"
2141 << " ===================================================\n"
2142 << " Solid type: G4ShiftedCone\n"
2143 << " Parameters: \n"
2144 << " inside -fDz radius: " << fRmin1/mm << " mm \n"
2145 << " outside -fDz radius: " << fRmax1/mm << " mm \n"
2146 << " inside +fDz radius: " << fRmin2/mm << " mm \n"
2147 << " outside +fDz radius: " << fRmax2/mm << " mm \n"
2148 << " Z1 : " << GetZ1()/mm << " mm \n"
2149 << " Z2 : " << GetZ2()/mm << " mm \n"
2150// << " starting angle of segment: " << fSPhi/degree << " degrees \n"
2151// << " delta angle of segment : " << fDPhi/degree << " degrees \n"
2152 << "-----------------------------------------------------------\n";
2153 os.precision(oldprc);
2154
2155 return os;
2156}
2157
2158
2159
2161//
2162// GetPointOnSurface
2163
2164G4ThreeVector G4ShiftedCone::GetPointOnSurface() const
2165{
2166 // declare working variables
2167 //
2168 G4double rone = (fRmax1-fRmax2)/(2.*fDz);
2169 G4double rtwo = (fRmin1-fRmin2)/(2.*fDz);
2170 G4double qone = (fRmax1 == fRmax2) ? 0. : fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2);
2171 G4double qtwo = (fRmin1 == fRmin2) ? 0. : fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2);
2172
2173 G4double slin = std::hypot(fRmin1-fRmin2, 2.*fDz);
2174 G4double slout = std::hypot(fRmax1-fRmax2, 2.*fDz);
2175 G4double Aone = 0.5*GetDeltaPhiAngle()*(fRmax2 + fRmax1)*slout; // outer surface
2176 G4double Atwo = 0.5*GetDeltaPhiAngle()*(fRmin2 + fRmin1)*slin; // inner surface
2177 G4double Athree = 0.5*GetDeltaPhiAngle()*(fRmax1*fRmax1-fRmin1*fRmin1); // base at -Dz
2178 G4double Afour = 0.5*GetDeltaPhiAngle()*(fRmax2*fRmax2-fRmin2*fRmin2); // base at +Dz
2179 G4double Afive = fDz*(fRmax1-fRmin1+fRmax2-fRmin2); // phi section
2180
2181 G4double phi = G4RandFlat::shoot(GetStartPhiAngle(),GetStartPhiAngle() + GetDeltaPhiAngle());
2182 G4double cosu = std::cos(phi);
2183 G4double sinu = std::sin(phi);
2184 G4double rRand1 = GetRadiusInRing(fRmin1, fRmax1);
2185 G4double rRand2 = GetRadiusInRing(fRmin2, fRmax2);
2186
2187 G4bool fPhiFullCone = true;
2188 if ( (GetStartPhiAngle() == 0.) && fPhiFullCone ) { Afive = 0.; }
2189 G4double chose = G4RandFlat::shoot(0.,Aone+Atwo+Athree+Afour+2.*Afive);
2190
2191 if( (chose >= 0.) && (chose < Aone) ) // outer surface
2192 {
2193 if(fRmax1 != fRmax2)
2194 {
2195 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2196 return G4ThreeVector (rone*cosu*(qone-zRand),
2197 rone*sinu*(qone-zRand), zRand + fZshift);
2198 }
2199 else
2200 {
2201 return G4ThreeVector(fRmax1*cosu, fRmax2*sinu,
2202 G4RandFlat::shoot(-1.*fDz,fDz) + fZshift);
2203 }
2204 }
2205 else if( (chose >= Aone) && (chose < Aone + Atwo) ) // inner surface
2206 {
2207 if(fRmin1 != fRmin2)
2208 {
2209 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2210 return G4ThreeVector (rtwo*cosu*(qtwo-zRand),
2211 rtwo*sinu*(qtwo-zRand), zRand + fZshift);
2212 }
2213 else
2214 {
2215 return G4ThreeVector(fRmin1*cosu, fRmin2*sinu,
2216 G4RandFlat::shoot(-1.*fDz,fDz) + fZshift);
2217 }
2218 }
2219 else if( (chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree) ) // base at -Dz
2220 {
2221 return G4ThreeVector (rRand1*cosu, rRand1*sinu, -1*fDz + fZshift);
2222 }
2223 else if( (chose >= Aone + Atwo + Athree)
2224 && (chose < Aone + Atwo + Athree + Afour) ) // base at +Dz
2225 {
2226 return G4ThreeVector (rRand2*cosu,rRand2*sinu,fDz + fZshift);
2227 }
2228 else if( (chose >= Aone + Atwo + Athree + Afour) // SPhi section
2229 && (chose < Aone + Atwo + Athree + Afour + Afive) )
2230 {
2231 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2232 rRand1 = G4RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
2233 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
2234 return G4ThreeVector (rRand1*GetCosStartPhi(),
2235 rRand1*GetSinStartPhi(), zRand + fZshift);
2236 }
2237 else // SPhi+DPhi section
2238 {
2239 G4double zRand = G4RandFlat::shoot(-1.*fDz,fDz);
2240 rRand1 = G4RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
2241 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
2242 return G4ThreeVector (rRand1*GetCosEndPhi(),
2243 rRand1*GetSinEndPhi(), zRand + fZshift);
2244 }
2245}
2246
2248//
2249// Methods for visualisation
2250
2251void G4ShiftedCone::DescribeYourselfTo (G4VGraphicsScene& scene) const
2252{
2253 scene.AddSolid (*this);
2254}
2255
2256G4Polyhedron* G4ShiftedCone::CreatePolyhedron () const
2257{
2258 G4double rmin[2] = { GetRmin1(), GetRmin2() };
2259 G4double rmax[2] = { GetRmax1(), GetRmax2() };
2260 G4double z[2] = { GetZ1(), GetZ2() };
2261 return new G4PolyhedronPcon(
2262 GetStartPhiAngle(), GetDeltaPhiAngle(), 2, z, rmin, rmax
2263 );
2264}
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:18
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:1412
G4ShiftedCone::halfRadTolerance
G4double halfRadTolerance
Definition G4ShiftedCone.h:240
G4ShiftedCone::DistanceToIn
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition G4ShiftedCone.cxx:679
G4ShiftedCone::GetRmin2
G4double GetRmin2() const
G4ShiftedCone::GetOuterRadiusMinusZ
G4double GetOuterRadiusMinusZ() const
G4ShiftedCone::DescribeYourselfTo
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition G4ShiftedCone.cxx:2251
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:2256
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:2118
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:2164
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:2127
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:2136
Generated on for ATLAS Offline Software by doxygen 1.16.1