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