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