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