| 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 // G4CutTubs implementation <<
27 // 26 //
>> 27 // $Id: G4CutTubs.cc 69788 2013-05-15 12:06:57Z gcosmo $
>> 28 // GEANT4 tag $Name: not supported by cvs2svn $
>> 29 //
>> 30 //
>> 31 // class G4CutTubs
>> 32 //
>> 33 // History:
>> 34 //
>> 35 // 05.04.12 M.Kelsey - GetPointOnSurface() throw flat in sqrt(r)
28 // 01.06.11 T.Nikitina - Derived from G4Tubs 36 // 01.06.11 T.Nikitina - Derived from G4Tubs
29 // 30.10.16 E.Tcherniaev - reimplemented Calcu << 37 //
30 // removed CreateRotat << 38 /////////////////////////////////////////////////////////////////////////
31 // ------------------------------------------- <<
32 39
33 #include "G4CutTubs.hh" 40 #include "G4CutTubs.hh"
34 41
35 #if !defined(G4GEOM_USE_UCTUBS) <<
36 <<
37 #include "G4GeomTools.hh" <<
38 #include "G4VoxelLimits.hh" 42 #include "G4VoxelLimits.hh"
39 #include "G4AffineTransform.hh" 43 #include "G4AffineTransform.hh"
40 #include "G4GeometryTolerance.hh" 44 #include "G4GeometryTolerance.hh"
41 #include "G4BoundingEnvelope.hh" <<
42 45
43 #include "G4VPVParameterisation.hh" 46 #include "G4VPVParameterisation.hh"
44 #include "G4QuickRand.hh" <<
45 47
46 #include "G4VGraphicsScene.hh" << 48 #include "Randomize.hh"
47 #include "G4Polyhedron.hh" <<
48 49
49 #include "G4AutoLock.hh" << 50 #include "meshdefs.hh"
50 51
51 namespace << 52 #include "G4VGraphicsScene.hh"
52 { << 53 #include "G4Polyhedron.hh"
53 G4Mutex zminmaxMutex = G4MUTEX_INITIALIZER; << 54 #include "G4NURBS.hh"
54 } << 55 #include "G4NURBStube.hh"
>> 56 #include "G4NURBScylinder.hh"
>> 57 #include "G4NURBStubesector.hh"
55 58
56 using namespace CLHEP; 59 using namespace CLHEP;
57 60
58 ////////////////////////////////////////////// 61 /////////////////////////////////////////////////////////////////////////
59 // 62 //
60 // Constructor - check parameters, convert ang 63 // Constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
61 // - note if pdphi>2PI then reset 64 // - note if pdphi>2PI then reset to 2PI
62 65
63 G4CutTubs::G4CutTubs( const G4String &pName, 66 G4CutTubs::G4CutTubs( const G4String &pName,
64 G4double pRMin, G4double 67 G4double pRMin, G4double pRMax,
65 G4double pDz, 68 G4double pDz,
66 G4double pSPhi, G4double 69 G4double pSPhi, G4double pDPhi,
67 G4ThreeVector pLowNorm,G 70 G4ThreeVector pLowNorm,G4ThreeVector pHighNorm )
68 : G4CSGSolid(pName), fRMin(pRMin), fRMax(pRM << 71 : G4Tubs(pName, pRMin, pRMax, pDz, pSPhi, pDPhi),
69 fSPhi(0.), fDPhi(0.), fZMin(0.), fZMax(0.) << 72 fPhiFullCutTube(true)
70 { 73 {
71 kRadTolerance = G4GeometryTolerance::GetInst 74 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
72 kAngTolerance = G4GeometryTolerance::GetInst 75 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
73 76
74 halfCarTolerance = kCarTolerance*0.5; <<
75 halfRadTolerance = kRadTolerance*0.5; <<
76 halfAngTolerance = kAngTolerance*0.5; <<
77 <<
78 if (pDz<=0) // Check z-len <<
79 { <<
80 std::ostringstream message; <<
81 message << "Negative Z half-length (" << p <<
82 G4Exception("G4CutTubs::G4CutTubs()", "Geo <<
83 } <<
84 if ( (pRMin >= pRMax) || (pRMin < 0) ) // Ch <<
85 { <<
86 std::ostringstream message; <<
87 message << "Invalid values for radii in so <<
88 << G4endl <<
89 << " pRMin = " << pRMin << <<
90 G4Exception("G4CutTubs::G4CutTubs()", "Geo <<
91 } <<
92 <<
93 // Check angles <<
94 // <<
95 CheckPhiAngles(pSPhi, pDPhi); <<
96 <<
97 // Check on Cutted Planes Normals 77 // Check on Cutted Planes Normals
98 // If there is NO CUT, propose to use G4Tubs 78 // If there is NO CUT, propose to use G4Tubs instead
99 // 79 //
100 if ( ( pLowNorm.x() == 0.0) && ( pLowNorm.y( << 80 if(pDPhi<twopi) { fPhiFullCutTube=false; }
101 && ( pHighNorm.x() == 0.0) && (pHighNorm.y << 81
>> 82 if ( ( !pLowNorm.x()) && ( !pLowNorm.y())
>> 83 && ( !pHighNorm.x()) && (!pHighNorm.y()) )
102 { 84 {
103 std::ostringstream message; 85 std::ostringstream message;
104 message << "Inexisting Low/High Normal to 86 message << "Inexisting Low/High Normal to Z plane or Parallel to Z."
105 << G4endl 87 << G4endl
106 << "Normals to Z plane are " << pL << 88 << "Normals to Z plane are (" << pLowNorm <<" and "
107 << pHighNorm << " in solid: " << G << 89 << pHighNorm << ") in solid: " << GetName();
108 G4Exception("G4CutTubs::G4CutTubs()", "Geo 90 G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids1001",
109 JustWarning, message, "Should 91 JustWarning, message, "Should use G4Tubs!");
110 } 92 }
111 93
112 // If Normal is (0,0,0),means parallel to R, 94 // If Normal is (0,0,0),means parallel to R, give it value of (0,0,+/-1)
113 // << 95 //
114 if (pLowNorm.mag2() == 0.) { pLowNorm.setZ( 96 if (pLowNorm.mag2() == 0.) { pLowNorm.setZ(-1.); }
115 if (pHighNorm.mag2()== 0.) { pHighNorm.setZ 97 if (pHighNorm.mag2()== 0.) { pHighNorm.setZ(1.); }
116 98
117 // Given Normals to Cut Planes have to be an << 99 // Given Normals to Cut Planes have to be an unit vectors.
118 // Normalize if it is needed. 100 // Normalize if it is needed.
119 // 101 //
120 if (pLowNorm.mag2() != 1.) { pLowNorm = pL 102 if (pLowNorm.mag2() != 1.) { pLowNorm = pLowNorm.unit(); }
121 if (pHighNorm.mag2()!= 1.) { pHighNorm = pH 103 if (pHighNorm.mag2()!= 1.) { pHighNorm = pHighNorm.unit(); }
122 104
123 // Normals to cutted planes have to point ou 105 // Normals to cutted planes have to point outside Solid
124 // 106 //
125 if( (pLowNorm.mag2() != 0.) && (pHighNorm.ma 107 if( (pLowNorm.mag2() != 0.) && (pHighNorm.mag2()!= 0. ) )
126 { 108 {
127 if( ( pLowNorm.z()>= 0. ) || ( pHighNorm.z 109 if( ( pLowNorm.z()>= 0. ) || ( pHighNorm.z() <= 0.))
128 { 110 {
129 std::ostringstream message; 111 std::ostringstream message;
130 message << "Invalid Low or High Normal t 112 message << "Invalid Low or High Normal to Z plane; "
131 "has to point outside Solid." 113 "has to point outside Solid." << G4endl
132 << "Invalid Norm to Z plane (" < 114 << "Invalid Norm to Z plane (" << pLowNorm << " or "
133 << pHighNorm << ") in solid: " < 115 << pHighNorm << ") in solid: " << GetName();
134 G4Exception("G4CutTubs::G4CutTubs()", "G 116 G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids0002",
135 FatalException, message); 117 FatalException, message);
136 } 118 }
137 } 119 }
138 fLowNorm = pLowNorm; 120 fLowNorm = pLowNorm;
139 fHighNorm = pHighNorm; 121 fHighNorm = pHighNorm;
140 122
141 // Check intersection of cut planes, they MU << 123 // Check Intersection of Cutted planes. They MUST NOT Intersect
142 // each other inside the lateral surface <<
143 // 124 //
144 if(IsCrossingCutPlanes()) 125 if(IsCrossingCutPlanes())
145 { 126 {
146 std::ostringstream message; 127 std::ostringstream message;
147 message << "Invalid normals to Z plane in << 128 message << "Invalid Low or High Normal to Z plane; "
148 << "Cut planes are crossing inside << 129 << "Crossing Cutted Planes." << G4endl
149 << " Solid type: G4CutTubs\n" << 130 << "Invalid Norm to Z plane (" << pLowNorm << " and "
150 << " Parameters: \n" << 131 << pHighNorm << ") in solid: " << GetName();
151 << " inner radius : " << fRMin/ <<
152 << " outer radius : " << fRMax/ <<
153 << " half length Z: " << fDz/mm <<
154 << " starting phi : " << fSPhi/ <<
155 << " delta phi : " << fDPhi/ <<
156 << " low Norm : " << fLowNo <<
157 << " high Norm : " << fHighN <<
158 G4Exception("G4CutTubs::G4CutTubs()", "Geo 132 G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids0002",
159 FatalException, message); 133 FatalException, message);
160 } 134 }
161 } 135 }
162 136
163 ////////////////////////////////////////////// 137 ///////////////////////////////////////////////////////////////////////
164 // 138 //
165 // Fake default constructor - sets only member 139 // Fake default constructor - sets only member data and allocates memory
166 // for usage restri 140 // for usage restricted to object persistency.
167 // 141 //
168 G4CutTubs::G4CutTubs( __void__& a ) 142 G4CutTubs::G4CutTubs( __void__& a )
169 : G4CSGSolid(a) << 143 : G4Tubs(a),
>> 144 fLowNorm(G4ThreeVector()), fHighNorm(G4ThreeVector()),
>> 145 fPhiFullCutTube(false)
170 { 146 {
171 } 147 }
172 148
173 ////////////////////////////////////////////// 149 //////////////////////////////////////////////////////////////////////////
174 // 150 //
175 // Destructor 151 // Destructor
176 152
177 G4CutTubs::~G4CutTubs() = default; << 153 G4CutTubs::~G4CutTubs()
>> 154 {
>> 155 }
178 156
179 ////////////////////////////////////////////// 157 //////////////////////////////////////////////////////////////////////////
180 // 158 //
181 // Copy constructor 159 // Copy constructor
182 160
183 G4CutTubs::G4CutTubs(const G4CutTubs&) = defau << 161 G4CutTubs::G4CutTubs(const G4CutTubs& rhs)
>> 162 : G4Tubs(rhs),
>> 163 fLowNorm(rhs.fLowNorm), fHighNorm(rhs.fHighNorm),
>> 164 fPhiFullCutTube(rhs.fPhiFullCutTube)
>> 165 {
>> 166 }
184 167
185 ////////////////////////////////////////////// 168 //////////////////////////////////////////////////////////////////////////
186 // 169 //
187 // Assignment operator 170 // Assignment operator
188 171
189 G4CutTubs& G4CutTubs::operator = (const G4CutT << 172 G4CutTubs& G4CutTubs::operator = (const G4CutTubs& rhs)
190 { 173 {
191 // Check assignment to self 174 // Check assignment to self
192 // 175 //
193 if (this == &rhs) { return *this; } 176 if (this == &rhs) { return *this; }
194 177
195 // Copy base class data 178 // Copy base class data
196 // 179 //
197 G4CSGSolid::operator=(rhs); << 180 G4Tubs::operator=(rhs);
198 181
199 // Copy data 182 // Copy data
200 // 183 //
201 kRadTolerance = rhs.kRadTolerance; kAngTole <<
202 fRMin = rhs.fRMin; fRMax = rhs.fRMax; fDz = <<
203 fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi; <<
204 fZMin = rhs.fZMin; fZMax = rhs.fZMax; <<
205 sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPh <<
206 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = r <<
207 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPh <<
208 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPh <<
209 fPhiFullCutTube = rhs.fPhiFullCutTube; <<
210 halfCarTolerance = rhs.halfCarTolerance; <<
211 halfRadTolerance = rhs.halfRadTolerance; <<
212 halfAngTolerance = rhs.halfAngTolerance; <<
213 fLowNorm = rhs.fLowNorm; fHighNorm = rhs.fH 184 fLowNorm = rhs.fLowNorm; fHighNorm = rhs.fHighNorm;
>> 185 fPhiFullCutTube = rhs.fPhiFullCutTube;
214 186
215 return *this; 187 return *this;
216 } 188 }
217 189
218 ////////////////////////////////////////////// << 190 ////////////////////////////////////////////////////////////////////////
219 // 191 //
220 // Get volume << 192 // Calculate extent under transform and specified limit
221 193
222 G4double G4CutTubs::GetCubicVolume() << 194 G4bool G4CutTubs::CalculateExtent( const EAxis pAxis,
>> 195 const G4VoxelLimits& pVoxelLimit,
>> 196 const G4AffineTransform& pTransform,
>> 197 G4double& pMin,
>> 198 G4double& pMax ) const
223 { 199 {
224 constexpr G4int nphi = 200, nrho = 100; << 200 if ( (!pTransform.IsRotated()) && (fDPhi == twopi) && (fRMin == 0) )
225 <<
226 if (fCubicVolume == 0.) <<
227 { 201 {
228 // get parameters << 202 // Special case handling for unrotated solid tubes
229 G4double rmin = GetInnerRadius(); << 203 // Compute x/y/z mins and maxs fro bounding box respecting limits,
230 G4double rmax = GetOuterRadius(); << 204 // with early returns if outside limits. Then switch() on pAxis,
231 G4double dz = GetZHalfLength(); << 205 // and compute exact x and y limit for x/y case
232 G4double sphi = GetStartPhiAngle(); << 206
233 G4double dphi = GetDeltaPhiAngle(); << 207 G4double xoffset, xMin, xMax;
>> 208 G4double yoffset, yMin, yMax;
>> 209 G4double zoffset, zMin, zMax;
>> 210
>> 211 G4double diff1, diff2, maxDiff, newMin, newMax;
>> 212 G4double xoff1, xoff2, yoff1, yoff2, delta;
>> 213
>> 214 xoffset = pTransform.NetTranslation().x();
>> 215 xMin = xoffset - fRMax;
>> 216 xMax = xoffset + fRMax;
234 217
235 // calculate volume << 218 if (pVoxelLimit.IsXLimited())
236 G4double volume = dz*dphi*(rmax*rmax - rmi <<
237 if (dphi < twopi) // make recalculation <<
238 { 219 {
239 // set values for calculation of h - dis << 220 if ( (xMin > pVoxelLimit.GetMaxXExtent())
240 // opposite points on bases << 221 || (xMax < pVoxelLimit.GetMinXExtent()) )
241 G4ThreeVector nbot = GetLowNorm(); << 222 {
242 G4ThreeVector ntop = GetHighNorm(); << 223 return false;
243 G4double nx = nbot.x()/nbot.z() - ntop.x << 224 }
244 G4double ny = nbot.y()/nbot.z() - ntop.y << 225 else
>> 226 {
>> 227 if (xMin < pVoxelLimit.GetMinXExtent())
>> 228 {
>> 229 xMin = pVoxelLimit.GetMinXExtent();
>> 230 }
>> 231 if (xMax > pVoxelLimit.GetMaxXExtent())
>> 232 {
>> 233 xMax = pVoxelLimit.GetMaxXExtent();
>> 234 }
>> 235 }
>> 236 }
>> 237 yoffset = pTransform.NetTranslation().y();
>> 238 yMin = yoffset - fRMax;
>> 239 yMax = yoffset + fRMax;
245 240
246 // compute volume by integration << 241 if ( pVoxelLimit.IsYLimited() )
247 G4double delrho = (rmax - rmin)/nrho; << 242 {
248 G4double delphi = dphi/nphi; << 243 if ( (yMin > pVoxelLimit.GetMaxYExtent())
249 volume = 0.; << 244 || (yMax < pVoxelLimit.GetMinYExtent()) )
250 for (G4int irho=0; irho<nrho; ++irho) <<
251 { 245 {
252 G4double r1 = rmin + delrho*irho; << 246 return false;
253 G4double r2 = rmin + delrho*(irho + 1 << 247 }
254 G4double rho = 0.5*(r1 + r2); << 248 else
255 G4double sector = 0.5*delphi*(r2*r2 - << 249 {
256 for (G4int iphi=0; iphi<nphi; ++iphi) << 250 if (yMin < pVoxelLimit.GetMinYExtent())
257 { 251 {
258 G4double phi = sphi + delphi*(iphi + << 252 yMin = pVoxelLimit.GetMinYExtent();
259 G4double h = nx*rho*std::cos(phi) + << 253 }
260 volume += sector*h; << 254 if (yMax > pVoxelLimit.GetMaxYExtent())
>> 255 {
>> 256 yMax=pVoxelLimit.GetMaxYExtent();
261 } 257 }
262 } 258 }
263 } 259 }
264 fCubicVolume = volume; << 260 zoffset = pTransform.NetTranslation().z();
265 } << 261 GetMaxMinZ(zMin,zMax);
266 return fCubicVolume; << 262 zMin += zoffset;
267 } << 263 zMax += zoffset;
268 264
269 ////////////////////////////////////////////// << 265 if ( pVoxelLimit.IsZLimited() )
270 // << 266 {
271 // Get surface area << 267 if ( (zMin > pVoxelLimit.GetMaxZExtent())
>> 268 || (zMax < pVoxelLimit.GetMinZExtent()) )
>> 269 {
>> 270 return false;
>> 271 }
>> 272 else
>> 273 {
>> 274 if (zMin < pVoxelLimit.GetMinZExtent())
>> 275 {
>> 276 zMin = pVoxelLimit.GetMinZExtent();
>> 277 }
>> 278 if (zMax > pVoxelLimit.GetMaxZExtent())
>> 279 {
>> 280 zMax = pVoxelLimit.GetMaxZExtent();
>> 281 }
>> 282 }
>> 283 }
>> 284 switch ( pAxis ) // Known to cut cylinder
>> 285 {
>> 286 case kXAxis :
>> 287 {
>> 288 yoff1 = yoffset - yMin;
>> 289 yoff2 = yMax - yoffset;
272 290
273 G4double G4CutTubs::GetSurfaceArea() << 291 if ( (yoff1 >= 0) && (yoff2 >= 0) ) // Y limits cross max/min x
274 { << 292 { // => no change
275 constexpr G4int nphi = 400; << 293 pMin = xMin;
>> 294 pMax = xMax;
>> 295 }
>> 296 else
>> 297 {
>> 298 // Y limits don't cross max/min x => compute max delta x,
>> 299 // hence new mins/maxs
276 300
277 if (fSurfaceArea == 0.) << 301 delta = fRMax*fRMax - yoff1*yoff1;
278 { << 302 diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
279 // get parameters << 303 delta = fRMax*fRMax - yoff2*yoff2;
280 G4double rmin = GetInnerRadius(); << 304 diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
281 G4double rmax = GetOuterRadius(); << 305 maxDiff = (diff1 > diff2) ? diff1:diff2;
282 G4double dz = GetZHalfLength(); << 306 newMin = xoffset - maxDiff;
283 G4double sphi = GetStartPhiAngle(); << 307 newMax = xoffset + maxDiff;
284 G4double dphi = GetDeltaPhiAngle(); << 308 pMin = (newMin < xMin) ? xMin : newMin;
285 G4ThreeVector nbot = GetLowNorm(); << 309 pMax = (newMax > xMax) ? xMax : newMax;
286 G4ThreeVector ntop = GetHighNorm(); << 310 }
287 << 311 break;
288 // calculate lateral surface area << 312 }
289 G4double sinner = 2.*dz*dphi*rmin; << 313 case kYAxis :
290 G4double souter = 2.*dz*dphi*rmax; << 314 {
291 if (dphi < twopi) // make recalculation << 315 xoff1 = xoffset - xMin;
292 { << 316 xoff2 = xMax - xoffset;
293 // set values for calculation of h - dis <<
294 // opposite points on bases <<
295 G4double nx = nbot.x()/nbot.z() - ntop.x <<
296 G4double ny = nbot.y()/nbot.z() - ntop.y <<
297 <<
298 // compute lateral surface area by integ <<
299 G4double delphi = dphi/nphi; <<
300 sinner = 0.; <<
301 souter = 0.; <<
302 for (G4int iphi=0; iphi<nphi; ++iphi) <<
303 { <<
304 G4double phi = sphi + delphi*(iphi + 0 <<
305 G4double cosphi = std::cos(phi); <<
306 G4double sinphi = std::sin(phi); <<
307 sinner += rmin*(nx*cosphi + ny*sinphi) <<
308 souter += rmax*(nx*cosphi + ny*sinphi) <<
309 } <<
310 sinner *= delphi*rmin; <<
311 souter *= delphi*rmax; <<
312 } <<
313 // set surface area <<
314 G4double scut = (dphi == twopi) ? 0. : 2. <<
315 G4double szero = 0.5*dphi*(rmax*rmax - rmi <<
316 G4double slow = szero/std::abs(nbot.z()); <<
317 G4double shigh = szero/std::abs(ntop.z()); <<
318 fSurfaceArea = sinner + souter + 2.*scut + <<
319 } <<
320 return fSurfaceArea; <<
321 } <<
322 317
323 ////////////////////////////////////////////// << 318 if ( (xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y
324 // << 319 { // => no change
325 // Get bounding box << 320 pMin = yMin;
>> 321 pMax = yMax;
>> 322 }
>> 323 else
>> 324 {
>> 325 // X limits don't cross max/min y => compute max delta y,
>> 326 // hence new mins/maxs
326 327
327 void G4CutTubs::BoundingLimits(G4ThreeVector& << 328 delta = fRMax*fRMax - xoff1*xoff1;
328 { << 329 diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
329 G4double rmin = GetInnerRadius(); << 330 delta = fRMax*fRMax - xoff2*xoff2;
330 G4double rmax = GetOuterRadius(); << 331 diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
331 G4double dz = GetZHalfLength(); << 332 maxDiff = (diff1 > diff2) ? diff1 : diff2;
332 G4double dphi = GetDeltaPhiAngle(); << 333 newMin = yoffset - maxDiff;
333 << 334 newMax = yoffset + maxDiff;
334 G4double sinSphi = GetSinStartPhi(); << 335 pMin = (newMin < yMin) ? yMin : newMin;
335 G4double cosSphi = GetCosStartPhi(); << 336 pMax = (newMax > yMax) ? yMax : newMax;
336 G4double sinEphi = GetSinEndPhi(); << 337 }
337 G4double cosEphi = GetCosEndPhi(); << 338 break;
338 << 339 }
339 G4ThreeVector norm; << 340 case kZAxis:
340 G4double mag, topx, topy, dists, diste; << 341 {
341 G4bool iftop; << 342 pMin = zMin;
342 << 343 pMax = zMax;
343 // Find Zmin << 344 break;
344 // << 345 }
345 G4double zmin; << 346 default:
346 norm = GetLowNorm(); << 347 break;
347 mag = std::sqrt(norm.x()*norm.x() + norm.y( << 348 }
348 topx = (mag == 0) ? 0 : -rmax*norm.x()/mag; << 349 pMin -= kCarTolerance;
349 topy = (mag == 0) ? 0 : -rmax*norm.y()/mag; << 350 pMax += kCarTolerance;
350 dists = sinSphi*topx - cosSphi*topy; << 351 return true;
351 diste = -sinEphi*topx + cosEphi*topy; <<
352 if (dphi > pi) <<
353 { <<
354 iftop = true; <<
355 if (dists > 0 && diste > 0)iftop = false; <<
356 } <<
357 else <<
358 { <<
359 iftop = false; <<
360 if (dists <= 0 && diste <= 0) iftop = true <<
361 } <<
362 if (iftop) <<
363 { <<
364 zmin = -(norm.x()*topx + norm.y()*topy)/no <<
365 } <<
366 else <<
367 { <<
368 G4double z1 = -rmin*(norm.x()*cosSphi + no <<
369 G4double z2 = -rmin*(norm.x()*cosEphi + no <<
370 G4double z3 = -rmax*(norm.x()*cosSphi + no <<
371 G4double z4 = -rmax*(norm.x()*cosEphi + no <<
372 zmin = std::min(std::min(std::min(z1,z2),z <<
373 } <<
374 <<
375 // Find Zmax <<
376 // <<
377 G4double zmax; <<
378 norm = GetHighNorm(); <<
379 mag = std::sqrt(norm.x()*norm.x() + norm.y( <<
380 topx = (mag == 0) ? 0 : -rmax*norm.x()/mag; <<
381 topy = (mag == 0) ? 0 : -rmax*norm.y()/mag; <<
382 dists = sinSphi*topx - cosSphi*topy; <<
383 diste = -sinEphi*topx + cosEphi*topy; <<
384 if (dphi > pi) <<
385 { <<
386 iftop = true; <<
387 if (dists > 0 && diste > 0) iftop = false; <<
388 } <<
389 else <<
390 { <<
391 iftop = false; <<
392 if (dists <= 0 && diste <= 0) iftop = true <<
393 } <<
394 if (iftop) <<
395 { <<
396 zmax = -(norm.x()*topx + norm.y()*topy)/no <<
397 } <<
398 else <<
399 { <<
400 G4double z1 = -rmin*(norm.x()*cosSphi + no <<
401 G4double z2 = -rmin*(norm.x()*cosEphi + no <<
402 G4double z3 = -rmax*(norm.x()*cosSphi + no <<
403 G4double z4 = -rmax*(norm.x()*cosEphi + no <<
404 zmax = std::max(std::max(std::max(z1,z2),z <<
405 } <<
406 <<
407 // Find bounding box <<
408 // <<
409 if (dphi < twopi) <<
410 { <<
411 G4TwoVector vmin,vmax; <<
412 G4GeomTools::DiskExtent(rmin,rmax, <<
413 GetSinStartPhi(),G <<
414 GetSinEndPhi(),Get <<
415 vmin,vmax); <<
416 pMin.set(vmin.x(),vmin.y(), zmin); <<
417 pMax.set(vmax.x(),vmax.y(), zmax); <<
418 } <<
419 else <<
420 { <<
421 pMin.set(-rmax,-rmax, zmin); <<
422 pMax.set( rmax, rmax, zmax); <<
423 } 352 }
424 << 353 else // Calculate rotated vertex coordinates
425 // Check correctness of the bounding box <<
426 // <<
427 if (pMin.x() >= pMax.x() || pMin.y() >= pMax <<
428 { 354 {
429 std::ostringstream message; << 355 G4int i, noEntries, noBetweenSections4;
430 message << "Bad bounding box (min >= max) << 356 G4bool existsAfterClip = false;
431 << GetName() << " !" << 357 G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform);
432 << "\npMin = " << pMin <<
433 << "\npMax = " << pMax; <<
434 G4Exception("G4CutTubs::BoundingLimits()", <<
435 JustWarning, message); <<
436 DumpInfo(); <<
437 } <<
438 } <<
439 <<
440 ////////////////////////////////////////////// <<
441 // <<
442 // Calculate extent under transform and specif <<
443 <<
444 G4bool G4CutTubs::CalculateExtent( const EAxis <<
445 const G4Vox <<
446 const G4Aff <<
447 G4dou <<
448 G4dou <<
449 { <<
450 G4ThreeVector bmin, bmax; <<
451 G4bool exist; <<
452 358
453 // Get bounding box << 359 pMin = kInfinity;
454 BoundingLimits(bmin,bmax); << 360 pMax = -kInfinity;
455 361
456 // Check bounding box << 362 noEntries = vertices->size();
457 G4BoundingEnvelope bbox(bmin,bmax); << 363 noBetweenSections4 = noEntries - 4;
458 #ifdef G4BBOX_EXTENT << 364
459 return bbox.CalculateExtent(pAxis,pVoxelLimi << 365 for ( i = 0 ; i < noEntries ; i += 4 )
460 #endif << 366 {
461 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 367 ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax);
462 { << 368 }
463 return exist = pMin < pMax; << 369 for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
464 } << 370 {
465 << 371 ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax);
466 // Get parameters of the solid << 372 }
467 G4double rmin = GetInnerRadius(); << 373 if ( (pMin != kInfinity) || (pMax != -kInfinity) )
468 G4double rmax = GetOuterRadius(); << 374 {
469 G4double dphi = GetDeltaPhiAngle(); << 375 existsAfterClip = true;
470 G4double zmin = bmin.z(); << 376 pMin -= kCarTolerance; // Add 2*tolerance to avoid precision troubles
471 G4double zmax = bmax.z(); << 377 pMax += kCarTolerance;
472 << 378 }
473 // Find bounding envelope and calculate exte << 379 else
474 // << 380 {
475 const G4int NSTEPS = 24; // numbe << 381 // Check for case where completely enveloping clipping volume
476 G4double astep = twopi/NSTEPS; // max a << 382 // If point inside then we are confident that the solid completely
477 G4int ksteps = (dphi <= astep) ? 1 : (G4i << 383 // envelopes the clipping volume. Hence set min/max extents according
478 G4double ang = dphi/ksteps; << 384 // to clipping volume extents along the specified axis.
479 << 385
480 G4double sinHalf = std::sin(0.5*ang); << 386 G4ThreeVector clipCentre(
481 G4double cosHalf = std::cos(0.5*ang); << 387 (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
482 G4double sinStep = 2.*sinHalf*cosHalf; << 388 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
483 G4double cosStep = 1. - 2.*sinHalf*sinHalf; << 389 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5 );
484 G4double rext = rmax/cosHalf; << 390
485 << 391 if ( Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside )
486 // bounding envelope for full cylinder consi << 392 {
487 // in other cases it is a sequence of quadri << 393 existsAfterClip = true;
488 if (rmin == 0 && dphi == twopi) << 394 pMin = pVoxelLimit.GetMinExtent(pAxis);
489 { << 395 pMax = pVoxelLimit.GetMaxExtent(pAxis);
490 G4double sinCur = sinHalf; << 396 }
491 G4double cosCur = cosHalf; << 397 }
492 << 398 delete vertices;
493 G4ThreeVectorList baseA(NSTEPS),baseB(NSTE << 399 return existsAfterClip;
494 for (G4int k=0; k<NSTEPS; ++k) <<
495 { <<
496 baseA[k].set(rext*cosCur,rext*sinCur,zmi <<
497 baseB[k].set(rext*cosCur,rext*sinCur,zma <<
498 <<
499 G4double sinTmp = sinCur; <<
500 sinCur = sinCur*cosStep + cosCur*sinStep <<
501 cosCur = cosCur*cosStep - sinTmp*sinStep <<
502 } <<
503 std::vector<const G4ThreeVectorList *> pol <<
504 polygons[0] = &baseA; <<
505 polygons[1] = &baseB; <<
506 G4BoundingEnvelope benv(bmin,bmax,polygons <<
507 exist = benv.CalculateExtent(pAxis,pVoxelL <<
508 } <<
509 else <<
510 { <<
511 G4double sinStart = GetSinStartPhi(); <<
512 G4double cosStart = GetCosStartPhi(); <<
513 G4double sinEnd = GetSinEndPhi(); <<
514 G4double cosEnd = GetCosEndPhi(); <<
515 G4double sinCur = sinStart*cosHalf + cos <<
516 G4double cosCur = cosStart*cosHalf - sin <<
517 <<
518 // set quadrilaterals <<
519 G4ThreeVectorList pols[NSTEPS+2]; <<
520 for (G4int k=0; k<ksteps+2; ++k) pols[k].r <<
521 pols[0][0].set(rmin*cosStart,rmin*sinStart <<
522 pols[0][1].set(rmin*cosStart,rmin*sinStart <<
523 pols[0][2].set(rmax*cosStart,rmax*sinStart <<
524 pols[0][3].set(rmax*cosStart,rmax*sinStart <<
525 for (G4int k=1; k<ksteps+1; ++k) <<
526 { <<
527 pols[k][0].set(rmin*cosCur,rmin*sinCur,z <<
528 pols[k][1].set(rmin*cosCur,rmin*sinCur,z <<
529 pols[k][2].set(rext*cosCur,rext*sinCur,z <<
530 pols[k][3].set(rext*cosCur,rext*sinCur,z <<
531 <<
532 G4double sinTmp = sinCur; <<
533 sinCur = sinCur*cosStep + cosCur*sinStep <<
534 cosCur = cosCur*cosStep - sinTmp*sinStep <<
535 } <<
536 pols[ksteps+1][0].set(rmin*cosEnd,rmin*sin <<
537 pols[ksteps+1][1].set(rmin*cosEnd,rmin*sin <<
538 pols[ksteps+1][2].set(rmax*cosEnd,rmax*sin <<
539 pols[ksteps+1][3].set(rmax*cosEnd,rmax*sin <<
540 <<
541 // set envelope and calculate extent <<
542 std::vector<const G4ThreeVectorList *> pol <<
543 polygons.resize(ksteps+2); <<
544 for (G4int k=0; k<ksteps+2; ++k) { polygon <<
545 G4BoundingEnvelope benv(bmin,bmax,polygons <<
546 exist = benv.CalculateExtent(pAxis,pVoxelL <<
547 } 400 }
548 return exist; <<
549 } 401 }
550 402
551 ////////////////////////////////////////////// << 403 ///////////////////////////////////////////////////////////////////////////
552 // 404 //
553 // Return whether point inside/outside/on surf 405 // Return whether point inside/outside/on surface
554 406
555 EInside G4CutTubs::Inside( const G4ThreeVector 407 EInside G4CutTubs::Inside( const G4ThreeVector& p ) const
556 { 408 {
557 G4ThreeVector vZ = G4ThreeVector(0,0,fDz); << 409 G4double zinLow,zinHigh,r2,pPhi=0.;
>> 410 G4double tolRMin,tolRMax;
>> 411 G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
558 EInside in = kInside; 412 EInside in = kInside;
>> 413 static const G4double halfCarTolerance=kCarTolerance*0.5;
>> 414 static const G4double halfRadTolerance=kRadTolerance*0.5;
>> 415 static const G4double halfAngTolerance=kAngTolerance*0.5;
>> 416
>> 417 // Check if point is contained in the cut plane in -/+ Z
559 418
560 // Check the lower cut plane 419 // Check the lower cut plane
561 // 420 //
562 G4double zinLow =(p+vZ).dot(fLowNorm); << 421 zinLow =(p+vZ).dot(fLowNorm);
563 if (zinLow > halfCarTolerance) { return kOu << 422 if (zinLow>halfCarTolerance) { return kOutside; }
564 423
565 // Check the higher cut plane 424 // Check the higher cut plane
566 // 425 //
567 G4double zinHigh = (p-vZ).dot(fHighNorm); << 426 zinHigh = (p-vZ).dot(fHighNorm);
568 if (zinHigh > halfCarTolerance) { return kO << 427 if (zinHigh>halfCarTolerance) { return kOutside; }
569 428
570 // Check radius 429 // Check radius
571 // 430 //
572 G4double r2 = p.x()*p.x() + p.y()*p.y() ; << 431 r2 = p.x()*p.x() + p.y()*p.y() ;
573 432
574 G4double tolRMin = fRMin - halfRadTolerance; << 433 // First check 'generous' boundaries R+tolerance
575 G4double tolRMax = fRMax + halfRadTolerance; << 434 //
>> 435 tolRMin = fRMin - halfRadTolerance ;
>> 436 tolRMax = fRMax + halfRadTolerance ;
576 if ( tolRMin < 0 ) { tolRMin = 0; } 437 if ( tolRMin < 0 ) { tolRMin = 0; }
577 << 438
578 if (r2 < tolRMin*tolRMin || r2 > tolRMax*tol << 439 if ( ((r2 < tolRMin*tolRMin) || (r2 > tolRMax*tolRMax))
579 << 440 && (r2 >=halfRadTolerance*halfRadTolerance) ) { return kOutside; }
580 // Check Phi cut << 441
>> 442 // Check Phi
581 // 443 //
582 if(!fPhiFullCutTube) 444 if(!fPhiFullCutTube)
583 { 445 {
584 if ((tolRMin == 0) && (std::fabs(p.x()) <= << 446 // Try outer tolerant phi boundaries only
585 && (std::fabs(p.y()) <= << 447
>> 448 if ( (tolRMin==0) && (std::fabs(p.x())<=halfCarTolerance)
>> 449 && (std::fabs(p.y())<=halfCarTolerance) )
586 { 450 {
587 return kSurface; 451 return kSurface;
588 } 452 }
>> 453
>> 454 pPhi = std::atan2(p.y(),p.x()) ;
589 455
590 G4double phi0 = std::atan2(p.y(),p.x()); << 456 if ( pPhi < -halfAngTolerance) { pPhi += twopi; } // 0<=pPhi<2pi
591 G4double phi1 = phi0 - twopi; << 457 if ( fSPhi >= 0 )
592 G4double phi2 = phi0 + twopi; << 458 {
593 << 459 if ( (std::fabs(pPhi) < halfAngTolerance)
594 in = kOutside; << 460 && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
595 G4double sphi = fSPhi - halfAngTolerance; << 461 {
596 G4double ephi = sphi + fDPhi + kAngToleran << 462 pPhi += twopi ; // 0 <= pPhi < 2pi
597 if ((phi0 >= sphi && phi0 <= ephi) || << 463 }
598 (phi1 >= sphi && phi1 <= ephi) || << 464 if ( (pPhi <= fSPhi - halfAngTolerance)
599 (phi2 >= sphi && phi2 <= ephi)) in = << 465 || (pPhi >= fSPhi + fDPhi + halfAngTolerance) )
600 if (in == kOutside) { return kOutside; } << 466 {
601 << 467 in = kOutside ;
602 sphi += kAngTolerance; << 468 }
603 ephi -= kAngTolerance; << 469 else if ( (pPhi <= fSPhi + halfAngTolerance)
604 if ((phi0 >= sphi && phi0 <= ephi) || << 470 || (pPhi >= fSPhi + fDPhi - halfAngTolerance) )
605 (phi1 >= sphi && phi1 <= ephi) || << 471 {
606 (phi2 >= sphi && phi2 <= ephi)) in = << 472 in=kSurface;
607 if (in == kSurface) { return kSurface; } << 473 }
>> 474 }
>> 475 else // fSPhi < 0
>> 476 {
>> 477 if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
>> 478 && (pPhi >= fSPhi + fDPhi + halfAngTolerance) )
>> 479 {
>> 480 in = kOutside;
>> 481 }
>> 482 else
>> 483 {
>> 484 in = kSurface ;
>> 485 }
>> 486 }
608 } 487 }
609 488
610 // Check on the Surface for Z 489 // Check on the Surface for Z
611 // 490 //
612 if ((zinLow >= -halfCarTolerance) || (zinHig << 491 if ((zinLow>=-halfCarTolerance)
>> 492 || (zinHigh>=-halfCarTolerance))
613 { 493 {
614 return kSurface; << 494 in=kSurface;
615 } 495 }
616 496
617 // Check on the Surface for R 497 // Check on the Surface for R
618 // 498 //
619 if (fRMin != 0.0) { tolRMin = fRMin + halfRa << 499 if (fRMin) { tolRMin = fRMin + halfRadTolerance ; }
620 else { tolRMin = 0; } << 500 else { tolRMin = 0 ; }
621 tolRMax = fRMax - halfRadTolerance; << 501 tolRMax = fRMax - halfRadTolerance ;
622 if (((r2 <= tolRMin*tolRMin) || (r2 >= tolRM << 502 if ( ((r2 <= tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax))&&
623 (r2 >= halfRadTolerance*halfRadToleranc << 503 (r2 >=halfRadTolerance*halfRadTolerance) )
624 { 504 {
625 return kSurface; 505 return kSurface;
626 } 506 }
627 507
628 return in; 508 return in;
629 } 509 }
630 510
631 ////////////////////////////////////////////// 511 ///////////////////////////////////////////////////////////////////////////
632 // 512 //
633 // Return unit normal of surface closest to p 513 // Return unit normal of surface closest to p
634 // - note if point on z axis, ignore phi divid 514 // - note if point on z axis, ignore phi divided sides
635 // - unsafe if point close to z axis a rmin=0 515 // - unsafe if point close to z axis a rmin=0 - no explicit checks
636 516
637 G4ThreeVector G4CutTubs::SurfaceNormal( const 517 G4ThreeVector G4CutTubs::SurfaceNormal( const G4ThreeVector& p ) const
638 { 518 {
639 G4int noSurfaces = 0; 519 G4int noSurfaces = 0;
640 G4double rho, pPhi; 520 G4double rho, pPhi;
641 G4double distZLow,distZHigh, distRMin, distR 521 G4double distZLow,distZHigh, distRMin, distRMax;
642 G4double distSPhi = kInfinity, distEPhi = kI 522 G4double distSPhi = kInfinity, distEPhi = kInfinity;
643 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 523 G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
644 524
>> 525 static const G4double halfCarTolerance = 0.5*kCarTolerance;
>> 526 static const G4double halfAngTolerance = 0.5*kAngTolerance;
>> 527
645 G4ThreeVector norm, sumnorm(0.,0.,0.); 528 G4ThreeVector norm, sumnorm(0.,0.,0.);
646 G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0); 529 G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0);
647 G4ThreeVector nR, nPs, nPe; 530 G4ThreeVector nR, nPs, nPe;
648 531
649 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()); 532 rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
650 533
651 distRMin = std::fabs(rho - fRMin); 534 distRMin = std::fabs(rho - fRMin);
652 distRMax = std::fabs(rho - fRMax); 535 distRMax = std::fabs(rho - fRMax);
653 536
654 // dist to Low Cut 537 // dist to Low Cut
655 // 538 //
656 distZLow =std::fabs((p+vZ).dot(fLowNorm)); 539 distZLow =std::fabs((p+vZ).dot(fLowNorm));
657 << 540
658 // dist to High Cut 541 // dist to High Cut
659 // 542 //
660 distZHigh = std::fabs((p-vZ).dot(fHighNorm)) 543 distZHigh = std::fabs((p-vZ).dot(fHighNorm));
661 544
662 if (!fPhiFullCutTube) // Protected agains << 545 if (!fPhiFullCutTube) // Protected against (0,0,z)
663 { 546 {
664 if ( rho > halfCarTolerance ) 547 if ( rho > halfCarTolerance )
665 { 548 {
666 pPhi = std::atan2(p.y(),p.x()); 549 pPhi = std::atan2(p.y(),p.x());
667 << 550
668 if(pPhi < fSPhi- halfCarTolerance) 551 if(pPhi < fSPhi- halfCarTolerance) { pPhi += twopi; }
669 else if(pPhi > fSPhi+fDPhi+ halfCarToler 552 else if(pPhi > fSPhi+fDPhi+ halfCarTolerance) { pPhi -= twopi; }
670 553
671 distSPhi = std::fabs(pPhi - fSPhi); << 554 distSPhi = std::fabs(pPhi - fSPhi);
672 distEPhi = std::fabs(pPhi - fSPhi - fDPh << 555 distEPhi = std::fabs(pPhi - fSPhi - fDPhi);
673 } 556 }
674 else if( fRMin == 0.0 ) << 557 else if( !fRMin )
675 { 558 {
676 distSPhi = 0.; << 559 distSPhi = 0.;
677 distEPhi = 0.; << 560 distEPhi = 0.;
678 } 561 }
679 nPs = G4ThreeVector( sinSPhi, -cosSPhi, 0 << 562 nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
680 nPe = G4ThreeVector( -sinEPhi, cosEPhi, 0 << 563 nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
681 } 564 }
682 if ( rho > halfCarTolerance ) { nR = G4Three 565 if ( rho > halfCarTolerance ) { nR = G4ThreeVector(p.x()/rho,p.y()/rho,0); }
683 566
684 if( distRMax <= halfCarTolerance ) << 567 if( distRMax <= halfCarTolerance )
685 { 568 {
686 ++noSurfaces; << 569 noSurfaces ++;
687 sumnorm += nR; 570 sumnorm += nR;
688 } 571 }
689 if( (fRMin != 0.0) && (distRMin <= halfCarTo << 572 if( fRMin && (distRMin <= halfCarTolerance) )
690 { 573 {
691 ++noSurfaces; << 574 noSurfaces ++;
692 sumnorm -= nR; 575 sumnorm -= nR;
693 } 576 }
694 if( fDPhi < twopi ) << 577 if( fDPhi < twopi )
695 { 578 {
696 if (distSPhi <= halfAngTolerance) << 579 if (distSPhi <= halfAngTolerance)
697 { 580 {
698 ++noSurfaces; << 581 noSurfaces ++;
699 sumnorm += nPs; 582 sumnorm += nPs;
700 } 583 }
701 if (distEPhi <= halfAngTolerance) << 584 if (distEPhi <= halfAngTolerance)
702 { 585 {
703 ++noSurfaces; << 586 noSurfaces ++;
704 sumnorm += nPe; 587 sumnorm += nPe;
705 } 588 }
706 } 589 }
707 if (distZLow <= halfCarTolerance) << 590 if (distZLow <= halfCarTolerance)
708 { 591 {
709 ++noSurfaces; << 592 noSurfaces ++;
710 sumnorm += fLowNorm; 593 sumnorm += fLowNorm;
711 } 594 }
712 if (distZHigh <= halfCarTolerance) << 595 if (distZHigh <= halfCarTolerance)
713 { 596 {
714 ++noSurfaces; << 597 noSurfaces ++;
715 sumnorm += fHighNorm; 598 sumnorm += fHighNorm;
716 } 599 }
717 if ( noSurfaces == 0 ) 600 if ( noSurfaces == 0 )
718 { 601 {
719 #ifdef G4CSGDEBUG 602 #ifdef G4CSGDEBUG
720 G4Exception("G4CutTubs::SurfaceNormal(p)", 603 G4Exception("G4CutTubs::SurfaceNormal(p)", "GeomSolids1002",
721 JustWarning, "Point p is not o 604 JustWarning, "Point p is not on surface !?" );
722 G4int oldprc = G4cout.precision(20); 605 G4int oldprc = G4cout.precision(20);
723 G4cout<< "G4CutTubs::SN ( "<<p.x()<<", "<< 606 G4cout<< "G4CutTubs::SN ( "<<p.x()<<", "<<p.y()<<", "<<p.z()<<" ); "
724 << G4endl << G4endl; 607 << G4endl << G4endl;
725 G4cout.precision(oldprc) ; 608 G4cout.precision(oldprc) ;
726 #endif << 609 #endif
727 norm = ApproxSurfaceNormal(p); 610 norm = ApproxSurfaceNormal(p);
728 } 611 }
729 else if ( noSurfaces == 1 ) { norm = sumnor 612 else if ( noSurfaces == 1 ) { norm = sumnorm; }
730 else { norm = sumnor 613 else { norm = sumnorm.unit(); }
731 614
732 return norm; 615 return norm;
733 } 616 }
734 617
735 ////////////////////////////////////////////// 618 /////////////////////////////////////////////////////////////////////////////
736 // 619 //
737 // Algorithm for SurfaceNormal() following the 620 // Algorithm for SurfaceNormal() following the original specification
738 // for points not on the surface 621 // for points not on the surface
739 622
740 G4ThreeVector G4CutTubs::ApproxSurfaceNormal( 623 G4ThreeVector G4CutTubs::ApproxSurfaceNormal( const G4ThreeVector& p ) const
741 { 624 {
742 enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ} <<
743 <<
744 ENorm side ; 625 ENorm side ;
745 G4ThreeVector norm ; 626 G4ThreeVector norm ;
746 G4double rho, phi ; 627 G4double rho, phi ;
747 G4double distZLow,distZHigh,distZ; 628 G4double distZLow,distZHigh,distZ;
748 G4double distRMin, distRMax, distSPhi, distE 629 G4double distRMin, distRMax, distSPhi, distEPhi, distMin ;
749 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 630 G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
750 631
751 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; 632 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
752 633
753 distRMin = std::fabs(rho - fRMin) ; 634 distRMin = std::fabs(rho - fRMin) ;
754 distRMax = std::fabs(rho - fRMax) ; 635 distRMax = std::fabs(rho - fRMax) ;
755 636
756 //dist to Low Cut 637 //dist to Low Cut
757 // 638 //
758 distZLow =std::fabs((p+vZ).dot(fLowNorm)); 639 distZLow =std::fabs((p+vZ).dot(fLowNorm));
759 640
760 //dist to High Cut 641 //dist to High Cut
761 // 642 //
762 distZHigh = std::fabs((p-vZ).dot(fHighNorm)) 643 distZHigh = std::fabs((p-vZ).dot(fHighNorm));
763 distZ=std::min(distZLow,distZHigh); 644 distZ=std::min(distZLow,distZHigh);
764 645
765 if (distRMin < distRMax) // First minimum 646 if (distRMin < distRMax) // First minimum
766 { 647 {
767 if ( distZ < distRMin ) 648 if ( distZ < distRMin )
768 { 649 {
769 distMin = distZ ; 650 distMin = distZ ;
770 side = kNZ ; 651 side = kNZ ;
771 } 652 }
772 else 653 else
773 { 654 {
774 distMin = distRMin ; 655 distMin = distRMin ;
775 side = kNRMin ; 656 side = kNRMin ;
776 } 657 }
777 } 658 }
778 else 659 else
779 { 660 {
780 if ( distZ < distRMax ) 661 if ( distZ < distRMax )
781 { 662 {
782 distMin = distZ ; 663 distMin = distZ ;
783 side = kNZ ; 664 side = kNZ ;
784 } 665 }
785 else 666 else
786 { 667 {
787 distMin = distRMax ; 668 distMin = distRMax ;
788 side = kNRMax ; 669 side = kNRMax ;
789 } 670 }
790 } << 671 }
791 if (!fPhiFullCutTube && (rho != 0.0) ) // << 672 if (!fPhiFullCutTube && rho ) // Protected against (0,0,z)
792 { 673 {
793 phi = std::atan2(p.y(),p.x()) ; 674 phi = std::atan2(p.y(),p.x()) ;
794 675
795 if ( phi < 0 ) { phi += twopi; } 676 if ( phi < 0 ) { phi += twopi; }
796 677
797 if ( fSPhi < 0 ) 678 if ( fSPhi < 0 )
798 { 679 {
799 distSPhi = std::fabs(phi - (fSPhi + twop 680 distSPhi = std::fabs(phi - (fSPhi + twopi))*rho ;
800 } 681 }
801 else 682 else
802 { 683 {
803 distSPhi = std::fabs(phi - fSPhi)*rho ; 684 distSPhi = std::fabs(phi - fSPhi)*rho ;
804 } 685 }
805 distEPhi = std::fabs(phi - fSPhi - fDPhi)* 686 distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
806 << 687
807 if (distSPhi < distEPhi) // Find new minim 688 if (distSPhi < distEPhi) // Find new minimum
808 { 689 {
809 if ( distSPhi < distMin ) 690 if ( distSPhi < distMin )
810 { 691 {
811 side = kNSPhi ; 692 side = kNSPhi ;
812 } 693 }
813 } 694 }
814 else 695 else
815 { 696 {
816 if ( distEPhi < distMin ) 697 if ( distEPhi < distMin )
817 { 698 {
818 side = kNEPhi ; 699 side = kNEPhi ;
819 } 700 }
820 } 701 }
821 } << 702 }
822 switch ( side ) 703 switch ( side )
823 { 704 {
824 case kNRMin : // Inner radius 705 case kNRMin : // Inner radius
825 { << 706 {
826 norm = G4ThreeVector(-p.x()/rho, -p.y()/ 707 norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, 0) ;
827 break ; 708 break ;
828 } 709 }
829 case kNRMax : // Outer radius 710 case kNRMax : // Outer radius
830 { << 711 {
831 norm = G4ThreeVector(p.x()/rho, p.y()/rh 712 norm = G4ThreeVector(p.x()/rho, p.y()/rho, 0) ;
832 break ; 713 break ;
833 } 714 }
834 case kNZ : // + or - dz 715 case kNZ : // + or - dz
835 { << 716 {
836 if ( distZHigh > distZLow ) { norm = fH 717 if ( distZHigh > distZLow ) { norm = fHighNorm ; }
837 else { norm = fL 718 else { norm = fLowNorm; }
838 break ; 719 break ;
839 } 720 }
840 case kNSPhi: 721 case kNSPhi:
841 { 722 {
842 norm = G4ThreeVector(sinSPhi, -cosSPhi, << 723 norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ;
843 break ; 724 break ;
844 } 725 }
845 case kNEPhi: 726 case kNEPhi:
846 { 727 {
847 norm = G4ThreeVector(-sinEPhi, cosEPhi, << 728 norm = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ;
848 break; 729 break;
849 } 730 }
850 default: // Should never reach this c 731 default: // Should never reach this case ...
851 { 732 {
852 DumpInfo(); 733 DumpInfo();
853 G4Exception("G4CutTubs::ApproxSurfaceNor 734 G4Exception("G4CutTubs::ApproxSurfaceNormal()",
854 "GeomSolids1002", JustWarnin 735 "GeomSolids1002", JustWarning,
855 "Undefined side for valid su 736 "Undefined side for valid surface normal to solid.");
856 break ; 737 break ;
857 } << 738 }
858 } << 739 }
859 return norm; 740 return norm;
860 } 741 }
861 742
862 ////////////////////////////////////////////// 743 ////////////////////////////////////////////////////////////////////
863 // 744 //
864 // 745 //
865 // Calculate distance to shape from outside, a 746 // Calculate distance to shape from outside, along normalised vector
866 // - return kInfinity if no intersection, or i 747 // - return kInfinity if no intersection, or intersection distance <= tolerance
867 // 748 //
868 // - Compute the intersection with the z plane << 749 // - Compute the intersection with the z planes
869 // - if at valid r, phi, return 750 // - if at valid r, phi, return
870 // 751 //
871 // -> If point is outer outer radius, compute 752 // -> If point is outer outer radius, compute intersection with rmax
872 // - if at valid phi,z return 753 // - if at valid phi,z return
873 // 754 //
874 // -> Compute intersection with inner radius, 755 // -> Compute intersection with inner radius, taking largest +ve root
875 // - if valid (in z,phi), save intersct 756 // - if valid (in z,phi), save intersction
876 // 757 //
877 // -> If phi segmented, compute intersectio 758 // -> If phi segmented, compute intersections with phi half planes
878 // - return smallest of valid phi inter 759 // - return smallest of valid phi intersections and
879 // inner radius intersection 760 // inner radius intersection
880 // 761 //
881 // NOTE: 762 // NOTE:
882 // - 'if valid' implies tolerant checking of i 763 // - 'if valid' implies tolerant checking of intersection points
883 764
884 G4double G4CutTubs::DistanceToIn( const G4Thre 765 G4double G4CutTubs::DistanceToIn( const G4ThreeVector& p,
885 const G4Thre 766 const G4ThreeVector& v ) const
886 { 767 {
887 G4double snxt = kInfinity ; // snxt = d 768 G4double snxt = kInfinity ; // snxt = default return value
888 G4double tolORMin2, tolIRMax2 ; // 'generou 769 G4double tolORMin2, tolIRMax2 ; // 'generous' radii squared
889 G4double tolORMax2, tolIRMin2; 770 G4double tolORMax2, tolIRMin2;
890 const G4double dRmax = 100.*fRMax; 771 const G4double dRmax = 100.*fRMax;
891 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 772 G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
>> 773
>> 774 static const G4double halfCarTolerance = 0.5*kCarTolerance;
>> 775 static const G4double halfRadTolerance = 0.5*kRadTolerance;
892 776
893 // Intersection point variables 777 // Intersection point variables
894 // 778 //
895 G4double Dist, sd=0, xi, yi, zi, rho2, inum, 779 G4double Dist, sd=0, xi, yi, zi, rho2, inum, iden, cosPsi, Comp,calf ;
896 G4double t1, t2, t3, b, c, d ; // Quadra << 780 G4double t1, t2, t3, b, c, d ; // Quadratic solver variables
897 G4double distZLow,distZHigh; 781 G4double distZLow,distZHigh;
898 // Calculate tolerant rmin and rmax 782 // Calculate tolerant rmin and rmax
899 783
900 if (fRMin > kRadTolerance) 784 if (fRMin > kRadTolerance)
901 { 785 {
902 tolORMin2 = (fRMin - halfRadTolerance)*(fR 786 tolORMin2 = (fRMin - halfRadTolerance)*(fRMin - halfRadTolerance) ;
903 tolIRMin2 = (fRMin + halfRadTolerance)*(fR 787 tolIRMin2 = (fRMin + halfRadTolerance)*(fRMin + halfRadTolerance) ;
904 } 788 }
905 else 789 else
906 { 790 {
907 tolORMin2 = 0.0 ; 791 tolORMin2 = 0.0 ;
908 tolIRMin2 = 0.0 ; 792 tolIRMin2 = 0.0 ;
909 } 793 }
910 tolORMax2 = (fRMax + halfRadTolerance)*(fRMa 794 tolORMax2 = (fRMax + halfRadTolerance)*(fRMax + halfRadTolerance) ;
911 tolIRMax2 = (fRMax - halfRadTolerance)*(fRMa 795 tolIRMax2 = (fRMax - halfRadTolerance)*(fRMax - halfRadTolerance) ;
912 796
913 // Intersection with ZCut surfaces 797 // Intersection with ZCut surfaces
914 798
915 // dist to Low Cut 799 // dist to Low Cut
916 // 800 //
917 distZLow =(p+vZ).dot(fLowNorm); 801 distZLow =(p+vZ).dot(fLowNorm);
918 802
919 // dist to High Cut 803 // dist to High Cut
920 // 804 //
921 distZHigh = (p-vZ).dot(fHighNorm); 805 distZHigh = (p-vZ).dot(fHighNorm);
922 806
923 if ( distZLow >= -halfCarTolerance ) 807 if ( distZLow >= -halfCarTolerance )
924 { 808 {
925 calf = v.dot(fLowNorm); 809 calf = v.dot(fLowNorm);
926 if (calf<0) 810 if (calf<0)
927 { 811 {
928 sd = -distZLow/calf; 812 sd = -distZLow/calf;
929 if(sd < 0.0) { sd = 0.0; } 813 if(sd < 0.0) { sd = 0.0; }
930 814
931 xi = p.x() + sd*v.x() ; 815 xi = p.x() + sd*v.x() ; // Intersection coords
932 yi = p.y() + sd*v.y() ; 816 yi = p.y() + sd*v.y() ;
933 rho2 = xi*xi + yi*yi ; 817 rho2 = xi*xi + yi*yi ;
934 818
935 // Check validity of intersection 819 // Check validity of intersection
936 820
937 if ((tolIRMin2 <= rho2) && (rho2 <= tolI 821 if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
938 { 822 {
939 if (!fPhiFullCutTube && (rho2 != 0.0)) << 823 if (!fPhiFullCutTube && rho2)
940 { 824 {
941 // Psi = angle made with central (av 825 // Psi = angle made with central (average) phi of shape
942 // 826 //
943 inum = xi*cosCPhi + yi*sinCPhi ; 827 inum = xi*cosCPhi + yi*sinCPhi ;
944 iden = std::sqrt(rho2) ; 828 iden = std::sqrt(rho2) ;
945 cosPsi = inum/iden ; 829 cosPsi = inum/iden ;
946 if (cosPsi >= cosHDPhiIT) { return 830 if (cosPsi >= cosHDPhiIT) { return sd ; }
947 } 831 }
948 else 832 else
949 { 833 {
950 return sd ; 834 return sd ;
951 } 835 }
952 } 836 }
953 } 837 }
954 else 838 else
955 { 839 {
956 if ( sd<halfCarTolerance ) 840 if ( sd<halfCarTolerance )
957 { 841 {
958 if(calf>=0) { sd=kInfinity; } 842 if(calf>=0) { sd=kInfinity; }
959 return sd ; // On/outside extent, and 843 return sd ; // On/outside extent, and heading away
960 } // -> cannot intersect 844 } // -> cannot intersect
961 } 845 }
962 } 846 }
963 847
964 if(distZHigh >= -halfCarTolerance ) 848 if(distZHigh >= -halfCarTolerance )
965 { 849 {
966 calf = v.dot(fHighNorm); 850 calf = v.dot(fHighNorm);
967 if (calf<0) 851 if (calf<0)
968 { 852 {
969 sd = -distZHigh/calf; 853 sd = -distZHigh/calf;
970 854
971 if(sd < 0.0) { sd = 0.0; } 855 if(sd < 0.0) { sd = 0.0; }
972 856
973 xi = p.x() + sd*v.x() ; 857 xi = p.x() + sd*v.x() ; // Intersection coords
974 yi = p.y() + sd*v.y() ; 858 yi = p.y() + sd*v.y() ;
975 rho2 = xi*xi + yi*yi ; 859 rho2 = xi*xi + yi*yi ;
976 860
977 // Check validity of intersection 861 // Check validity of intersection
978 862
979 if ((tolIRMin2 <= rho2) && (rho2 <= tolI 863 if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
980 { 864 {
981 if (!fPhiFullCutTube && (rho2 != 0.0)) << 865 if (!fPhiFullCutTube && rho2)
982 { 866 {
983 // Psi = angle made with central (av 867 // Psi = angle made with central (average) phi of shape
984 // 868 //
985 inum = xi*cosCPhi + yi*sinCPhi ; 869 inum = xi*cosCPhi + yi*sinCPhi ;
986 iden = std::sqrt(rho2) ; 870 iden = std::sqrt(rho2) ;
987 cosPsi = inum/iden ; 871 cosPsi = inum/iden ;
988 if (cosPsi >= cosHDPhiIT) { return 872 if (cosPsi >= cosHDPhiIT) { return sd ; }
989 } 873 }
990 else 874 else
991 { 875 {
992 return sd ; 876 return sd ;
993 } 877 }
994 } 878 }
995 } 879 }
996 else 880 else
997 { 881 {
998 if ( sd<halfCarTolerance ) 882 if ( sd<halfCarTolerance )
999 { << 883 {
1000 if(calf>=0) { sd=kInfinity; } 884 if(calf>=0) { sd=kInfinity; }
1001 return sd ; // On/outside extent, an 885 return sd ; // On/outside extent, and heading away
1002 } // -> cannot intersect 886 } // -> cannot intersect
1003 } 887 }
1004 } 888 }
1005 889
1006 // -> Can not intersect z surfaces 890 // -> Can not intersect z surfaces
1007 // 891 //
1008 // Intersection with rmax (possible return) 892 // Intersection with rmax (possible return) and rmin (must also check phi)
1009 // 893 //
1010 // Intersection point (xi,yi,zi) on line x= 894 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1011 // 895 //
1012 // Intersects with x^2+y^2=R^2 896 // Intersects with x^2+y^2=R^2
1013 // 897 //
1014 // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v 898 // 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
1015 // t1 t2 899 // t1 t2 t3
1016 900
1017 t1 = 1.0 - v.z()*v.z() ; 901 t1 = 1.0 - v.z()*v.z() ;
1018 t2 = p.x()*v.x() + p.y()*v.y() ; 902 t2 = p.x()*v.x() + p.y()*v.y() ;
1019 t3 = p.x()*p.x() + p.y()*p.y() ; 903 t3 = p.x()*p.x() + p.y()*p.y() ;
1020 if ( t1 > 0 ) // Check not || to z a 904 if ( t1 > 0 ) // Check not || to z axis
1021 { 905 {
1022 b = t2/t1 ; 906 b = t2/t1 ;
1023 c = t3 - fRMax*fRMax ; 907 c = t3 - fRMax*fRMax ;
1024 << 908
1025 if ((t3 >= tolORMax2) && (t2<0)) // Thi 909 if ((t3 >= tolORMax2) && (t2<0)) // This also handles the tangent case
1026 { 910 {
1027 // Try outer cylinder intersection, c=( 911 // Try outer cylinder intersection, c=(t3-fRMax*fRMax)/t1;
1028 912
1029 c /= t1 ; 913 c /= t1 ;
1030 d = b*b - c ; 914 d = b*b - c ;
1031 915
1032 if (d >= 0) // If real root 916 if (d >= 0) // If real root
1033 { 917 {
1034 sd = c/(-b+std::sqrt(d)); 918 sd = c/(-b+std::sqrt(d));
1035 if (sd >= 0) // If 'forwards' 919 if (sd >= 0) // If 'forwards'
1036 { 920 {
1037 if ( sd>dRmax ) // Avoid rounding e 921 if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
1038 { // 64 bits systems. 922 { // 64 bits systems. Split long distances and recompute
1039 G4double fTerm = sd-std::fmod(sd, 923 G4double fTerm = sd-std::fmod(sd,dRmax);
1040 sd = fTerm + DistanceToIn(p+fTerm 924 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1041 } << 925 }
1042 // Check z intersection 926 // Check z intersection
1043 // 927 //
1044 zi = p.z() + sd*v.z() ; 928 zi = p.z() + sd*v.z() ;
1045 xi = p.x() + sd*v.x() ; 929 xi = p.x() + sd*v.x() ;
1046 yi = p.y() + sd*v.y() ; 930 yi = p.y() + sd*v.y() ;
1047 if ((-xi*fLowNorm.x()-yi*fLowNorm.y 931 if ((-xi*fLowNorm.x()-yi*fLowNorm.y()
1048 -(zi+fDz)*fLowNorm.z())>-halfC 932 -(zi+fDz)*fLowNorm.z())>-halfCarTolerance)
1049 { 933 {
1050 if ((-xi*fHighNorm.x()-yi*fHighNo 934 if ((-xi*fHighNorm.x()-yi*fHighNorm.y()
1051 +(fDz-zi)*fHighNorm.z())>-ha 935 +(fDz-zi)*fHighNorm.z())>-halfCarTolerance)
1052 { 936 {
1053 // Z ok. Check phi intersection 937 // Z ok. Check phi intersection if reqd
1054 // 938 //
1055 if (fPhiFullCutTube) 939 if (fPhiFullCutTube)
1056 { 940 {
1057 return sd ; 941 return sd ;
1058 } 942 }
1059 else 943 else
1060 { 944 {
1061 xi = p.x() + sd*v.x() ; 945 xi = p.x() + sd*v.x() ;
1062 yi = p.y() + sd*v.y() ; 946 yi = p.y() + sd*v.y() ;
1063 cosPsi = (xi*cosCPhi + yi*sin 947 cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMax ;
1064 if (cosPsi >= cosHDPhiIT) { 948 if (cosPsi >= cosHDPhiIT) { return sd ; }
1065 } 949 }
1066 } // end if std::fabs(zi) 950 } // end if std::fabs(zi)
1067 } 951 }
1068 } // end if (sd>=0) 952 } // end if (sd>=0)
1069 } // end if (d>=0) 953 } // end if (d>=0)
1070 } // end if (r>=fRMax) 954 } // end if (r>=fRMax)
1071 else << 955 else
1072 { 956 {
1073 // Inside outer radius : 957 // Inside outer radius :
1074 // check not inside, and heading throug 958 // check not inside, and heading through tubs (-> 0 to in)
1075 if ((t3 > tolIRMin2) && (t2 < 0) 959 if ((t3 > tolIRMin2) && (t2 < 0)
1076 && (std::fabs(p.z()) <= std::fabs(GetC 960 && (std::fabs(p.z()) <= std::fabs(GetCutZ(p))-halfCarTolerance ))
1077 { 961 {
1078 // Inside both radii, delta r -ve, in 962 // Inside both radii, delta r -ve, inside z extent
1079 963
1080 if (!fPhiFullCutTube) 964 if (!fPhiFullCutTube)
1081 { 965 {
1082 inum = p.x()*cosCPhi + p.y()*sinC 966 inum = p.x()*cosCPhi + p.y()*sinCPhi ;
1083 iden = std::sqrt(t3) ; 967 iden = std::sqrt(t3) ;
1084 cosPsi = inum/iden ; 968 cosPsi = inum/iden ;
1085 if (cosPsi >= cosHDPhiIT) 969 if (cosPsi >= cosHDPhiIT)
1086 { 970 {
1087 // In the old version, the small 971 // In the old version, the small negative tangent for the point
1088 // on surface was not taken in ac 972 // on surface was not taken in account, and returning 0.0 ...
1089 // New version: check the tangent << 973 // New version: check the tangent for the point on surface and
1090 // if no intersection, return kIn 974 // if no intersection, return kInfinity, if intersection instead
1091 // return sd. 975 // return sd.
1092 // 976 //
1093 c = t3-fRMax*fRMax; << 977 c = t3-fRMax*fRMax;
1094 if ( c<=0.0 ) 978 if ( c<=0.0 )
1095 { 979 {
1096 return 0.0; 980 return 0.0;
1097 } 981 }
1098 else 982 else
1099 { 983 {
1100 c = c/t1 ; 984 c = c/t1 ;
1101 d = b*b-c; 985 d = b*b-c;
1102 if ( d>=0.0 ) 986 if ( d>=0.0 )
1103 { 987 {
1104 snxt = c/(-b+std::sqrt(d)); / 988 snxt = c/(-b+std::sqrt(d)); // using safe solution
1105 / << 989 // for quadratic equation
1106 if ( snxt < halfCarTolerance 990 if ( snxt < halfCarTolerance ) { snxt=0; }
1107 return snxt ; 991 return snxt ;
1108 } << 992 }
1109 else 993 else
1110 { 994 {
1111 return kInfinity; 995 return kInfinity;
1112 } 996 }
1113 } 997 }
1114 } << 998 }
1115 } 999 }
1116 else 1000 else
1117 { << 1001 {
1118 // In the old version, the small ne 1002 // In the old version, the small negative tangent for the point
1119 // on surface was not taken in acco 1003 // on surface was not taken in account, and returning 0.0 ...
1120 // New version: check the tangent f << 1004 // New version: check the tangent for the point on surface and
1121 // if no intersection, return kInfi 1005 // if no intersection, return kInfinity, if intersection instead
1122 // return sd. 1006 // return sd.
1123 // 1007 //
1124 c = t3 - fRMax*fRMax; << 1008 c = t3 - fRMax*fRMax;
1125 if ( c<=0.0 ) 1009 if ( c<=0.0 )
1126 { 1010 {
1127 return 0.0; 1011 return 0.0;
1128 } 1012 }
1129 else 1013 else
1130 { 1014 {
1131 c = c/t1 ; 1015 c = c/t1 ;
1132 d = b*b-c; 1016 d = b*b-c;
1133 if ( d>=0.0 ) 1017 if ( d>=0.0 )
1134 { 1018 {
1135 snxt= c/(-b+std::sqrt(d)); // u 1019 snxt= c/(-b+std::sqrt(d)); // using safe solution
1136 // f << 1020 // for quadratic equation
1137 if ( snxt < halfCarTolerance ) 1021 if ( snxt < halfCarTolerance ) { snxt=0; }
1138 return snxt ; 1022 return snxt ;
1139 } << 1023 }
1140 else 1024 else
1141 { 1025 {
1142 return kInfinity; 1026 return kInfinity;
1143 } 1027 }
1144 } 1028 }
1145 } // end if (!fPhiFullCutTube) 1029 } // end if (!fPhiFullCutTube)
1146 } // end if (t3>tolIRMin2) 1030 } // end if (t3>tolIRMin2)
1147 } // end if (Inside Outer Radius) << 1031 } // end if (Inside Outer Radius)
1148 << 1032
1149 if ( fRMin != 0.0 ) // Try inner cylin << 1033 if ( fRMin ) // Try inner cylinder intersection
1150 { 1034 {
1151 c = (t3 - fRMin*fRMin)/t1 ; 1035 c = (t3 - fRMin*fRMin)/t1 ;
1152 d = b*b - c ; 1036 d = b*b - c ;
1153 if ( d >= 0.0 ) // If real root 1037 if ( d >= 0.0 ) // If real root
1154 { 1038 {
1155 // Always want 2nd root - we are outs 1039 // Always want 2nd root - we are outside and know rmax Hit was bad
1156 // - If on surface of rmin also need 1040 // - If on surface of rmin also need farthest root
1157 << 1041
1158 sd =( b > 0. )? c/(-b - std::sqrt(d)) 1042 sd =( b > 0. )? c/(-b - std::sqrt(d)) : (-b + std::sqrt(d));
1159 if (sd >= -10*halfCarTolerance) // c 1043 if (sd >= -10*halfCarTolerance) // check forwards
1160 { 1044 {
1161 // Check z intersection 1045 // Check z intersection
1162 // 1046 //
1163 if (sd < 0.0) { sd = 0.0; } 1047 if (sd < 0.0) { sd = 0.0; }
1164 if (sd>dRmax) // Avoid rounding err 1048 if (sd>dRmax) // Avoid rounding errors due to precision issues seen
1165 { // 64 bits systems. S 1049 { // 64 bits systems. Split long distances and recompute
1166 G4double fTerm = sd-std::fmod(sd, 1050 G4double fTerm = sd-std::fmod(sd,dRmax);
1167 sd = fTerm + DistanceToIn(p+fTerm 1051 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1168 } << 1052 }
1169 zi = p.z() + sd*v.z() ; 1053 zi = p.z() + sd*v.z() ;
1170 xi = p.x() + sd*v.x() ; 1054 xi = p.x() + sd*v.x() ;
1171 yi = p.y() + sd*v.y() ; 1055 yi = p.y() + sd*v.y() ;
1172 if ((-xi*fLowNorm.x()-yi*fLowNorm.y 1056 if ((-xi*fLowNorm.x()-yi*fLowNorm.y()
1173 -(zi+fDz)*fLowNorm.z())>-halfC 1057 -(zi+fDz)*fLowNorm.z())>-halfCarTolerance)
1174 { 1058 {
1175 if ((-xi*fHighNorm.x()-yi*fHighNo 1059 if ((-xi*fHighNorm.x()-yi*fHighNorm.y()
1176 +(fDz-zi)*fHighNorm.z())>-ha 1060 +(fDz-zi)*fHighNorm.z())>-halfCarTolerance)
1177 { 1061 {
1178 // Z ok. Check phi 1062 // Z ok. Check phi
1179 // 1063 //
1180 if ( fPhiFullCutTube ) 1064 if ( fPhiFullCutTube )
1181 { 1065 {
1182 return sd ; << 1066 return sd ;
1183 } 1067 }
1184 else 1068 else
1185 { 1069 {
1186 cosPsi = (xi*cosCPhi + yi*sin 1070 cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMin ;
1187 if (cosPsi >= cosHDPhiIT) 1071 if (cosPsi >= cosHDPhiIT)
1188 { 1072 {
1189 // Good inner radius isect 1073 // Good inner radius isect
1190 // - but earlier phi isect 1074 // - but earlier phi isect still possible
1191 // 1075 //
1192 snxt = sd ; 1076 snxt = sd ;
1193 } 1077 }
1194 } 1078 }
1195 } // end if std::fabs(zi) 1079 } // end if std::fabs(zi)
1196 } 1080 }
1197 } // end if (sd>=0) 1081 } // end if (sd>=0)
1198 } // end if (d>=0) 1082 } // end if (d>=0)
1199 } // end if (fRMin) 1083 } // end if (fRMin)
1200 } 1084 }
1201 1085
1202 // Phi segment intersection 1086 // Phi segment intersection
1203 // 1087 //
1204 // o Tolerant of points inside phi planes b 1088 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1205 // 1089 //
1206 // o NOTE: Large duplication of code betwee 1090 // o NOTE: Large duplication of code between sphi & ephi checks
1207 // -> only diffs: sphi -> ephi, Com 1091 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1208 // intersection check <=0 -> >=0 1092 // intersection check <=0 -> >=0
1209 // -> use some form of loop Constru 1093 // -> use some form of loop Construct ?
1210 // 1094 //
1211 if ( !fPhiFullCutTube ) 1095 if ( !fPhiFullCutTube )
1212 { 1096 {
1213 // First phi surface (Starting phi) 1097 // First phi surface (Starting phi)
1214 // 1098 //
1215 Comp = v.x()*sinSPhi - v.y()*cosSPhi ; 1099 Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
1216 << 1100
1217 if ( Comp < 0 ) // Component in outwards 1101 if ( Comp < 0 ) // Component in outwards normal dirn
1218 { 1102 {
1219 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) 1103 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1220 1104
1221 if ( Dist < halfCarTolerance ) 1105 if ( Dist < halfCarTolerance )
1222 { 1106 {
1223 sd = Dist/Comp ; 1107 sd = Dist/Comp ;
1224 1108
1225 if (sd < snxt) 1109 if (sd < snxt)
1226 { 1110 {
1227 if ( sd < 0 ) { sd = 0.0; } 1111 if ( sd < 0 ) { sd = 0.0; }
1228 zi = p.z() + sd*v.z() ; 1112 zi = p.z() + sd*v.z() ;
1229 xi = p.x() + sd*v.x() ; 1113 xi = p.x() + sd*v.x() ;
1230 yi = p.y() + sd*v.y() ; 1114 yi = p.y() + sd*v.y() ;
1231 if ((-xi*fLowNorm.x()-yi*fLowNorm.y 1115 if ((-xi*fLowNorm.x()-yi*fLowNorm.y()
1232 -(zi+fDz)*fLowNorm.z())>-halfC 1116 -(zi+fDz)*fLowNorm.z())>-halfCarTolerance)
1233 { 1117 {
1234 if ((-xi*fHighNorm.x()-yi*fHighNo 1118 if ((-xi*fHighNorm.x()-yi*fHighNorm.y()
1235 +(fDz-zi)*fHighNorm.z())>-ha << 1119 +(fDz-zi)*fHighNorm.z())>-halfCarTolerance)
1236 { 1120 {
1237 rho2 = xi*xi + yi*yi ; 1121 rho2 = xi*xi + yi*yi ;
1238 if ( ( (rho2 >= tolIRMin2) && ( 1122 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1239 || ( (rho2 > tolORMin2) && ( 1123 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1240 && ( v.y()*cosSPhi - v.x()* 1124 && ( v.y()*cosSPhi - v.x()*sinSPhi > 0 )
1241 && ( v.x()*cosSPhi + v.y()* 1125 && ( v.x()*cosSPhi + v.y()*sinSPhi >= 0 ) )
1242 || ( (rho2 > tolIRMax2) && (r 1126 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1243 && (v.y()*cosSPhi - v.x()*s 1127 && (v.y()*cosSPhi - v.x()*sinSPhi > 0)
1244 && (v.x()*cosSPhi + v.y()*s 1128 && (v.x()*cosSPhi + v.y()*sinSPhi < 0) ) )
1245 { 1129 {
1246 // z and r intersections good 1130 // z and r intersections good
1247 // - check intersecting with 1131 // - check intersecting with correct half-plane
1248 // 1132 //
1249 if ((yi*cosCPhi-xi*sinCPhi) < 1133 if ((yi*cosCPhi-xi*sinCPhi) <= halfCarTolerance) { snxt = sd; }
1250 } 1134 }
1251 } //two Z conditions 1135 } //two Z conditions
1252 } 1136 }
1253 } 1137 }
1254 } << 1138 }
1255 } 1139 }
1256 << 1140
1257 // Second phi surface (Ending phi) 1141 // Second phi surface (Ending phi)
1258 // 1142 //
1259 Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ; 1143 Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
1260 << 1144
1261 if (Comp < 0 ) // Component in outwards 1145 if (Comp < 0 ) // Component in outwards normal dirn
1262 { 1146 {
1263 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) 1147 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1264 1148
1265 if ( Dist < halfCarTolerance ) 1149 if ( Dist < halfCarTolerance )
1266 { 1150 {
1267 sd = Dist/Comp ; 1151 sd = Dist/Comp ;
1268 1152
1269 if (sd < snxt) 1153 if (sd < snxt)
1270 { 1154 {
1271 if ( sd < 0 ) { sd = 0; } 1155 if ( sd < 0 ) { sd = 0; }
1272 zi = p.z() + sd*v.z() ; 1156 zi = p.z() + sd*v.z() ;
1273 xi = p.x() + sd*v.x() ; 1157 xi = p.x() + sd*v.x() ;
1274 yi = p.y() + sd*v.y() ; 1158 yi = p.y() + sd*v.y() ;
1275 if ((-xi*fLowNorm.x()-yi*fLowNorm.y 1159 if ((-xi*fLowNorm.x()-yi*fLowNorm.y()
1276 -(zi+fDz)*fLowNorm.z())>-halfC 1160 -(zi+fDz)*fLowNorm.z())>-halfCarTolerance)
1277 { 1161 {
1278 if ((-xi*fHighNorm.x()-yi*fHighNo 1162 if ((-xi*fHighNorm.x()-yi*fHighNorm.y()
1279 +(fDz-zi)*fHighNorm.z())>-ha 1163 +(fDz-zi)*fHighNorm.z())>-halfCarTolerance)
1280 { 1164 {
1281 xi = p.x() + sd*v.x() ; 1165 xi = p.x() + sd*v.x() ;
1282 yi = p.y() + sd*v.y() ; 1166 yi = p.y() + sd*v.y() ;
1283 rho2 = xi*xi + yi*yi ; 1167 rho2 = xi*xi + yi*yi ;
1284 if ( ( (rho2 >= tolIRMin2) && ( 1168 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1285 || ( (rho2 > tolORMin2) && 1169 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1286 && (v.x()*sinEPhi - v.y() 1170 && (v.x()*sinEPhi - v.y()*cosEPhi > 0)
1287 && (v.x()*cosEPhi + v.y() 1171 && (v.x()*cosEPhi + v.y()*sinEPhi >= 0) )
1288 || ( (rho2 > tolIRMax2) && 1172 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1289 && (v.x()*sinEPhi - v.y() 1173 && (v.x()*sinEPhi - v.y()*cosEPhi > 0)
1290 && (v.x()*cosEPhi + v.y() 1174 && (v.x()*cosEPhi + v.y()*sinEPhi < 0) ) )
1291 { 1175 {
1292 // z and r intersections good 1176 // z and r intersections good
1293 // - check intersecting with 1177 // - check intersecting with correct half-plane
1294 // 1178 //
1295 if ( (yi*cosCPhi-xi*sinCPhi) 1179 if ( (yi*cosCPhi-xi*sinCPhi) >= -halfCarTolerance )
1296 { 1180 {
1297 snxt = sd; 1181 snxt = sd;
1298 } 1182 }
1299 } //?? >=-halfCarTolerance 1183 } //?? >=-halfCarTolerance
1300 } 1184 }
1301 } // two Z conditions 1185 } // two Z conditions
1302 } 1186 }
1303 } 1187 }
1304 } // Comp < 0 1188 } // Comp < 0
1305 } // !fPhiFullTube << 1189 } // !fPhiFullTube
1306 if ( snxt<halfCarTolerance ) { snxt=0; } 1190 if ( snxt<halfCarTolerance ) { snxt=0; }
1307 1191
1308 return snxt ; 1192 return snxt ;
1309 } 1193 }
1310 << 1194
1311 ///////////////////////////////////////////// 1195 //////////////////////////////////////////////////////////////////
1312 // 1196 //
1313 // Calculate distance to shape from outside, 1197 // Calculate distance to shape from outside, along normalised vector
1314 // - return kInfinity if no intersection, or 1198 // - return kInfinity if no intersection, or intersection distance <= tolerance
1315 // 1199 //
1316 // - Compute the intersection with the z plan << 1200 // - Compute the intersection with the z planes
1317 // - if at valid r, phi, return 1201 // - if at valid r, phi, return
1318 // 1202 //
1319 // -> If point is outer outer radius, compute 1203 // -> If point is outer outer radius, compute intersection with rmax
1320 // - if at valid phi,z return 1204 // - if at valid phi,z return
1321 // 1205 //
1322 // -> Compute intersection with inner radius, 1206 // -> Compute intersection with inner radius, taking largest +ve root
1323 // - if valid (in z,phi), save intersc 1207 // - if valid (in z,phi), save intersction
1324 // 1208 //
1325 // -> If phi segmented, compute intersecti 1209 // -> If phi segmented, compute intersections with phi half planes
1326 // - return smallest of valid phi inte 1210 // - return smallest of valid phi intersections and
1327 // inner radius intersection 1211 // inner radius intersection
1328 // 1212 //
1329 // NOTE: 1213 // NOTE:
1330 // - Precalculations for phi trigonometry are 1214 // - Precalculations for phi trigonometry are Done `just in time'
1331 // - `if valid' implies tolerant checking of 1215 // - `if valid' implies tolerant checking of intersection points
1332 // Calculate distance (<= actual) to closes 1216 // Calculate distance (<= actual) to closest surface of shape from outside
1333 // - Calculate distance to z, radial planes 1217 // - Calculate distance to z, radial planes
1334 // - Only to phi planes if outside phi extent 1218 // - Only to phi planes if outside phi extent
1335 // - Return 0 if point inside 1219 // - Return 0 if point inside
1336 1220
1337 G4double G4CutTubs::DistanceToIn( const G4Thr 1221 G4double G4CutTubs::DistanceToIn( const G4ThreeVector& p ) const
1338 { 1222 {
1339 G4double safRMin,safRMax,safZLow,safZHigh,s 1223 G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho,cosPsi;
1340 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 1224 G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
1341 1225
1342 // Distance to R 1226 // Distance to R
1343 // 1227 //
1344 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) 1228 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1345 1229
1346 safRMin = fRMin- rho ; 1230 safRMin = fRMin- rho ;
1347 safRMax = rho - fRMax ; 1231 safRMax = rho - fRMax ;
1348 1232
1349 // Distances to ZCut(Low/High) 1233 // Distances to ZCut(Low/High)
1350 1234
1351 // Dist to Low Cut 1235 // Dist to Low Cut
1352 // 1236 //
1353 safZLow = (p+vZ).dot(fLowNorm); 1237 safZLow = (p+vZ).dot(fLowNorm);
1354 1238
1355 // Dist to High Cut 1239 // Dist to High Cut
1356 // 1240 //
1357 safZHigh = (p-vZ).dot(fHighNorm); 1241 safZHigh = (p-vZ).dot(fHighNorm);
1358 1242
1359 safe = std::max(safZLow,safZHigh); 1243 safe = std::max(safZLow,safZHigh);
1360 1244
1361 if ( safRMin > safe ) { safe = safRMin; } 1245 if ( safRMin > safe ) { safe = safRMin; }
1362 if ( safRMax> safe ) { safe = safRMax; } 1246 if ( safRMax> safe ) { safe = safRMax; }
1363 1247
1364 // Distance to Phi 1248 // Distance to Phi
1365 // 1249 //
1366 if ( (!fPhiFullCutTube) && ((rho) != 0.0) ) << 1250 if ( (!fPhiFullCutTube) && (rho) )
1367 { 1251 {
1368 // Psi=angle from central phi to point 1252 // Psi=angle from central phi to point
1369 // 1253 //
1370 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi) 1254 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
1371 << 1255
1372 if ( cosPsi < cosHDPhi ) << 1256 if ( cosPsi < std::cos(fDPhi*0.5) )
1373 { 1257 {
1374 // Point lies outside phi range 1258 // Point lies outside phi range
1375 << 1259
1376 if ( (p.y()*cosCPhi - p.x()*sinCPhi) < 1260 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
1377 { 1261 {
1378 safePhi = std::fabs(p.x()*sinSPhi - 1262 safePhi = std::fabs(p.x()*sinSPhi - p.y()*cosSPhi) ;
1379 } 1263 }
1380 else 1264 else
1381 { 1265 {
1382 safePhi = std::fabs(p.x()*sinEPhi - 1266 safePhi = std::fabs(p.x()*sinEPhi - p.y()*cosEPhi) ;
1383 } 1267 }
1384 if ( safePhi > safe ) { safe = safePh 1268 if ( safePhi > safe ) { safe = safePhi; }
1385 } 1269 }
1386 } 1270 }
1387 if ( safe < 0 ) { safe = 0; } 1271 if ( safe < 0 ) { safe = 0; }
1388 1272
1389 return safe ; 1273 return safe ;
1390 } 1274 }
1391 1275
1392 ///////////////////////////////////////////// 1276 //////////////////////////////////////////////////////////////////////////////
1393 // 1277 //
1394 // Calculate distance to surface of shape fro 1278 // Calculate distance to surface of shape from `inside', allowing for tolerance
1395 // - Only Calc rmax intersection if no valid 1279 // - Only Calc rmax intersection if no valid rmin intersection
1396 1280
1397 G4double G4CutTubs::DistanceToOut( const G4Th 1281 G4double G4CutTubs::DistanceToOut( const G4ThreeVector& p,
1398 const G4Th 1282 const G4ThreeVector& v,
1399 const G4bo 1283 const G4bool calcNorm,
1400 G4bo << 1284 G4bool *validNorm,
1401 G4Th << 1285 G4ThreeVector *n ) const
1402 { << 1286 {
1403 enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,k <<
1404 <<
1405 ESide side=kNull , sider=kNull, sidephi=kNu 1287 ESide side=kNull , sider=kNull, sidephi=kNull ;
1406 G4double snxt=kInfinity, srd=kInfinity,sz=k 1288 G4double snxt=kInfinity, srd=kInfinity,sz=kInfinity, sphi=kInfinity ;
1407 G4double deltaR, t1, t2, t3, b, c, d2, roMi 1289 G4double deltaR, t1, t2, t3, b, c, d2, roMin2 ;
1408 G4double distZLow,distZHigh,calfH,calfL; 1290 G4double distZLow,distZHigh,calfH,calfL;
1409 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 1291 G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
1410 << 1292 static const G4double halfCarTolerance = kCarTolerance*0.5;
>> 1293 static const G4double halfAngTolerance = kAngTolerance*0.5;
>> 1294
1411 // Vars for phi intersection: 1295 // Vars for phi intersection:
1412 // 1296 //
1413 G4double pDistS, compS, pDistE, compE, sphi 1297 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ;
1414 << 1298
1415 // Z plane intersection 1299 // Z plane intersection
1416 // Distances to ZCut(Low/High) 1300 // Distances to ZCut(Low/High)
1417 1301
1418 // dist to Low Cut 1302 // dist to Low Cut
1419 // 1303 //
1420 distZLow =(p+vZ).dot(fLowNorm); 1304 distZLow =(p+vZ).dot(fLowNorm);
1421 1305
1422 // dist to High Cut 1306 // dist to High Cut
1423 // 1307 //
1424 distZHigh = (p-vZ).dot(fHighNorm); 1308 distZHigh = (p-vZ).dot(fHighNorm);
1425 1309
1426 calfH = v.dot(fHighNorm); 1310 calfH = v.dot(fHighNorm);
1427 calfL = v.dot(fLowNorm); 1311 calfL = v.dot(fLowNorm);
1428 1312
1429 if (calfH > 0 ) 1313 if (calfH > 0 )
1430 { 1314 {
1431 if ( distZHigh < halfCarTolerance ) 1315 if ( distZHigh < halfCarTolerance )
1432 { 1316 {
1433 snxt = -distZHigh/calfH ; 1317 snxt = -distZHigh/calfH ;
1434 side = kPZ ; 1318 side = kPZ ;
1435 } 1319 }
1436 else 1320 else
1437 { 1321 {
1438 if (calcNorm) 1322 if (calcNorm)
1439 { 1323 {
1440 *n = G4ThreeVector(0,0,1) ; 1324 *n = G4ThreeVector(0,0,1) ;
1441 *validNorm = true ; 1325 *validNorm = true ;
1442 } 1326 }
1443 return snxt = 0 ; 1327 return snxt = 0 ;
1444 } 1328 }
1445 } 1329 }
1446 if ( calfL>0) 1330 if ( calfL>0)
1447 { 1331 {
1448 << 1332
1449 if ( distZLow < halfCarTolerance ) 1333 if ( distZLow < halfCarTolerance )
1450 { 1334 {
1451 sz = -distZLow/calfL ; 1335 sz = -distZLow/calfL ;
1452 if(sz<snxt){ 1336 if(sz<snxt){
1453 snxt=sz; 1337 snxt=sz;
1454 side = kMZ ; 1338 side = kMZ ;
1455 } 1339 }
1456 << 1340
1457 } 1341 }
1458 else 1342 else
1459 { 1343 {
1460 if (calcNorm) 1344 if (calcNorm)
1461 { 1345 {
1462 *n = G4ThreeVector(0,0,-1) ; 1346 *n = G4ThreeVector(0,0,-1) ;
1463 *validNorm = true ; 1347 *validNorm = true ;
1464 } 1348 }
1465 return snxt = 0.0 ; 1349 return snxt = 0.0 ;
1466 } 1350 }
1467 } 1351 }
1468 if((calfH<=0)&&(calfL<=0)) 1352 if((calfH<=0)&&(calfL<=0))
1469 { 1353 {
1470 snxt = kInfinity ; // Travel perpendic 1354 snxt = kInfinity ; // Travel perpendicular to z axis
1471 side = kNull; 1355 side = kNull;
1472 } 1356 }
1473 // Radial Intersections 1357 // Radial Intersections
1474 // 1358 //
1475 // Find intersection with cylinders at rmax 1359 // Find intersection with cylinders at rmax/rmin
1476 // Intersection point (xi,yi,zi) on line x= 1360 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1477 // 1361 //
1478 // Intersects with x^2+y^2=R^2 1362 // Intersects with x^2+y^2=R^2
1479 // 1363 //
1480 // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v 1364 // 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
1481 // 1365 //
1482 // t1 t2 1366 // t1 t2 t3
1483 1367
1484 t1 = 1.0 - v.z()*v.z() ; // since v 1368 t1 = 1.0 - v.z()*v.z() ; // since v normalised
1485 t2 = p.x()*v.x() + p.y()*v.y() ; 1369 t2 = p.x()*v.x() + p.y()*v.y() ;
1486 t3 = p.x()*p.x() + p.y()*p.y() ; 1370 t3 = p.x()*p.x() + p.y()*p.y() ;
1487 1371
1488 if ( snxt > 10*(fDz+fRMax) ) { roi2 = 2*fR 1372 if ( snxt > 10*(fDz+fRMax) ) { roi2 = 2*fRMax*fRMax; }
1489 else { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t 1373 else { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t3; } // radius^2 on +-fDz
1490 1374
1491 if ( t1 > 0 ) // Check not parallel 1375 if ( t1 > 0 ) // Check not parallel
1492 { 1376 {
1493 // Calculate srd, r exit distance 1377 // Calculate srd, r exit distance
1494 << 1378
1495 if ( (t2 >= 0.0) && (roi2 > fRMax*(fRMax 1379 if ( (t2 >= 0.0) && (roi2 > fRMax*(fRMax + kRadTolerance)) )
1496 { 1380 {
1497 // Delta r not negative => leaving via 1381 // Delta r not negative => leaving via rmax
1498 1382
1499 deltaR = t3 - fRMax*fRMax ; 1383 deltaR = t3 - fRMax*fRMax ;
1500 1384
1501 // NOTE: Should use rho-fRMax<-kRadTole 1385 // NOTE: Should use rho-fRMax<-kRadTolerance*0.5
1502 // - avoid sqrt for efficiency 1386 // - avoid sqrt for efficiency
1503 1387
1504 if ( deltaR < -kRadTolerance*fRMax ) 1388 if ( deltaR < -kRadTolerance*fRMax )
1505 { 1389 {
1506 b = t2/t1 ; 1390 b = t2/t1 ;
1507 c = deltaR/t1 ; 1391 c = deltaR/t1 ;
1508 d2 = b*b-c; 1392 d2 = b*b-c;
1509 if( d2 >= 0 ) { srd = c/( -b - std::s 1393 if( d2 >= 0 ) { srd = c/( -b - std::sqrt(d2)); }
1510 else { srd = 0.; } 1394 else { srd = 0.; }
1511 sider = kRMax ; 1395 sider = kRMax ;
1512 } 1396 }
1513 else 1397 else
1514 { 1398 {
1515 // On tolerant boundary & heading out 1399 // On tolerant boundary & heading outwards (or perpendicular to)
1516 // outer radial surface -> leaving im 1400 // outer radial surface -> leaving immediately
1517 1401
1518 if ( calcNorm ) << 1402 if ( calcNorm )
1519 { 1403 {
1520 *n = G4ThreeVector(p.x()/fR 1404 *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1521 *validNorm = true ; 1405 *validNorm = true ;
1522 } 1406 }
1523 return snxt = 0 ; // Leaving by rmax 1407 return snxt = 0 ; // Leaving by rmax immediately
1524 } 1408 }
1525 } << 1409 }
1526 else if ( t2 < 0. ) // i.e. t2 < 0; Poss 1410 else if ( t2 < 0. ) // i.e. t2 < 0; Possible rmin intersection
1527 { 1411 {
1528 roMin2 = t3 - t2*t2/t1 ; // min ro2 of << 1412 roMin2 = t3 - t2*t2/t1 ; // min ro2 of the plane of movement
1529 1413
1530 if ( (fRMin != 0.0) && (roMin2 < fRMin* << 1414 if ( fRMin && (roMin2 < fRMin*(fRMin - kRadTolerance)) )
1531 { 1415 {
1532 deltaR = t3 - fRMin*fRMin ; 1416 deltaR = t3 - fRMin*fRMin ;
1533 b = t2/t1 ; 1417 b = t2/t1 ;
1534 c = deltaR/t1 ; 1418 c = deltaR/t1 ;
1535 d2 = b*b - c ; 1419 d2 = b*b - c ;
1536 1420
1537 if ( d2 >= 0 ) // Leaving via rmin 1421 if ( d2 >= 0 ) // Leaving via rmin
1538 { 1422 {
1539 // NOTE: SHould use rho-rmin>kRadTo 1423 // NOTE: SHould use rho-rmin>kRadTolerance*0.5
1540 // - avoid sqrt for efficiency 1424 // - avoid sqrt for efficiency
1541 1425
1542 if (deltaR > kRadTolerance*fRMin) 1426 if (deltaR > kRadTolerance*fRMin)
1543 { 1427 {
1544 srd = c/(-b+std::sqrt(d2)); << 1428 srd = c/(-b+std::sqrt(d2));
1545 sider = kRMin ; 1429 sider = kRMin ;
1546 } 1430 }
1547 else 1431 else
1548 { 1432 {
1549 if ( calcNorm ) { *validNorm = fa 1433 if ( calcNorm ) { *validNorm = false; } // Concave side
1550 return snxt = 0.0; 1434 return snxt = 0.0;
1551 } 1435 }
1552 } 1436 }
1553 else // No rmin intersect -> must 1437 else // No rmin intersect -> must be rmax intersect
1554 { 1438 {
1555 deltaR = t3 - fRMax*fRMax ; 1439 deltaR = t3 - fRMax*fRMax ;
1556 c = deltaR/t1 ; 1440 c = deltaR/t1 ;
1557 d2 = b*b-c; 1441 d2 = b*b-c;
1558 if( d2 >=0. ) 1442 if( d2 >=0. )
1559 { 1443 {
1560 srd = -b + std::sqrt(d2) ; 1444 srd = -b + std::sqrt(d2) ;
1561 sider = kRMax ; 1445 sider = kRMax ;
1562 } 1446 }
1563 else // Case: On the border+t2<kRad 1447 else // Case: On the border+t2<kRadTolerance
1564 // (v is perpendicular t 1448 // (v is perpendicular to the surface)
1565 { 1449 {
1566 if (calcNorm) 1450 if (calcNorm)
1567 { 1451 {
1568 *n = G4ThreeVector(p.x()/fRMax, 1452 *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1569 *validNorm = true ; 1453 *validNorm = true ;
1570 } 1454 }
1571 return snxt = 0.0; 1455 return snxt = 0.0;
1572 } 1456 }
1573 } 1457 }
1574 } 1458 }
1575 else if ( roi2 > fRMax*(fRMax + kRadTol 1459 else if ( roi2 > fRMax*(fRMax + kRadTolerance) )
1576 // No rmin intersect -> must be rm 1460 // No rmin intersect -> must be rmax intersect
1577 { 1461 {
1578 deltaR = t3 - fRMax*fRMax ; 1462 deltaR = t3 - fRMax*fRMax ;
1579 b = t2/t1 ; 1463 b = t2/t1 ;
1580 c = deltaR/t1; 1464 c = deltaR/t1;
1581 d2 = b*b-c; 1465 d2 = b*b-c;
1582 if( d2 >= 0 ) 1466 if( d2 >= 0 )
1583 { 1467 {
1584 srd = -b + std::sqrt(d2) ; 1468 srd = -b + std::sqrt(d2) ;
1585 sider = kRMax ; 1469 sider = kRMax ;
1586 } 1470 }
1587 else // Case: On the border+t2<kRadTo 1471 else // Case: On the border+t2<kRadTolerance
1588 // (v is perpendicular to 1472 // (v is perpendicular to the surface)
1589 { 1473 {
1590 if (calcNorm) 1474 if (calcNorm)
1591 { 1475 {
1592 *n = G4ThreeVector(p.x()/fRMax,p. 1476 *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1593 *validNorm = true ; 1477 *validNorm = true ;
1594 } 1478 }
1595 return snxt = 0.0; 1479 return snxt = 0.0;
1596 } 1480 }
1597 } 1481 }
1598 } 1482 }
1599 // Phi Intersection 1483 // Phi Intersection
1600 1484
1601 if ( !fPhiFullCutTube ) 1485 if ( !fPhiFullCutTube )
1602 { 1486 {
1603 // add angle calculation with correctio << 1487 // add angle calculation with correction
1604 // of the difference in domain of atan2 1488 // of the difference in domain of atan2 and Sphi
1605 // 1489 //
1606 vphi = std::atan2(v.y(),v.x()) ; 1490 vphi = std::atan2(v.y(),v.x()) ;
1607 << 1491
1608 if ( vphi < fSPhi - halfAngTolerance ) 1492 if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; }
1609 else if ( vphi > fSPhi + fDPhi + halfAn 1493 else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; }
1610 1494
1611 1495
1612 if ( (p.x() != 0.0) || (p.y() != 0.0) ) << 1496 if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
1613 { 1497 {
1614 // pDist -ve when inside 1498 // pDist -ve when inside
1615 1499
1616 pDistS = p.x()*sinSPhi - p.y()*cosSPh 1500 pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
1617 pDistE = -p.x()*sinEPhi + p.y()*cosEP 1501 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1618 1502
1619 // Comp -ve when in direction of outw 1503 // Comp -ve when in direction of outwards normal
1620 1504
1621 compS = -sinSPhi*v.x() + cosSPhi*v. 1505 compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
1622 compE = sinEPhi*v.x() - cosEPhi*v. 1506 compE = sinEPhi*v.x() - cosEPhi*v.y() ;
1623 << 1507
1624 sidephi = kNull; 1508 sidephi = kNull;
1625 << 1509
1626 if( ( (fDPhi <= pi) && ( (pDistS <= h 1510 if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
1627 && (pDistE <= h 1511 && (pDistE <= halfCarTolerance) ) )
1628 || ( (fDPhi > pi) && ((pDistS <= h << 1512 || ( (fDPhi > pi) && !((pDistS > halfCarTolerance)
1629 || (pDistE <= << 1513 && (pDistE > halfCarTolerance) ) ) )
1630 { 1514 {
1631 // Inside both phi *full* planes 1515 // Inside both phi *full* planes
1632 << 1516
1633 if ( compS < 0 ) 1517 if ( compS < 0 )
1634 { 1518 {
1635 sphi = pDistS/compS ; 1519 sphi = pDistS/compS ;
1636 << 1520
1637 if (sphi >= -halfCarTolerance) 1521 if (sphi >= -halfCarTolerance)
1638 { 1522 {
1639 xi = p.x() + sphi*v.x() ; 1523 xi = p.x() + sphi*v.x() ;
1640 yi = p.y() + sphi*v.y() ; 1524 yi = p.y() + sphi*v.y() ;
1641 << 1525
1642 // Check intersecting with corr 1526 // Check intersecting with correct half-plane
1643 // (if not -> no intersect) 1527 // (if not -> no intersect)
1644 // 1528 //
1645 if( (std::fabs(xi)<=kCarToleran 1529 if( (std::fabs(xi)<=kCarTolerance)
1646 && (std::fabs(yi)<=kCarToleran 1530 && (std::fabs(yi)<=kCarTolerance) )
1647 { 1531 {
1648 sidephi = kSPhi; 1532 sidephi = kSPhi;
1649 if (((fSPhi-halfAngTolerance) 1533 if (((fSPhi-halfAngTolerance)<=vphi)
1650 &&((fSPhi+fDPhi+halfAngTol 1534 &&((fSPhi+fDPhi+halfAngTolerance)>=vphi))
1651 { 1535 {
1652 sphi = kInfinity; 1536 sphi = kInfinity;
1653 } 1537 }
1654 } 1538 }
1655 else if ( yi*cosCPhi-xi*sinCPhi 1539 else if ( yi*cosCPhi-xi*sinCPhi >=0 )
1656 { 1540 {
1657 sphi = kInfinity ; 1541 sphi = kInfinity ;
1658 } 1542 }
1659 else 1543 else
1660 { 1544 {
1661 sidephi = kSPhi ; 1545 sidephi = kSPhi ;
1662 if ( pDistS > -halfCarToleran 1546 if ( pDistS > -halfCarTolerance )
1663 { 1547 {
1664 sphi = 0.0 ; // Leave by sp 1548 sphi = 0.0 ; // Leave by sphi immediately
1665 } << 1549 }
1666 } << 1550 }
1667 } 1551 }
1668 else 1552 else
1669 { 1553 {
1670 sphi = kInfinity ; 1554 sphi = kInfinity ;
1671 } 1555 }
1672 } 1556 }
1673 else 1557 else
1674 { 1558 {
1675 sphi = kInfinity ; 1559 sphi = kInfinity ;
1676 } 1560 }
1677 1561
1678 if ( compE < 0 ) 1562 if ( compE < 0 )
1679 { 1563 {
1680 sphi2 = pDistE/compE ; 1564 sphi2 = pDistE/compE ;
1681 << 1565
1682 // Only check further if < starti 1566 // Only check further if < starting phi intersection
1683 // 1567 //
1684 if ( (sphi2 > -halfCarTolerance) 1568 if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
1685 { 1569 {
1686 xi = p.x() + sphi2*v.x() ; 1570 xi = p.x() + sphi2*v.x() ;
1687 yi = p.y() + sphi2*v.y() ; 1571 yi = p.y() + sphi2*v.y() ;
1688 << 1572
1689 if ((std::fabs(xi)<=kCarToleran 1573 if ((std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance))
1690 { 1574 {
1691 // Leaving via ending phi 1575 // Leaving via ending phi
1692 // 1576 //
1693 if( (fSPhi-halfAngTolerance > << 1577 if( !((fSPhi-halfAngTolerance <= vphi)
1694 ||(fSPhi+fDPhi+halfAngTo << 1578 &&(fSPhi+fDPhi+halfAngTolerance >= vphi)) )
1695 { 1579 {
1696 sidephi = kEPhi ; 1580 sidephi = kEPhi ;
1697 if ( pDistE <= -halfCarTole 1581 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; }
1698 else 1582 else { sphi = 0.0 ; }
1699 } 1583 }
1700 } << 1584 }
1701 else // Check intersecting w << 1585 else // Check intersecting with correct half-plane
1702 1586
1703 if ( (yi*cosCPhi-xi*sinCPhi) >= 1587 if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
1704 { 1588 {
1705 // Leaving via ending phi 1589 // Leaving via ending phi
1706 // 1590 //
1707 sidephi = kEPhi ; 1591 sidephi = kEPhi ;
1708 if ( pDistE <= -halfCarTolera 1592 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; }
1709 else 1593 else { sphi = 0.0 ; }
1710 } 1594 }
1711 } 1595 }
1712 } 1596 }
1713 } 1597 }
1714 else 1598 else
1715 { 1599 {
1716 sphi = kInfinity ; 1600 sphi = kInfinity ;
1717 } 1601 }
1718 } 1602 }
1719 else 1603 else
1720 { 1604 {
1721 // On z axis + travel not || to z axi 1605 // On z axis + travel not || to z axis -> if phi of vector direction
1722 // within phi of shape, Step limited 1606 // within phi of shape, Step limited by rmax, else Step =0
1723 << 1607
1724 if ( (fSPhi - halfAngTolerance <= vph 1608 if ( (fSPhi - halfAngTolerance <= vphi)
1725 && (vphi <= fSPhi + fDPhi + halfAn 1609 && (vphi <= fSPhi + fDPhi + halfAngTolerance ) )
1726 { 1610 {
1727 sphi = kInfinity ; 1611 sphi = kInfinity ;
1728 } 1612 }
1729 else 1613 else
1730 { 1614 {
1731 sidephi = kSPhi ; // arbitrary << 1615 sidephi = kSPhi ; // arbitrary
1732 sphi = 0.0 ; 1616 sphi = 0.0 ;
1733 } 1617 }
1734 } 1618 }
1735 if (sphi < snxt) // Order intersecttio 1619 if (sphi < snxt) // Order intersecttions
1736 { 1620 {
1737 snxt = sphi ; 1621 snxt = sphi ;
1738 side = sidephi ; 1622 side = sidephi ;
1739 } 1623 }
1740 } 1624 }
1741 if (srd < snxt) // Order intersections 1625 if (srd < snxt) // Order intersections
1742 { 1626 {
1743 snxt = srd ; 1627 snxt = srd ;
1744 side = sider ; 1628 side = sider ;
1745 } 1629 }
1746 } 1630 }
1747 if (calcNorm) 1631 if (calcNorm)
1748 { 1632 {
1749 switch(side) 1633 switch(side)
1750 { 1634 {
1751 case kRMax: 1635 case kRMax:
1752 // Note: returned vector not normalis 1636 // Note: returned vector not normalised
1753 // (divide by fRMax for unit vector) 1637 // (divide by fRMax for unit vector)
1754 // 1638 //
1755 xi = p.x() + snxt*v.x() ; 1639 xi = p.x() + snxt*v.x() ;
1756 yi = p.y() + snxt*v.y() ; 1640 yi = p.y() + snxt*v.y() ;
1757 *n = G4ThreeVector(xi/fRMax,yi/fRMax, 1641 *n = G4ThreeVector(xi/fRMax,yi/fRMax,0) ;
1758 *validNorm = true ; 1642 *validNorm = true ;
1759 break ; 1643 break ;
1760 1644
1761 case kRMin: 1645 case kRMin:
1762 *validNorm = false ; // Rmin is inco 1646 *validNorm = false ; // Rmin is inconvex
1763 break ; 1647 break ;
1764 1648
1765 case kSPhi: 1649 case kSPhi:
1766 if ( fDPhi <= pi ) 1650 if ( fDPhi <= pi )
1767 { 1651 {
1768 *n = G4ThreeVector(sinSPhi, 1652 *n = G4ThreeVector(sinSPhi,-cosSPhi,0) ;
1769 *validNorm = true ; 1653 *validNorm = true ;
1770 } 1654 }
1771 else 1655 else
1772 { 1656 {
1773 *validNorm = false ; 1657 *validNorm = false ;
1774 } 1658 }
1775 break ; 1659 break ;
1776 1660
1777 case kEPhi: 1661 case kEPhi:
1778 if (fDPhi <= pi) 1662 if (fDPhi <= pi)
1779 { 1663 {
1780 *n = G4ThreeVector(-sinEPhi,cosEPhi 1664 *n = G4ThreeVector(-sinEPhi,cosEPhi,0) ;
1781 *validNorm = true ; 1665 *validNorm = true ;
1782 } 1666 }
1783 else 1667 else
1784 { 1668 {
1785 *validNorm = false ; 1669 *validNorm = false ;
1786 } 1670 }
1787 break ; 1671 break ;
1788 1672
1789 case kPZ: 1673 case kPZ:
1790 *n = fHighNorm ; 1674 *n = fHighNorm ;
1791 *validNorm = true ; 1675 *validNorm = true ;
1792 break ; 1676 break ;
1793 1677
1794 case kMZ: 1678 case kMZ:
1795 *n = fLowNorm ; 1679 *n = fLowNorm ;
1796 *validNorm = true ; 1680 *validNorm = true ;
1797 break ; 1681 break ;
1798 1682
1799 default: 1683 default:
1800 G4cout << G4endl ; 1684 G4cout << G4endl ;
1801 DumpInfo(); 1685 DumpInfo();
1802 std::ostringstream message; 1686 std::ostringstream message;
1803 G4long oldprc = message.precision(16) << 1687 G4int oldprc = message.precision(16);
1804 message << "Undefined side for valid 1688 message << "Undefined side for valid surface normal to solid."
1805 << G4endl 1689 << G4endl
1806 << "Position:" << G4endl << 1690 << "Position:" << G4endl << G4endl
1807 << "p.x() = " << p.x()/mm < 1691 << "p.x() = " << p.x()/mm << " mm" << G4endl
1808 << "p.y() = " << p.y()/mm < 1692 << "p.y() = " << p.y()/mm << " mm" << G4endl
1809 << "p.z() = " << p.z()/mm < 1693 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
1810 << "Direction:" << G4endl << 1694 << "Direction:" << G4endl << G4endl
1811 << "v.x() = " << v.x() << G 1695 << "v.x() = " << v.x() << G4endl
1812 << "v.y() = " << v.y() << G 1696 << "v.y() = " << v.y() << G4endl
1813 << "v.z() = " << v.z() << G 1697 << "v.z() = " << v.z() << G4endl << G4endl
1814 << "Proposed distance :" << G 1698 << "Proposed distance :" << G4endl << G4endl
1815 << "snxt = " << snxt/mm << 1699 << "snxt = " << snxt/mm << " mm" << G4endl ;
1816 message.precision(oldprc) ; 1700 message.precision(oldprc) ;
1817 G4Exception("G4CutTubs::DistanceToOut 1701 G4Exception("G4CutTubs::DistanceToOut(p,v,..)", "GeomSolids1002",
1818 JustWarning, message); 1702 JustWarning, message);
1819 break ; 1703 break ;
1820 } 1704 }
1821 } 1705 }
1822 if ( snxt<halfCarTolerance ) { snxt=0 ; } 1706 if ( snxt<halfCarTolerance ) { snxt=0 ; }
1823 return snxt ; 1707 return snxt ;
1824 } 1708 }
1825 1709
1826 ///////////////////////////////////////////// 1710 //////////////////////////////////////////////////////////////////////////
1827 // 1711 //
1828 // Calculate distance (<=actual) to closest s 1712 // Calculate distance (<=actual) to closest surface of shape from inside
1829 1713
1830 G4double G4CutTubs::DistanceToOut( const G4Th 1714 G4double G4CutTubs::DistanceToOut( const G4ThreeVector& p ) const
1831 { 1715 {
1832 G4double safRMin,safRMax,safZLow,safZHigh,s 1716 G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho;
1833 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 1717 G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
1834 1718
1835 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) 1719 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; // Distance to R
1836 1720
1837 safRMin = rho - fRMin ; 1721 safRMin = rho - fRMin ;
1838 safRMax = fRMax - rho ; 1722 safRMax = fRMax - rho ;
1839 1723
1840 // Distances to ZCut(Low/High) 1724 // Distances to ZCut(Low/High)
1841 1725
1842 // Dist to Low Cut 1726 // Dist to Low Cut
1843 // 1727 //
1844 safZLow = std::fabs((p+vZ).dot(fLowNorm)); 1728 safZLow = std::fabs((p+vZ).dot(fLowNorm));
1845 1729
1846 // Dist to High Cut 1730 // Dist to High Cut
1847 // 1731 //
1848 safZHigh = std::fabs((p-vZ).dot(fHighNorm)) 1732 safZHigh = std::fabs((p-vZ).dot(fHighNorm));
1849 safe = std::min(safZLow,safZHigh); 1733 safe = std::min(safZLow,safZHigh);
1850 1734
1851 if ( safRMin < safe ) { safe = safRMin; } 1735 if ( safRMin < safe ) { safe = safRMin; }
1852 if ( safRMax< safe ) { safe = safRMax; } 1736 if ( safRMax< safe ) { safe = safRMax; }
1853 1737
1854 // Check if phi divided, Calc distances clo 1738 // Check if phi divided, Calc distances closest phi plane
1855 // 1739 //
1856 if ( !fPhiFullCutTube ) 1740 if ( !fPhiFullCutTube )
1857 { 1741 {
1858 if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 ) 1742 if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 )
1859 { 1743 {
1860 safePhi = -(p.x()*sinSPhi - p.y()*cosSP 1744 safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
1861 } 1745 }
1862 else 1746 else
1863 { 1747 {
1864 safePhi = (p.x()*sinEPhi - p.y()*cosEPh 1748 safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
1865 } 1749 }
1866 if (safePhi < safe) { safe = safePhi ; } 1750 if (safePhi < safe) { safe = safePhi ; }
1867 } 1751 }
1868 if ( safe < 0 ) { safe = 0; } 1752 if ( safe < 0 ) { safe = 0; }
1869 1753
1870 return safe ; 1754 return safe ;
1871 } 1755 }
1872 1756
>> 1757 /////////////////////////////////////////////////////////////////////////
>> 1758 //
>> 1759 // Create a List containing the transformed vertices
>> 1760 // Ordering [0-3] -fDz cross section
>> 1761 // [4-7] +fDz cross section such that [0] is below [4],
>> 1762 // [1] below [5] etc.
>> 1763 // Note:
>> 1764 // Caller has deletion resposibility
>> 1765
>> 1766 G4ThreeVectorList*
>> 1767 G4CutTubs::CreateRotatedVertices( const G4AffineTransform& pTransform ) const
>> 1768 {
>> 1769 G4ThreeVectorList* vertices ;
>> 1770 G4ThreeVector vertex0, vertex1, vertex2, vertex3 ;
>> 1771 G4double meshAngle, meshRMax, crossAngle,
>> 1772 cosCrossAngle, sinCrossAngle, sAngle;
>> 1773 G4double rMaxX, rMaxY, rMinX, rMinY, meshRMin ;
>> 1774 G4int crossSection, noCrossSections;
>> 1775
>> 1776 // Compute no of cross-sections necessary to mesh tube
>> 1777 //
>> 1778 noCrossSections = G4int(fDPhi/kMeshAngleDefault) + 1 ;
>> 1779
>> 1780 if ( noCrossSections < kMinMeshSections )
>> 1781 {
>> 1782 noCrossSections = kMinMeshSections ;
>> 1783 }
>> 1784 else if (noCrossSections>kMaxMeshSections)
>> 1785 {
>> 1786 noCrossSections = kMaxMeshSections ;
>> 1787 }
>> 1788 // noCrossSections = 4 ;
>> 1789
>> 1790 meshAngle = fDPhi/(noCrossSections - 1) ;
>> 1791 // meshAngle = fDPhi/(noCrossSections) ;
>> 1792
>> 1793 meshRMax = (fRMax+100*kCarTolerance)/std::cos(meshAngle*0.5) ;
>> 1794 meshRMin = fRMin - 100*kCarTolerance ;
>> 1795
>> 1796 // If complete in phi, set start angle such that mesh will be at fRMax
>> 1797 // on the x axis. Will give better extent calculations when not rotated.
>> 1798
>> 1799 if (fPhiFullCutTube && (fSPhi == 0) ) { sAngle = -meshAngle*0.5 ; }
>> 1800 else { sAngle = fSPhi ; }
>> 1801
>> 1802 vertices = new G4ThreeVectorList();
>> 1803
>> 1804 if ( vertices )
>> 1805 {
>> 1806 vertices->reserve(noCrossSections*4);
>> 1807 for (crossSection = 0 ; crossSection < noCrossSections ; crossSection++ )
>> 1808 {
>> 1809 // Compute coordinates of cross section at section crossSection
>> 1810
>> 1811 crossAngle = sAngle + crossSection*meshAngle ;
>> 1812 cosCrossAngle = std::cos(crossAngle) ;
>> 1813 sinCrossAngle = std::sin(crossAngle) ;
>> 1814
>> 1815 rMaxX = meshRMax*cosCrossAngle ;
>> 1816 rMaxY = meshRMax*sinCrossAngle ;
>> 1817
>> 1818 if(meshRMin <= 0.0)
>> 1819 {
>> 1820 rMinX = 0.0 ;
>> 1821 rMinY = 0.0 ;
>> 1822 }
>> 1823 else
>> 1824 {
>> 1825 rMinX = meshRMin*cosCrossAngle ;
>> 1826 rMinY = meshRMin*sinCrossAngle ;
>> 1827 }
>> 1828 vertex0 = G4ThreeVector(rMinX,rMinY,GetCutZ(G4ThreeVector(rMinX,rMinY,-fDz))) ;
>> 1829 vertex1 = G4ThreeVector(rMaxX,rMaxY,GetCutZ(G4ThreeVector(rMaxX,rMaxY,-fDz))) ;
>> 1830 vertex2 = G4ThreeVector(rMaxX,rMaxY,GetCutZ(G4ThreeVector(rMaxX,rMaxY,+fDz))) ;
>> 1831 vertex3 = G4ThreeVector(rMinX,rMinY,GetCutZ(G4ThreeVector(rMinX,rMinY,+fDz))) ;
>> 1832
>> 1833 vertices->push_back(pTransform.TransformPoint(vertex0)) ;
>> 1834 vertices->push_back(pTransform.TransformPoint(vertex1)) ;
>> 1835 vertices->push_back(pTransform.TransformPoint(vertex2)) ;
>> 1836 vertices->push_back(pTransform.TransformPoint(vertex3)) ;
>> 1837 }
>> 1838 }
>> 1839 else
>> 1840 {
>> 1841 DumpInfo();
>> 1842 G4Exception("G4CutTubs::CreateRotatedVertices()",
>> 1843 "GeomSolids0003", FatalException,
>> 1844 "Error in allocation of vertices. Out of memory !");
>> 1845 }
>> 1846 return vertices ;
>> 1847 }
>> 1848
1873 ///////////////////////////////////////////// 1849 //////////////////////////////////////////////////////////////////////////
1874 // 1850 //
1875 // Stream object contents to an output stream 1851 // Stream object contents to an output stream
1876 1852
1877 G4GeometryType G4CutTubs::GetEntityType() con 1853 G4GeometryType G4CutTubs::GetEntityType() const
1878 { 1854 {
1879 return {"G4CutTubs"}; << 1855 return G4String("G4CutTubs");
1880 } 1856 }
1881 1857
1882 ///////////////////////////////////////////// 1858 //////////////////////////////////////////////////////////////////////////
1883 // 1859 //
1884 // Make a clone of the object 1860 // Make a clone of the object
1885 // 1861 //
1886 G4VSolid* G4CutTubs::Clone() const 1862 G4VSolid* G4CutTubs::Clone() const
1887 { 1863 {
1888 return new G4CutTubs(*this); 1864 return new G4CutTubs(*this);
1889 } 1865 }
1890 1866
1891 ///////////////////////////////////////////// 1867 //////////////////////////////////////////////////////////////////////////
1892 // 1868 //
1893 // Stream object contents to an output stream 1869 // Stream object contents to an output stream
1894 1870
1895 std::ostream& G4CutTubs::StreamInfo( std::ost 1871 std::ostream& G4CutTubs::StreamInfo( std::ostream& os ) const
1896 { 1872 {
1897 G4long oldprc = os.precision(16); << 1873 G4int oldprc = os.precision(16);
1898 os << "------------------------------------ 1874 os << "-----------------------------------------------------------\n"
1899 << " *** Dump for solid - " << GetNam 1875 << " *** Dump for solid - " << GetName() << " ***\n"
1900 << " ================================ 1876 << " ===================================================\n"
1901 << " Solid type: G4CutTubs\n" 1877 << " Solid type: G4CutTubs\n"
1902 << " Parameters: \n" 1878 << " Parameters: \n"
1903 << " inner radius : " << fRMin/mm << 1879 << " inner radius : " << fRMin/mm << " mm \n"
1904 << " outer radius : " << fRMax/mm << 1880 << " outer radius : " << fRMax/mm << " mm \n"
1905 << " half length Z: " << fDz/mm << " 1881 << " half length Z: " << fDz/mm << " mm \n"
1906 << " starting phi : " << fSPhi/degree 1882 << " starting phi : " << fSPhi/degree << " degrees \n"
1907 << " delta phi : " << fDPhi/degree 1883 << " delta phi : " << fDPhi/degree << " degrees \n"
1908 << " low Norm : " << fLowNorm << 1884 << " low Norm : " << fLowNorm << " \n"
1909 << " high Norm : " <<fHighNorm 1885 << " high Norm : " <<fHighNorm << " \n"
1910 << "------------------------------------ 1886 << "-----------------------------------------------------------\n";
1911 os.precision(oldprc); 1887 os.precision(oldprc);
1912 1888
1913 return os; 1889 return os;
1914 } 1890 }
1915 1891
1916 ///////////////////////////////////////////// 1892 /////////////////////////////////////////////////////////////////////////
1917 // 1893 //
1918 // GetPointOnSurface 1894 // GetPointOnSurface
1919 1895
1920 G4ThreeVector G4CutTubs::GetPointOnSurface() 1896 G4ThreeVector G4CutTubs::GetPointOnSurface() const
1921 { 1897 {
1922 // Set min and max z << 1898 G4double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1923 if (fZMin == 0. && fZMax == 0.) << 1899 aOne, aTwo, aThr, aFou;
>> 1900 G4double rRand;
>> 1901
>> 1902 aOne = 2.*fDz*fDPhi*fRMax;
>> 1903 aTwo = 2.*fDz*fDPhi*fRMin;
>> 1904 aThr = 0.5*fDPhi*(fRMax*fRMax-fRMin*fRMin);
>> 1905 aFou = 2.*fDz*(fRMax-fRMin);
>> 1906
>> 1907 phi = RandFlat::shoot(fSPhi, fSPhi+fDPhi);
>> 1908 cosphi = std::cos(phi);
>> 1909 sinphi = std::sin(phi);
>> 1910
>> 1911 rRand = GetRadiusInRing(fRMin,fRMax);
>> 1912
>> 1913 if( (fSPhi == 0) && (fDPhi == twopi) ) { aFou = 0; }
>> 1914
>> 1915 chose = RandFlat::shoot(0.,aOne+aTwo+2.*aThr+2.*aFou);
>> 1916
>> 1917 if( (chose >=0) && (chose < aOne) )
>> 1918 {
>> 1919 xRand = fRMax*cosphi;
>> 1920 yRand = fRMax*sinphi;
>> 1921 zRand = RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
>> 1922 GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
>> 1923 return G4ThreeVector (xRand, yRand, zRand);
>> 1924 }
>> 1925 else if( (chose >= aOne) && (chose < aOne + aTwo) )
>> 1926 {
>> 1927 xRand = fRMin*cosphi;
>> 1928 yRand = fRMin*sinphi;
>> 1929 zRand = RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
>> 1930 GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
>> 1931 return G4ThreeVector (xRand, yRand, zRand);
>> 1932 }
>> 1933 else if( (chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr) )
>> 1934 {
>> 1935 xRand = rRand*cosphi;
>> 1936 yRand = rRand*sinphi;
>> 1937 zRand = GetCutZ(G4ThreeVector(xRand,yRand,fDz));
>> 1938 return G4ThreeVector (xRand, yRand, zRand);
>> 1939 }
>> 1940 else if( (chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr) )
>> 1941 {
>> 1942 xRand = rRand*cosphi;
>> 1943 yRand = rRand*sinphi;
>> 1944 zRand = GetCutZ(G4ThreeVector(xRand,yRand,-fDz));
>> 1945 return G4ThreeVector (xRand, yRand, zRand);
>> 1946 }
>> 1947 else if( (chose >= aOne + aTwo + 2.*aThr)
>> 1948 && (chose < aOne + aTwo + 2.*aThr + aFou) )
>> 1949 {
>> 1950 xRand = rRand*std::cos(fSPhi);
>> 1951 yRand = rRand*std::sin(fSPhi);
>> 1952 zRand = RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
>> 1953 GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
>> 1954 return G4ThreeVector (xRand, yRand, zRand);
>> 1955 }
>> 1956 else
1924 { 1957 {
1925 G4AutoLock l(&zminmaxMutex); << 1958 xRand = rRand*std::cos(fSPhi+fDPhi);
1926 G4ThreeVector bmin, bmax; << 1959 yRand = rRand*std::sin(fSPhi+fDPhi);
1927 BoundingLimits(bmin,bmax); << 1960 zRand = RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)),
1928 fZMin = bmin.z(); << 1961 GetCutZ(G4ThreeVector(xRand,yRand,fDz)));
1929 fZMax = bmax.z(); << 1962 return G4ThreeVector (xRand, yRand, zRand);
1930 l.unlock(); << 1963 }
1931 } <<
1932 <<
1933 // Set parameters <<
1934 G4double hmax = fZMax - fZMin; <<
1935 G4double sphi = fSPhi; <<
1936 G4double dphi = fDPhi; <<
1937 G4double rmin = fRMin; <<
1938 G4double rmax = fRMax; <<
1939 G4double rrmax = rmax*rmax; <<
1940 G4double rrmin = rmin*rmin; <<
1941 <<
1942 G4ThreeVector nbot = GetLowNorm(); <<
1943 G4ThreeVector ntop = GetHighNorm(); <<
1944 <<
1945 // Set array of surface areas <<
1946 G4double sbase = 0.5*dphi*(rrmax - rrmin); <<
1947 G4double sbot = sbase/std::abs(nbot.z()); <<
1948 G4double stop = sbase/std::abs(ntop.z()); <<
1949 G4double scut = (dphi == twopi) ? 0. : hmax <<
1950 G4double ssurf[6] = { scut, scut, sbot, sto <<
1951 ssurf[1] += ssurf[0]; <<
1952 ssurf[2] += ssurf[1]; <<
1953 ssurf[3] += ssurf[2]; <<
1954 ssurf[4] += ssurf[3]; <<
1955 ssurf[5] += ssurf[4]; <<
1956 <<
1957 constexpr G4int ntry = 100000; <<
1958 for (G4int i=0; i<ntry; ++i) <<
1959 { <<
1960 // Select surface <<
1961 G4double select = ssurf[5]*G4QuickRand(); <<
1962 G4int k = 5; <<
1963 k -= (G4int)(select <= ssurf[4]); <<
1964 k -= (G4int)(select <= ssurf[3]); <<
1965 k -= (G4int)(select <= ssurf[2]); <<
1966 k -= (G4int)(select <= ssurf[1]); <<
1967 k -= (G4int)(select <= ssurf[0]); <<
1968 <<
1969 // Generate point on selected surface (re <<
1970 G4ThreeVector p(0,0,0); <<
1971 switch(k) <<
1972 { <<
1973 case 0: // cut at start phi <<
1974 { <<
1975 G4double r = rmin + (rmax - rmin)*G4Q <<
1976 p.set(r*cosSPhi, r*sinSPhi, fZMin + h <<
1977 break; <<
1978 } <<
1979 case 1: // cut at end phi <<
1980 { <<
1981 G4double r = rmin + (rmax - rmin)*G4Q <<
1982 p.set(r*cosEPhi, r*sinEPhi, fZMin + h <<
1983 break; <<
1984 } <<
1985 case 2: // base at low z <<
1986 { <<
1987 G4double r = std::sqrt(rrmin + (rrmax <<
1988 G4double phi = sphi + dphi*G4QuickRan <<
1989 G4double x = r*std::cos(phi); <<
1990 G4double y = r*std::sin(phi); <<
1991 G4double z = -fDz - (x*nbot.x() + y*n <<
1992 return {x, y, z}; <<
1993 } <<
1994 case 3: // base at high z <<
1995 { <<
1996 G4double r = std::sqrt(rrmin + (rrmax <<
1997 G4double phi = sphi + dphi*G4QuickRan <<
1998 G4double x = r*std::cos(phi); <<
1999 G4double y = r*std::sin(phi); <<
2000 G4double z = fDz - (x*ntop.x() + y*nt <<
2001 return {x, y, z}; <<
2002 } <<
2003 case 4: // external lateral surface <<
2004 { <<
2005 G4double phi = sphi + dphi*G4QuickRan <<
2006 G4double z = fZMin + hmax*G4QuickRand <<
2007 G4double x = rmax*std::cos(phi); <<
2008 G4double y = rmax*std::sin(phi); <<
2009 p.set(x, y, z); <<
2010 break; <<
2011 } <<
2012 case 5: // internal lateral surface <<
2013 { <<
2014 G4double phi = sphi + dphi*G4QuickRan <<
2015 G4double z = fZMin + hmax*G4QuickRand <<
2016 G4double x = rmin*std::cos(phi); <<
2017 G4double y = rmin*std::sin(phi); <<
2018 p.set(x, y, z); <<
2019 break; <<
2020 } <<
2021 } <<
2022 if ((ntop.dot(p) - fDz*ntop.z()) > 0.) co <<
2023 if ((nbot.dot(p) + fDz*nbot.z()) > 0.) co <<
2024 return p; <<
2025 } <<
2026 // Just in case, if all attempts to generat <<
2027 // Normally should never happen <<
2028 G4double x = rmax*std::cos(sphi + 0.5*dphi) <<
2029 G4double y = rmax*std::sin(sphi + 0.5*dphi) <<
2030 G4double z = fDz - (x*ntop.x() + y*ntop.y() <<
2031 return {x, y, z}; <<
2032 } 1964 }
2033 1965
2034 ///////////////////////////////////////////// 1966 ///////////////////////////////////////////////////////////////////////////
2035 // 1967 //
2036 // Methods for visualisation 1968 // Methods for visualisation
2037 1969
2038 void G4CutTubs::DescribeYourselfTo ( G4VGraph << 1970 void G4CutTubs::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
2039 { 1971 {
2040 scene.AddSolid (*this) ; 1972 scene.AddSolid (*this) ;
2041 } 1973 }
2042 1974
2043 G4Polyhedron* G4CutTubs::CreatePolyhedron () << 1975 G4Polyhedron* G4CutTubs::CreatePolyhedron () const
2044 { 1976 {
2045 typedef G4double G4double3[3]; 1977 typedef G4double G4double3[3];
2046 typedef G4int G4int4[4]; 1978 typedef G4int G4int4[4];
2047 1979
2048 auto ph = new G4Polyhedron; << 1980 G4Polyhedron *ph = new G4Polyhedron;
2049 G4Polyhedron *ph1 = new G4PolyhedronTubs (f << 1981 G4Polyhedron *ph1 = G4Tubs::CreatePolyhedron();
2050 G4int nn=ph1->GetNoVertices(); 1982 G4int nn=ph1->GetNoVertices();
2051 G4int nf=ph1->GetNoFacets(); 1983 G4int nf=ph1->GetNoFacets();
2052 auto xyz = new G4double3[nn]; // number of << 1984 G4double3* xyz = new G4double3[nn]; // number of nodes
2053 auto faces = new G4int4[nf] ; // number of << 1985 G4int4* faces = new G4int4[nf] ; // number of faces
2054 1986
2055 for(G4int i=0; i<nn; ++i) << 1987 for(G4int i=0;i<nn;i++)
2056 { 1988 {
2057 xyz[i][0]=ph1->GetVertex(i+1).x(); 1989 xyz[i][0]=ph1->GetVertex(i+1).x();
2058 xyz[i][1]=ph1->GetVertex(i+1).y(); 1990 xyz[i][1]=ph1->GetVertex(i+1).y();
2059 G4double tmpZ=ph1->GetVertex(i+1).z(); 1991 G4double tmpZ=ph1->GetVertex(i+1).z();
2060 if(tmpZ>=fDz-kCarTolerance) 1992 if(tmpZ>=fDz-kCarTolerance)
2061 { 1993 {
2062 xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][ 1994 xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][0],xyz[i][1],fDz));
2063 } 1995 }
2064 else if(tmpZ<=-fDz+kCarTolerance) 1996 else if(tmpZ<=-fDz+kCarTolerance)
2065 { 1997 {
2066 xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][ 1998 xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][0],xyz[i][1],-fDz));
2067 } 1999 }
2068 else 2000 else
2069 { 2001 {
2070 xyz[i][2]=tmpZ; 2002 xyz[i][2]=tmpZ;
2071 } 2003 }
2072 } 2004 }
2073 G4int iNodes[4]; 2005 G4int iNodes[4];
2074 G4int* iEdge = nullptr; << 2006 G4int *iEdge=0;
2075 G4int n; 2007 G4int n;
2076 for(G4int i=0; i<nf ; ++i) << 2008 for(G4int i=0;i<nf;i++)
2077 { 2009 {
2078 ph1->GetFacet(i+1,n,iNodes,iEdge); 2010 ph1->GetFacet(i+1,n,iNodes,iEdge);
2079 for(G4int k=0; k<n; ++k) << 2011 for(G4int k=0;k<n;k++)
2080 { 2012 {
2081 faces[i][k]=iNodes[k]; 2013 faces[i][k]=iNodes[k];
2082 } 2014 }
2083 for(G4int k=n; k<4; ++k) << 2015 for(G4int k=n;k<4;k++)
2084 { 2016 {
2085 faces[i][k]=0; 2017 faces[i][k]=0;
2086 } 2018 }
2087 } 2019 }
2088 ph->createPolyhedron(nn,nf,xyz,faces); 2020 ph->createPolyhedron(nn,nf,xyz,faces);
2089 2021
2090 delete [] xyz; 2022 delete [] xyz;
2091 delete [] faces; 2023 delete [] faces;
2092 delete ph1; 2024 delete ph1;
2093 2025
2094 return ph; 2026 return ph;
2095 } 2027 }
2096 2028
2097 // Auxilary Methods for Solid 2029 // Auxilary Methods for Solid
2098 << 2030
2099 ///////////////////////////////////////////// << 2031 ///////////////////////////////////////////////////////////////////////////
2100 // << 2032 // Return true if Cutted planes are crossing
2101 // Check set of points on the outer lateral s << 2033 // Check Intersection Points on OX and OY axes
2102 // if the cut planes are crossing inside the <<
2103 // <<
2104 2034
2105 G4bool G4CutTubs::IsCrossingCutPlanes() const 2035 G4bool G4CutTubs::IsCrossingCutPlanes() const
2106 { 2036 {
2107 constexpr G4int npoints = 30; << 2037 G4double zXLow1,zXLow2,zYLow1,zYLow2;
>> 2038 G4double zXHigh1,zXHigh2,zYHigh1,zYHigh2;
>> 2039
>> 2040 zXLow1 = GetCutZ(G4ThreeVector(-fRMax, 0,-fDz));
>> 2041 zXLow2 = GetCutZ(G4ThreeVector( fRMax, 0,-fDz));
>> 2042 zYLow1 = GetCutZ(G4ThreeVector( 0,-fRMax,-fDz));
>> 2043 zYLow2 = GetCutZ(G4ThreeVector( 0, fRMax,-fDz));
>> 2044 zXHigh1 = GetCutZ(G4ThreeVector(-fRMax, 0, fDz));
>> 2045 zXHigh2 = GetCutZ(G4ThreeVector( fRMax, 0, fDz));
>> 2046 zYHigh1 = GetCutZ(G4ThreeVector( 0,-fRMax, fDz));
>> 2047 zYHigh2 = GetCutZ(G4ThreeVector( 0, fRMax, fDz));
>> 2048 if ( (zXLow1>zXHigh1) ||(zXLow2>zXHigh2)
>> 2049 || (zYLow1>zYHigh1) ||(zYLow2>zYHigh2)) { return true; }
2108 2050
2109 // set values for calculation of h - distan <<
2110 // opposite points on bases <<
2111 G4ThreeVector nbot = GetLowNorm(); <<
2112 G4ThreeVector ntop = GetHighNorm(); <<
2113 if (std::abs(nbot.z()) < kCarTolerance) ret <<
2114 if (std::abs(ntop.z()) < kCarTolerance) ret <<
2115 G4double nx = nbot.x()/nbot.z() - ntop.x()/ <<
2116 G4double ny = nbot.y()/nbot.z() - ntop.y()/ <<
2117 <<
2118 // check points <<
2119 G4double cosphi = GetCosStartPhi(); <<
2120 G4double sinphi = GetSinStartPhi(); <<
2121 G4double delphi = GetDeltaPhiAngle()/npoint <<
2122 G4double cosdel = std::cos(delphi); <<
2123 G4double sindel = std::sin(delphi); <<
2124 G4double hzero = 2.*GetZHalfLength()/GetOut <<
2125 for (G4int i=0; i<npoints+1; ++i) <<
2126 { <<
2127 G4double h = nx*cosphi + ny*sinphi + hzer <<
2128 if (h < 0.) return true; <<
2129 G4double sintmp = sinphi; <<
2130 sinphi = sintmp*cosdel + cosphi*sindel; <<
2131 cosphi = cosphi*cosdel - sintmp*sindel; <<
2132 } <<
2133 return false; 2051 return false;
2134 } 2052 }
2135 2053
2136 ///////////////////////////////////////////// 2054 ///////////////////////////////////////////////////////////////////////////
2137 // 2055 //
2138 // Return real Z coordinate of point on Cutte 2056 // Return real Z coordinate of point on Cutted +/- fDZ plane
2139 2057
2140 G4double G4CutTubs::GetCutZ(const G4ThreeVect 2058 G4double G4CutTubs::GetCutZ(const G4ThreeVector& p) const
2141 { 2059 {
2142 G4double newz = p.z(); // p.z() should be 2060 G4double newz = p.z(); // p.z() should be either +fDz or -fDz
2143 if (p.z()<0) 2061 if (p.z()<0)
2144 { 2062 {
2145 if(fLowNorm.z()!=0.) 2063 if(fLowNorm.z()!=0.)
2146 { 2064 {
2147 newz = -fDz-(p.x()*fLowNorm.x()+p.y()* 2065 newz = -fDz-(p.x()*fLowNorm.x()+p.y()*fLowNorm.y())/fLowNorm.z();
2148 } 2066 }
2149 } 2067 }
2150 else 2068 else
2151 { 2069 {
2152 if(fHighNorm.z()!=0.) 2070 if(fHighNorm.z()!=0.)
2153 { 2071 {
2154 newz = fDz-(p.x()*fHighNorm.x()+p.y()* 2072 newz = fDz-(p.x()*fHighNorm.x()+p.y()*fHighNorm.y())/fHighNorm.z();
2155 } 2073 }
2156 } 2074 }
2157 return newz; 2075 return newz;
2158 } 2076 }
2159 #endif << 2077
>> 2078 ///////////////////////////////////////////////////////////////////////////
>> 2079 //
>> 2080 // Calculate Min and Max Z for CutZ
>> 2081
>> 2082 void G4CutTubs::GetMaxMinZ(G4double& zmin,G4double& zmax)const
>> 2083
>> 2084 {
>> 2085 G4double phiLow = std::atan2(fLowNorm.y(),fLowNorm.x());
>> 2086 G4double phiHigh= std::atan2(fHighNorm.y(),fHighNorm.x());
>> 2087
>> 2088 G4double xc=0, yc=0,z1;
>> 2089 G4double z[8];
>> 2090 G4bool in_range_low = false;
>> 2091 G4bool in_range_hi = false;
>> 2092
>> 2093 G4int i;
>> 2094 for (i=0; i<2; i++)
>> 2095 {
>> 2096 if (phiLow<0) { phiLow+=twopi; }
>> 2097 G4double ddp = phiLow-fSPhi;
>> 2098 if (ddp<0) { ddp += twopi; }
>> 2099 if (ddp <= fDPhi)
>> 2100 {
>> 2101 xc = fRMin*std::cos(phiLow);
>> 2102 yc = fRMin*std::sin(phiLow);
>> 2103 z1 = GetCutZ(G4ThreeVector(xc, yc, -fDz));
>> 2104 xc = fRMax*std::cos(phiLow);
>> 2105 yc = fRMax*std::sin(phiLow);
>> 2106 z1 = std::min(z1, GetCutZ(G4ThreeVector(xc, yc, -fDz)));
>> 2107 if (in_range_low) { zmin = std::min(zmin, z1); }
>> 2108 else { zmin = z1; }
>> 2109 in_range_low = true;
>> 2110 }
>> 2111 phiLow += pi;
>> 2112 if (phiLow>twopi) { phiLow-=twopi; }
>> 2113 }
>> 2114 for (i=0; i<2; i++)
>> 2115 {
>> 2116 if (phiHigh<0) { phiHigh+=twopi; }
>> 2117 G4double ddp = phiHigh-fSPhi;
>> 2118 if (ddp<0) { ddp += twopi; }
>> 2119 if (ddp <= fDPhi)
>> 2120 {
>> 2121 xc = fRMin*std::cos(phiHigh);
>> 2122 yc = fRMin*std::sin(phiHigh);
>> 2123 z1 = GetCutZ(G4ThreeVector(xc, yc, fDz));
>> 2124 xc = fRMax*std::cos(phiHigh);
>> 2125 yc = fRMax*std::sin(phiHigh);
>> 2126 z1 = std::min(z1, GetCutZ(G4ThreeVector(xc, yc, fDz)));
>> 2127 if (in_range_hi) { zmax = std::min(zmax, z1); }
>> 2128 else { zmax = z1; }
>> 2129 in_range_hi = true;
>> 2130 }
>> 2131 phiHigh += pi;
>> 2132 if (phiLow>twopi) { phiHigh-=twopi; }
>> 2133 }
>> 2134
>> 2135 xc = fRMin*std::cos(fSPhi);
>> 2136 yc = fRMin*std::sin(fSPhi);
>> 2137 z[0] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
>> 2138 z[4] = GetCutZ(G4ThreeVector(xc, yc, fDz));
>> 2139
>> 2140 xc = fRMin*std::cos(fSPhi+fDPhi);
>> 2141 yc = fRMin*std::sin(fSPhi+fDPhi);
>> 2142 z[1] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
>> 2143 z[5] = GetCutZ(G4ThreeVector(xc, yc, fDz));
>> 2144
>> 2145 xc = fRMax*std::cos(fSPhi);
>> 2146 yc = fRMax*std::sin(fSPhi);
>> 2147 z[2] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
>> 2148 z[6] = GetCutZ(G4ThreeVector(xc, yc, fDz));
>> 2149
>> 2150 xc = fRMax*std::cos(fSPhi+fDPhi);
>> 2151 yc = fRMax*std::sin(fSPhi+fDPhi);
>> 2152 z[3] = GetCutZ(G4ThreeVector(xc, yc, -fDz));
>> 2153 z[7] = GetCutZ(G4ThreeVector(xc, yc, fDz));
>> 2154
>> 2155 // Find min/max
>> 2156
>> 2157 z1=z[0];
>> 2158 for (i = 1; i < 4; i++)
>> 2159 {
>> 2160 if(z[i] < z[i-1])z1=z[i];
>> 2161 }
>> 2162
>> 2163 if (in_range_low)
>> 2164 {
>> 2165 zmin = std::min(zmin, z1);
>> 2166 }
>> 2167 else
>> 2168 {
>> 2169 zmin = z1;
>> 2170 }
>> 2171 z1=z[4];
>> 2172 for (i = 1; i < 4; i++)
>> 2173 {
>> 2174 if(z[4+i] > z[4+i-1]) { z1=z[4+i]; }
>> 2175 }
>> 2176
>> 2177 if (in_range_hi) { zmax = std::max(zmax, z1); }
>> 2178 else { zmax = z1; }
>> 2179 }
2160 2180