ATLAS Offline Software
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 
64 using namespace CLHEP;
65 
67 //
68 // Private enum: Not for external use
69 namespace {
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 
82 G4ShiftedCone::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 //
134 G4ShiftedCone::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 
149 G4ShiftedCone::~G4ShiftedCone()
150 {
151 }
152 
154 //
155 // Copy constructor
156 
157 G4ShiftedCone::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 
176 G4ShiftedCone& 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 
210 EInside 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 
268 void 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 
283 void 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 
326 G4bool 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 
443 G4ThreeVector 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 
549 G4ThreeVector 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 
678 G4double 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 
1351 G4double 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 
1415 G4double 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 
2041 G4double 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 
2121 G4GeometryType G4ShiftedCone::GetEntityType() const
2122 {
2123  return G4String("G4ShiftedCone");
2124 }
2125 
2127 //
2128 // Make a clone of the object
2129 //
2130 G4VSolid* G4ShiftedCone::Clone() const
2131 {
2132  return new G4ShiftedCone(*this);
2133 }
2134 
2136 //
2137 // Stream object contents to an output stream
2138 
2139 std::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 
2167 G4ThreeVector 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 
2254 void G4ShiftedCone::DescribeYourselfTo (G4VGraphicsScene& scene) const
2255 {
2256  scene.AddSolid (*this);
2257 }
2258 
2259 G4Polyhedron* 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 }
G4ShiftedCone::fDz
G4double fDz
Definition: G4ShiftedCone.h:227
G4ShiftedCone::DescribeYourselfTo
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition: G4ShiftedCone.cxx:2254
G4ShiftedCone::fZshift
G4double fZshift
Definition: G4ShiftedCone.h:227
G4ShiftedCone::~G4ShiftedCone
~G4ShiftedCone()
Definition: G4ShiftedCone.cxx:149
PlotCalibFromCool.norm
norm
Definition: PlotCalibFromCool.py:100
G4ShiftedCone::GetRmin1
G4double GetRmin1() const
G4ShiftedCone::GetCosStartPhi
G4double GetCosStartPhi() const
G4ShiftedCone::Clone
G4VSolid * Clone() const
Definition: G4ShiftedCone.cxx:2130
python.SystemOfUnits.mm
float mm
Definition: SystemOfUnits.py:98
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:67
G4ShiftedCone::operator=
G4ShiftedCone & operator=(const G4ShiftedCone &rhs)
Definition: G4ShiftedCone.cxx:176
G4ShiftedCone::kNRMax
@ kNRMax
Definition: G4ShiftedCone.h:220
G4ShiftedCone::GetInnerRadiusPlusZ
G4double GetInnerRadiusPlusZ() const
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
hist_file_dump.d
d
Definition: hist_file_dump.py:142
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
LUCID_EventTPCnv_Dict::t3
std::vector< LUCID_RawData_p1 > t3
Definition: LUCID_EventTPCnvDict.h:28
G4ShiftedCone::GetRmax2
G4double GetRmax2() const
ORAlgo::Normal
@ Normal
ALFA_EventTPCnv_Dict::t1
std::vector< ALFA_RawDataCollection_p1 > t1
Definition: ALFA_EventTPCnvDict.h:43
G4ShiftedCone::GetEntityType
G4GeometryType GetEntityType() const
Definition: G4ShiftedCone.cxx:2121
deg
#define deg
Definition: SbPolyhedron.cxx:17
G4ShiftedCone::GetZ1
G4double GetZ1() const
G4ShiftedCone::kMZ
@ kMZ
Definition: G4ShiftedCone.h:216
G4ShiftedCone::Inside
EInside Inside(const G4ThreeVector &p) const
Definition: G4ShiftedCone.cxx:210
G4ShiftedCone::CalculateExtent
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
Definition: G4ShiftedCone.cxx:326
columnar::operator=
AccessorTemplate & operator=(AccessorTemplate &&that)
Definition: VectorColumn.h:88
G4ShiftedCone::SurfaceNormal
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4ShiftedCone.cxx:443
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
G4ShiftedCone::CreatePolyhedron
G4Polyhedron * CreatePolyhedron() const
Definition: G4ShiftedCone.cxx:2259
G4ShiftedCone::GetInnerRadiusMinusZ
G4double GetInnerRadiusMinusZ() const
G4ShiftedCone
Definition: G4ShiftedCone.h:87
G4ShiftedCone::GetSinEndPhi
G4double GetSinEndPhi() const
ReweightUtils.message
message
Definition: ReweightUtils.py:15
G4ShiftedCone::StreamInfo
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4ShiftedCone.cxx:2139
G4ShiftedCone::ESide
ESide
Definition: G4ShiftedCone.h:216
G4ShiftedCone::halfCarTolerance
G4double halfCarTolerance
Definition: G4ShiftedCone.h:240
G4ShiftedCone::GetRmax1
G4double GetRmax1() const
G4ShiftedCone::kRMin
@ kRMin
Definition: G4ShiftedCone.h:216
TRT::Hit::side
@ side
Definition: HitInfo.h:83
python.selector.AtlRunQuerySelectorLhcOlc.sd
sd
Definition: AtlRunQuerySelectorLhcOlc.py:611
xAOD::TauJetParameters::dRmax
@ dRmax
Get maximal dR of tracks associated to calo-seeded tau.
Definition: TauDefs.h:226
G4ShiftedCone::DistanceToIn
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4ShiftedCone.cxx:678
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
G4ShiftedCone::ComputeDimensions
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4ShiftedCone.cxx:268
WriteCellNoiseToCool.exist
exist
Definition: WriteCellNoiseToCool.py:543
CheckAppliedSFs.e3
e3
Definition: CheckAppliedSFs.py:264
G4ShiftedCone::GetOuterRadiusMinusZ
G4double GetOuterRadiusMinusZ() const
G4ShiftedCone::GetSinStartPhi
G4double GetSinStartPhi() const
G4ShiftedCone::halfAngTolerance
G4double halfAngTolerance
Definition: G4ShiftedCone.h:240
z
#define z
G4ShiftedCone::fRmax1
G4double fRmax1
Definition: G4ShiftedCone.h:226
CheckAppliedSFs.bmin
bmin
Definition: CheckAppliedSFs.py:242
beamspotman.n
n
Definition: beamspotman.py:727
G4ShiftedCone::GetCosEndPhi
G4double GetCosEndPhi() const
CLHEP
STD'S.
Definition: CaloNoiseCompCondAlg.h:58
G4ShiftedCone::fRmin2
G4double fRmin2
Definition: G4ShiftedCone.h:226
G4ShiftedCone::fRmin1
G4double fRmin1
Definition: G4ShiftedCone.h:226
G4ShiftedCone::G4ShiftedCone
G4ShiftedCone(const G4String &pName, G4double pZ1, G4double pZ2, G4double pRmin1, G4double pRmax1, G4double pRmin2, G4double pRmax2)
Definition: G4ShiftedCone.cxx:82
ReadFromCoolCompare.os
os
Definition: ReadFromCoolCompare.py:231
G4ShiftedCone.h
twopi
constexpr double twopi
Definition: VertexPointEstimator.cxx:16
CheckAppliedSFs.bmax
bmax
Definition: CheckAppliedSFs.py:264
G4ShiftedCone::BoundingLimits
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
Definition: G4ShiftedCone.cxx:283
G4ShiftedCone::kRadTolerance
G4double kRadTolerance
Definition: G4ShiftedCone.h:222
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:76
G4ShiftedCone::GetZ2
G4double GetZ2() const
G4ShiftedCone::kNRMin
@ kNRMin
Definition: G4ShiftedCone.h:220
G4ShiftedCone::ApproxSurfaceNormal
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
Definition: G4ShiftedCone.cxx:549
G4ShiftedCone::kNull
@ kNull
Definition: G4ShiftedCone.h:216
PlotCalibFromCool.rl
rl
Definition: PlotCalibFromCool.py:529
G4ShiftedCone::GetRmin2
G4double GetRmin2() const
G4ShiftedCone::GetPointOnSurface
G4ThreeVector GetPointOnSurface() const
Definition: G4ShiftedCone.cxx:2167
python.PyAthena.v
v
Definition: PyAthena.py:154
G4ShiftedCone::kPZ
@ kPZ
Definition: G4ShiftedCone.h:216
ALFA_EventTPCnv_Dict::t2
std::vector< ALFA_RawDataContainer_p1 > t2
Definition: ALFA_EventTPCnvDict.h:44
a
TList * a
Definition: liststreamerinfos.cxx:10
G4ShiftedCone::ENorm
ENorm
Definition: G4ShiftedCone.h:220
G4ShiftedCone::kRMax
@ kRMax
Definition: G4ShiftedCone.h:216
G4ShiftedCone::kAngTolerance
G4double kAngTolerance
Definition: G4ShiftedCone.h:222
G4ShiftedCone::kNZ
@ kNZ
Definition: G4ShiftedCone.h:220
G4ShiftedCone::GetDeltaPhiAngle
G4double GetDeltaPhiAngle() const
G4ShiftedCone::GetStartPhiAngle
G4double GetStartPhiAngle() const
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
G4ShiftedCone::fRmax2
G4double fRmax2
Definition: G4ShiftedCone.h:226
python.compressB64.c
def c
Definition: compressB64.py:93
python.SystemOfUnits.degree
tuple degree
Definition: SystemOfUnits.py:121
G4ShiftedCone::halfRadTolerance
G4double halfRadTolerance
Definition: G4ShiftedCone.h:240
fitman.rho
rho
Definition: fitman.py:532
fitman.k
k
Definition: fitman.py:528
G4ShiftedCone::GetOuterRadiusPlusZ
G4double GetOuterRadiusPlusZ() const
Generated on Thu Oct 30 2025 21:10:27 for ATLAS Offline Software by   doxygen 1.8.18