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