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