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