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