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