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