| 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: G4Tubs.cc,v 1.11.2.1 1999/12/07 20:48:33 gunter Exp $
7 // * conditions of the Geant4 Software License << 9 // GEANT4 tag $Name: geant4-01-00 $
8 // * LICENSE and available at http://cern.ch/ << 10 //
9 // * include a list of copyright holders. << 11 //
10 // * << 12 // class G4Tubs
11 // * Neither the authors of this software syst << 13 //
12 // * institutes,nor the agencies providing fin << 14 // History:
13 // * work make any representation or warran << 15 //
14 // * regarding this software system or assum << 16 // 1994-95 P.Kent, implementation
15 // * use. Please see the license in the file << 17 //
16 // * for the full disclaimer and the limitatio << 18 // 18.06.98 V.Grichine, n-normalisation in DistanceToOut(p.v)
17 // * << 19 // 09.10.98 V.Grichine, modifications in Distance ToOut(p,v,...)
18 // * This code implementation is the result << 20 // 23.03.99 V.Grichine, bug fixed in DistanceToIn(p,v)
19 // * technical work of the GEANT4 collaboratio << 21 // 25.05.99 V.Grichine, bugs fixed in DistanceToIn(p,v)
20 // * By using, copying, modifying or distri << 22 // 28.05.99 V.Grichine, bugs fixed in Distance ToOut(p,v,...)
21 // * any work based on the software) you ag << 23 // 13.10.99 V.Grichine, bugs fixed in DistanceToIn(p,v)
22 // * use in resulting scientific publicati << 24 // 19.11.99 V. Grichine, side = kNull in Distance ToOut(p,v,...)
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* <<
25 // <<
26 // G4Tubs implementation <<
27 // <<
28 // 1994-95 P.Kent: first implementation <<
29 // 08.08.00 V.Grichine: more stable roots of 2 <<
30 // 07.12.00 V.Grichine: phi-section algorithm <<
31 // 03.05.05 V.Grichine: SurfaceNormal(p) accor <<
32 // 24.08.16 E.Tcherniaev: reimplemented Calcul <<
33 // ------------------------------------------- <<
34 25
35 #include "G4Tubs.hh" <<
36 26
37 #if !defined(G4GEOM_USE_UTUBS) << 27 #include "G4Tubs.hh"
38 28
39 #include "G4GeomTools.hh" <<
40 #include "G4VoxelLimits.hh" 29 #include "G4VoxelLimits.hh"
41 #include "G4AffineTransform.hh" 30 #include "G4AffineTransform.hh"
42 #include "G4GeometryTolerance.hh" <<
43 #include "G4BoundingEnvelope.hh" <<
44 31
45 #include "G4VPVParameterisation.hh" 32 #include "G4VPVParameterisation.hh"
46 #include "G4QuickRand.hh" << 33
>> 34 #include "meshdefs.hh"
47 35
48 #include "G4VGraphicsScene.hh" 36 #include "G4VGraphicsScene.hh"
49 #include "G4Polyhedron.hh" 37 #include "G4Polyhedron.hh"
50 << 38 #include "G4NURBS.hh"
51 using namespace CLHEP; << 39 #include "G4NURBStube.hh"
>> 40 #include "G4NURBScylinder.hh"
>> 41 #include "G4NURBStubesector.hh"
>> 42 #include "G4VisExtent.hh"
52 43
53 ////////////////////////////////////////////// 44 /////////////////////////////////////////////////////////////////////////
54 // 45 //
55 // Constructor - check parameters, convert ang 46 // Constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
56 // - note if pdphi>2PI then reset 47 // - note if pdphi>2PI then reset to 2PI
57 48
58 G4Tubs::G4Tubs( const G4String& pName, << 49 G4Tubs::G4Tubs(const G4String &pName,
59 G4double pRMin, G4double << 50 G4double pRMin,
60 G4double pDz, << 51 G4double pRMax,
61 G4double pSPhi, G4double << 52 G4double pDz,
62 : G4CSGSolid(pName), fRMin(pRMin), fRMax(pR << 53 G4double pSPhi,
63 fSPhi(0), fDPhi(0), << 54 G4double pDPhi)
64 fInvRmax( pRMax > 0.0 ? 1.0/pRMax : 0.0 ) << 55 : G4CSGSolid(pName)
65 fInvRmin( pRMin > 0.0 ? 1.0/pRMin : 0.0 ) <<
66 { 56 {
67 kRadTolerance = G4GeometryTolerance::GetInst << 57 // Check z-len
68 kAngTolerance = G4GeometryTolerance::GetInst << 58 if (pDz>0)
>> 59 {
>> 60 fDz=pDz;
>> 61 }
>> 62 else
>> 63 {
>> 64 G4Exception("Error in G4Tubs::G4Tubs - invalid z half-length");
>> 65 }
69 66
70 halfCarTolerance=kCarTolerance*0.5; << 67 // Check radii
71 halfRadTolerance=kRadTolerance*0.5; <<
72 halfAngTolerance=kAngTolerance*0.5; <<
73 68
74 if (pDz<=0) // Check z-len << 69 if (pRMin<pRMax&&pRMin>=0)
75 { << 70 {
76 std::ostringstream message; << 71 fRMin=pRMin; fRMax=pRMax;
77 message << "Negative Z half-length (" << p << 72 }
78 G4Exception("G4Tubs::G4Tubs()", "GeomSolid << 73 else
79 } << 74 {
80 if ( (pRMin >= pRMax) || (pRMin < 0) ) // Ch << 75 G4Exception("Error in G4Tubs::G4Tubs - invalid radii");
81 { << 76 }
82 std::ostringstream message; <<
83 message << "Invalid values for radii in so <<
84 << G4endl <<
85 << " pRMin = " << pRMin << <<
86 G4Exception("G4Tubs::G4Tubs()", "GeomSolid <<
87 } <<
88 <<
89 // Check angles <<
90 // <<
91 CheckPhiAngles(pSPhi, pDPhi); <<
92 } <<
93 77
94 ////////////////////////////////////////////// << 78 // Check angles
95 // <<
96 // Fake default constructor - sets only member <<
97 // for usage restri <<
98 // <<
99 G4Tubs::G4Tubs( __void__& a ) <<
100 : G4CSGSolid(a) <<
101 { <<
102 } <<
103 79
104 ////////////////////////////////////////////// << 80 if (pDPhi>=2.0*M_PI)
105 // << 81 {
106 // Destructor << 82 fDPhi=2*M_PI;
>> 83 }
>> 84 else
>> 85 {
>> 86 if (pDPhi>0)
>> 87 {
>> 88 fDPhi = pDPhi;
>> 89 }
>> 90 else
>> 91 {
>> 92 G4Exception("Error in G4Tubs::G4Tubs - invalid dphi");
>> 93 }
>> 94 }
>> 95
>> 96 // Ensure psphi in 0-2PI or -2PI-0 range if shape crosses 0
107 97
108 G4Tubs::~G4Tubs() = default; << 98 fSPhi = pSPhi;
109 99
110 ////////////////////////////////////////////// << 100 if ( fSPhi < 0 )
111 // << 101 {
112 // Copy constructor << 102 fSPhi=2.0*M_PI-fmod(fabs(fSPhi),2.0*M_PI);
>> 103 }
>> 104 else
>> 105 {
>> 106 fSPhi=fmod(fSPhi,2.0*M_PI);
>> 107 }
113 108
114 G4Tubs::G4Tubs(const G4Tubs&) = default; << 109 if (fSPhi+fDPhi>2.0*M_PI)
>> 110 {
>> 111 fSPhi -= 2.0*M_PI ;
>> 112 }
>> 113 }
115 114
116 ////////////////////////////////////////////// 115 //////////////////////////////////////////////////////////////////////////
117 // 116 //
118 // Assignment operator << 117 // Destructor
119 <<
120 G4Tubs& G4Tubs::operator = (const G4Tubs& rhs) <<
121 { <<
122 // Check assignment to self <<
123 // <<
124 if (this == &rhs) { return *this; } <<
125 <<
126 // Copy base class data <<
127 // <<
128 G4CSGSolid::operator=(rhs); <<
129 <<
130 // Copy data <<
131 // <<
132 kRadTolerance = rhs.kRadTolerance; kAngTole <<
133 fRMin = rhs.fRMin; fRMax = rhs.fRMax; fDz = <<
134 fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi; <<
135 sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPh <<
136 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = r <<
137 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPh <<
138 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPh <<
139 fPhiFullTube = rhs.fPhiFullTube; <<
140 fInvRmax = rhs.fInvRmax; <<
141 fInvRmin = rhs.fInvRmin; <<
142 halfCarTolerance = rhs.halfCarTolerance; <<
143 halfRadTolerance = rhs.halfRadTolerance; <<
144 halfAngTolerance = rhs.halfAngTolerance; <<
145 118
146 return *this; << 119 G4Tubs::~G4Tubs()
147 } << 120 {;}
148 121
149 ////////////////////////////////////////////// 122 /////////////////////////////////////////////////////////////////////////
150 // 123 //
151 // Dispatch to parameterisation for replicatio 124 // Dispatch to parameterisation for replication mechanism dimension
152 // computation & modification. 125 // computation & modification.
153 126
154 void G4Tubs::ComputeDimensions( G4VPVPar << 127 void G4Tubs::ComputeDimensions(G4VPVParameterisation* p,
155 const G4int n, << 128 const G4int n,
156 const G4VPhysi << 129 const G4VPhysicalVolume* pRep)
157 { 130 {
158 p->ComputeDimensions(*this,n,pRep) ; << 131 p->ComputeDimensions(*this,n,pRep);
159 } 132 }
160 133
161 ////////////////////////////////////////////// << 134 ////////////////////////////////////////////////////////////////////////
162 // <<
163 // Get bounding box <<
164 <<
165 void G4Tubs::BoundingLimits(G4ThreeVector& pMi <<
166 { <<
167 G4double rmin = GetInnerRadius(); <<
168 G4double rmax = GetOuterRadius(); <<
169 G4double dz = GetZHalfLength(); <<
170 <<
171 // Find bounding box <<
172 // <<
173 if (GetDeltaPhiAngle() < twopi) <<
174 { <<
175 G4TwoVector vmin,vmax; <<
176 G4GeomTools::DiskExtent(rmin,rmax, <<
177 GetSinStartPhi(),G <<
178 GetSinEndPhi(),Get <<
179 vmin,vmax); <<
180 pMin.set(vmin.x(),vmin.y(),-dz); <<
181 pMax.set(vmax.x(),vmax.y(), dz); <<
182 } <<
183 else <<
184 { <<
185 pMin.set(-rmax,-rmax,-dz); <<
186 pMax.set( rmax, rmax, dz); <<
187 } <<
188 <<
189 // Check correctness of the bounding box <<
190 // <<
191 if (pMin.x() >= pMax.x() || pMin.y() >= pMax <<
192 { <<
193 std::ostringstream message; <<
194 message << "Bad bounding box (min >= max) <<
195 << GetName() << " !" <<
196 << "\npMin = " << pMin <<
197 << "\npMax = " << pMax; <<
198 G4Exception("G4Tubs::BoundingLimits()", "G <<
199 JustWarning, message); <<
200 DumpInfo(); <<
201 } <<
202 } <<
203 <<
204 ////////////////////////////////////////////// <<
205 // 135 //
206 // Calculate extent under transform and specif 136 // Calculate extent under transform and specified limit
207 137
208 G4bool G4Tubs::CalculateExtent( const EAxis << 138 G4bool G4Tubs::CalculateExtent(const EAxis pAxis,
209 const G4VoxelL << 139 const G4VoxelLimits& pVoxelLimit,
210 const G4Affine << 140 const G4AffineTransform& pTransform,
211 G4double << 141 G4double& pMin, G4double& pMax) const
212 G4double <<
213 { 142 {
214 G4ThreeVector bmin, bmax; << 143 if (!pTransform.IsRotated()&&fDPhi==2.0*M_PI&&fRMin==0)
215 G4bool exist; << 144 {
216 << 145 // Special case handling for unrotated solid tubes
217 // Get bounding box << 146 // Compute x/y/z mins and maxs fro bounding box respecting limits,
218 BoundingLimits(bmin,bmax); << 147 // with early returns if outside limits. Then switch() on pAxis,
219 << 148 // and compute exact x and y limit for x/y case
220 // Check bounding box << 149
221 G4BoundingEnvelope bbox(bmin,bmax); << 150 G4double xoffset,xMin,xMax;
222 #ifdef G4BBOX_EXTENT << 151 G4double yoffset,yMin,yMax;
223 return bbox.CalculateExtent(pAxis,pVoxelLimi << 152 G4double zoffset,zMin,zMax;
224 #endif << 153
225 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 154 G4double diff1,diff2,maxDiff,newMin,newMax;
226 { << 155 G4double xoff1,xoff2,yoff1,yoff2;
227 return exist = pMin < pMax; << 156
228 } << 157 xoffset=pTransform.NetTranslation().x();
229 << 158 xMin=xoffset-fRMax;
230 // Get parameters of the solid << 159 xMax=xoffset+fRMax;
231 G4double rmin = GetInnerRadius(); << 160 if (pVoxelLimit.IsXLimited())
232 G4double rmax = GetOuterRadius(); << 161 {
233 G4double dz = GetZHalfLength(); << 162 if (xMin>pVoxelLimit.GetMaxXExtent()
234 G4double dphi = GetDeltaPhiAngle(); << 163 ||xMax<pVoxelLimit.GetMinXExtent())
235 << 164 {
236 // Find bounding envelope and calculate exte << 165 return false;
237 // << 166 }
238 const G4int NSTEPS = 24; // numbe << 167 else
239 G4double astep = twopi/NSTEPS; // max a << 168 {
240 G4int ksteps = (dphi <= astep) ? 1 : (G4i << 169 if (xMin<pVoxelLimit.GetMinXExtent())
241 G4double ang = dphi/ksteps; << 170 {
242 << 171 xMin=pVoxelLimit.GetMinXExtent();
243 G4double sinHalf = std::sin(0.5*ang); << 172 }
244 G4double cosHalf = std::cos(0.5*ang); << 173 if (xMax>pVoxelLimit.GetMaxXExtent())
245 G4double sinStep = 2.*sinHalf*cosHalf; << 174 {
246 G4double cosStep = 1. - 2.*sinHalf*sinHalf; << 175 xMax=pVoxelLimit.GetMaxXExtent();
247 G4double rext = rmax/cosHalf; << 176 }
248 << 177 }
249 // bounding envelope for full cylinder consi << 178 }
250 // in other cases it is a sequence of quadri << 179
251 if (rmin == 0 && dphi == twopi) << 180 yoffset=pTransform.NetTranslation().y();
252 { << 181 yMin=yoffset-fRMax;
253 G4double sinCur = sinHalf; << 182 yMax=yoffset+fRMax;
254 G4double cosCur = cosHalf; << 183 if (pVoxelLimit.IsYLimited())
255 << 184 {
256 G4ThreeVectorList baseA(NSTEPS),baseB(NSTE << 185 if (yMin>pVoxelLimit.GetMaxYExtent()
257 for (G4int k=0; k<NSTEPS; ++k) << 186 ||yMax<pVoxelLimit.GetMinYExtent())
258 { << 187 {
259 baseA[k].set(rext*cosCur,rext*sinCur,-dz << 188 return false;
260 baseB[k].set(rext*cosCur,rext*sinCur, dz << 189 }
261 << 190 else
262 G4double sinTmp = sinCur; << 191 {
263 sinCur = sinCur*cosStep + cosCur*sinStep << 192 if (yMin<pVoxelLimit.GetMinYExtent())
264 cosCur = cosCur*cosStep - sinTmp*sinStep << 193 {
265 } << 194 yMin=pVoxelLimit.GetMinYExtent();
266 std::vector<const G4ThreeVectorList *> pol << 195 }
267 polygons[0] = &baseA; << 196 if (yMax>pVoxelLimit.GetMaxYExtent())
268 polygons[1] = &baseB; << 197 {
269 G4BoundingEnvelope benv(bmin,bmax,polygons << 198 yMax=pVoxelLimit.GetMaxYExtent();
270 exist = benv.CalculateExtent(pAxis,pVoxelL << 199 }
271 } << 200 }
272 else << 201 }
273 { << 202
274 G4double sinStart = GetSinStartPhi(); << 203
275 G4double cosStart = GetCosStartPhi(); << 204 zoffset=pTransform.NetTranslation().z();
276 G4double sinEnd = GetSinEndPhi(); << 205 zMin=zoffset-fDz;
277 G4double cosEnd = GetCosEndPhi(); << 206 zMax=zoffset+fDz;
278 G4double sinCur = sinStart*cosHalf + cos << 207 if (pVoxelLimit.IsZLimited())
279 G4double cosCur = cosStart*cosHalf - sin << 208 {
280 << 209 if (zMin>pVoxelLimit.GetMaxZExtent()
281 // set quadrilaterals << 210 ||zMax<pVoxelLimit.GetMinZExtent())
282 G4ThreeVectorList pols[NSTEPS+2]; << 211 {
283 for (G4int k=0; k<ksteps+2; ++k) pols[k].r << 212 return false;
284 pols[0][0].set(rmin*cosStart,rmin*sinStart << 213 }
285 pols[0][1].set(rmin*cosStart,rmin*sinStart << 214 else
286 pols[0][2].set(rmax*cosStart,rmax*sinStart << 215 {
287 pols[0][3].set(rmax*cosStart,rmax*sinStart << 216 if (zMin<pVoxelLimit.GetMinZExtent())
288 for (G4int k=1; k<ksteps+1; ++k) << 217 {
289 { << 218 zMin=pVoxelLimit.GetMinZExtent();
290 pols[k][0].set(rmin*cosCur,rmin*sinCur, << 219 }
291 pols[k][1].set(rmin*cosCur,rmin*sinCur,- << 220 if (zMax>pVoxelLimit.GetMaxZExtent())
292 pols[k][2].set(rext*cosCur,rext*sinCur,- << 221 {
293 pols[k][3].set(rext*cosCur,rext*sinCur, << 222 zMax=pVoxelLimit.GetMaxZExtent();
294 << 223 }
295 G4double sinTmp = sinCur; << 224 }
296 sinCur = sinCur*cosStep + cosCur*sinStep << 225 }
297 cosCur = cosCur*cosStep - sinTmp*sinStep << 226
298 } << 227 // Known to cut cylinder
299 pols[ksteps+1][0].set(rmin*cosEnd,rmin*sin << 228 switch (pAxis)
300 pols[ksteps+1][1].set(rmin*cosEnd,rmin*sin << 229 {
301 pols[ksteps+1][2].set(rmax*cosEnd,rmax*sin << 230 case kXAxis:
302 pols[ksteps+1][3].set(rmax*cosEnd,rmax*sin << 231 yoff1=yoffset-yMin;
303 << 232 yoff2=yMax-yoffset;
304 // set envelope and calculate extent << 233 if (yoff1>=0&&yoff2>=0)
305 std::vector<const G4ThreeVectorList *> pol << 234 {
306 polygons.resize(ksteps+2); << 235 // Y limits cross max/min x => no change
307 for (G4int k=0; k<ksteps+2; ++k) polygons[ << 236 pMin=xMin;
308 G4BoundingEnvelope benv(bmin,bmax,polygons << 237 pMax=xMax;
309 exist = benv.CalculateExtent(pAxis,pVoxelL << 238 }
310 } << 239 else
311 return exist; << 240 {
>> 241 // Y limits don't cross max/min x => compute max delta x, hence new mins/maxs
>> 242 diff1=sqrt(fRMax*fRMax-yoff1*yoff1);
>> 243 diff2=sqrt(fRMax*fRMax-yoff2*yoff2);
>> 244 maxDiff=(diff1>diff2) ? diff1:diff2;
>> 245 newMin=xoffset-maxDiff;
>> 246 newMax=xoffset+maxDiff;
>> 247 pMin=(newMin<xMin) ? xMin : newMin;
>> 248 pMax=(newMax>xMax) ? xMax : newMax;
>> 249 }
>> 250
>> 251 break;
>> 252 case kYAxis:
>> 253 xoff1=xoffset-xMin;
>> 254 xoff2=xMax-xoffset;
>> 255 if (xoff1>=0&&xoff2>=0)
>> 256 {
>> 257 // X limits cross max/min y => no change
>> 258 pMin=yMin;
>> 259 pMax=yMax;
>> 260 }
>> 261 else
>> 262 {
>> 263 // X limits don't cross max/min y => compute max delta y, hence new mins/maxs
>> 264 diff1=sqrt(fRMax*fRMax-xoff1*xoff1);
>> 265 diff2=sqrt(fRMax*fRMax-xoff2*xoff2);
>> 266 maxDiff=(diff1>diff2) ? diff1:diff2;
>> 267 newMin=yoffset-maxDiff;
>> 268 newMax=yoffset+maxDiff;
>> 269 pMin=(newMin<yMin) ? yMin : newMin;
>> 270 pMax=(newMax>yMax) ? yMax : newMax;
>> 271 }
>> 272 break;
>> 273 case kZAxis:
>> 274 pMin=zMin;
>> 275 pMax=zMax;
>> 276 break;
>> 277 }
>> 278
>> 279 pMin-=kCarTolerance;
>> 280 pMax+=kCarTolerance;
>> 281
>> 282 return true;
>> 283
>> 284 }
>> 285 else
>> 286 {
>> 287 G4int i,noEntries,noBetweenSections4;
>> 288 G4bool existsAfterClip=false;
>> 289
>> 290 // Calculate rotated vertex coordinates
>> 291 G4ThreeVectorList *vertices;
>> 292
>> 293 vertices=CreateRotatedVertices(pTransform);
>> 294
>> 295 pMin=+kInfinity;
>> 296 pMax=-kInfinity;
>> 297
>> 298 noEntries=vertices->entries();
>> 299 noBetweenSections4=noEntries-4;
>> 300
>> 301 for (i=0;i<noEntries;i+=4)
>> 302 {
>> 303 ClipCrossSection(vertices,i,pVoxelLimit,pAxis,pMin,pMax);
>> 304 }
>> 305
>> 306 for (i=0;i<noBetweenSections4;i+=4)
>> 307 {
>> 308 ClipBetweenSections(vertices,i,pVoxelLimit,pAxis,pMin,pMax);
>> 309 }
>> 310
>> 311 if (pMin!=kInfinity||pMax!=-kInfinity)
>> 312 {
>> 313 existsAfterClip=true;
>> 314
>> 315 // Add 2*tolerance to avoid precision troubles
>> 316 pMin-=kCarTolerance;
>> 317 pMax+=kCarTolerance;
>> 318
>> 319 }
>> 320 else
>> 321 {
>> 322 // Check for case where completely enveloping clipping volume
>> 323 // If point inside then we are confident that the solid completely
>> 324 // envelopes the clipping volume. Hence set min/max extents according
>> 325 // to clipping volume extents along the specified axis.
>> 326 G4ThreeVector clipCentre(
>> 327 (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
>> 328 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
>> 329 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
>> 330
>> 331 if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
>> 332 {
>> 333 existsAfterClip=true;
>> 334 pMin=pVoxelLimit.GetMinExtent(pAxis);
>> 335 pMax=pVoxelLimit.GetMaxExtent(pAxis);
>> 336 }
>> 337 }
>> 338 delete vertices;
>> 339 return existsAfterClip;
>> 340 }
312 } 341 }
313 342
314 ////////////////////////////////////////////// 343 ///////////////////////////////////////////////////////////////////////////
315 // 344 //
316 // Return whether point inside/outside/on surf 345 // Return whether point inside/outside/on surface
317 346
318 EInside G4Tubs::Inside( const G4ThreeVector& p << 347 EInside G4Tubs::Inside(const G4ThreeVector& p) const
319 { 348 {
320 G4double r2,pPhi,tolRMin,tolRMax; << 349 G4double r2,pPhi,tolRMin,tolRMax;
321 EInside in = kOutside ; << 350 EInside in=kOutside;
322 << 351 if (fabs(p.z())<=fDz-kCarTolerance*0.5)
323 if (std::fabs(p.z()) <= fDz - halfCarToleran << 352 {
324 { << 353 r2=p.x()*p.x()+p.y()*p.y();
325 r2 = p.x()*p.x() + p.y()*p.y() ; << 354 if (fRMin) tolRMin=fRMin+kRadTolerance*0.5;
326 << 355 else tolRMin=0;
327 if (fRMin != 0.0) { tolRMin = fRMin + half << 356 tolRMax=fRMax-kRadTolerance*0.5;
328 else { tolRMin = 0 ; } << 357
329 << 358 if (r2>=tolRMin*tolRMin && r2<=tolRMax*tolRMax)
330 tolRMax = fRMax - halfRadTolerance ; << 359 {
331 << 360 if (fDPhi==2*M_PI||r2==0)
332 if ((r2 >= tolRMin*tolRMin) && (r2 <= tolR << 361 {
333 { << 362 in=kInside;
334 if ( fPhiFullTube ) << 363 }
335 { << 364 else
336 in = kInside ; << 365 {
337 } << 366 // Try inner tolerant phi boundaries (=>inside)
338 else << 367 // if not inside, try outer tolerant phi boundaries
339 { << 368 pPhi=atan2(p.y(),p.x());
340 // Try inner tolerant phi boundaries ( << 369 if (pPhi<0) pPhi+=2*M_PI; // 0<=pPhi<2pi
341 // if not inside, try outer tolerant p << 370 if (fSPhi>=0)
342 << 371 {
343 if ( (tolRMin==0) && (std::fabs(p.x()) << 372 if (pPhi>=fSPhi+kAngTolerance*0.5 &&
344 && (std::fabs(p.y()) << 373 pPhi<=fSPhi+fDPhi-kAngTolerance*0.5)
345 { << 374 {
346 in=kSurface; << 375 in=kInside;
347 } << 376 }
348 else << 377 else if (pPhi>=fSPhi-kAngTolerance*0.5 &&
349 { << 378 pPhi<=fSPhi+fDPhi+kAngTolerance*0.5)
350 pPhi = std::atan2(p.y(),p.x()) ; << 379 {
351 if ( pPhi < -halfAngTolerance ) { p << 380 in=kSurface;
352 << 381 }
353 if ( fSPhi >= 0 ) << 382 }
354 { << 383 else
355 if ( (std::fabs(pPhi) < halfAngTol << 384 {
356 && (std::fabs(fSPhi + fDPhi - tw << 385 if (pPhi<fSPhi+2*M_PI) pPhi+=2*M_PI;
357 { << 386 if (pPhi>=fSPhi+2*M_PI+kAngTolerance*0.5 &&
358 pPhi += twopi ; // 0 <= pPhi < 2 << 387 pPhi<=fSPhi+fDPhi+2*M_PI-kAngTolerance*0.5)
359 } << 388 {
360 if ( (pPhi >= fSPhi + halfAngToler << 389 in=kInside;
361 && (pPhi <= fSPhi + fDPhi - half << 390 }
362 { << 391 else if (pPhi>=fSPhi+2*M_PI-kAngTolerance*0.5 &&
363 in = kInside ; << 392 pPhi<=fSPhi+fDPhi+2*M_PI+kAngTolerance*0.5)
364 } << 393 {
365 else if ( (pPhi >= fSPhi - halfAng << 394 in=kSurface;
366 && (pPhi <= fSPhi + fDPhi + << 395 }
367 { << 396 }
368 in = kSurface ; << 397
369 } << 398
370 } << 399 }
371 else // fSPhi < 0 << 400 }
372 { << 401 else
373 if ( (pPhi <= fSPhi + twopi - half << 402 {
374 && (pPhi >= fSPhi + fDPhi + hal << 403 // Try generous boundaries
375 else if ( (pPhi <= fSPhi + twopi + << 404 tolRMin=fRMin-kRadTolerance*0.5;
376 && (pPhi >= fSPhi + fDPhi << 405 tolRMax=fRMax+kRadTolerance*0.5;
377 { << 406 if (tolRMin<0) tolRMin=0;
378 in = kSurface ; << 407 if (r2>=tolRMin*tolRMin && r2 <= tolRMax*tolRMax)
379 } << 408 {
380 else << 409 if (fDPhi==2*M_PI||r2==0)
381 { << 410 {
382 in = kInside ; << 411 // Continuous in phi or on z-axis
383 } << 412 in=kSurface;
384 } << 413 }
385 } << 414 else
386 } << 415 {
387 } << 416 // Try outer tolerant phi boundaries only
388 else // Try generous boundaries << 417 pPhi=atan2(p.y(),p.x());
389 { << 418 if (pPhi<0) pPhi+=2*M_PI; // 0<=pPhi<2pi
390 tolRMin = fRMin - halfRadTolerance ; << 419 if (fSPhi>=0)
391 tolRMax = fRMax + halfRadTolerance ; << 420 {
392 << 421 if (pPhi>=fSPhi-kAngTolerance*0.5 &&
393 if ( tolRMin < 0 ) { tolRMin = 0; } << 422 pPhi<=fSPhi+fDPhi+kAngTolerance*0.5)
394 << 423 {
395 if ( (r2 >= tolRMin*tolRMin) && (r2 <= t << 424 in=kSurface;
396 { << 425 }
397 if (fPhiFullTube || (r2 <=halfRadToler << 426 }
398 { // Continuous << 427 else
399 in = kSurface ; << 428 {
400 } << 429 if (pPhi<fSPhi+2*M_PI) pPhi+=2*M_PI;
401 else // Try outer tolerant phi boundar << 430 if (pPhi>=fSPhi+2*M_PI-kAngTolerance*0.5 &&
402 { << 431 pPhi<=fSPhi+fDPhi+2*M_PI+kAngTolerance*0.5)
403 pPhi = std::atan2(p.y(),p.x()) ; << 432 {
404 << 433 in=kSurface;
405 if ( pPhi < -halfAngTolerance) { pP << 434 }
406 if ( fSPhi >= 0 ) << 435 }
407 { << 436
408 if ( (std::fabs(pPhi) < halfAngTol << 437
409 && (std::fabs(fSPhi + fDPhi - tw << 438 }
410 { << 439 }
411 pPhi += twopi ; // 0 <= pPhi < 2 << 440 }
412 } << 441 }
413 if ( (pPhi >= fSPhi - halfAngToler << 442 else if (fabs(p.z())<=fDz+kCarTolerance*0.5)
414 && (pPhi <= fSPhi + fDPhi + half << 443 {
415 { << 444 // Check within tolerant r limits
416 in = kSurface ; << 445 r2=p.x()*p.x()+p.y()*p.y();
417 } << 446 tolRMin=fRMin-kRadTolerance*0.5;
418 } << 447 tolRMax=fRMax+kRadTolerance*0.5;
419 else // fSPhi < 0 << 448 if (tolRMin<0) tolRMin=0;
420 { << 449 if (r2>=tolRMin*tolRMin && r2 <= tolRMax*tolRMax)
421 if ( (pPhi <= fSPhi + twopi - half << 450 {
422 && (pPhi >= fSPhi + fDPhi + half << 451 if (fDPhi==2*M_PI||r2==0)
423 else << 452 {
424 { << 453 // Continuous in phi or on z-axis
425 in = kSurface ; << 454 in=kSurface;
426 } << 455 }
427 } << 456 else
428 } << 457 {
429 } << 458 // Try outer tolerant phi boundaries
430 } << 459 pPhi=atan2(p.y(),p.x());
431 } << 460 if (pPhi<0) pPhi+=2*M_PI; // 0<=pPhi<2pi
432 else if (std::fabs(p.z()) <= fDz + halfCarTo << 461 if (fSPhi>=0)
433 { / << 462 {
434 r2 = p.x()*p.x() + p.y()*p.y() ; << 463 if (pPhi>=fSPhi-kAngTolerance*0.5 &&
435 tolRMin = fRMin - halfRadTolerance ; << 464 pPhi<=fSPhi+fDPhi+kAngTolerance*0.5)
436 tolRMax = fRMax + halfRadTolerance ; << 465 {
437 << 466 in=kSurface;
438 if ( tolRMin < 0 ) { tolRMin = 0; } << 467 }
439 << 468 }
440 if ( (r2 >= tolRMin*tolRMin) && (r2 <= tol << 469 else
441 { << 470 {
442 if (fPhiFullTube || (r2 <=halfRadToleran << 471 if (pPhi<fSPhi+2*M_PI) pPhi+=2*M_PI;
443 { // Continuous i << 472 if (pPhi>=fSPhi+2*M_PI-kAngTolerance*0.5 &&
444 in = kSurface ; << 473 pPhi<=fSPhi+fDPhi+2*M_PI+kAngTolerance*0.5)
445 } << 474 {
446 else // Try outer tolerant phi boundarie << 475 in=kSurface;
447 { << 476 }
448 pPhi = std::atan2(p.y(),p.x()) ; << 477 }
449 << 478
450 if ( pPhi < -halfAngTolerance ) { pPh << 479 }
451 if ( fSPhi >= 0 ) << 480 }
452 { << 481 }
453 if ( (std::fabs(pPhi) < halfAngToler << 482 return in;
454 && (std::fabs(fSPhi + fDPhi - twop <<
455 { <<
456 pPhi += twopi ; // 0 <= pPhi < 2pi <<
457 } <<
458 if ( (pPhi >= fSPhi - halfAngToleran <<
459 && (pPhi <= fSPhi + fDPhi + halfAn <<
460 { <<
461 in = kSurface; <<
462 } <<
463 } <<
464 else // fSPhi < 0 <<
465 { <<
466 if ( (pPhi <= fSPhi + twopi - halfAn <<
467 && (pPhi >= fSPhi + fDPhi + halfA <<
468 else <<
469 { <<
470 in = kSurface ; <<
471 } <<
472 } <<
473 } <<
474 } <<
475 } <<
476 return in; <<
477 } 483 }
478 484
479 ////////////////////////////////////////////// 485 ///////////////////////////////////////////////////////////////////////////
480 // 486 //
481 // Return unit normal of surface closest to p 487 // Return unit normal of surface closest to p
482 // - note if point on z axis, ignore phi divid 488 // - note if point on z axis, ignore phi divided sides
483 // - unsafe if point close to z axis a rmin=0 489 // - unsafe if point close to z axis a rmin=0 - no explicit checks
484 490
485 G4ThreeVector G4Tubs::SurfaceNormal( const G4T << 491 G4ThreeVector G4Tubs::SurfaceNormal( const G4ThreeVector& p) const
486 { <<
487 G4int noSurfaces = 0; <<
488 G4double rho, pPhi; <<
489 G4double distZ, distRMin, distRMax; <<
490 G4double distSPhi = kInfinity, distEPhi = kI <<
491 <<
492 G4ThreeVector norm, sumnorm(0.,0.,0.); <<
493 G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0); <<
494 G4ThreeVector nR, nPs, nPe; <<
495 <<
496 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()); <<
497 <<
498 distRMin = std::fabs(rho - fRMin); <<
499 distRMax = std::fabs(rho - fRMax); <<
500 distZ = std::fabs(std::fabs(p.z()) - fDz) <<
501 <<
502 if (!fPhiFullTube) // Protected against ( <<
503 { <<
504 if ( rho > halfCarTolerance ) <<
505 { <<
506 pPhi = std::atan2(p.y(),p.x()); <<
507 <<
508 if (pPhi < fSPhi-halfCarTolerance) <<
509 else if (pPhi > fSPhi+fDPhi+halfCarToler <<
510 <<
511 distSPhi = std::fabs( pPhi - fSPhi ); <<
512 distEPhi = std::fabs( pPhi - fSPhi - fDP <<
513 } <<
514 else if ( fRMin == 0.0 ) <<
515 { <<
516 distSPhi = 0.; <<
517 distEPhi = 0.; <<
518 } <<
519 nPs = G4ThreeVector( sinSPhi, -cosSPhi, 0 <<
520 nPe = G4ThreeVector( -sinEPhi, cosEPhi, 0 <<
521 } <<
522 if ( rho > halfCarTolerance ) { nR = G4Three <<
523 <<
524 if( distRMax <= halfCarTolerance ) <<
525 { <<
526 ++noSurfaces; <<
527 sumnorm += nR; <<
528 } <<
529 if( (fRMin != 0.0) && (distRMin <= halfCarTo <<
530 { <<
531 ++noSurfaces; <<
532 sumnorm -= nR; <<
533 } <<
534 if( fDPhi < twopi ) <<
535 { <<
536 if (distSPhi <= halfAngTolerance) <<
537 { <<
538 ++noSurfaces; <<
539 sumnorm += nPs; <<
540 } <<
541 if (distEPhi <= halfAngTolerance) <<
542 { <<
543 ++noSurfaces; <<
544 sumnorm += nPe; <<
545 } <<
546 } <<
547 if (distZ <= halfCarTolerance) <<
548 { <<
549 ++noSurfaces; <<
550 if ( p.z() >= 0.) { sumnorm += nZ; } <<
551 else { sumnorm -= nZ; } <<
552 } <<
553 if ( noSurfaces == 0 ) <<
554 { <<
555 #ifdef G4CSGDEBUG <<
556 G4Exception("G4Tubs::SurfaceNormal(p)", "G <<
557 JustWarning, "Point p is not o <<
558 G4long oldprc = G4cout.precision(20); <<
559 G4cout<< "G4Tubs::SN ( "<<p.x()<<", "<<p.y <<
560 << G4endl << G4endl; <<
561 G4cout.precision(oldprc) ; <<
562 #endif <<
563 norm = ApproxSurfaceNormal(p); <<
564 } <<
565 else if ( noSurfaces == 1 ) { norm = sumnor <<
566 else { norm = sumnor <<
567 <<
568 return norm; <<
569 } <<
570 <<
571 ////////////////////////////////////////////// <<
572 // <<
573 // Algorithm for SurfaceNormal() following the <<
574 // for points not on the surface <<
575 <<
576 G4ThreeVector G4Tubs::ApproxSurfaceNormal( con <<
577 { 492 {
578 ENorm side ; << 493 ENorm side;
579 G4ThreeVector norm ; << 494 G4ThreeVector norm;
580 G4double rho, phi ; << 495 G4double rho,phi;
581 G4double distZ, distRMin, distRMax, distSPhi << 496 G4double distZ,distRMin,distRMax,distSPhi,distEPhi,distMin;
582 << 497
583 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; << 498 rho=sqrt(p.x()*p.x()+p.y()*p.y());
584 << 499
585 distRMin = std::fabs(rho - fRMin) ; << 500 distRMin=fabs(rho-fRMin);
586 distRMax = std::fabs(rho - fRMax) ; << 501 distRMax=fabs(rho-fRMax);
587 distZ = std::fabs(std::fabs(p.z()) - fDz) << 502 distZ=fabs(fabs(p.z())-fDz);
588 << 503
589 if (distRMin < distRMax) // First minimum << 504 // First minimum
590 { << 505
591 if ( distZ < distRMin ) << 506 if (distRMin < distRMax)
592 { << 507 {
593 distMin = distZ ; << 508 if (distZ<distRMin)
594 side = kNZ ; << 509 {
595 } << 510 distMin=distZ;
596 else << 511 side=kNZ;
597 { << 512 }
598 distMin = distRMin ; << 513 else
599 side = kNRMin ; << 514 {
600 } << 515 distMin=distRMin;
>> 516 side=kNRMin;
>> 517 }
601 } 518 }
602 else 519 else
603 { 520 {
604 if ( distZ < distRMax ) << 521 if (distZ<distRMax)
605 { << 522 {
606 distMin = distZ ; << 523 distMin=distZ;
607 side = kNZ ; << 524 side=kNZ;
608 } << 525 }
609 else << 526 else
610 { << 527 {
611 distMin = distRMax ; << 528 distMin=distRMax;
612 side = kNRMax ; << 529 side=kNRMax;
613 } << 530 }
614 } << 531 }
615 if (!fPhiFullTube && (rho != 0.0) ) // Pro << 532
616 { << 533 if (fDPhi < 2.0*M_PI && rho )
617 phi = std::atan2(p.y(),p.x()) ; << 534 {
618 << 535 // Protected against (0,0,z) (above)
619 if ( phi < 0 ) { phi += twopi; } << 536 phi=atan2(p.y(),p.x());
620 << 537
621 if ( fSPhi < 0 ) << 538 if (phi<0) phi+=2*M_PI;
622 { << 539
623 distSPhi = std::fabs(phi - (fSPhi + twop << 540 if (fSPhi<0)
624 } << 541 {
625 else << 542 distSPhi=fabs(phi-(fSPhi+2.0*M_PI))*rho;
626 { << 543 }
627 distSPhi = std::fabs(phi - fSPhi)*rho ; << 544 else
628 } << 545 {
629 distEPhi = std::fabs(phi - fSPhi - fDPhi)* << 546 distSPhi=fabs(phi-fSPhi)*rho;
630 << 547 }
631 if (distSPhi < distEPhi) // Find new minim << 548 distEPhi=fabs(phi-fSPhi-fDPhi)*rho;
632 { << 549
633 if ( distSPhi < distMin ) << 550 // Find new minimum
634 { << 551 if (distSPhi < distEPhi)
635 side = kNSPhi ; << 552 {
636 } << 553 if (distSPhi<distMin)
637 } << 554 {
638 else << 555 side = kNSPhi ;
639 { << 556 }
640 if ( distEPhi < distMin ) << 557 }
641 { << 558 else
642 side = kNEPhi ; << 559 {
643 } << 560 if (distEPhi<distMin)
644 } << 561 {
645 } << 562 side=kNEPhi;
646 switch ( side ) << 563 }
647 { << 564 }
648 case kNRMin : // Inner radius << 565 }
649 { << 566
650 norm = G4ThreeVector(-p.x()/rho, -p.y()/ << 567 switch (side)
651 break ; << 568 {
652 } << 569 case kNRMin: // Inner radius
653 case kNRMax : // Outer radius << 570 norm=G4ThreeVector(-p.x()/rho,-p.y()/rho,0);
654 { << 571 break;
655 norm = G4ThreeVector(p.x()/rho, p.y()/rh << 572 case kNRMax: // Outer radius
656 break ; << 573 norm=G4ThreeVector(p.x()/rho,p.y()/rho,0);
657 } << 574 break;
658 case kNZ : // + or - dz << 575 case kNZ: // + or - dz
659 { << 576 if (p.z()>0)
660 if ( p.z() > 0 ) { norm = G4ThreeVector << 577 { norm=G4ThreeVector(0,0,1); }
661 else { norm = G4ThreeVector << 578 else
662 break ; << 579 { norm=G4ThreeVector(0,0,-1); }
663 } << 580 break;
664 case kNSPhi: << 581 case kNSPhi:
665 { << 582 norm=G4ThreeVector(sin(fSPhi),-cos(fSPhi),0);
666 norm = G4ThreeVector(sinSPhi, -cosSPhi, << 583 break;
667 break ; << 584 case kNEPhi:
668 } << 585 norm=G4ThreeVector(-sin(fSPhi+fDPhi),cos(fSPhi+fDPhi),0);
669 case kNEPhi: << 586 break;
670 { << 587 default:
671 norm = G4ThreeVector(-sinEPhi, cosEPhi, << 588 G4Exception("Logic error in G4Tubs::SurfaceNormal");
672 break; << 589 break;
673 } << 590
674 default: // Should never reach this c << 591 } // end case
675 { <<
676 DumpInfo(); <<
677 G4Exception("G4Tubs::ApproxSurfaceNormal <<
678 "GeomSolids1002", JustWarnin <<
679 "Undefined side for valid su <<
680 break ; <<
681 } <<
682 } <<
683 return norm; 592 return norm;
684 } 593 }
685 594
686 ////////////////////////////////////////////// 595 ////////////////////////////////////////////////////////////////////
687 // 596 //
688 // 597 //
689 // Calculate distance to shape from outside, a 598 // Calculate distance to shape from outside, along normalised vector
690 // - return kInfinity if no intersection, or i 599 // - return kInfinity if no intersection, or intersection distance <= tolerance
691 // 600 //
692 // - Compute the intersection with the z plane << 601 // - Compute the intersection with the z planes
693 // - if at valid r, phi, return 602 // - if at valid r, phi, return
694 // 603 //
695 // -> If point is outer outer radius, compute 604 // -> If point is outer outer radius, compute intersection with rmax
696 // - if at valid phi,z return 605 // - if at valid phi,z return
697 // 606 //
698 // -> Compute intersection with inner radius, 607 // -> Compute intersection with inner radius, taking largest +ve root
699 // - if valid (in z,phi), save intersct 608 // - if valid (in z,phi), save intersction
700 // 609 //
701 // -> If phi segmented, compute intersectio 610 // -> If phi segmented, compute intersections with phi half planes
702 // - return smallest of valid phi inter 611 // - return smallest of valid phi intersections and
703 // inner radius intersection 612 // inner radius intersection
704 // 613 //
705 // NOTE: 614 // NOTE:
706 // - 'if valid' implies tolerant checking of i << 615 // - Precalculations for phi trigonometry are Done `just in time'
>> 616 // - `if valid' implies tolerant checking of intersection points
707 617
708 G4double G4Tubs::DistanceToIn( const G4ThreeVe << 618 G4double G4Tubs::DistanceToIn(const G4ThreeVector& p,
709 const G4ThreeVe << 619 const G4ThreeVector& v ) const
710 { 620 {
711 G4double snxt = kInfinity ; // snxt = d << 621 G4double snxt=kInfinity; // snxt = default return value
712 G4double tolORMin2, tolIRMax2 ; // 'generou <<
713 G4double tolORMax2, tolIRMin2, tolODz, tolID <<
714 const G4double dRmax = 100.*fRMax; <<
715 <<
716 // Intersection point variables <<
717 // <<
718 G4double Dist, sd, xi, yi, zi, rho2, inum, i <<
719 G4double t1, t2, t3, b, c, d ; // Quadra <<
720 <<
721 // Calculate tolerant rmin and rmax <<
722 <<
723 if (fRMin > kRadTolerance) <<
724 { <<
725 tolORMin2 = (fRMin - halfRadTolerance)*(fR <<
726 tolIRMin2 = (fRMin + halfRadTolerance)*(fR <<
727 } <<
728 else <<
729 { <<
730 tolORMin2 = 0.0 ; <<
731 tolIRMin2 = 0.0 ; <<
732 } <<
733 tolORMax2 = (fRMax + halfRadTolerance)*(fRMa <<
734 tolIRMax2 = (fRMax - halfRadTolerance)*(fRMa <<
735 622
736 // Intersection with Z surfaces << 623 // Precalculated trig for phi intersections - used by r,z intersections to
>> 624 // check validity
737 625
738 tolIDz = fDz - halfCarTolerance ; << 626 G4bool seg; // true if segmented
739 tolODz = fDz + halfCarTolerance ; <<
740 627
741 if (std::fabs(p.z()) >= tolIDz) << 628 G4double hDPhi,hDPhiOT,hDPhiIT,cosHDPhiOT,cosHDPhiIT;
742 { << 629 // half dphi + outer tolerance
743 if ( p.z()*v.z() < 0 ) // at +Z going i <<
744 { <<
745 sd = (std::fabs(p.z()) - fDz)/std::fabs( <<
746 630
747 if(sd < 0.0) { sd = 0.0; } << 631 G4double cPhi,sinCPhi,cosCPhi; // central phi
748 632
749 xi = p.x() + sd*v.x() ; << 633 G4double tolORMin2,tolIRMax2; // `generous' radii squared
750 yi = p.y() + sd*v.y() ; <<
751 rho2 = xi*xi + yi*yi ; <<
752 634
753 // Check validity of intersection << 635 G4double tolORMax2,tolIRMin2,tolODz,tolIDz;
754 636
755 if ((tolIRMin2 <= rho2) && (rho2 <= tolI << 637 // Intersection point variables
756 { << 638 G4double Dist,s,xi,yi,zi,rho2,inum,iden,cosPsi;
757 if (!fPhiFullTube && (rho2 != 0.0)) <<
758 { <<
759 // Psi = angle made with central (av <<
760 // <<
761 inum = xi*cosCPhi + yi*sinCPhi ; <<
762 iden = std::sqrt(rho2) ; <<
763 cosPsi = inum/iden ; <<
764 if (cosPsi >= cosHDPhiIT) { return <<
765 } <<
766 else <<
767 { <<
768 return sd ; <<
769 } <<
770 } <<
771 } <<
772 else <<
773 { <<
774 if ( snxt<halfCarTolerance ) { snxt=0; <<
775 return snxt ; // On/outside extent, and <<
776 // -> cannot intersect <<
777 } <<
778 } <<
779 639
780 // -> Can not intersect z surfaces <<
781 // <<
782 // Intersection with rmax (possible return) <<
783 // <<
784 // Intersection point (xi,yi,zi) on line x=p <<
785 // <<
786 // Intersects with x^2+y^2=R^2 <<
787 // <<
788 // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v. <<
789 // t1 t2 <<
790 <<
791 t1 = 1.0 - v.z()*v.z() ; <<
792 t2 = p.x()*v.x() + p.y()*v.y() ; <<
793 t3 = p.x()*p.x() + p.y()*p.y() ; <<
794 640
795 if ( t1 > 0 ) // Check not || to z ax << 641 G4double t1,t2,t3,b,c,d; // Quadratic solver variables
796 { <<
797 b = t2/t1 ; <<
798 c = t3 - fRMax*fRMax ; <<
799 if ((t3 >= tolORMax2) && (t2<0)) // This <<
800 { <<
801 // Try outer cylinder intersection <<
802 // c=(t3-fRMax*fRMax)/t1; <<
803 642
804 c /= t1 ; << 643 G4double Comp;
805 d = b*b - c ; << 644 G4double cosSPhi,sinSPhi; // Trig for phi start intersect
806 645
807 if (d >= 0) // If real root << 646 G4double ePhi,cosEPhi,sinEPhi; // for phi end intersect
808 { <<
809 sd = c/(-b+std::sqrt(d)); <<
810 if (sd >= 0) // If 'forwards' <<
811 { <<
812 if ( sd>dRmax ) // Avoid rounding er <<
813 { // 64 bits systems. <<
814 G4double fTerm = sd-std::fmod(sd,d <<
815 sd = fTerm + DistanceToIn(p+fTerm* <<
816 } <<
817 // Check z intersection <<
818 // <<
819 zi = p.z() + sd*v.z() ; <<
820 if (std::fabs(zi)<=tolODz) <<
821 { <<
822 // Z ok. Check phi intersection if <<
823 // <<
824 if (fPhiFullTube) <<
825 { <<
826 return sd ; <<
827 } <<
828 else <<
829 { <<
830 xi = p.x() + sd*v.x() ; <<
831 yi = p.y() + sd*v.y() ; <<
832 cosPsi = (xi*cosCPhi + yi*sinCPh <<
833 if (cosPsi >= cosHDPhiIT) { ret <<
834 } <<
835 } // end if std::fabs(zi) <<
836 } // end if (sd>=0) <<
837 } // end if (d>=0) <<
838 } // end if (r>=fRMax) <<
839 else <<
840 { <<
841 // Inside outer radius : <<
842 // check not inside, and heading through <<
843 647
844 if ((t3 > tolIRMin2) && (t2 < 0) && (std << 648 // Set phi divided flag and precalcs
845 { <<
846 // Inside both radii, delta r -ve, ins <<
847 649
848 if (!fPhiFullTube) << 650 if (fDPhi<2.0*M_PI)
849 { << 651 {
850 inum = p.x()*cosCPhi + p.y()*sinCP << 652 seg=true;
851 iden = std::sqrt(t3) ; << 653 hDPhi=0.5*fDPhi; // half delta phi
852 cosPsi = inum/iden ; << 654 cPhi=fSPhi+hDPhi;;
853 if (cosPsi >= cosHDPhiIT) << 655 hDPhiOT=hDPhi+0.5*kAngTolerance; // outers tol' half delta phi
854 { << 656 hDPhiIT=hDPhi-0.5*kAngTolerance;
855 // In the old version, the small n << 657 sinCPhi=sin(cPhi);
856 // on surface was not taken in acc << 658 cosCPhi=cos(cPhi);
857 // New version: check the tangent << 659 cosHDPhiOT=cos(hDPhiOT);
858 // if no intersection, return kInf << 660 cosHDPhiIT=cos(hDPhiIT);
859 // return sd. <<
860 // <<
861 c = t3-fRMax*fRMax; <<
862 if ( c<=0.0 ) <<
863 { <<
864 return 0.0; <<
865 } <<
866 else <<
867 { <<
868 c = c/t1 ; <<
869 d = b*b-c; <<
870 if ( d>=0.0 ) <<
871 { <<
872 snxt = c/(-b+std::sqrt(d)); // <<
873 // <<
874 if ( snxt < halfCarTolerance ) <<
875 return snxt ; <<
876 } <<
877 else <<
878 { <<
879 return kInfinity; <<
880 } <<
881 } <<
882 } <<
883 } <<
884 else <<
885 { <<
886 // In the old version, the small neg <<
887 // on surface was not taken in accou <<
888 // New version: check the tangent fo <<
889 // if no intersection, return kInfin <<
890 // return sd. <<
891 // <<
892 c = t3 - fRMax*fRMax; <<
893 if ( c<=0.0 ) <<
894 { <<
895 return 0.0; <<
896 } <<
897 else <<
898 { <<
899 c = c/t1 ; <<
900 d = b*b-c; <<
901 if ( d>=0.0 ) <<
902 { <<
903 snxt= c/(-b+std::sqrt(d)); // us <<
904 // fo <<
905 if ( snxt < halfCarTolerance ) { <<
906 return snxt ; <<
907 } <<
908 else <<
909 { <<
910 return kInfinity; <<
911 } <<
912 } <<
913 } // end if (!fPhiFullTube) <<
914 } // end if (t3>tolIRMin2) <<
915 } // end if (Inside Outer Radius) <<
916 if ( fRMin != 0.0 ) // Try inner cylind <<
917 { <<
918 c = (t3 - fRMin*fRMin)/t1 ; <<
919 d = b*b - c ; <<
920 if ( d >= 0.0 ) // If real root <<
921 { <<
922 // Always want 2nd root - we are outsi <<
923 // - If on surface of rmin also need f <<
924 <<
925 sd =( b > 0. )? c/(-b - std::sqrt(d)) <<
926 if (sd >= -halfCarTolerance) // check <<
927 { <<
928 // Check z intersection <<
929 // <<
930 if(sd < 0.0) { sd = 0.0; } <<
931 if ( sd>dRmax ) // Avoid rounding er <<
932 { // 64 bits systems. <<
933 G4double fTerm = sd-std::fmod(sd,d <<
934 sd = fTerm + DistanceToIn(p+fTerm* <<
935 } <<
936 zi = p.z() + sd*v.z() ; <<
937 if (std::fabs(zi) <= tolODz) <<
938 { <<
939 // Z ok. Check phi <<
940 // <<
941 if ( fPhiFullTube ) <<
942 { <<
943 return sd ; <<
944 } <<
945 else <<
946 { <<
947 xi = p.x() + sd*v.x() ; <<
948 yi = p.y() + sd*v.y() ; <<
949 cosPsi = (xi*cosCPhi + yi*sinCPh <<
950 if (cosPsi >= cosHDPhiIT) <<
951 { <<
952 // Good inner radius isect <<
953 // - but earlier phi isect sti <<
954 <<
955 snxt = sd ; <<
956 } <<
957 } <<
958 } // end if std::fabs(zi) <<
959 } // end if (sd>=0) <<
960 } // end if (d>=0) <<
961 } // end if (fRMin) <<
962 } 661 }
963 << 662 else
964 // Phi segment intersection <<
965 // <<
966 // o Tolerant of points inside phi planes by <<
967 // <<
968 // o NOTE: Large duplication of code between <<
969 // -> only diffs: sphi -> ephi, Comp <<
970 // intersection check <=0 -> >=0 <<
971 // -> use some form of loop Construc <<
972 // <<
973 if ( !fPhiFullTube ) <<
974 { 663 {
975 // First phi surface (Starting phi) << 664 seg=false;
976 // << 665 }
977 Comp = v.x()*sinSPhi - v.y()*cosSPhi ; <<
978 <<
979 if ( Comp < 0 ) // Component in outwards <<
980 { <<
981 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ; <<
982 <<
983 if ( Dist < halfCarTolerance ) <<
984 { <<
985 sd = Dist/Comp ; <<
986 <<
987 if (sd < snxt) <<
988 { <<
989 if ( sd < 0 ) { sd = 0.0; } <<
990 zi = p.z() + sd*v.z() ; <<
991 if ( std::fabs(zi) <= tolODz ) <<
992 { <<
993 xi = p.x() + sd*v.x() ; <<
994 yi = p.y() + sd*v.y() ; <<
995 rho2 = xi*xi + yi*yi ; <<
996 <<
997 if ( ( (rho2 >= tolIRMin2) && (rho <<
998 || ( (rho2 > tolORMin2) && (rho <<
999 && ( v.y()*cosSPhi - v.x()*sin <<
1000 && ( v.x()*cosSPhi + v.y()*si <<
1001 || ( (rho2 > tolIRMax2) && (rho <<
1002 && (v.y()*cosSPhi - v.x()*sin <<
1003 && (v.x()*cosSPhi + v.y()*sin <<
1004 { <<
1005 // z and r intersections good <<
1006 // - check intersecting with co <<
1007 // <<
1008 if ((yi*cosCPhi-xi*sinCPhi) <= <<
1009 } <<
1010 } <<
1011 } <<
1012 } <<
1013 } <<
1014 <<
1015 // Second phi surface (Ending phi) <<
1016 666
1017 Comp = -(v.x()*sinEPhi - v.y()*cosEPhi << 667 // Calculate tolerant rmin and rmax
1018 668
1019 if (Comp < 0 ) // Component in outwards << 669 if (fRMin > kRadTolerance)
1020 { << 670 {
1021 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) << 671 tolORMin2=(fRMin-0.5*kRadTolerance)*(fRMin-0.5*kRadTolerance);
>> 672 tolIRMin2=(fRMin+0.5*kRadTolerance)*(fRMin+0.5*kRadTolerance);
>> 673 }
>> 674 else
>> 675 {
>> 676 tolORMin2=0;
>> 677 tolIRMin2=0;
>> 678 }
>> 679 tolORMax2=(fRMax+0.5*kRadTolerance)*(fRMax+0.5*kRadTolerance);
>> 680 tolIRMax2=(fRMax-0.5*kRadTolerance)*(fRMax-0.5*kRadTolerance);
1022 681
1023 if ( Dist < halfCarTolerance ) << 682 // Intersection with Z surfaces
1024 { <<
1025 sd = Dist/Comp ; <<
1026 683
1027 if (sd < snxt) << 684 tolIDz=fDz-kCarTolerance*0.5;
1028 { << 685 tolODz=fDz+kCarTolerance*0.5;
1029 if ( sd < 0 ) { sd = 0; } << 686
1030 zi = p.z() + sd*v.z() ; << 687 if (fabs(p.z()) >= tolIDz)
1031 if ( std::fabs(zi) <= tolODz ) << 688 {
1032 { << 689 if (p.z()*v.z()<0) // at +Z going in -Z or visa versa
1033 xi = p.x() + sd*v.x() ; << 690 {
1034 yi = p.y() + sd*v.y() ; << 691 s=(fabs(p.z())-fDz)/fabs(v.z()); // Z intersect distance
1035 rho2 = xi*xi + yi*yi ; << 692 if(s < 0.0) s = 0.0 ;
1036 if ( ( (rho2 >= tolIRMin2) && (rh << 693 xi=p.x()+s*v.x(); // Intersection coords
1037 || ( (rho2 > tolORMin2) && ( << 694 yi=p.y()+s*v.y();
1038 && (v.x()*sinEPhi - v.y()*c << 695 rho2=xi*xi+yi*yi;
1039 && (v.x()*cosEPhi + v.y()*s << 696
1040 || ( (rho2 > tolIRMax2) && (r << 697 // Check validity of intersection
1041 && (v.x()*sinEPhi - v.y()*c << 698
1042 && (v.x()*cosEPhi + v.y()*s << 699 if (tolIRMin2 <= rho2 && rho2 <= tolIRMax2)
1043 { << 700 {
1044 // z and r intersections good << 701 if (seg && rho2)
1045 // - check intersecting with co << 702 {
1046 // << 703 // Psi = angle made with central (average) phi of shape
1047 if ( (yi*cosCPhi-xi*sinCPhi) >= << 704
1048 } //?? >= << 705 inum = xi*cosCPhi+yi*sinCPhi;
1049 } << 706 iden = sqrt(rho2);
1050 } << 707 cosPsi = inum/iden;
1051 } << 708
1052 } // Comp < 0 << 709 if (cosPsi >= cosHDPhiIT)
1053 } // !fPhiFullTube << 710 {
1054 if ( snxt<halfCarTolerance ) { snxt=0; } << 711 return s;
1055 return snxt ; << 712 }
>> 713 }
>> 714 else
>> 715 {
>> 716 return s;
>> 717 }
>> 718 }
>> 719 }
>> 720 else
>> 721 {
>> 722 return snxt ; // On/outside extent, and heading away
>> 723 // -> cannot intersect
>> 724 }
>> 725 }
>> 726
>> 727 // -> Can not intersect z surfaces
>> 728 //
>> 729 // Intersection with rmax (possible return) and rmin (must also check phi)
>> 730 //
>> 731 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
>> 732 //
>> 733 // Intersects with x^2+y^2=R^2
>> 734 //
>> 735 // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0
>> 736 // t1 t2 t3
>> 737
>> 738 t1=1.0-v.z()*v.z();
>> 739 t2=p.x()*v.x()+p.y()*v.y();
>> 740 t3=p.x()*p.x()+p.y()*p.y();
>> 741
>> 742 if (t1>0) // Check not || to z axis
>> 743 {
>> 744 b=t2/t1;
>> 745 c=t3-fRMax*fRMax;
>> 746
>> 747 if (t3 >= tolORMax2 && t2<0) // This also handles the tangent case
>> 748 {
>> 749 // Try outer cylinder intersection
>> 750 // c=(t3-fRMax*fRMax)/t1;
>> 751
>> 752 c /= t1 ;
>> 753 d=b*b-c;
>> 754
>> 755 if (d >= 0) // If real root
>> 756 {
>> 757 s = -b-sqrt(d) ;
>> 758
>> 759 if (s >= 0) // If 'forwards'
>> 760 {
>> 761 // Check z intersection
>> 762 zi=p.z()+s*v.z();
>> 763
>> 764 if (fabs(zi)<=tolODz)
>> 765 {
>> 766 // Z ok. Check phi intersection if reqd
>> 767
>> 768 if (!seg)
>> 769 {
>> 770 return s;
>> 771 }
>> 772 else
>> 773 {
>> 774 xi=p.x()+s*v.x();
>> 775 yi=p.y()+s*v.y();
>> 776 cosPsi=(xi*cosCPhi+yi*sinCPhi)/fRMax;
>> 777
>> 778 if (cosPsi >= cosHDPhiIT)
>> 779 {
>> 780 return s;
>> 781 }
>> 782 }
>> 783 } // end if fabs(zi)
>> 784 } // end if (s>=0)
>> 785 } // end if (d>=0)
>> 786 } // end if (r>=fRMax)
>> 787
>> 788 /* *************************************************************************
>> 789 // This code is shown here just to make explicit the logic contained
>> 790 // in the condition if (t3>tolIRMin2 && t2<0 && fabs(p.z())<=tolIDz),
>> 791 // which is placed below.
>> 792
>> 793 // else if (t3>=tolIRMax2 && t2<0 && fabs(p.z())<=tolIDz)
>> 794 // { // Point on Rmax surface
>> 795 // if (!seg)
>> 796 // {
>> 797 // return s = 0 ; // No Phi cut and move inside
>> 798 // }
>> 799 // else
>> 800 // {
>> 801 // cosPsi=(p.x()*cosCPhi+
>> 802 // p.y()*sinCPhi)/fRMax;
>> 803 // if (cosPsi>=cosHDPhiOT)
>> 804 // {
>> 805 // return s = 0 ; // On real Rmax surface and move inside
>> 806 // }
>> 807 // }
>> 808 // }
>> 809
>> 810 ************************************************************************* */
>> 811
>> 812 else
>> 813 {
>> 814 // Inside outer radius :
>> 815 // check not inside, and heading through tubs (-> 0 to in)
>> 816
>> 817 if (t3>tolIRMin2 && t2<0 && fabs(p.z())<=tolIDz)
>> 818 {
>> 819 // Inside both radii, delta r -ve, inside z extent
>> 820
>> 821 if (seg)
>> 822 {
>> 823 inum = p.x()*cosCPhi+p.y()*sinCPhi;
>> 824 iden = sqrt(t3);
>> 825 cosPsi = inum/iden;
>> 826
>> 827 if (cosPsi>=cosHDPhiIT)
>> 828 {
>> 829 return 0;
>> 830 }
>> 831 }
>> 832 else
>> 833 {
>> 834 return 0;
>> 835 }
>> 836 }
>> 837 }
>> 838 if (fRMin) // Try inner cylinder intersection
>> 839 {
>> 840 c=(t3-fRMin*fRMin)/t1;
>> 841 d=b*b-c;
>> 842
>> 843 if (d>=0) // If real root
>> 844 {
>> 845 // Always want 2nd root - we are outside and know rmax Hit was bad
>> 846 // - If on surface of rmin also need farthest root
>> 847
>> 848 s=-b+sqrt(d);
>> 849
>> 850 if (s >= -0.5*kCarTolerance) // check forwards
>> 851 {
>> 852 // Check z intersection
>> 853 if(s < 0.0) s = 0.0 ;
>> 854 zi=p.z()+s*v.z();
>> 855
>> 856 if (fabs(zi)<=tolODz)
>> 857 {
>> 858 // Z ok. Check phi
>> 859 if (!seg)
>> 860 {
>> 861 return s;
>> 862 }
>> 863 else
>> 864 {
>> 865 xi=p.x()+s*v.x();
>> 866 yi=p.y()+s*v.y();
>> 867 cosPsi=(xi*cosCPhi+yi*sinCPhi)/fRMin;
>> 868
>> 869 if (cosPsi >= cosHDPhiIT)
>> 870 {
>> 871 // Good inner radius isect - but earlier phi isect still possible
>> 872
>> 873 snxt=s;
>> 874 }
>> 875 }
>> 876 } // end if fabs(zi)
>> 877 } // end if (s>=0)
>> 878 } // end if (d>=0)
>> 879 } // end if (fRMin)
>> 880 }
>> 881
>> 882 // Phi segment intersection
>> 883 //
>> 884 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
>> 885 //
>> 886 // o NOTE: Large duplication of code between sphi & ephi checks
>> 887 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
>> 888 // intersection check <=0 -> >=0
>> 889 // -> use some form of loop Construct ?
>> 890 //
>> 891 if (seg)
>> 892 {
>> 893 // First phi surface (`S'tarting phi)
>> 894
>> 895 sinSPhi=sin(fSPhi);
>> 896 cosSPhi=cos(fSPhi);
>> 897 Comp=v.x()*sinSPhi-v.y()*cosSPhi;
>> 898
>> 899 if (Comp<0) // Component in outwards normal dirn
>> 900 {
>> 901 Dist=(p.y()*cosSPhi-p.x()*sinSPhi);
>> 902
>> 903 if (Dist<kCarTolerance*0.5)
>> 904 {
>> 905 s=Dist/Comp;
>> 906
>> 907 if (s<snxt)
>> 908 {
>> 909 if (s<0)
>> 910 {
>> 911 s=0;
>> 912 }
>> 913 zi=p.z()+s*v.z();
>> 914
>> 915 if (fabs(zi)<=tolODz)
>> 916 {
>> 917 xi=p.x()+s*v.x();
>> 918 yi=p.y()+s*v.y();
>> 919 rho2=xi*xi+yi*yi;
>> 920
>> 921 if ( ( rho2 >= tolIRMin2 && rho2 <= tolIRMax2 )
>> 922
>> 923 || ( rho2 > tolORMin2 && rho2 < tolIRMin2 &&
>> 924 ( v.y()*cosSPhi - v.x()*sinSPhi > 0 ) &&
>> 925 ( v.x()*cosSPhi + v.y()*sinSPhi >= 0 ) )
>> 926
>> 927 || ( rho2 > tolIRMax2 && rho2 < tolORMax2 &&
>> 928 ( v.y()*cosSPhi - v.x()*sinSPhi > 0 ) &&
>> 929 ( v.x()*cosSPhi + v.y()*sinSPhi < 0 ) ) )
>> 930 {
>> 931 // z and r intersections good - check intersecting with correct half-plane
>> 932
>> 933 if ((yi*cosCPhi-xi*sinCPhi) <= 0)
>> 934 {
>> 935 snxt=s;
>> 936 }
>> 937 }
>> 938 }
>> 939 }
>> 940 }
>> 941 }
>> 942
>> 943 // Second phi surface (`E'nding phi)
>> 944
>> 945 ePhi=fSPhi+fDPhi;
>> 946 sinEPhi=sin(ePhi);
>> 947 cosEPhi=cos(ePhi);
>> 948 Comp=-(v.x()*sinEPhi-v.y()*cosEPhi);
>> 949
>> 950 if (Comp<0) // Component in outwards normal dirn
>> 951 {
>> 952 Dist=-(p.y()*cosEPhi-p.x()*sinEPhi);
>> 953
>> 954 if (Dist<kCarTolerance*0.5)
>> 955 {
>> 956 s=Dist/Comp;
>> 957
>> 958 if (s<snxt)
>> 959 {
>> 960 if (s<0)
>> 961 {
>> 962 s=0;
>> 963 }
>> 964 zi=p.z()+s*v.z();
>> 965
>> 966 if (fabs(zi)<=tolODz)
>> 967 {
>> 968 xi=p.x()+s*v.x();
>> 969 yi=p.y()+s*v.y();
>> 970 rho2=xi*xi+yi*yi;
>> 971
>> 972 if ( ( rho2 >= tolIRMin2 && rho2 <= tolIRMax2 )
>> 973
>> 974 || ( rho2 > tolORMin2 && rho2 < tolIRMin2 &&
>> 975 ( v.x()*sinEPhi - v.y()*cosEPhi > 0 ) &&
>> 976 ( v.x()*cosEPhi + v.y()*sinEPhi >= 0 ) )
>> 977
>> 978 || ( rho2 > tolIRMax2 && rho2 < tolORMax2 &&
>> 979 ( v.x()*sinEPhi - v.y()*cosEPhi > 0 ) &&
>> 980 ( v.x()*cosEPhi + v.y()*sinEPhi < 0 ) ) )
>> 981 {
>> 982 // z and r intersections good - check intersecting with correct half-plane
>> 983
>> 984 if ((yi*cosCPhi-xi*sinCPhi)>=0)
>> 985 {
>> 986 snxt=s;
>> 987 }
>> 988 }
>> 989 }
>> 990 }
>> 991 }
>> 992 } // Comp < 0
>> 993 } // seg != 0
>> 994 return snxt;
1056 } 995 }
1057 << 996
1058 ///////////////////////////////////////////// 997 //////////////////////////////////////////////////////////////////
1059 // 998 //
1060 // Calculate distance to shape from outside, 999 // Calculate distance to shape from outside, along normalised vector
1061 // - return kInfinity if no intersection, or 1000 // - return kInfinity if no intersection, or intersection distance <= tolerance
1062 // 1001 //
1063 // - Compute the intersection with the z plan << 1002 // - Compute the intersection with the z planes
1064 // - if at valid r, phi, return 1003 // - if at valid r, phi, return
1065 // 1004 //
1066 // -> If point is outer outer radius, compute 1005 // -> If point is outer outer radius, compute intersection with rmax
1067 // - if at valid phi,z return 1006 // - if at valid phi,z return
1068 // 1007 //
1069 // -> Compute intersection with inner radius, 1008 // -> Compute intersection with inner radius, taking largest +ve root
1070 // - if valid (in z,phi), save intersc 1009 // - if valid (in z,phi), save intersction
1071 // 1010 //
1072 // -> If phi segmented, compute intersecti 1011 // -> If phi segmented, compute intersections with phi half planes
1073 // - return smallest of valid phi inte 1012 // - return smallest of valid phi intersections and
1074 // inner radius intersection 1013 // inner radius intersection
1075 // 1014 //
1076 // NOTE: 1015 // NOTE:
1077 // - Precalculations for phi trigonometry are 1016 // - Precalculations for phi trigonometry are Done `just in time'
1078 // - `if valid' implies tolerant checking of 1017 // - `if valid' implies tolerant checking of intersection points
1079 // Calculate distance (<= actual) to closes << 1018 // Calculate distance (<= actual) to closest surface of shape from outside
1080 // - Calculate distance to z, radial planes 1019 // - Calculate distance to z, radial planes
1081 // - Only to phi planes if outside phi extent 1020 // - Only to phi planes if outside phi extent
1082 // - Return 0 if point inside 1021 // - Return 0 if point inside
1083 1022
1084 G4double G4Tubs::DistanceToIn( const G4ThreeV << 1023 G4double G4Tubs::DistanceToIn(const G4ThreeVector& p) const
1085 { 1024 {
1086 G4double safe=0.0, rho, safe1, safe2, safe3 << 1025 G4double safe,rho,safe1,safe2,safe3;
1087 G4double safePhi, cosPsi ; << 1026 G4double phiC,cosPhiC,sinPhiC,safePhi,ePhi,cosPsi;
1088 <<
1089 rho = std::sqrt(p.x()*p.x() + p.y()*p.y() <<
1090 safe1 = fRMin - rho ; <<
1091 safe2 = rho - fRMax ; <<
1092 safe3 = std::fabs(p.z()) - fDz ; <<
1093 <<
1094 if ( safe1 > safe2 ) { safe = safe1; } <<
1095 else { safe = safe2; } <<
1096 if ( safe3 > safe ) { safe = safe3; } <<
1097 <<
1098 if ( (!fPhiFullTube) && ((rho) != 0.0) ) <<
1099 { <<
1100 // Psi=angle from central phi to point <<
1101 // <<
1102 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/ <<
1103 <<
1104 if ( cosPsi < cosHDPhi ) <<
1105 { <<
1106 // Point lies outside phi range <<
1107 1027
1108 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= << 1028 rho=sqrt(p.x()*p.x()+p.y()*p.y());
1109 { << 1029 safe1=fRMin-rho;
1110 safePhi = std::fabs(p.x()*sinSPhi - p << 1030 safe2=rho-fRMax;
1111 } << 1031 safe3=fabs(p.z())-fDz;
1112 else << 1032
1113 { << 1033 if (safe1>safe2) safe=safe1;
1114 safePhi = std::fabs(p.x()*sinEPhi - p << 1034 else safe=safe2;
1115 } << 1035 if (safe3>safe) safe=safe3;
1116 if ( safePhi > safe ) { safe = safePhi << 1036
1117 } << 1037 if (fDPhi < 2.0*M_PI && rho)
1118 } << 1038 {
1119 if ( safe < 0 ) { safe = 0; } << 1039 phiC=fSPhi+fDPhi*0.5;
1120 return safe ; << 1040 cosPhiC=cos(phiC);
>> 1041 sinPhiC=sin(phiC);
>> 1042 // Psi=angle from central phi to point
>> 1043 cosPsi=(p.x()*cosPhiC+p.y()*sinPhiC)/rho;
>> 1044 if (cosPsi<cos(fDPhi*0.5))
>> 1045 {
>> 1046 // Point lies outside phi range
>> 1047 if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)
>> 1048 {
>> 1049 safePhi=fabs(p.x()*sin(fSPhi)-p.y()*cos(fSPhi));
>> 1050 }
>> 1051 else
>> 1052 {
>> 1053 ePhi=fSPhi+fDPhi;
>> 1054 safePhi=fabs(p.x()*sin(ePhi)-p.y()*cos(ePhi));
>> 1055 }
>> 1056 if (safePhi>safe) safe=safePhi;
>> 1057 }
>> 1058 }
>> 1059 if (safe<0) safe=0;
>> 1060 return safe;
1121 } 1061 }
1122 1062
1123 ///////////////////////////////////////////// 1063 //////////////////////////////////////////////////////////////////////////////
1124 // 1064 //
1125 // Calculate distance to surface of shape fro 1065 // Calculate distance to surface of shape from `inside', allowing for tolerance
1126 // - Only Calc rmax intersection if no valid 1066 // - Only Calc rmax intersection if no valid rmin intersection
1127 1067
1128 G4double G4Tubs::DistanceToOut( const G4Three 1068 G4double G4Tubs::DistanceToOut( const G4ThreeVector& p,
1129 const G4Three 1069 const G4ThreeVector& v,
1130 const G4bool << 1070 const G4bool calcNorm,
1131 G4bool* << 1071 G4bool *validNorm,
1132 G4Three << 1072 G4ThreeVector *n ) const
1133 { 1073 {
1134 ESide side=kNull , sider=kNull, sidephi=kNu << 1074 ESide side = kNull ,sider,sidephi;
1135 G4double snxt, srd=kInfinity, sphi=kInfinit <<
1136 G4double deltaR, t1, t2, t3, b, c, d2, roMi <<
1137 <<
1138 // Vars for phi intersection: <<
1139 <<
1140 G4double pDistS, compS, pDistE, compE, sphi <<
1141 <<
1142 // Z plane intersection <<
1143 <<
1144 if (v.z() > 0 ) <<
1145 { <<
1146 pdist = fDz - p.z() ; <<
1147 if ( pdist > halfCarTolerance ) <<
1148 { <<
1149 snxt = pdist/v.z() ; <<
1150 side = kPZ ; <<
1151 } <<
1152 else <<
1153 { <<
1154 if (calcNorm) <<
1155 { <<
1156 *n = G4ThreeVector(0,0,1) ; <<
1157 *validNorm = true ; <<
1158 } <<
1159 return snxt = 0 ; <<
1160 } <<
1161 } <<
1162 else if ( v.z() < 0 ) <<
1163 { <<
1164 pdist = fDz + p.z() ; <<
1165 <<
1166 if ( pdist > halfCarTolerance ) <<
1167 { <<
1168 snxt = -pdist/v.z() ; <<
1169 side = kMZ ; <<
1170 } <<
1171 else <<
1172 { <<
1173 if (calcNorm) <<
1174 { <<
1175 *n = G4ThreeVector(0,0,-1) ; <<
1176 *validNorm = true ; <<
1177 } <<
1178 return snxt = 0.0 ; <<
1179 } <<
1180 } <<
1181 else <<
1182 { <<
1183 snxt = kInfinity ; // Travel perpendic <<
1184 side = kNull; <<
1185 } <<
1186 1075
1187 // Radial Intersections << 1076 G4double snxt,sr,sphi,pdist;
1188 // <<
1189 // Find intersection with cylinders at rmax <<
1190 // Intersection point (xi,yi,zi) on line x= <<
1191 // <<
1192 // Intersects with x^2+y^2=R^2 <<
1193 // <<
1194 // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v <<
1195 // <<
1196 // t1 t2 <<
1197 <<
1198 t1 = 1.0 - v.z()*v.z() ; // since v <<
1199 t2 = p.x()*v.x() + p.y()*v.y() ; <<
1200 t3 = p.x()*p.x() + p.y()*p.y() ; <<
1201 1077
1202 if ( snxt > 10*(fDz+fRMax) ) { roi2 = 2*fR << 1078 G4double deltaR,t1,t2,t3,b,c,d2;
1203 else { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t <<
1204 1079
1205 if ( t1 > 0 ) // Check not parallel << 1080 // Vars for phi intersection:
1206 { <<
1207 // Calculate srd, r exit distance <<
1208 1081
1209 if ( (t2 >= 0.0) && (roi2 > fRMax*(fRMax << 1082 G4double sinSPhi,cosSPhi,ePhi,sinEPhi,cosEPhi;
1210 { << 1083 G4double cPhi,sinCPhi,cosCPhi;
1211 // Delta r not negative => leaving via << 1084 G4double pDistS,compS,pDistE,compE,sphi2,xi,yi,vphi;
1212 1085
1213 deltaR = t3 - fRMax*fRMax ; << 1086 // Z plane intersection
1214 1087
1215 // NOTE: Should use rho-fRMax<-kRadTole << 1088 if (v.z()>0)
1216 // - avoid sqrt for efficiency << 1089 {
>> 1090 pdist=fDz-p.z();
1217 1091
1218 if ( deltaR < -kRadTolerance*fRMax ) << 1092 if (pdist>kCarTolerance*0.5)
1219 { 1093 {
1220 b = t2/t1 ; << 1094 snxt=pdist/v.z();
1221 c = deltaR/t1 ; << 1095 side=kPZ;
1222 d2 = b*b-c; <<
1223 if( d2 >= 0 ) { srd = c/( -b - std::s <<
1224 else { srd = 0.; } <<
1225 sider = kRMax ; <<
1226 } 1096 }
1227 else 1097 else
1228 { 1098 {
1229 // On tolerant boundary & heading out << 1099 if (calcNorm)
1230 // outer radial surface -> leaving im << 1100 {
1231 << 1101 *n=G4ThreeVector(0,0,1);
1232 if ( calcNorm ) << 1102 *validNorm=true;
1233 { << 1103 }
1234 G4double invRho = FastInverseRxy( p << 1104 return snxt=0;
1235 *n = G4ThreeVector(p.x()*in << 1105 }
1236 *validNorm = true ; << 1106 }
1237 } << 1107 else if (v.z()<0)
1238 return snxt = 0 ; // Leaving by rmax << 1108 {
1239 } << 1109 pdist=fDz+p.z();
1240 } <<
1241 else if ( t2 < 0. ) // i.e. t2 < 0; Poss <<
1242 { <<
1243 roMin2 = t3 - t2*t2/t1 ; // min ro2 of <<
1244 1110
1245 if ( (fRMin != 0.0) && (roMin2 < fRMin* << 1111 if (pdist>kCarTolerance*0.5)
1246 { 1112 {
1247 deltaR = t3 - fRMin*fRMin ; << 1113 snxt=-pdist/v.z();
1248 b = t2/t1 ; << 1114 side=kMZ;
1249 c = deltaR/t1 ; <<
1250 d2 = b*b - c ; <<
1251 <<
1252 if ( d2 >= 0 ) // Leaving via rmin <<
1253 { <<
1254 // NOTE: SHould use rho-rmin>kRadTo <<
1255 // - avoid sqrt for efficiency <<
1256 <<
1257 if (deltaR > kRadTolerance*fRMin) <<
1258 { <<
1259 srd = c/(-b+std::sqrt(d2)); <<
1260 sider = kRMin ; <<
1261 } <<
1262 else <<
1263 { <<
1264 if ( calcNorm ) { <<
1265 *validNorm = false; <<
1266 } // Concave side <<
1267 return snxt = 0.0; <<
1268 } <<
1269 } <<
1270 else // No rmin intersect -> must <<
1271 { <<
1272 deltaR = t3 - fRMax*fRMax ; <<
1273 c = deltaR/t1 ; <<
1274 d2 = b*b-c; <<
1275 if( d2 >=0. ) <<
1276 { <<
1277 srd = -b + std::sqrt(d2) ; <<
1278 sider = kRMax ; <<
1279 } <<
1280 else // Case: On the border+t2<kRad <<
1281 // (v is perpendicular t <<
1282 { <<
1283 if (calcNorm) <<
1284 { <<
1285 G4double invRho = FastInverseRx <<
1286 *n = G4ThreeVector(p.x()*invRho <<
1287 *validNorm = true ; <<
1288 } <<
1289 return snxt = 0.0; <<
1290 } <<
1291 } <<
1292 } <<
1293 else if ( roi2 > fRMax*(fRMax + kRadTol <<
1294 // No rmin intersect -> must be rm <<
1295 { <<
1296 deltaR = t3 - fRMax*fRMax ; <<
1297 b = t2/t1 ; <<
1298 c = deltaR/t1; <<
1299 d2 = b*b-c; <<
1300 if( d2 >= 0 ) <<
1301 { <<
1302 srd = -b + std::sqrt(d2) ; <<
1303 sider = kRMax ; <<
1304 } <<
1305 else // Case: On the border+t2<kRadTo <<
1306 // (v is perpendicular to <<
1307 { <<
1308 if (calcNorm) <<
1309 { <<
1310 G4double invRho = FastInverseRxy( <<
1311 *n = G4ThreeVector(p.x()*invRho,p <<
1312 *validNorm = true ; <<
1313 } <<
1314 return snxt = 0.0; <<
1315 } <<
1316 } <<
1317 } <<
1318 <<
1319 // Phi Intersection <<
1320 <<
1321 if ( !fPhiFullTube ) <<
1322 { <<
1323 // add angle calculation with correctio <<
1324 // of the difference in domain of atan2 <<
1325 // <<
1326 vphi = std::atan2(v.y(),v.x()) ; <<
1327 <<
1328 if ( vphi < fSPhi - halfAngTolerance ) <<
1329 else if ( vphi > fSPhi + fDPhi + halfAn <<
1330 <<
1331 <<
1332 if ( (p.x() != 0.0) || (p.y() != 0.0) ) <<
1333 { <<
1334 // pDist -ve when inside <<
1335 <<
1336 pDistS = p.x()*sinSPhi - p.y()*cosSPh <<
1337 pDistE = -p.x()*sinEPhi + p.y()*cosEP <<
1338 <<
1339 // Comp -ve when in direction of outw <<
1340 <<
1341 compS = -sinSPhi*v.x() + cosSPhi*v.y( <<
1342 compE = sinEPhi*v.x() - cosEPhi*v.y( <<
1343 <<
1344 sidephi = kNull; <<
1345 <<
1346 if( ( (fDPhi <= pi) && ( (pDistS <= h <<
1347 && (pDistE <= h <<
1348 || ( (fDPhi > pi) && ((pDistS <= h <<
1349 || (pDistE <= <<
1350 { <<
1351 // Inside both phi *full* planes <<
1352 <<
1353 if ( compS < 0 ) <<
1354 { <<
1355 sphi = pDistS/compS ; <<
1356 <<
1357 if (sphi >= -halfCarTolerance) <<
1358 { <<
1359 xi = p.x() + sphi*v.x() ; <<
1360 yi = p.y() + sphi*v.y() ; <<
1361 <<
1362 // Check intersecting with corr <<
1363 // (if not -> no intersect) <<
1364 // <<
1365 if((std::fabs(xi)<=kCarToleranc <<
1366 { <<
1367 sidephi = kSPhi; <<
1368 if (((fSPhi-halfAngTolerance) <<
1369 &&((fSPhi+fDPhi+halfAngTol <<
1370 { <<
1371 sphi = kInfinity; <<
1372 } <<
1373 } <<
1374 else if ( yi*cosCPhi-xi*sinCPhi <<
1375 { <<
1376 sphi = kInfinity ; <<
1377 } <<
1378 else <<
1379 { <<
1380 sidephi = kSPhi ; <<
1381 if ( pDistS > -halfCarToleran <<
1382 { <<
1383 sphi = 0.0 ; // Leave by sp <<
1384 } <<
1385 } <<
1386 } <<
1387 else <<
1388 { <<
1389 sphi = kInfinity ; <<
1390 } <<
1391 } <<
1392 else <<
1393 { <<
1394 sphi = kInfinity ; <<
1395 } <<
1396 <<
1397 if ( compE < 0 ) <<
1398 { <<
1399 sphi2 = pDistE/compE ; <<
1400 <<
1401 // Only check further if < starti <<
1402 // <<
1403 if ( (sphi2 > -halfCarTolerance) <<
1404 { <<
1405 xi = p.x() + sphi2*v.x() ; <<
1406 yi = p.y() + sphi2*v.y() ; <<
1407 <<
1408 if((std::fabs(xi)<=kCarToleranc <<
1409 { <<
1410 // Leaving via ending phi <<
1411 // <<
1412 if( (fSPhi-halfAngTolerance > <<
1413 ||(fSPhi+fDPhi+halfAngTo <<
1414 { <<
1415 sidephi = kEPhi ; <<
1416 if ( pDistE <= -halfCarTole <<
1417 else <<
1418 } <<
1419 } <<
1420 else // Check intersecting w <<
1421 <<
1422 if ( (yi*cosCPhi-xi*sinCPhi) >= <<
1423 { <<
1424 // Leaving via ending phi <<
1425 // <<
1426 sidephi = kEPhi ; <<
1427 if ( pDistE <= -halfCarTolera <<
1428 else <<
1429 } <<
1430 } <<
1431 } <<
1432 } <<
1433 else <<
1434 { <<
1435 sphi = kInfinity ; <<
1436 } <<
1437 } 1115 }
1438 else 1116 else
1439 { 1117 {
1440 // On z axis + travel not || to z axi << 1118 if (calcNorm)
1441 // within phi of shape, Step limited << 1119 {
1442 << 1120 *n=G4ThreeVector(0,0,-1);
1443 if ( (fSPhi - halfAngTolerance <= vph << 1121 *validNorm=true;
1444 && (vphi <= fSPhi + fDPhi + halfAn << 1122 }
1445 { << 1123 return snxt=0;
1446 sphi = kInfinity ; << 1124 }
1447 } << 1125 }
1448 else << 1126 else
1449 { << 1127 {
1450 sidephi = kSPhi ; // arbitrary << 1128 snxt=kInfinity; // Travel perpendicular to z axis
1451 sphi = 0.0 ; << 1129 side=kNull;
1452 } << 1130 }
1453 } << 1131
1454 if (sphi < snxt) // Order intersecttio << 1132 // Radial Intersections
1455 { << 1133 //
1456 snxt = sphi ; << 1134 // Find intersction with cylinders at rmax/rmin
1457 side = sidephi ; << 1135 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1458 } << 1136 //
>> 1137 // Intersects with x^2+y^2=R^2
>> 1138 //
>> 1139 // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0
>> 1140 //
>> 1141 // t1 t2 t3
>> 1142
>> 1143 t1=1.0-v.z()*v.z(); // since v normalised
>> 1144 t2=p.x()*v.x()+p.y()*v.y();
>> 1145 t3=p.x()*p.x()+p.y()*p.y();
>> 1146
>> 1147 if (t1>0) // Check not parallel
>> 1148 {
>> 1149 // Calculate sr, r exit distance
>> 1150
>> 1151 if (t2>=0)
>> 1152 {
>> 1153 // Delta r not negative => leaving via rmax
>> 1154
>> 1155 deltaR=t3-fRMax*fRMax;
>> 1156
>> 1157 // NOTE: Should use rho-fRMax<-kRadTolerance*0.5 - avoid sqrt for efficiency
>> 1158
>> 1159 if (deltaR<-kRadTolerance*fRMax)
>> 1160 {
>> 1161 b=t2/t1;
>> 1162 c=deltaR/t1;
>> 1163 sr=-b+sqrt(b*b-c);
>> 1164 sider=kRMax;
>> 1165 }
>> 1166 else
>> 1167 {
>> 1168 // On tolerant boundary & heading outwards (or perpendicular to) outer
>> 1169 // radial surface -> leaving immediately
>> 1170
>> 1171 if (calcNorm)
>> 1172 {
>> 1173 // if ( p.x() || p.y() )
>> 1174 // {
>> 1175 // *n=G4ThreeVector(p.x(),p.y(),0);
>> 1176 // }
>> 1177 // else
>> 1178 // {
>> 1179 // *n=v;
>> 1180 // }
>> 1181 *n=G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0);
>> 1182 *validNorm=true;
>> 1183 }
>> 1184 return snxt=0; // Leaving by rmax immediately
>> 1185 }
>> 1186 }
>> 1187 else // i.e. t2 < 0
>> 1188 {
>> 1189 // Possible rmin intersection
>> 1190
>> 1191 if ( fRMin )
>> 1192 {
>> 1193 deltaR=t3-fRMin*fRMin;
>> 1194 b=t2/t1;
>> 1195 c=deltaR/t1;
>> 1196 d2=b*b-c;
>> 1197
>> 1198 if (d2>=0) // Leaving via rmin
>> 1199 {
>> 1200 // NOTE: SHould use rho-rmin>kRadTolerance*0.5 - avoid sqrt for efficiency
>> 1201
>> 1202 if (deltaR > kRadTolerance*fRMin)
>> 1203 {
>> 1204 sr = -b-sqrt(d2) ;
>> 1205 sider = kRMin ;
>> 1206 }
>> 1207 else
>> 1208 {
>> 1209 if (calcNorm)
>> 1210 {
>> 1211 *validNorm = false ; // Concave side
>> 1212 }
>> 1213 return snxt=0;
>> 1214 }
>> 1215 }
>> 1216 else // No rmin intersect -> must be rmax intersect
>> 1217 {
>> 1218 deltaR=t3-fRMax*fRMax;
>> 1219 c=deltaR/t1;
>> 1220 sr=-b+sqrt(b*b-c);
>> 1221 sider=kRMax;
>> 1222 }
>> 1223 }
>> 1224 else // No rmin intersect -> must be rmax intersect
>> 1225 {
>> 1226 deltaR=t3-fRMax*fRMax;
>> 1227 b=t2/t1;
>> 1228 c=deltaR/t1;
>> 1229 sr=-b+sqrt(b*b-c);
>> 1230 sider=kRMax;
>> 1231 }
>> 1232 }
>> 1233
>> 1234 // Phi Intersection
>> 1235
>> 1236
>> 1237 if (fDPhi<2.0*M_PI)
>> 1238 {
>> 1239 sinSPhi=sin(fSPhi);
>> 1240 cosSPhi=cos(fSPhi);
>> 1241 ePhi=fSPhi+fDPhi;
>> 1242 sinEPhi=sin(ePhi);
>> 1243 cosEPhi=cos(ePhi);
>> 1244 cPhi=fSPhi+fDPhi*0.5;
>> 1245 sinCPhi=sin(cPhi);
>> 1246 cosCPhi=cos(cPhi);
>> 1247
>> 1248
>> 1249 if (p.x()||p.y()) // Check if on z axis (rho not needed later)
>> 1250 {
>> 1251 // pDist -ve when inside
>> 1252
>> 1253 pDistS=p.x()*sinSPhi-p.y()*cosSPhi;
>> 1254 pDistE=-p.x()*sinEPhi+p.y()*cosEPhi;
>> 1255
>> 1256 // Comp -ve when in direction of outwards normal
>> 1257
>> 1258 compS=-sinSPhi*v.x()+cosSPhi*v.y();
>> 1259 compE=sinEPhi*v.x()-cosEPhi*v.y();
>> 1260 sidephi=kNull;
>> 1261
>> 1262 // if ( pDistS <= 0 && pDistE <= 0 )
>> 1263
>> 1264 if( ( fDPhi <= pi && ( pDistS <= 0.5*kCarTolerance &&
>> 1265 pDistE <= 0.5*kCarTolerance ) ) ||
>> 1266 ( fDPhi > pi && !( pDistS > 0.5*kCarTolerance &&
>> 1267 pDistE > 0.5*kCarTolerance ) ) )
>> 1268 {
>> 1269 // Inside both phi *full* planes
>> 1270
>> 1271 if (compS<0)
>> 1272 {
>> 1273 sphi=pDistS/compS;
>> 1274 xi=p.x()+sphi*v.x();
>> 1275 yi=p.y()+sphi*v.y();
>> 1276
>> 1277 // Check intersecting with correct half-plane (if not -> no intersect)
>> 1278
>> 1279 if ((yi*cosCPhi-xi*sinCPhi)>=0)
>> 1280 {
>> 1281 sphi=kInfinity;
>> 1282 }
>> 1283 else
>> 1284 {
>> 1285 sidephi=kSPhi;
>> 1286
>> 1287 if (pDistS>-kCarTolerance*0.5)
>> 1288 {
>> 1289 sphi=0; // Leave by sphi immediately
>> 1290 }
>> 1291 }
>> 1292 }
>> 1293 else
>> 1294 {
>> 1295 sphi=kInfinity ;
>> 1296 }
>> 1297 if (compE<0)
>> 1298 {
>> 1299 sphi2=pDistE/compE;
>> 1300
>> 1301 // Only check further if < starting phi intersection
>> 1302
>> 1303 if (sphi2 > -0.5*kCarTolerance && sphi2 < sphi )
>> 1304 {
>> 1305 xi=p.x()+sphi2*v.x();
>> 1306 yi=p.y()+sphi2*v.y();
>> 1307
>> 1308 // Check intersecting with correct half-plane
>> 1309
>> 1310 if ((yi*cosCPhi-xi*sinCPhi) >= 0)
>> 1311 {
>> 1312 // Leaving via ending phi
>> 1313 sidephi=kEPhi;
>> 1314
>> 1315 if (pDistE<=-kCarTolerance*0.5)
>> 1316 {
>> 1317 sphi=sphi2;
>> 1318 }
>> 1319 else
>> 1320 {
>> 1321 sphi=0;
>> 1322 }
>> 1323 }
>> 1324 }
>> 1325 }
>> 1326 }
>> 1327 else sphi = kInfinity ;
>> 1328
>> 1329 /* *******************************************
>> 1330
>> 1331 else if ( pDistS >= 0 && pDistE >= 0 )
>> 1332 {
>> 1333 // Outside both *full* phi planes
>> 1334 if (pDistS <= pDistE)
>> 1335 {
>> 1336 sidephi = kSPhi ;
>> 1337 }
>> 1338 else
>> 1339 {
>> 1340 sidephi = kEPhi ;
>> 1341 }
>> 1342
>> 1343 if (fDPhi>M_PI)
>> 1344 {
>> 1345 if (compS<0&&compE<0)
>> 1346 {
>> 1347 sphi=0;
>> 1348 }
>> 1349 else
>> 1350 {
>> 1351 sphi=kInfinity;
>> 1352 }
>> 1353 }
>> 1354 else
>> 1355 {
>> 1356 // if towards both >=0 then once inside (after error) will remain inside
>> 1357
>> 1358 if (compS>=0&&compE>=0)
>> 1359 {
>> 1360 sphi=kInfinity;
>> 1361 }
>> 1362 else
>> 1363 {
>> 1364 sphi=0;
>> 1365 }
>> 1366 }
>> 1367 }
>> 1368 else if ( pDistS > 0 && pDistE < 0 )
>> 1369 {
>> 1370 // Outside full starting plane, inside full ending plane
>> 1371
>> 1372 if (fDPhi>M_PI)
>> 1373 {
>> 1374 if (compE<0)
>> 1375 {
>> 1376 sphi=pDistE/compE;
>> 1377 xi=p.x()+sphi*v.x();
>> 1378 yi=p.y()+sphi*v.y();
>> 1379
>> 1380 // Check intersection in correct half-plane (if not -> not leaving phi extent)
>> 1381
>> 1382 if ( (yi*cosCPhi-xi*sinCPhi) <= 0 )
>> 1383 {
>> 1384 sphi=kInfinity;
>> 1385 }
>> 1386 else // Leaving via Ending phi
>> 1387 {
>> 1388 sidephi = kEPhi ;
>> 1389 if (pDistE>-kCarTolerance*0.5)
>> 1390 {
>> 1391 sphi=0;
>> 1392 }
>> 1393 }
>> 1394 }
>> 1395 else
>> 1396 {
>> 1397 sphi=kInfinity;
>> 1398 }
>> 1399 }
>> 1400 else
>> 1401 {
>> 1402 if (compS>=0)
>> 1403 {
>> 1404 if (compE<0)
>> 1405 {
>> 1406 sphi=pDistE/compE;
>> 1407 xi=p.x()+sphi*v.x();
>> 1408 yi=p.y()+sphi*v.y();
>> 1409
>> 1410 // Check intersection in correct half-plane (if not -> remain in extent)
>> 1411
>> 1412 if ( (yi*cosCPhi-xi*sinCPhi) <= 0 )
>> 1413 {
>> 1414 sphi=kInfinity;
>> 1415 }
>> 1416 else // otherwise leaving via Ending phi
>> 1417 {
>> 1418 sidephi=kEPhi;
>> 1419 }
>> 1420 }
>> 1421 else
>> 1422 {
>> 1423 sphi=kInfinity;
>> 1424 }
>> 1425 }
>> 1426 else // leaving immediately by starting phi
>> 1427 {
>> 1428 sidephi=kSPhi;
>> 1429 sphi=0;
>> 1430 }
>> 1431 }
>> 1432 }
>> 1433 else
>> 1434 {
>> 1435 // Must be pDistS < 0 && pDistE > 0
>> 1436 // Inside full starting plane, outside full ending plane
>> 1437
>> 1438 if (fDPhi>M_PI)
>> 1439 {
>> 1440 if (compS<0)
>> 1441 {
>> 1442 sphi=pDistS/compS;
>> 1443 xi=p.x()+sphi*v.x();
>> 1444 yi=p.y()+sphi*v.y();
>> 1445
>> 1446 // Check intersection in correct half-plane (if not -> not leaving phi extent)
>> 1447
>> 1448 if ((yi*cosCPhi-xi*sinCPhi)>=0)
>> 1449 {
>> 1450 sphi=kInfinity;
>> 1451 }
>> 1452 else // Leaving via Starting phi
>> 1453 {
>> 1454 sidephi = kSPhi ;
>> 1455 if ( pDistS >- kCarTolerance*0.5 )
>> 1456 {
>> 1457 sphi=0;
>> 1458 }
>> 1459 }
>> 1460 }
>> 1461 else
>> 1462 {
>> 1463 sphi=kInfinity;
>> 1464 }
>> 1465 }
>> 1466 else
>> 1467 {
>> 1468 if (compE>=0)
>> 1469 {
>> 1470 if (compS<0)
>> 1471 {
>> 1472 sphi=pDistS/compS;
>> 1473 xi=p.x()+sphi*v.x();
>> 1474 yi=p.y()+sphi*v.y();
>> 1475
>> 1476 // Check intersection in correct half-plane (if not -> remain in extent)
>> 1477
>> 1478 if ((yi*cosCPhi-xi*sinCPhi)>=0)
>> 1479 {
>> 1480 sphi=kInfinity;
>> 1481 }
>> 1482 else // otherwise leaving via Starting phi
>> 1483 {
>> 1484 sidephi=kSPhi;
>> 1485 }
>> 1486 }
>> 1487 else
>> 1488 {
>> 1489 sphi=kInfinity;
>> 1490 }
>> 1491 }
>> 1492 else // leaving immediately by ending
>> 1493 {
>> 1494 sidephi=kEPhi;
>> 1495 sphi=0;
>> 1496 }
>> 1497 }
>> 1498 }
>> 1499
>> 1500 ****************************** */
>> 1501
>> 1502 }
>> 1503 else
>> 1504 {
>> 1505 // On z axis + travel not || to z axis -> if phi of vector direction
>> 1506 // within phi of shape, Step limited by rmax, else Step =0
>> 1507
>> 1508 vphi=atan2(v.y(),v.x());
>> 1509
>> 1510 if (fSPhi<vphi&&vphi<fSPhi+fDPhi)
>> 1511 {
>> 1512 sphi=kInfinity;
>> 1513 }
>> 1514 else
>> 1515 {
>> 1516 sidephi = kSPhi ; // arbitrary
>> 1517 sphi=0;
>> 1518 }
>> 1519 }
>> 1520 if (sphi<snxt) // Order intersecttions
>> 1521 {
>> 1522 snxt=sphi;
>> 1523 side=sidephi;
>> 1524 }
>> 1525 }
>> 1526 if (sr<snxt) // Order intersections
>> 1527 {
>> 1528 snxt=sr;
>> 1529 side=sider;
>> 1530 }
>> 1531 }
>> 1532 if (calcNorm)
>> 1533 {
>> 1534 switch(side)
>> 1535 {
>> 1536 case kRMax:
>> 1537 // Note: returned vector not normalised
>> 1538 // (divide by fRMax for unit vector)
>> 1539 xi=p.x()+snxt*v.x();
>> 1540 yi=p.y()+snxt*v.y();
>> 1541 *n=G4ThreeVector(xi/fRMax,yi/fRMax,0);
>> 1542 *validNorm=true;
>> 1543 break;
>> 1544
>> 1545 case kRMin:
>> 1546 *validNorm=false; // Rmin is inconvex
>> 1547 break;
>> 1548
>> 1549 case kSPhi:
>> 1550 if (fDPhi<=M_PI)
>> 1551 {
>> 1552 *n=G4ThreeVector(sin(fSPhi),-cos(fSPhi),0);
>> 1553 *validNorm=true;
>> 1554 }
>> 1555 else
>> 1556 {
>> 1557 *validNorm=false;
>> 1558 }
>> 1559 break;
>> 1560
>> 1561 case kEPhi:
>> 1562 if (fDPhi<=M_PI)
>> 1563 {
>> 1564 *n=G4ThreeVector(-sin(fSPhi+fDPhi),cos(fSPhi+fDPhi),0);
>> 1565 *validNorm=true;
>> 1566 }
>> 1567 else
>> 1568 {
>> 1569 *validNorm=false;
>> 1570 }
>> 1571 break;
>> 1572
>> 1573 case kPZ:
>> 1574 *n=G4ThreeVector(0,0,1);
>> 1575 *validNorm=true;
>> 1576 break;
>> 1577
>> 1578 case kMZ:
>> 1579 *n=G4ThreeVector(0,0,-1);
>> 1580 *validNorm=true;
>> 1581 break;
>> 1582
>> 1583 default:
>> 1584 G4Exception("Invalid enum in G4Tubs::DistanceToOut");
>> 1585 break;
>> 1586 }
1459 } 1587 }
1460 if (srd < snxt) // Order intersections << 1588 return snxt;
1461 { <<
1462 snxt = srd ; <<
1463 side = sider ; <<
1464 } <<
1465 } <<
1466 if (calcNorm) <<
1467 { <<
1468 switch(side) <<
1469 { <<
1470 case kRMax: <<
1471 // Note: returned vector not normalis <<
1472 // (divide by fRMax for unit vector) <<
1473 // <<
1474 xi = p.x() + snxt*v.x() ; <<
1475 yi = p.y() + snxt*v.y() ; <<
1476 *n = G4ThreeVector(xi/fRMax,yi/fRMax, <<
1477 *validNorm = true ; <<
1478 break ; <<
1479 <<
1480 case kRMin: <<
1481 *validNorm = false ; // Rmin is inco <<
1482 break ; <<
1483 <<
1484 case kSPhi: <<
1485 if ( fDPhi <= pi ) <<
1486 { <<
1487 *n = G4ThreeVector(sinSPhi, <<
1488 *validNorm = true ; <<
1489 } <<
1490 else <<
1491 { <<
1492 *validNorm = false ; <<
1493 } <<
1494 break ; <<
1495 <<
1496 case kEPhi: <<
1497 if (fDPhi <= pi) <<
1498 { <<
1499 *n = G4ThreeVector(-sinEPhi,cosEPhi <<
1500 *validNorm = true ; <<
1501 } <<
1502 else <<
1503 { <<
1504 *validNorm = false ; <<
1505 } <<
1506 break ; <<
1507 <<
1508 case kPZ: <<
1509 *n = G4ThreeVector(0,0,1) ; <<
1510 *validNorm = true ; <<
1511 break ; <<
1512 <<
1513 case kMZ: <<
1514 *n = G4ThreeVector(0,0,-1) ; <<
1515 *validNorm = true ; <<
1516 break ; <<
1517 <<
1518 default: <<
1519 G4cout << G4endl ; <<
1520 DumpInfo(); <<
1521 std::ostringstream message; <<
1522 G4long oldprc = message.precision(16) <<
1523 message << "Undefined side for valid <<
1524 << G4endl <<
1525 << "Position:" << G4endl << <<
1526 << "p.x() = " << p.x()/mm < <<
1527 << "p.y() = " << p.y()/mm < <<
1528 << "p.z() = " << p.z()/mm < <<
1529 << "Direction:" << G4endl << <<
1530 << "v.x() = " << v.x() << G <<
1531 << "v.y() = " << v.y() << G <<
1532 << "v.z() = " << v.z() << G <<
1533 << "Proposed distance :" << G <<
1534 << "snxt = " << snxt/mm << <<
1535 message.precision(oldprc) ; <<
1536 G4Exception("G4Tubs::DistanceToOut(p, <<
1537 JustWarning, message); <<
1538 break ; <<
1539 } <<
1540 } <<
1541 if ( snxt<halfCarTolerance ) { snxt=0 ; } <<
1542 <<
1543 return snxt ; <<
1544 } 1589 }
1545 1590
1546 ///////////////////////////////////////////// 1591 //////////////////////////////////////////////////////////////////////////
1547 // 1592 //
1548 // Calculate distance (<=actual) to closest s << 1593 // Calcluate distance (<=actual) to closest surface of shape from inside
1549 1594
1550 G4double G4Tubs::DistanceToOut( const G4Three << 1595 G4double G4Tubs::DistanceToOut(const G4ThreeVector& p) const
1551 { 1596 {
1552 G4double safe=0.0, rho, safeR1, safeR2, saf << 1597 G4double safe,rho,safeR1,safeR2,safeZ;
1553 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) << 1598 G4double safePhi,phiC,cosPhiC,sinPhiC,ePhi;
1554 << 1599 rho=sqrt(p.x()*p.x()+p.y()*p.y());
1555 #ifdef G4CSGDEBUG << 1600
1556 if( Inside(p) == kOutside ) << 1601 if (fRMin)
1557 { << 1602 {
1558 G4long oldprc = G4cout.precision(16) ; << 1603 safeR1=rho-fRMin;
1559 G4cout << G4endl ; << 1604 safeR2=fRMax-rho;
1560 DumpInfo(); << 1605 if (safeR1<safeR2)
1561 G4cout << "Position:" << G4endl << G4end << 1606 {
1562 G4cout << "p.x() = " << p.x()/mm << " m << 1607 safe=safeR1;
1563 G4cout << "p.y() = " << p.y()/mm << " m << 1608 }
1564 G4cout << "p.z() = " << p.z()/mm << " m << 1609 else
1565 G4cout.precision(oldprc) ; << 1610 {
1566 G4Exception("G4Tubs::DistanceToOut(p)", " << 1611 safe=safeR2;
1567 JustWarning, "Point p is outs << 1612 }
1568 } << 1613 }
1569 #endif <<
1570 <<
1571 if ( fRMin != 0.0 ) <<
1572 { <<
1573 safeR1 = rho - fRMin ; <<
1574 safeR2 = fRMax - rho ; <<
1575 <<
1576 if ( safeR1 < safeR2 ) { safe = safeR1 ; <<
1577 else { safe = safeR2 ; <<
1578 } <<
1579 else <<
1580 { <<
1581 safe = fRMax - rho ; <<
1582 } <<
1583 safeZ = fDz - std::fabs(p.z()) ; <<
1584 <<
1585 if ( safeZ < safe ) { safe = safeZ ; } <<
1586 <<
1587 // Check if phi divided, Calc distances clo <<
1588 // <<
1589 if ( !fPhiFullTube ) <<
1590 { <<
1591 if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 ) <<
1592 { <<
1593 safePhi = -(p.x()*sinSPhi - p.y()*cosSP <<
1594 } <<
1595 else 1614 else
1596 { << 1615 {
1597 safePhi = (p.x()*sinEPhi - p.y()*cosEPh << 1616 safe=fRMax-rho;
1598 } << 1617 }
1599 if (safePhi < safe) { safe = safePhi ; } << 1618
1600 } << 1619 safeZ=fDz-fabs(p.z());
1601 if ( safe < 0 ) { safe = 0 ; } << 1620 if (safeZ<safe) safe=safeZ;
1602 << 1621
1603 return safe ; << 1622 // Check if phi divided, Calc distances closest phi plane
>> 1623 if (fDPhi<2.0*M_PI)
>> 1624 {
>> 1625 // Above/below central phi of Tubs?
>> 1626 phiC=fSPhi+fDPhi*0.5;
>> 1627 cosPhiC=cos(phiC);
>> 1628 sinPhiC=sin(phiC);
>> 1629 if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)
>> 1630 {
>> 1631 safePhi=-(p.x()*sin(fSPhi)-p.y()*cos(fSPhi));
>> 1632 }
>> 1633 else
>> 1634 {
>> 1635 ePhi=fSPhi+fDPhi;
>> 1636 safePhi=(p.x()*sin(ePhi)-p.y()*cos(ePhi));
>> 1637 }
>> 1638 if (safePhi<safe) safe=safePhi;
>> 1639 }
>> 1640 if (safe<0) safe=0;
>> 1641 return safe;
1604 } 1642 }
1605 1643
1606 ///////////////////////////////////////////// << 1644 /////////////////////////////////////////////////////////////////////////
1607 // 1645 //
1608 // Stream object contents to an output stream << 1646 // Create a List containing the transformed vertices
>> 1647 // Ordering [0-3] -fDz cross section
>> 1648 // [4-7] +fDz cross section such that [0] is below [4],
>> 1649 // [1] below [5] etc.
>> 1650 // Note:
>> 1651 // Caller has deletion resposibility
>> 1652 // Potential improvement: For last slice, use actual ending angle
>> 1653 // to avoid rounding error problems.
1609 1654
1610 G4GeometryType G4Tubs::GetEntityType() const << 1655 G4ThreeVectorList*
>> 1656 G4Tubs::CreateRotatedVertices(const G4AffineTransform& pTransform) const
1611 { 1657 {
1612 return {"G4Tubs"}; << 1658 G4ThreeVectorList *vertices;
>> 1659 G4ThreeVector vertex0,vertex1,vertex2,vertex3;
>> 1660 G4double meshAngle,meshRMax,crossAngle,cosCrossAngle,sinCrossAngle,sAngle;
>> 1661 G4double rMaxX,rMaxY,rMinX,rMinY;
>> 1662 G4int crossSection,noCrossSections;
>> 1663
>> 1664 // Compute no of cross-sections necessary to mesh tube
>> 1665 noCrossSections=G4int (fDPhi/kMeshAngleDefault)+1;
>> 1666 if (noCrossSections<kMinMeshSections)
>> 1667 {
>> 1668 noCrossSections=kMinMeshSections;
>> 1669 }
>> 1670 else if (noCrossSections>kMaxMeshSections)
>> 1671 {
>> 1672 noCrossSections=kMaxMeshSections;
>> 1673 }
>> 1674
>> 1675 meshAngle=fDPhi/(noCrossSections-1);
>> 1676 meshRMax=fRMax/cos(meshAngle*0.5);
>> 1677
>> 1678 // If complete in phi, set start angle such that mesh will be at fRMax
>> 1679 // on the x axis. Will give better extent calculations when not rotated.
>> 1680 if (fDPhi==M_PI*2.0&&fSPhi==0)
>> 1681 {
>> 1682 sAngle=-meshAngle*0.5;
>> 1683 }
>> 1684 else
>> 1685 {
>> 1686 sAngle=fSPhi;
>> 1687 }
>> 1688
>> 1689 vertices=new G4ThreeVectorList(noCrossSections*4);
>> 1690 if (vertices)
>> 1691 {
>> 1692 for (crossSection=0;crossSection<noCrossSections;crossSection++)
>> 1693 {
>> 1694 // Compute coordinates of cross section at section crossSection
>> 1695 crossAngle=sAngle+crossSection*meshAngle;
>> 1696 cosCrossAngle=cos(crossAngle);
>> 1697 sinCrossAngle=sin(crossAngle);
>> 1698
>> 1699 rMaxX=meshRMax*cosCrossAngle;
>> 1700 rMaxY=meshRMax*sinCrossAngle;
>> 1701 rMinX=fRMin*cosCrossAngle;
>> 1702 rMinY=fRMin*sinCrossAngle;
>> 1703 vertex0=G4ThreeVector(rMinX,rMinY,-fDz);
>> 1704 vertex1=G4ThreeVector(rMaxX,rMaxY,-fDz);
>> 1705 vertex2=G4ThreeVector(rMaxX,rMaxY,+fDz);
>> 1706 vertex3=G4ThreeVector(rMinX,rMinY,+fDz);
>> 1707
>> 1708 vertices->insert(pTransform.TransformPoint(vertex0));
>> 1709 vertices->insert(pTransform.TransformPoint(vertex1));
>> 1710 vertices->insert(pTransform.TransformPoint(vertex2));
>> 1711 vertices->insert(pTransform.TransformPoint(vertex3));
>> 1712 }
>> 1713 }
>> 1714 else
>> 1715 {
>> 1716 G4Exception("G4Tubs::CreateRotatedVertices Out of memory - Cannot alloc vertices");
>> 1717 }
>> 1718 return vertices;
1613 } 1719 }
1614 1720
1615 ///////////////////////////////////////////// << 1721 ///////////////////////////////////////////////////////////////////////////
1616 // <<
1617 // Make a clone of the object <<
1618 // 1722 //
1619 G4VSolid* G4Tubs::Clone() const << 1723 // Methods for visualisation
1620 { <<
1621 return new G4Tubs(*this); <<
1622 } <<
1623 1724
1624 ///////////////////////////////////////////// <<
1625 // <<
1626 // Stream object contents to an output stream <<
1627 1725
1628 std::ostream& G4Tubs::StreamInfo( std::ostrea << 1726 void G4Tubs::DescribeYourselfTo (G4VGraphicsScene& scene) const {
1629 { << 1727 scene.AddThis (*this);
1630 G4long oldprc = os.precision(16); << 1728 }
1631 os << "------------------------------------ <<
1632 << " *** Dump for solid - " << GetNam <<
1633 << " ================================ <<
1634 << " Solid type: G4Tubs\n" <<
1635 << " Parameters: \n" <<
1636 << " inner radius : " << fRMin/mm << <<
1637 << " outer radius : " << fRMax/mm << <<
1638 << " half length Z: " << fDz/mm << " <<
1639 << " starting phi : " << fSPhi/degree <<
1640 << " delta phi : " << fDPhi/degree <<
1641 << "------------------------------------ <<
1642 os.precision(oldprc); <<
1643 1729
1644 return os; << 1730 G4VisExtent G4Tubs::GetExtent() const {
>> 1731 // Define the sides of the box into which the G4Tubs instance would fit.
>> 1732 return G4VisExtent (-fRMax, fRMax, -fRMax, fRMax, -fDz, fDz);
1645 } 1733 }
1646 1734
1647 ///////////////////////////////////////////// << 1735 G4Polyhedron* G4Tubs::CreatePolyhedron () const {
1648 // << 1736 return new G4PolyhedronTubs (fRMin, fRMax, fDz, fSPhi, fDPhi);
1649 // GetPointOnSurface << 1737 }
1650 1738
1651 G4ThreeVector G4Tubs::GetPointOnSurface() con << 1739 G4NURBS* G4Tubs::CreateNURBS () const {
1652 { << 1740 G4NURBS* pNURBS;
1653 G4double Rmax = fRMax; << 1741 if (fRMin != 0) {
1654 G4double Rmin = fRMin; << 1742 if (fDPhi >= 2.0 * M_PI) {
1655 G4double hz = 2.*fDz; // height << 1743 pNURBS = new G4NURBStube (fRMin, fRMax, fDz);
1656 G4double lext = fDPhi*Rmax; // length of ex <<
1657 G4double lint = fDPhi*Rmin; // length of in <<
1658 <<
1659 // Set array of surface areas <<
1660 // <<
1661 G4double RRmax = Rmax * Rmax; <<
1662 G4double RRmin = Rmin * Rmin; <<
1663 G4double sbase = 0.5*fDPhi*(RRmax - RRmin); <<
1664 G4double scut = (fDPhi == twopi) ? 0. : hz* <<
1665 G4double ssurf[6] = { scut, scut, sbase, sb <<
1666 ssurf[1] += ssurf[0]; <<
1667 ssurf[2] += ssurf[1]; <<
1668 ssurf[3] += ssurf[2]; <<
1669 ssurf[4] += ssurf[3]; <<
1670 ssurf[5] += ssurf[4]; <<
1671 <<
1672 // Select surface <<
1673 // <<
1674 G4double select = ssurf[5]*G4QuickRand(); <<
1675 G4int k = 5; <<
1676 k -= (G4int)(select <= ssurf[4]); <<
1677 k -= (G4int)(select <= ssurf[3]); <<
1678 k -= (G4int)(select <= ssurf[2]); <<
1679 k -= (G4int)(select <= ssurf[1]); <<
1680 k -= (G4int)(select <= ssurf[0]); <<
1681 <<
1682 // Generate point on selected surface <<
1683 // <<
1684 switch(k) <<
1685 { <<
1686 case 0: // start phi cut <<
1687 { <<
1688 G4double r = Rmin + (Rmax - Rmin)*G4Qui <<
1689 return { r*cosSPhi, r*sinSPhi, hz*G4Qui <<
1690 } <<
1691 case 1: // end phi cut <<
1692 { <<
1693 G4double r = Rmin + (Rmax - Rmin)*G4Qui <<
1694 return { r*cosEPhi, r*sinEPhi, hz*G4Qui <<
1695 } 1744 }
1696 case 2: // base at -dz << 1745 else {
1697 { << 1746 pNURBS = new G4NURBStubesector (fRMin, fRMax, fDz, fSPhi, fSPhi + fDPhi);
1698 G4double r = std::sqrt(RRmin + (RRmax - <<
1699 G4double phi = fSPhi + fDPhi*G4QuickRan <<
1700 return { r*std::cos(phi), r*std::sin(ph <<
1701 } <<
1702 case 3: // base at +dz <<
1703 { <<
1704 G4double r = std::sqrt(RRmin + (RRmax - <<
1705 G4double phi = fSPhi + fDPhi*G4QuickRan <<
1706 return { r*std::cos(phi), r*std::sin(ph <<
1707 } 1747 }
1708 case 4: // external lateral surface << 1748 }
1709 { << 1749 else {
1710 G4double phi = fSPhi + fDPhi*G4QuickRan << 1750 if (fDPhi >= 2.0 * M_PI) {
1711 G4double z = hz*G4QuickRand() - fDz; << 1751 pNURBS = new G4NURBScylinder (fRMax, fDz);
1712 G4double x = Rmax*std::cos(phi); <<
1713 G4double y = Rmax*std::sin(phi); <<
1714 return { x,y,z }; <<
1715 } 1752 }
1716 case 5: // internal lateral surface << 1753 else {
1717 { << 1754 const G4double epsilon = 1.e-4; // Cylinder sector not yet available!
1718 G4double phi = fSPhi + fDPhi*G4QuickRan << 1755 pNURBS = new G4NURBStubesector (epsilon, fRMax, fDz,
1719 G4double z = hz*G4QuickRand() - fDz; << 1756 fSPhi, fSPhi + fDPhi);
1720 G4double x = Rmin*std::cos(phi); <<
1721 G4double y = Rmin*std::sin(phi); <<
1722 return { x,y,z }; <<
1723 } 1757 }
1724 } 1758 }
1725 return {0., 0., 0.}; << 1759 return pNURBS;
1726 } 1760 }
1727 1761
1728 ///////////////////////////////////////////// <<
1729 // 1762 //
1730 // Methods for visualisation << 1763 //
1731 << 1764 /////////////////////////////////// End of G4Tubs.cc ////////////////////////
1732 void G4Tubs::DescribeYourselfTo ( G4VGraphics <<
1733 { <<
1734 scene.AddSolid (*this) ; <<
1735 } <<
1736 <<
1737 G4Polyhedron* G4Tubs::CreatePolyhedron () con <<
1738 { <<
1739 return new G4PolyhedronTubs (fRMin, fRMax, <<
1740 } <<
1741 1765
1742 #endif <<
1743 1766