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