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