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