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