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