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