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