| 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 // G4TwistTubsHypeSide implementation 26 // G4TwistTubsHypeSide implementation
27 // 27 //
28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepbu 28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp), created.
29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), 29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
30 // from original version in Jupi 30 // from original version in Jupiter-2.5.02 application.
31 // ------------------------------------------- 31 // --------------------------------------------------------------------
32 32
33 #include "G4TwistTubsHypeSide.hh" 33 #include "G4TwistTubsHypeSide.hh"
34 #include "G4PhysicalConstants.hh" 34 #include "G4PhysicalConstants.hh"
35 #include "G4GeometryTolerance.hh" 35 #include "G4GeometryTolerance.hh"
36 36
37 //============================================ 37 //=====================================================================
38 //* constructors ----------------------------- 38 //* constructors ------------------------------------------------------
39 39
40 G4TwistTubsHypeSide::G4TwistTubsHypeSide(const 40 G4TwistTubsHypeSide::G4TwistTubsHypeSide(const G4String& name,
41 const 41 const G4RotationMatrix& rot,
42 const 42 const G4ThreeVector& tlate,
43 const 43 const G4int handedness,
44 const 44 const G4double kappa,
45 const 45 const G4double tanstereo,
46 const 46 const G4double r0,
47 const 47 const EAxis axis0,
48 const 48 const EAxis axis1,
49 49 G4double axis0min,
50 50 G4double axis1min,
51 51 G4double axis0max,
52 52 G4double axis1max )
53 : G4VTwistSurface(name, rot, tlate, handedne 53 : G4VTwistSurface(name, rot, tlate, handedness, axis0, axis1,
54 axis0min, axis1min, axis0m 54 axis0min, axis1min, axis0max, axis1max),
55 fKappa(kappa), fTanStereo(tanstereo), 55 fKappa(kappa), fTanStereo(tanstereo),
56 fTan2Stereo(tanstereo*tanstereo), fR0(r0), 56 fTan2Stereo(tanstereo*tanstereo), fR0(r0), fR02(r0*r0), fDPhi(twopi)
57 { 57 {
58 if ( (axis0 == kZAxis) && (axis1 == kPhi) ) 58 if ( (axis0 == kZAxis) && (axis1 == kPhi) )
59 { 59 {
60 G4Exception("G4TwistTubsHypeSide::G4Twis 60 G4Exception("G4TwistTubsHypeSide::G4TwistTubsHypeSide()",
61 "GeomSolids0002", FatalError 61 "GeomSolids0002", FatalErrorInArgument,
62 "Should swap axis0 and axis1 62 "Should swap axis0 and axis1!");
63 } 63 }
64 fInside.gp.set(kInfinity, kInfinity, kInfin 64 fInside.gp.set(kInfinity, kInfinity, kInfinity);
65 fInside.inside = kOutside; 65 fInside.inside = kOutside;
66 fIsValidNorm = false; 66 fIsValidNorm = false;
67 SetCorners(); 67 SetCorners();
68 SetBoundaries(); 68 SetBoundaries();
69 } 69 }
70 70
71 G4TwistTubsHypeSide::G4TwistTubsHypeSide(const 71 G4TwistTubsHypeSide::G4TwistTubsHypeSide(const G4String& name,
72 72 G4double EndInnerRadius[2],
73 73 G4double EndOuterRadius[2],
74 74 G4double DPhi,
75 75 G4double EndPhi[2],
76 76 G4double EndZ[2],
77 77 G4double InnerRadius,
78 78 G4double OuterRadius,
79 79 G4double Kappa,
80 80 G4double TanInnerStereo,
81 81 G4double TanOuterStereo,
82 82 G4int handedness)
83 : G4VTwistSurface(name) 83 : G4VTwistSurface(name)
84 { 84 {
85 85
86 fHandedness = handedness; // +z = +ve, -z 86 fHandedness = handedness; // +z = +ve, -z = -ve
87 fAxis[0] = kPhi; 87 fAxis[0] = kPhi;
88 fAxis[1] = kZAxis; 88 fAxis[1] = kZAxis;
89 fAxisMin[0] = kInfinity; // we cann 89 fAxisMin[0] = kInfinity; // we cannot fix boundary min of Phi,
90 fAxisMax[0] = kInfinity; // because 90 fAxisMax[0] = kInfinity; // because it depends on z.
91 fAxisMin[1] = EndZ[0]; 91 fAxisMin[1] = EndZ[0];
92 fAxisMax[1] = EndZ[1]; 92 fAxisMax[1] = EndZ[1];
93 fKappa = Kappa; 93 fKappa = Kappa;
94 fDPhi = DPhi ; 94 fDPhi = DPhi ;
95 95
96 if (handedness < 0) // inner hyperbolic su 96 if (handedness < 0) // inner hyperbolic surface
97 { 97 {
98 fTanStereo = TanInnerStereo; 98 fTanStereo = TanInnerStereo;
99 fR0 = InnerRadius; 99 fR0 = InnerRadius;
100 } 100 }
101 else // outer hyperbolic su 101 else // outer hyperbolic surface
102 { 102 {
103 fTanStereo = TanOuterStereo; 103 fTanStereo = TanOuterStereo;
104 fR0 = OuterRadius; 104 fR0 = OuterRadius;
105 } 105 }
106 fTan2Stereo = fTanStereo * fTanStereo; 106 fTan2Stereo = fTanStereo * fTanStereo;
107 fR02 = fR0 * fR0; 107 fR02 = fR0 * fR0;
108 108
109 fTrans.set(0, 0, 0); 109 fTrans.set(0, 0, 0);
110 fIsValidNorm = false; 110 fIsValidNorm = false;
111 111
112 fInside.gp.set(kInfinity, kInfinity, kInfin 112 fInside.gp.set(kInfinity, kInfinity, kInfinity);
113 fInside.inside = kOutside; 113 fInside.inside = kOutside;
114 114
115 SetCorners(EndInnerRadius, EndOuterRadius, 115 SetCorners(EndInnerRadius, EndOuterRadius, DPhi, EndPhi, EndZ) ;
116 116
117 SetBoundaries(); 117 SetBoundaries();
118 } 118 }
119 119
120 //============================================ 120 //=====================================================================
121 //* Fake default constructor ----------------- 121 //* Fake default constructor ------------------------------------------
122 122
123 G4TwistTubsHypeSide::G4TwistTubsHypeSide( __vo 123 G4TwistTubsHypeSide::G4TwistTubsHypeSide( __void__& a )
124 : G4VTwistSurface(a) << 124 : G4VTwistSurface(a), fKappa(0.), fTanStereo(0.), fTan2Stereo(0.),
>> 125 fR0(0.), fR02(0.), fDPhi(0.)
125 { 126 {
126 } 127 }
127 128
128 //============================================ 129 //=====================================================================
129 //* destructor ------------------------------- 130 //* destructor --------------------------------------------------------
130 131
131 G4TwistTubsHypeSide::~G4TwistTubsHypeSide() = << 132 G4TwistTubsHypeSide::~G4TwistTubsHypeSide()
>> 133 {
>> 134 }
132 135
133 //============================================ 136 //=====================================================================
134 //* GetNormal -------------------------------- 137 //* GetNormal ---------------------------------------------------------
135 138
136 G4ThreeVector G4TwistTubsHypeSide::GetNormal(c 139 G4ThreeVector G4TwistTubsHypeSide::GetNormal(const G4ThreeVector& tmpxx,
137 140 G4bool isGlobal)
138 { 141 {
139 // GetNormal returns a normal vector at a s 142 // GetNormal returns a normal vector at a surface (or very close
140 // to surface) point at tmpxx. 143 // to surface) point at tmpxx.
141 // If isGlobal=true, it returns the normal 144 // If isGlobal=true, it returns the normal in global coordinate.
142 145
143 G4ThreeVector xx; 146 G4ThreeVector xx;
144 if (isGlobal) 147 if (isGlobal)
145 { 148 {
146 xx = ComputeLocalPoint(tmpxx); 149 xx = ComputeLocalPoint(tmpxx);
147 if ((xx - fCurrentNormal.p).mag() < 0.5 150 if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance)
148 { 151 {
149 return ComputeGlobalDirection(fCurren 152 return ComputeGlobalDirection(fCurrentNormal.normal);
150 } 153 }
151 } 154 }
152 else 155 else
153 { 156 {
154 xx = tmpxx; 157 xx = tmpxx;
155 if (xx == fCurrentNormal.p) 158 if (xx == fCurrentNormal.p)
156 { 159 {
157 return fCurrentNormal.normal; 160 return fCurrentNormal.normal;
158 } 161 }
159 } 162 }
160 163
161 fCurrentNormal.p = xx; 164 fCurrentNormal.p = xx;
162 165
163 G4ThreeVector normal( xx.x(), xx.y(), -xx.z 166 G4ThreeVector normal( xx.x(), xx.y(), -xx.z() * fTan2Stereo);
164 normal *= fHandedness; 167 normal *= fHandedness;
165 normal = normal.unit(); 168 normal = normal.unit();
166 169
167 if (isGlobal) 170 if (isGlobal)
168 { 171 {
169 fCurrentNormal.normal = ComputeLocalDire 172 fCurrentNormal.normal = ComputeLocalDirection(normal);
170 } 173 }
171 else 174 else
172 { 175 {
173 fCurrentNormal.normal = normal; 176 fCurrentNormal.normal = normal;
174 } 177 }
175 return fCurrentNormal.normal; 178 return fCurrentNormal.normal;
176 } 179 }
177 180
178 //============================================ 181 //=====================================================================
179 //* Inside() --------------------------------- 182 //* Inside() ----------------------------------------------------------
180 183
181 EInside G4TwistTubsHypeSide::Inside(const G4Th 184 EInside G4TwistTubsHypeSide::Inside(const G4ThreeVector& gp)
182 { 185 {
183 // Inside returns 186 // Inside returns
184 const G4double halftol 187 const G4double halftol
185 = 0.5 * G4GeometryTolerance::GetInstance( 188 = 0.5 * G4GeometryTolerance::GetInstance()->GetRadialTolerance();
186 189
187 if (fInside.gp == gp) 190 if (fInside.gp == gp)
188 { 191 {
189 return fInside.inside; 192 return fInside.inside;
190 } 193 }
191 fInside.gp = gp; 194 fInside.gp = gp;
192 195
193 G4ThreeVector p = ComputeLocalPoint(gp); 196 G4ThreeVector p = ComputeLocalPoint(gp);
194 197
195 198
196 if (p.mag() < DBL_MIN) 199 if (p.mag() < DBL_MIN)
197 { 200 {
198 fInside.inside = kOutside; 201 fInside.inside = kOutside;
199 return fInside.inside; 202 return fInside.inside;
200 } 203 }
201 204
202 G4double rhohype = GetRhoAtPZ(p); 205 G4double rhohype = GetRhoAtPZ(p);
203 G4double distanceToOut = fHandedness * (rho 206 G4double distanceToOut = fHandedness * (rhohype - p.getRho());
204 // +ve : inside 207 // +ve : inside
205 208
206 if (distanceToOut < -halftol) 209 if (distanceToOut < -halftol)
207 { 210 {
208 fInside.inside = kOutside; 211 fInside.inside = kOutside;
209 } 212 }
210 else 213 else
211 { 214 {
212 G4int areacode = GetAreaCode(p); 215 G4int areacode = GetAreaCode(p);
213 if (IsOutside(areacode)) 216 if (IsOutside(areacode))
214 { 217 {
215 fInside.inside = kOutside; 218 fInside.inside = kOutside;
216 } 219 }
217 else if (IsBoundary(areacode)) 220 else if (IsBoundary(areacode))
218 { 221 {
219 fInside.inside = kSurface; 222 fInside.inside = kSurface;
220 } 223 }
221 else if (IsInside(areacode)) 224 else if (IsInside(areacode))
222 { 225 {
223 if (distanceToOut <= halftol) 226 if (distanceToOut <= halftol)
224 { 227 {
225 fInside.inside = kSurface; 228 fInside.inside = kSurface;
226 } 229 }
227 else 230 else
228 { 231 {
229 fInside.inside = kInside; 232 fInside.inside = kInside;
230 } 233 }
231 } 234 }
232 else 235 else
233 { 236 {
234 G4cout << "WARNING - G4TwistTubsHypeS 237 G4cout << "WARNING - G4TwistTubsHypeSide::Inside()" << G4endl
235 << " Invalid option ! 238 << " Invalid option !" << G4endl
236 << " name, areacode, 239 << " name, areacode, distanceToOut = "
237 << GetName() << ", " << std::h 240 << GetName() << ", " << std::hex << areacode
238 << std::dec << ", " << distanc 241 << std::dec << ", " << distanceToOut << G4endl;
239 } 242 }
240 } 243 }
241 244
242 return fInside.inside; 245 return fInside.inside;
243 } 246 }
244 247
245 //============================================ 248 //=====================================================================
246 //* DistanceToSurface ------------------------ 249 //* DistanceToSurface -------------------------------------------------
247 250
248 G4int G4TwistTubsHypeSide::DistanceToSurface(c 251 G4int G4TwistTubsHypeSide::DistanceToSurface(const G4ThreeVector& gp,
249 c 252 const G4ThreeVector& gv,
250 253 G4ThreeVector gxx[],
251 254 G4double distance[],
252 255 G4int areacode[],
253 256 G4bool isvalid[],
254 257 EValidate validate)
255 { 258 {
256 // Decide if and where a line intersects wi 259 // Decide if and where a line intersects with a hyperbolic
257 // surface (of infinite extent) 260 // surface (of infinite extent)
258 // 261 //
259 // Arguments: 262 // Arguments:
260 // p - (in) Point on trajectory 263 // p - (in) Point on trajectory
261 // v - (in) Vector along trajecto 264 // v - (in) Vector along trajectory
262 // r2 - (in) Square of radius at z 265 // r2 - (in) Square of radius at z = 0
263 // tan2phi - (in) std::tan(stereo)**2 266 // tan2phi - (in) std::tan(stereo)**2
264 // s - (out) Up to two points of 267 // s - (out) Up to two points of intersection, where the
265 // intersection point i 268 // intersection point is p + s*v, and if there are
266 // two intersections, s 269 // two intersections, s[0] < s[1]. May be negative.
267 // Returns: 270 // Returns:
268 // The number of intersections. If 0, t 271 // The number of intersections. If 0, the trajectory misses.
269 // 272 //
270 // 273 //
271 // Equation of a line: 274 // Equation of a line:
272 // 275 //
273 // x = x0 + s*tx y = y0 + s*ty 276 // x = x0 + s*tx y = y0 + s*ty z = z0 + s*tz
274 // 277 //
275 // Equation of a hyperbolic surface: 278 // Equation of a hyperbolic surface:
276 // 279 //
277 // x**2 + y**2 = r**2 + (z*tanPhi)**2 280 // x**2 + y**2 = r**2 + (z*tanPhi)**2
278 // 281 //
279 // Solution is quadratic: 282 // Solution is quadratic:
280 // 283 //
281 // a*s**2 + b*s + c = 0 284 // a*s**2 + b*s + c = 0
282 // 285 //
283 // where: 286 // where:
284 // 287 //
285 // a = tx**2 + ty**2 - (tz*tanPhi)**2 288 // a = tx**2 + ty**2 - (tz*tanPhi)**2
286 // 289 //
287 // b = 2*( x0*tx + y0*ty - z0*tz*tanPhi**2 290 // b = 2*( x0*tx + y0*ty - z0*tz*tanPhi**2 )
288 // 291 //
289 // c = x0**2 + y0**2 - r**2 - (z0*tanPhi)* 292 // c = x0**2 + y0**2 - r**2 - (z0*tanPhi)**2
290 // 293 //
291 294
292 fCurStatWithV.ResetfDone(validate, &gp, &gv 295 fCurStatWithV.ResetfDone(validate, &gp, &gv);
293 296
294 if (fCurStatWithV.IsDone()) 297 if (fCurStatWithV.IsDone())
295 { 298 {
296 for (G4int i=0; i<fCurStatWithV.GetNXX() 299 for (G4int i=0; i<fCurStatWithV.GetNXX(); ++i)
297 { 300 {
298 gxx[i] = fCurStatWithV.GetXX(i); 301 gxx[i] = fCurStatWithV.GetXX(i);
299 distance[i] = fCurStatWithV.GetDistan 302 distance[i] = fCurStatWithV.GetDistance(i);
300 areacode[i] = fCurStatWithV.GetAreaco 303 areacode[i] = fCurStatWithV.GetAreacode(i);
301 isvalid[i] = fCurStatWithV.IsValid(i 304 isvalid[i] = fCurStatWithV.IsValid(i);
302 } 305 }
303 return fCurStatWithV.GetNXX(); 306 return fCurStatWithV.GetNXX();
304 } 307 }
305 else // initialize 308 else // initialize
306 { 309 {
307 for (auto i=0; i<2; ++i) 310 for (auto i=0; i<2; ++i)
308 { 311 {
309 distance[i] = kInfinity; 312 distance[i] = kInfinity;
310 areacode[i] = sOutside; 313 areacode[i] = sOutside;
311 isvalid[i] = false; 314 isvalid[i] = false;
312 gxx[i].set(kInfinity, kInfinity, kInf 315 gxx[i].set(kInfinity, kInfinity, kInfinity);
313 } 316 }
314 } 317 }
315 318
316 G4ThreeVector p = ComputeLocalPoint(gp); 319 G4ThreeVector p = ComputeLocalPoint(gp);
317 G4ThreeVector v = ComputeLocalDirection(gv) 320 G4ThreeVector v = ComputeLocalDirection(gv);
318 G4ThreeVector xx[2]; 321 G4ThreeVector xx[2];
319 322
320 // 323 //
321 // special case! p is on origin. 324 // special case! p is on origin.
322 // 325 //
323 326
324 if (p.mag() == 0) 327 if (p.mag() == 0)
325 { 328 {
326 // p is origin. 329 // p is origin.
327 // unique solution of 2-dimension questi 330 // unique solution of 2-dimension question in r-z plane
328 // Equations: 331 // Equations:
329 // r^2 = fR02 + z^2*fTan2Stere0 332 // r^2 = fR02 + z^2*fTan2Stere0
330 // r = beta*z 333 // r = beta*z
331 // where 334 // where
332 // beta = vrho / vz 335 // beta = vrho / vz
333 // Solution (z value of intersection poi 336 // Solution (z value of intersection point):
334 // xxz = +- std::sqrt (fR02 / (beta^2 337 // xxz = +- std::sqrt (fR02 / (beta^2 - fTan2Stereo))
335 // 338 //
336 339
337 G4double vz = v.z(); 340 G4double vz = v.z();
338 G4double absvz = std::fabs(vz); 341 G4double absvz = std::fabs(vz);
339 G4double vrho = v.getRho(); 342 G4double vrho = v.getRho();
340 G4double vslope = vrho/vz; 343 G4double vslope = vrho/vz;
341 G4double vslope2 = vslope * vslope; 344 G4double vslope2 = vslope * vslope;
342 if (vrho == 0 || (vrho/absvz) <= (absvz* 345 if (vrho == 0 || (vrho/absvz) <= (absvz*std::fabs(fTanStereo)/absvz))
343 { 346 {
344 // vz/vrho is bigger than slope of as 347 // vz/vrho is bigger than slope of asymptonic line
345 distance[0] = kInfinity; 348 distance[0] = kInfinity;
346 fCurStatWithV.SetCurrentStatus(0, gxx 349 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
347 isvali 350 isvalid[0], 0, validate, &gp, &gv);
348 return 0; 351 return 0;
349 } 352 }
350 353
351 if (vz != 0.0) << 354 if (vz)
352 { 355 {
353 G4double xxz = std::sqrt(fR02 / (vsl 356 G4double xxz = std::sqrt(fR02 / (vslope2 - fTan2Stereo))
354 * (vz / std::fabs(vz)) 357 * (vz / std::fabs(vz)) ;
355 G4double t = xxz / vz; 358 G4double t = xxz / vz;
356 xx[0].set(t*v.x(), t*v.y(), xxz); 359 xx[0].set(t*v.x(), t*v.y(), xxz);
357 } 360 }
358 else 361 else
359 { 362 {
360 // p.z = 0 && v.z =0 363 // p.z = 0 && v.z =0
361 xx[0].set(v.x()*fR0, v.y()*fR0, 0); 364 xx[0].set(v.x()*fR0, v.y()*fR0, 0); // v is a unit vector.
362 } 365 }
363 distance[0] = xx[0].mag(); 366 distance[0] = xx[0].mag();
364 gxx[0] = ComputeGlobalPoint(xx[0]); 367 gxx[0] = ComputeGlobalPoint(xx[0]);
365 368
366 if (validate == kValidateWithTol) 369 if (validate == kValidateWithTol)
367 { 370 {
368 areacode[0] = GetAreaCode(xx[0]); 371 areacode[0] = GetAreaCode(xx[0]);
369 if (!IsOutside(areacode[0])) 372 if (!IsOutside(areacode[0]))
370 { 373 {
371 if (distance[0] >= 0) isvalid[0] = 374 if (distance[0] >= 0) isvalid[0] = true;
372 } 375 }
373 } 376 }
374 else if (validate == kValidateWithoutTol 377 else if (validate == kValidateWithoutTol)
375 { 378 {
376 areacode[0] = GetAreaCode(xx[0], fals 379 areacode[0] = GetAreaCode(xx[0], false);
377 if (IsInside(areacode[0])) 380 if (IsInside(areacode[0]))
378 { 381 {
379 if (distance[0] >= 0) isvalid[0] = 382 if (distance[0] >= 0) isvalid[0] = true;
380 } 383 }
381 } 384 }
382 else // kDontValidate 385 else // kDontValidate
383 { 386 {
384 areacode[0] = sInside; 387 areacode[0] = sInside;
385 if (distance[0] >= 0) isvalid[0] = 388 if (distance[0] >= 0) isvalid[0] = true;
386 } 389 }
387 390
388 fCurStatWithV.SetCurrentStatus(0, gxx[0] 391 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
389 isvalid[0 392 isvalid[0], 1, validate, &gp, &gv);
390 return 1; 393 return 1;
391 } 394 }
392 395
393 // 396 //
394 // special case end. 397 // special case end.
395 // 398 //
396 399
397 G4double a = v.x()*v.x() + v.y()*v.y() - v. 400 G4double a = v.x()*v.x() + v.y()*v.y() - v.z()*v.z()*fTan2Stereo;
398 G4double b = 2.0 401 G4double b = 2.0
399 * ( p.x() * v.x() + p.y() * v.y( 402 * ( p.x() * v.x() + p.y() * v.y() - p.z() * v.z() * fTan2Stereo );
400 G4double c = p.x()*p.x() + p.y()*p.y() - fR 403 G4double c = p.x()*p.x() + p.y()*p.y() - fR02 - p.z()*p.z()*fTan2Stereo;
401 G4double D = b*b - 4*a*c; //discri 404 G4double D = b*b - 4*a*c; //discriminant
402 G4int vout = 0; 405 G4int vout = 0;
403 406
404 if (std::fabs(a) < DBL_MIN) 407 if (std::fabs(a) < DBL_MIN)
405 { 408 {
406 if (std::fabs(b) > DBL_MIN) / 409 if (std::fabs(b) > DBL_MIN) // single solution
407 { 410 {
408 distance[0] = -c/b; 411 distance[0] = -c/b;
409 xx[0] = p + distance[0]*v; 412 xx[0] = p + distance[0]*v;
410 gxx[0] = ComputeGlobalPoint(xx[0]); 413 gxx[0] = ComputeGlobalPoint(xx[0]);
411 414
412 if (validate == kValidateWithTol) 415 if (validate == kValidateWithTol)
413 { 416 {
414 areacode[0] = GetAreaCode(xx[0]); 417 areacode[0] = GetAreaCode(xx[0]);
415 if (!IsOutside(areacode[0])) 418 if (!IsOutside(areacode[0]))
416 { 419 {
417 if (distance[0] >= 0) isvalid[0 420 if (distance[0] >= 0) isvalid[0] = true;
418 } 421 }
419 } 422 }
420 else if (validate == kValidateWithout 423 else if (validate == kValidateWithoutTol)
421 { 424 {
422 areacode[0] = GetAreaCode(xx[0], f 425 areacode[0] = GetAreaCode(xx[0], false);
423 if (IsInside(areacode[0])) 426 if (IsInside(areacode[0]))
424 { 427 {
425 if (distance[0] >= 0) isvalid[0 428 if (distance[0] >= 0) isvalid[0] = true;
426 } 429 }
427 } 430 }
428 else // kDontValidate 431 else // kDontValidate
429 { 432 {
430 areacode[0] = sInside; 433 areacode[0] = sInside;
431 if (distance[0] >= 0) isvalid[0 434 if (distance[0] >= 0) isvalid[0] = true;
432 } 435 }
433 436
434 fCurStatWithV.SetCurrentStatus(0, gxx 437 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
435 isvali 438 isvalid[0], 1, validate, &gp, &gv);
436 vout = 1; 439 vout = 1;
437 } 440 }
438 else 441 else
439 { 442 {
440 // if a=b=0 and c != 0, p is origin an 443 // if a=b=0 and c != 0, p is origin and v is parallel to asymptotic line
441 // if a=b=c=0, p is on surface and v i 444 // if a=b=c=0, p is on surface and v is paralell to stereo wire.
442 // return distance = infinity 445 // return distance = infinity
443 446
444 fCurStatWithV.SetCurrentStatus(0, gxx 447 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
445 isvali 448 isvalid[0], 0, validate, &gp, &gv);
446 vout = 0; 449 vout = 0;
447 } 450 }
448 } 451 }
449 else if (D > DBL_MIN) // double so 452 else if (D > DBL_MIN) // double solutions
450 { 453 {
451 D = std::sqrt(D); 454 D = std::sqrt(D);
452 G4double factor = 0.5/a; 455 G4double factor = 0.5/a;
453 G4double tmpdist[2] = {kInfinity, k 456 G4double tmpdist[2] = {kInfinity, kInfinity};
454 G4ThreeVector tmpxx[2] ; 457 G4ThreeVector tmpxx[2] ;
455 G4int tmpareacode[2] = {sOutside 458 G4int tmpareacode[2] = {sOutside, sOutside};
456 G4bool tmpisvalid[2] = {false, f 459 G4bool tmpisvalid[2] = {false, false};
457 460
458 for (auto i=0; i<2; ++i) 461 for (auto i=0; i<2; ++i)
459 { 462 {
460 tmpdist[i] = factor*(-b - D); 463 tmpdist[i] = factor*(-b - D);
461 D = -D; 464 D = -D;
462 tmpxx[i] = p + tmpdist[i]*v; 465 tmpxx[i] = p + tmpdist[i]*v;
463 466
464 if (validate == kValidateWithTol) 467 if (validate == kValidateWithTol)
465 { 468 {
466 tmpareacode[i] = GetAreaCode(tmpxx 469 tmpareacode[i] = GetAreaCode(tmpxx[i]);
467 if (!IsOutside(tmpareacode[i])) 470 if (!IsOutside(tmpareacode[i]))
468 { 471 {
469 if (tmpdist[i] >= 0) tmpisvalid 472 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
470 continue; 473 continue;
471 } 474 }
472 } 475 }
473 else if (validate == kValidateWithout 476 else if (validate == kValidateWithoutTol)
474 { 477 {
475 tmpareacode[i] = GetAreaCode(tmpxx 478 tmpareacode[i] = GetAreaCode(tmpxx[i], false);
476 if (IsInside(tmpareacode[i])) 479 if (IsInside(tmpareacode[i]))
477 { 480 {
478 if (tmpdist[i] >= 0) tmpisvalid 481 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
479 continue; 482 continue;
480 } 483 }
481 } 484 }
482 else // kDontValidate 485 else // kDontValidate
483 { 486 {
484 tmpareacode[i] = sInside; 487 tmpareacode[i] = sInside;
485 if (tmpdist[i] >= 0) tmpisvalid 488 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
486 continue; 489 continue;
487 } 490 }
488 } 491 }
489 492
490 if (tmpdist[0] <= tmpdist[1]) 493 if (tmpdist[0] <= tmpdist[1])
491 { 494 {
492 distance[0] = tmpdist[0]; 495 distance[0] = tmpdist[0];
493 distance[1] = tmpdist[1]; 496 distance[1] = tmpdist[1];
494 xx[0] = tmpxx[0]; 497 xx[0] = tmpxx[0];
495 xx[1] = tmpxx[1]; 498 xx[1] = tmpxx[1];
496 gxx[0] = ComputeGlobalPoint(tmp 499 gxx[0] = ComputeGlobalPoint(tmpxx[0]);
497 gxx[1] = ComputeGlobalPoint(tmp 500 gxx[1] = ComputeGlobalPoint(tmpxx[1]);
498 areacode[0] = tmpareacode[0]; 501 areacode[0] = tmpareacode[0];
499 areacode[1] = tmpareacode[1]; 502 areacode[1] = tmpareacode[1];
500 isvalid[0] = tmpisvalid[0]; 503 isvalid[0] = tmpisvalid[0];
501 isvalid[1] = tmpisvalid[1]; 504 isvalid[1] = tmpisvalid[1];
502 } 505 }
503 else 506 else
504 { 507 {
505 distance[0] = tmpdist[1]; 508 distance[0] = tmpdist[1];
506 distance[1] = tmpdist[0]; 509 distance[1] = tmpdist[0];
507 xx[0] = tmpxx[1]; 510 xx[0] = tmpxx[1];
508 xx[1] = tmpxx[0]; 511 xx[1] = tmpxx[0];
509 gxx[0] = ComputeGlobalPoint(tmp 512 gxx[0] = ComputeGlobalPoint(tmpxx[1]);
510 gxx[1] = ComputeGlobalPoint(tmp 513 gxx[1] = ComputeGlobalPoint(tmpxx[0]);
511 areacode[0] = tmpareacode[1]; 514 areacode[0] = tmpareacode[1];
512 areacode[1] = tmpareacode[0]; 515 areacode[1] = tmpareacode[0];
513 isvalid[0] = tmpisvalid[1]; 516 isvalid[0] = tmpisvalid[1];
514 isvalid[1] = tmpisvalid[0]; 517 isvalid[1] = tmpisvalid[0];
515 } 518 }
516 519
517 fCurStatWithV.SetCurrentStatus(0, gxx[0] 520 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
518 isvalid[0 521 isvalid[0], 2, validate, &gp, &gv);
519 fCurStatWithV.SetCurrentStatus(1, gxx[1] 522 fCurStatWithV.SetCurrentStatus(1, gxx[1], distance[1], areacode[1],
520 isvalid[1 523 isvalid[1], 2, validate, &gp, &gv);
521 vout = 2; 524 vout = 2;
522 } 525 }
523 else 526 else
524 { 527 {
525 // if D<0, no solution 528 // if D<0, no solution
526 // if D=0, just grazing the surfaces, re 529 // if D=0, just grazing the surfaces, return kInfinity
527 530
528 fCurStatWithV.SetCurrentStatus(0, gxx[0] 531 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
529 isvalid[0 532 isvalid[0], 0, validate, &gp, &gv);
530 vout = 0; 533 vout = 0;
531 } 534 }
532 return vout; 535 return vout;
533 } 536 }
534 537
535 //============================================ 538 //=====================================================================
536 //* DistanceToSurface ------------------------ 539 //* DistanceToSurface -------------------------------------------------
537 540
538 G4int G4TwistTubsHypeSide::DistanceToSurface(c 541 G4int G4TwistTubsHypeSide::DistanceToSurface(const G4ThreeVector& gp,
539 542 G4ThreeVector gxx[],
540 543 G4double distance[],
541 544 G4int areacode[])
542 { 545 {
543 // Find the approximate distance of a poin 546 // Find the approximate distance of a point of a hyperbolic surface.
544 // The distance must be an underestimate. 547 // The distance must be an underestimate.
545 // It will also be nice (although not nece 548 // It will also be nice (although not necessary) that the estimate is
546 // always finite no matter how close the p 549 // always finite no matter how close the point is.
547 // 550 //
548 // We arranged G4Hype::ApproxDistOutside a 551 // We arranged G4Hype::ApproxDistOutside and G4Hype::ApproxDistInside
549 // for this function. See these discriptio 552 // for this function. See these discriptions.
550 553
551 const G4double halftol 554 const G4double halftol
552 = 0.5 * G4GeometryTolerance::GetInstance( 555 = 0.5 * G4GeometryTolerance::GetInstance()->GetRadialTolerance();
553 556
554 fCurStat.ResetfDone(kDontValidate, &gp); 557 fCurStat.ResetfDone(kDontValidate, &gp);
555 558
556 if (fCurStat.IsDone()) 559 if (fCurStat.IsDone())
557 { 560 {
558 for (G4int i=0; i<fCurStat.GetNXX(); ++i 561 for (G4int i=0; i<fCurStat.GetNXX(); ++i)
559 { 562 {
560 gxx[i] = fCurStat.GetXX(i); 563 gxx[i] = fCurStat.GetXX(i);
561 distance[i] = fCurStat.GetDistance(i) 564 distance[i] = fCurStat.GetDistance(i);
562 areacode[i] = fCurStat.GetAreacode(i) 565 areacode[i] = fCurStat.GetAreacode(i);
563 } 566 }
564 return fCurStat.GetNXX(); 567 return fCurStat.GetNXX();
565 } 568 }
566 else // initialize 569 else // initialize
567 { 570 {
568 for (auto i=0; i<2; ++i) 571 for (auto i=0; i<2; ++i)
569 { 572 {
570 distance[i] = kInfinity; 573 distance[i] = kInfinity;
571 areacode[i] = sOutside; 574 areacode[i] = sOutside;
572 gxx[i].set(kInfinity, kInfinity, kInf 575 gxx[i].set(kInfinity, kInfinity, kInfinity);
573 } 576 }
574 } 577 }
575 578
576 579
577 G4ThreeVector p = ComputeLocalPoint(gp); 580 G4ThreeVector p = ComputeLocalPoint(gp);
578 G4ThreeVector xx; 581 G4ThreeVector xx;
579 582
580 // 583 //
581 // special case! 584 // special case!
582 // If p is on surface, return distance = 0 585 // If p is on surface, return distance = 0 immediatery .
583 // 586 //
584 G4ThreeVector lastgxx[2]; 587 G4ThreeVector lastgxx[2];
585 for (auto i=0; i<2; ++i) 588 for (auto i=0; i<2; ++i)
586 { 589 {
587 lastgxx[i] = fCurStatWithV.GetXX(i); 590 lastgxx[i] = fCurStatWithV.GetXX(i);
588 } 591 }
589 592
590 if ((gp - lastgxx[0]).mag() < halftol || (g 593 if ((gp - lastgxx[0]).mag() < halftol || (gp - lastgxx[1]).mag() < halftol)
591 { 594 {
592 // last winner, or last poststep point i 595 // last winner, or last poststep point is on the surface.
593 xx = p; 596 xx = p;
594 gxx[0] = gp; 597 gxx[0] = gp;
595 distance[0] = 0; 598 distance[0] = 0;
596 599
597 G4bool isvalid = true; 600 G4bool isvalid = true;
598 fCurStat.SetCurrentStatus(0, gxx[0], dis 601 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
599 isvalid, 1, kD 602 isvalid, 1, kDontValidate, &gp);
600 603
601 return 1; 604 return 1;
602 } 605 }
603 // 606 //
604 // special case end 607 // special case end
605 // 608 //
606 609
607 G4double prho = p.getRho(); 610 G4double prho = p.getRho();
608 G4double pz = std::fabs(p.z()); 611 G4double pz = std::fabs(p.z()); // use symmetry
609 G4double r1 = std::sqrt(fR02 + pz * 612 G4double r1 = std::sqrt(fR02 + pz * pz * fTan2Stereo);
610 613
611 G4ThreeVector pabsz(p.x(), p.y(), pz); 614 G4ThreeVector pabsz(p.x(), p.y(), pz);
612 615
613 if (prho > r1 + halftol) // p is outside 616 if (prho > r1 + halftol) // p is outside of Hyperbolic surface
614 { 617 {
615 // First point xx1 618 // First point xx1
616 G4double t = r1 / prho; 619 G4double t = r1 / prho;
617 G4ThreeVector xx1(t * pabsz.x(), t * pab 620 G4ThreeVector xx1(t * pabsz.x(), t * pabsz.y() , pz);
618 621
619 // Second point xx2 622 // Second point xx2
620 G4double z2 = (prho * fTanStereo + pz) / 623 G4double z2 = (prho * fTanStereo + pz) / (1 + fTan2Stereo);
621 G4double r2 = std::sqrt(fR02 + z2 * z2 * 624 G4double r2 = std::sqrt(fR02 + z2 * z2 * fTan2Stereo);
622 t = r2 / prho; 625 t = r2 / prho;
623 G4ThreeVector xx2(t * pabsz.x(), t * pab 626 G4ThreeVector xx2(t * pabsz.x(), t * pabsz.y() , z2);
624 627
625 G4double len = (xx2 - xx1).mag(); 628 G4double len = (xx2 - xx1).mag();
626 if (len < DBL_MIN) 629 if (len < DBL_MIN)
627 { 630 {
628 // xx2 = xx1?? I guess we 631 // xx2 = xx1?? I guess we
629 // must have really bracketed the nor 632 // must have really bracketed the normal
630 distance[0] = (pabsz - xx1).mag(); 633 distance[0] = (pabsz - xx1).mag();
631 xx = xx1; 634 xx = xx1;
632 } 635 }
633 else 636 else
634 { 637 {
635 distance[0] = DistanceToLine(pabsz, x 638 distance[0] = DistanceToLine(pabsz, xx1, (xx2 - xx1) , xx);
636 } 639 }
637 640
638 } 641 }
639 else if (prho < r1 - halftol) // p is insi 642 else if (prho < r1 - halftol) // p is inside of Hyperbolic surface.
640 { 643 {
641 // First point xx1 644 // First point xx1
642 G4double t; 645 G4double t;
643 G4ThreeVector xx1; 646 G4ThreeVector xx1;
644 if (prho < DBL_MIN) 647 if (prho < DBL_MIN)
645 { 648 {
646 xx1.set(r1, 0. , pz); 649 xx1.set(r1, 0. , pz);
647 } 650 }
648 else 651 else
649 { 652 {
650 t = r1 / prho; 653 t = r1 / prho;
651 xx1.set(t * pabsz.x(), t * pabsz.y() 654 xx1.set(t * pabsz.x(), t * pabsz.y() , pz);
652 } 655 }
653 656
654 // dr, dz is tangential vector of Hyparb 657 // dr, dz is tangential vector of Hyparbolic surface at xx1
655 // dr = r, dz = z*tan2stereo 658 // dr = r, dz = z*tan2stereo
656 G4double dr = pz * fTan2Stereo; 659 G4double dr = pz * fTan2Stereo;
657 G4double dz = r1; 660 G4double dz = r1;
658 G4double tanbeta = dr / dz; 661 G4double tanbeta = dr / dz;
659 G4double pztanbeta = pz * tanbeta; 662 G4double pztanbeta = pz * tanbeta;
660 663
661 // Second point xx2 664 // Second point xx2
662 // xx2 is intersection between x-axis an 665 // xx2 is intersection between x-axis and tangential vector
663 G4double r2 = r1 - pztanbeta; 666 G4double r2 = r1 - pztanbeta;
664 G4ThreeVector xx2; 667 G4ThreeVector xx2;
665 if (prho < DBL_MIN) 668 if (prho < DBL_MIN)
666 { 669 {
667 xx2.set(r2, 0. , 0.); 670 xx2.set(r2, 0. , 0.);
668 } 671 }
669 else 672 else
670 { 673 {
671 t = r2 / prho; 674 t = r2 / prho;
672 xx2.set(t * pabsz.x(), t * pabsz.y() 675 xx2.set(t * pabsz.x(), t * pabsz.y() , 0.);
673 } 676 }
674 677
675 G4ThreeVector d = xx2 - xx1; 678 G4ThreeVector d = xx2 - xx1;
676 distance[0] = DistanceToLine(pabsz, xx1, 679 distance[0] = DistanceToLine(pabsz, xx1, d, xx);
677 680
678 } 681 }
679 else // p is on Hyperbolic surface. 682 else // p is on Hyperbolic surface.
680 { 683 {
681 distance[0] = 0; 684 distance[0] = 0;
682 xx.set(p.x(), p.y(), pz); 685 xx.set(p.x(), p.y(), pz);
683 } 686 }
684 687
685 if (p.z() < 0) 688 if (p.z() < 0)
686 { 689 {
687 G4ThreeVector tmpxx(xx.x(), xx.y(), -xx. 690 G4ThreeVector tmpxx(xx.x(), xx.y(), -xx.z());
688 xx = tmpxx; 691 xx = tmpxx;
689 } 692 }
690 693
691 gxx[0] = ComputeGlobalPoint(xx); 694 gxx[0] = ComputeGlobalPoint(xx);
692 areacode[0] = sInside; 695 areacode[0] = sInside;
693 G4bool isvalid = true; 696 G4bool isvalid = true;
694 fCurStat.SetCurrentStatus(0, gxx[0], distan 697 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
695 isvalid, 1, kDont 698 isvalid, 1, kDontValidate, &gp);
696 return 1; 699 return 1;
697 } 700 }
698 701
699 //============================================ 702 //=====================================================================
700 //* GetAreaCode ------------------------------ 703 //* GetAreaCode -------------------------------------------------------
701 704
702 G4int G4TwistTubsHypeSide::GetAreaCode(const G 705 G4int G4TwistTubsHypeSide::GetAreaCode(const G4ThreeVector& xx,
703 G 706 G4bool withTol)
704 { 707 {
705 const G4double ctol = 0.5 * kCarTolerance; 708 const G4double ctol = 0.5 * kCarTolerance;
706 G4int areacode = sInside; 709 G4int areacode = sInside;
707 710
708 if ((fAxis[0] == kPhi && fAxis[1] == kZAxis 711 if ((fAxis[0] == kPhi && fAxis[1] == kZAxis))
709 { 712 {
710 G4int zaxis = 1; 713 G4int zaxis = 1;
711 714
712 if (withTol) 715 if (withTol)
713 { 716 {
714 G4bool isoutside = false; 717 G4bool isoutside = false;
715 G4int phiareacode = GetAreaCodeIn 718 G4int phiareacode = GetAreaCodeInPhi(xx);
716 G4bool isoutsideinphi = IsOutside(phi 719 G4bool isoutsideinphi = IsOutside(phiareacode);
717 720
718 // test boundary of phiaxis 721 // test boundary of phiaxis
719 722
720 if ((phiareacode & sAxisMin) == sAxis 723 if ((phiareacode & sAxisMin) == sAxisMin)
721 { 724 {
722 areacode |= (sAxis0 & (sAxisPhi | 725 areacode |= (sAxis0 & (sAxisPhi | sAxisMin)) | sBoundary;
723 if (isoutsideinphi) isoutside = tr 726 if (isoutsideinphi) isoutside = true;
724 } 727 }
725 else if ((phiareacode & sAxisMax) == 728 else if ((phiareacode & sAxisMax) == sAxisMax)
726 { 729 {
727 areacode |= (sAxis0 & (sAxisPhi | 730 areacode |= (sAxis0 & (sAxisPhi | sAxisMax)) | sBoundary;
728 if (isoutsideinphi) isoutside = tr 731 if (isoutsideinphi) isoutside = true;
729 } 732 }
730 733
731 // test boundary of zaxis 734 // test boundary of zaxis
732 735
733 if (xx.z() < fAxisMin[zaxis] + ctol) 736 if (xx.z() < fAxisMin[zaxis] + ctol)
734 { 737 {
735 areacode |= (sAxis1 & (sAxisZ | sA 738 areacode |= (sAxis1 & (sAxisZ | sAxisMin));
736 if ((areacode & sBoundary) != 0) << 739 if (areacode & sBoundary) areacode |= sCorner; // xx is on corner.
737 else areaco 740 else areacode |= sBoundary;
738 741
739 if (xx.z() <= fAxisMin[zaxis] - ct 742 if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true;
740 743
741 } 744 }
742 else if (xx.z() > fAxisMax[zaxis] - c 745 else if (xx.z() > fAxisMax[zaxis] - ctol)
743 { 746 {
744 areacode |= (sAxis1 & (sAxisZ | sA 747 areacode |= (sAxis1 & (sAxisZ | sAxisMax));
745 if ((areacode & sBoundary) != 0) << 748 if (areacode & sBoundary) areacode |= sCorner; // xx is on corner.
746 else areaco 749 else areacode |= sBoundary;
747 750
748 if (xx.z() >= fAxisMax[zaxis] + ct 751 if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true;
749 } 752 }
750 753
751 // if isoutside = true, clear sInside 754 // if isoutside = true, clear sInside bit.
752 // if not on boundary, add boundary i 755 // if not on boundary, add boundary information.
753 756
754 if (isoutside) 757 if (isoutside)
755 { 758 {
756 G4int tmpareacode = areacode & (~s 759 G4int tmpareacode = areacode & (~sInside);
757 areacode = tmpareacode; 760 areacode = tmpareacode;
758 } 761 }
759 else if ((areacode & sBoundary) != sB 762 else if ((areacode & sBoundary) != sBoundary)
760 { 763 {
761 areacode |= (sAxis0 & sAxisPhi) | 764 areacode |= (sAxis0 & sAxisPhi) | (sAxis1 & sAxisZ);
762 } 765 }
763 766
764 return areacode; 767 return areacode;
765 } 768 }
766 else 769 else
767 { 770 {
768 G4int phiareacode = GetAreaCodeInPhi( 771 G4int phiareacode = GetAreaCodeInPhi(xx, false);
769 772
770 // test boundary of z-axis 773 // test boundary of z-axis
771 774
772 if (xx.z() < fAxisMin[zaxis]) 775 if (xx.z() < fAxisMin[zaxis])
773 { 776 {
774 areacode |= (sAxis1 & (sAxisZ | sA 777 areacode |= (sAxis1 & (sAxisZ | sAxisMin)) | sBoundary;
775 778
776 } 779 }
777 else if (xx.z() > fAxisMax[zaxis]) 780 else if (xx.z() > fAxisMax[zaxis])
778 { 781 {
779 areacode |= (sAxis1 & (sAxisZ | sA 782 areacode |= (sAxis1 & (sAxisZ | sAxisMax)) | sBoundary;
780 } 783 }
781 784
782 // boundary of phi-axis 785 // boundary of phi-axis
783 786
784 if (phiareacode == sAxisMin) 787 if (phiareacode == sAxisMin)
785 { 788 {
786 areacode |= (sAxis0 & (sAxisPhi | 789 areacode |= (sAxis0 & (sAxisPhi | sAxisMin));
787 if ((areacode & sBoundary) != 0) << 790 if (areacode & sBoundary) areacode |= sCorner; // xx is on corner.
788 else areaco 791 else areacode |= sBoundary;
789 792
790 } 793 }
791 else if (phiareacode == sAxisMax) 794 else if (phiareacode == sAxisMax)
792 { 795 {
793 areacode |= (sAxis0 & (sAxisPhi | 796 areacode |= (sAxis0 & (sAxisPhi | sAxisMax));
794 if ((areacode & sBoundary) != 0) << 797 if (areacode & sBoundary) areacode |= sCorner; // xx is on corner.
795 else areaco 798 else areacode |= sBoundary;
796 } 799 }
797 800
798 // if not on boundary, add boundary i 801 // if not on boundary, add boundary information.
799 802
800 if ((areacode & sBoundary) != sBounda 803 if ((areacode & sBoundary) != sBoundary)
801 { 804 {
802 areacode |= (sAxis0 & sAxisPhi) | 805 areacode |= (sAxis0 & sAxisPhi) | (sAxis1 & sAxisZ);
803 } 806 }
804 return areacode; 807 return areacode;
805 } 808 }
806 } 809 }
807 else 810 else
808 { 811 {
809 std::ostringstream message; 812 std::ostringstream message;
810 message << "Feature NOT implemented !" < 813 message << "Feature NOT implemented !" << G4endl
811 << " fAxis[0] = " << fAxi 814 << " fAxis[0] = " << fAxis[0] << G4endl
812 << " fAxis[1] = " << fAxi 815 << " fAxis[1] = " << fAxis[1];
813 G4Exception("G4TwistTubsHypeSide::GetAre 816 G4Exception("G4TwistTubsHypeSide::GetAreaCode()",
814 "GeomSolids0001", FatalExcep 817 "GeomSolids0001", FatalException, message);
815 } 818 }
816 return areacode; 819 return areacode;
817 } 820 }
818 821
819 //============================================ 822 //=====================================================================
820 //* GetAreaCodeInPhi ------------------------- 823 //* GetAreaCodeInPhi --------------------------------------------------
821 824
822 G4int G4TwistTubsHypeSide::GetAreaCodeInPhi(co 825 G4int G4TwistTubsHypeSide::GetAreaCodeInPhi(const G4ThreeVector& xx,
823 826 G4bool withTol)
824 { 827 {
825 828
826 G4ThreeVector lowerlimit; // lower phi-boun 829 G4ThreeVector lowerlimit; // lower phi-boundary limit at z = xx.z()
827 G4ThreeVector upperlimit; // upper phi-boun 830 G4ThreeVector upperlimit; // upper phi-boundary limit at z = xx.z()
828 lowerlimit = GetBoundaryAtPZ(sAxis0 & sAxis 831 lowerlimit = GetBoundaryAtPZ(sAxis0 & sAxisMin, xx);
829 upperlimit = GetBoundaryAtPZ(sAxis0 & sAxis 832 upperlimit = GetBoundaryAtPZ(sAxis0 & sAxisMax, xx);
830 833
831 G4int areacode = sInside; 834 G4int areacode = sInside;
832 G4bool isoutside = false; 835 G4bool isoutside = false;
833 836
834 if (withTol) 837 if (withTol)
835 { 838 {
836 if (AmIOnLeftSide(xx, lowerlimit) >= 0) 839 if (AmIOnLeftSide(xx, lowerlimit) >= 0) // xx is on lowerlimit
837 { 840 {
838 areacode |= (sAxisMin | sBoundary); 841 areacode |= (sAxisMin | sBoundary);
839 if (AmIOnLeftSide(xx, lowerlimit) > 0 842 if (AmIOnLeftSide(xx, lowerlimit) > 0) isoutside = true;
840 843
841 } 844 }
842 else if (AmIOnLeftSide(xx, upperlimit) < 845 else if (AmIOnLeftSide(xx, upperlimit) <= 0) // xx is on upperlimit
843 { 846 {
844 areacode |= (sAxisMax | sBoundary); 847 areacode |= (sAxisMax | sBoundary);
845 if (AmIOnLeftSide(xx, upperlimit) < 0 848 if (AmIOnLeftSide(xx, upperlimit) < 0) isoutside = true;
846 } 849 }
847 850
848 // if isoutside = true, clear inside bit 851 // if isoutside = true, clear inside bit.
849 852
850 if (isoutside) 853 if (isoutside)
851 { 854 {
852 G4int tmpareacode = areacode & (~sIns 855 G4int tmpareacode = areacode & (~sInside);
853 areacode = tmpareacode; 856 areacode = tmpareacode;
854 } 857 }
855 } 858 }
856 else 859 else
857 { 860 {
858 if (AmIOnLeftSide(xx, lowerlimit, false) 861 if (AmIOnLeftSide(xx, lowerlimit, false) >= 0)
859 { 862 {
860 areacode |= (sAxisMin | sBoundary); 863 areacode |= (sAxisMin | sBoundary);
861 } 864 }
862 else if (AmIOnLeftSide(xx, upperlimit, f 865 else if (AmIOnLeftSide(xx, upperlimit, false) <= 0)
863 { 866 {
864 areacode |= (sAxisMax | sBoundary); 867 areacode |= (sAxisMax | sBoundary);
865 } 868 }
866 } 869 }
867 870
868 return areacode; 871 return areacode;
869 } 872 }
870 873
871 //============================================ 874 //=====================================================================
872 //* SetCorners(EndInnerRadius, EndOuterRadius, 875 //* SetCorners(EndInnerRadius, EndOuterRadius,DPhi,EndPhi,EndZ) -------
873 876
874 void G4TwistTubsHypeSide::SetCorners( G4double 877 void G4TwistTubsHypeSide::SetCorners( G4double EndInnerRadius[2],
875 G4double 878 G4double EndOuterRadius[2],
876 G4double 879 G4double DPhi,
877 G4double 880 G4double endPhi[2],
878 G4double 881 G4double endZ[2] )
879 { 882 {
880 // Set Corner points in local coodinate. 883 // Set Corner points in local coodinate.
881 884
882 if (fAxis[0] == kPhi && fAxis[1] == kZAxis) 885 if (fAxis[0] == kPhi && fAxis[1] == kZAxis) {
883 886
884 G4double endRad[2]; 887 G4double endRad[2];
885 G4double halfdphi = 0.5*DPhi; 888 G4double halfdphi = 0.5*DPhi;
886 889
887 for (auto i=0; i<2; ++i) // i=0,1 : -ve 890 for (auto i=0; i<2; ++i) // i=0,1 : -ve z, +ve z
888 { 891 {
889 endRad[i] = (fHandedness == 1 ? EndOut 892 endRad[i] = (fHandedness == 1 ? EndOuterRadius[i] : EndInnerRadius[i]);
890 } 893 }
891 894
892 G4int zmin = 0 ; // at -ve z 895 G4int zmin = 0 ; // at -ve z
893 G4int zmax = 1 ; // at +ve z 896 G4int zmax = 1 ; // at +ve z
894 897
895 G4double x, y, z; 898 G4double x, y, z;
896 899
897 // corner of Axis0min and Axis1min 900 // corner of Axis0min and Axis1min
898 x = endRad[zmin]*std::cos(endPhi[zmin] - 901 x = endRad[zmin]*std::cos(endPhi[zmin] - halfdphi);
899 y = endRad[zmin]*std::sin(endPhi[zmin] - 902 y = endRad[zmin]*std::sin(endPhi[zmin] - halfdphi);
900 z = endZ[zmin]; 903 z = endZ[zmin];
901 SetCorner(sC0Min1Min, x, y, z); 904 SetCorner(sC0Min1Min, x, y, z);
902 905
903 // corner of Axis0max and Axis1min 906 // corner of Axis0max and Axis1min
904 x = endRad[zmin]*std::cos(endPhi[zmin] + 907 x = endRad[zmin]*std::cos(endPhi[zmin] + halfdphi);
905 y = endRad[zmin]*std::sin(endPhi[zmin] + 908 y = endRad[zmin]*std::sin(endPhi[zmin] + halfdphi);
906 z = endZ[zmin]; 909 z = endZ[zmin];
907 SetCorner(sC0Max1Min, x, y, z); 910 SetCorner(sC0Max1Min, x, y, z);
908 911
909 // corner of Axis0max and Axis1max 912 // corner of Axis0max and Axis1max
910 x = endRad[zmax]*std::cos(endPhi[zmax] + 913 x = endRad[zmax]*std::cos(endPhi[zmax] + halfdphi);
911 y = endRad[zmax]*std::sin(endPhi[zmax] + 914 y = endRad[zmax]*std::sin(endPhi[zmax] + halfdphi);
912 z = endZ[zmax]; 915 z = endZ[zmax];
913 SetCorner(sC0Max1Max, x, y, z); 916 SetCorner(sC0Max1Max, x, y, z);
914 917
915 // corner of Axis0min and Axis1max 918 // corner of Axis0min and Axis1max
916 x = endRad[zmax]*std::cos(endPhi[zmax] - 919 x = endRad[zmax]*std::cos(endPhi[zmax] - halfdphi);
917 y = endRad[zmax]*std::sin(endPhi[zmax] - 920 y = endRad[zmax]*std::sin(endPhi[zmax] - halfdphi);
918 z = endZ[zmax]; 921 z = endZ[zmax];
919 SetCorner(sC0Min1Max, x, y, z); 922 SetCorner(sC0Min1Max, x, y, z);
920 923
921 } 924 }
922 else 925 else
923 { 926 {
924 std::ostringstream message; 927 std::ostringstream message;
925 message << "Feature NOT implemented !" < 928 message << "Feature NOT implemented !" << G4endl
926 << " fAxis[0] = " << fAxi 929 << " fAxis[0] = " << fAxis[0] << G4endl
927 << " fAxis[1] = " << fAxi 930 << " fAxis[1] = " << fAxis[1];
928 G4Exception("G4TwistTubsHypeSide::SetCor 931 G4Exception("G4TwistTubsHypeSide::SetCorners()",
929 "GeomSolids0001", FatalExcep 932 "GeomSolids0001", FatalException, message);
930 } 933 }
931 } 934 }
932 935
933 //============================================ 936 //=====================================================================
934 //* SetCorners() ----------------------------- 937 //* SetCorners() ------------------------------------------------------
935 938
936 void G4TwistTubsHypeSide::SetCorners() 939 void G4TwistTubsHypeSide::SetCorners()
937 { 940 {
938 G4Exception("G4TwistTubsHypeSide::SetCorner 941 G4Exception("G4TwistTubsHypeSide::SetCorners()",
939 "GeomSolids0001", FatalExceptio 942 "GeomSolids0001", FatalException,
940 "Method NOT implemented !"); 943 "Method NOT implemented !");
941 } 944 }
942 945
943 //============================================ 946 //=====================================================================
944 //* SetBoundaries() -------------------------- 947 //* SetBoundaries() ---------------------------------------------------
945 948
946 void G4TwistTubsHypeSide::SetBoundaries() 949 void G4TwistTubsHypeSide::SetBoundaries()
947 { 950 {
948 // Set direction-unit vector of phi-boundar 951 // Set direction-unit vector of phi-boundary-lines in local coodinate.
949 // sAxis0 must be kPhi. 952 // sAxis0 must be kPhi.
950 // This fanction set lower phi-boundary and 953 // This fanction set lower phi-boundary and upper phi-boundary.
951 954
952 if (fAxis[0] == kPhi && fAxis[1] == kZAxis) 955 if (fAxis[0] == kPhi && fAxis[1] == kZAxis)
953 { 956 {
954 G4ThreeVector direction; 957 G4ThreeVector direction;
955 // sAxis0 & sAxisMin 958 // sAxis0 & sAxisMin
956 direction = GetCorner(sC0Min1Max) - GetC 959 direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
957 direction = direction.unit(); 960 direction = direction.unit();
958 SetBoundary(sAxis0 & (sAxisPhi | sAxisMi 961 SetBoundary(sAxis0 & (sAxisPhi | sAxisMin), direction,
959 GetCorner(sC0Min1Min), sAxi 962 GetCorner(sC0Min1Min), sAxisZ);
960 963
961 // sAxis0 & sAxisMax 964 // sAxis0 & sAxisMax
962 direction = GetCorner(sC0Max1Max) - GetC 965 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
963 direction = direction.unit(); 966 direction = direction.unit();
964 SetBoundary(sAxis0 & (sAxisPhi | sAxisMa 967 SetBoundary(sAxis0 & (sAxisPhi | sAxisMax), direction,
965 GetCorner(sC0Max1Min), sAxis 968 GetCorner(sC0Max1Min), sAxisZ);
966 969
967 // sAxis1 & sAxisMin 970 // sAxis1 & sAxisMin
968 direction = GetCorner(sC0Max1Min) - GetC 971 direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
969 direction = direction.unit(); 972 direction = direction.unit();
970 SetBoundary(sAxis1 & (sAxisZ | sAxisMin) 973 SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction,
971 GetCorner(sC0Min1Min), sAxi 974 GetCorner(sC0Min1Min), sAxisPhi);
972 975
973 // sAxis1 & sAxisMax 976 // sAxis1 & sAxisMax
974 direction = GetCorner(sC0Max1Max) - GetC 977 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
975 direction = direction.unit(); 978 direction = direction.unit();
976 SetBoundary(sAxis1 & (sAxisZ | sAxisMax) 979 SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction,
977 GetCorner(sC0Min1Max), sAxis 980 GetCorner(sC0Min1Max), sAxisPhi);
978 } 981 }
979 else 982 else
980 { 983 {
981 std::ostringstream message; 984 std::ostringstream message;
982 message << "Feature NOT implemented !" < 985 message << "Feature NOT implemented !" << G4endl
983 << " fAxis[0] = " << fAxi 986 << " fAxis[0] = " << fAxis[0] << G4endl
984 << " fAxis[1] = " << fAxi 987 << " fAxis[1] = " << fAxis[1];
985 G4Exception("G4TwistTubsHypeSide::SetBou 988 G4Exception("G4TwistTubsHypeSide::SetBoundaries()",
986 "GeomSolids0001", FatalExcep 989 "GeomSolids0001", FatalException, message);
987 } 990 }
988 } 991 }
989 992
990 //============================================ 993 //=====================================================================
991 //* GetFacets() ------------------------------ 994 //* GetFacets() -------------------------------------------------------
992 995
993 void G4TwistTubsHypeSide::GetFacets( G4int k, 996 void G4TwistTubsHypeSide::GetFacets( G4int k, G4int n, G4double xyz[][3],
994 G4int fac 997 G4int faces[][4], G4int iside )
995 { 998 {
996 G4double z ; // the two parameters for t 999 G4double z ; // the two parameters for the surface equation
997 G4double x,xmin,xmax ; 1000 G4double x,xmin,xmax ;
998 1001
999 G4ThreeVector p ; // a point on the surface 1002 G4ThreeVector p ; // a point on the surface, given by (z,u)
1000 1003
1001 G4int nnode ; 1004 G4int nnode ;
1002 G4int nface ; 1005 G4int nface ;
1003 1006
1004 // calculate the (n-1)*(k-1) vertices 1007 // calculate the (n-1)*(k-1) vertices
1005 1008
1006 for ( G4int i = 0 ; i<n ; ++i ) 1009 for ( G4int i = 0 ; i<n ; ++i )
1007 { 1010 {
1008 z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin 1011 z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin[1])/(n-1) ;
1009 1012
1010 for ( G4int j = 0 ; j<k ; ++j ) 1013 for ( G4int j = 0 ; j<k ; ++j )
1011 { 1014 {
1012 nnode = GetNode(i,j,k,n,iside) ; 1015 nnode = GetNode(i,j,k,n,iside) ;
1013 1016
1014 xmin = GetBoundaryMin(z) ; 1017 xmin = GetBoundaryMin(z) ;
1015 xmax = GetBoundaryMax(z) ; 1018 xmax = GetBoundaryMax(z) ;
1016 1019
1017 if (fHandedness < 0) // inner hyperbol 1020 if (fHandedness < 0) // inner hyperbolic surface
1018 { 1021 {
1019 x = xmin + j*(xmax-xmin)/(k-1) ; 1022 x = xmin + j*(xmax-xmin)/(k-1) ;
1020 } 1023 }
1021 else // outer hyperbol 1024 else // outer hyperbolic surface
1022 { 1025 {
1023 x = xmax - j*(xmax-xmin)/(k-1) ; 1026 x = xmax - j*(xmax-xmin)/(k-1) ;
1024 } 1027 }
1025 1028
1026 p = SurfacePoint(x,z,true) ; // surfac 1029 p = SurfacePoint(x,z,true) ; // surface point in global coord.system
1027 1030
1028 xyz[nnode][0] = p.x() ; 1031 xyz[nnode][0] = p.x() ;
1029 xyz[nnode][1] = p.y() ; 1032 xyz[nnode][1] = p.y() ;
1030 xyz[nnode][2] = p.z() ; 1033 xyz[nnode][2] = p.z() ;
1031 1034
1032 if ( i<n-1 && j<k-1 ) // clock wise 1035 if ( i<n-1 && j<k-1 ) // clock wise filling
1033 { 1036 {
1034 nface = GetFace(i,j,k,n,iside) ; 1037 nface = GetFace(i,j,k,n,iside) ;
1035 faces[nface][0] = GetEdgeVisibility(i,j,k,n 1038 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1)
1036 * ( GetNode(i ,j ,k 1039 * ( GetNode(i ,j ,k,n,iside)+1) ;
1037 faces[nface][1] = GetEdgeVisibility(i,j,k,n 1040 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1)
1038 * ( GetNode(i+1,j ,k 1041 * ( GetNode(i+1,j ,k,n,iside)+1) ;
1039 faces[nface][2] = GetEdgeVisibility(i,j,k,n 1042 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1)
1040 * ( GetNode(i+1,j+1,k 1043 * ( GetNode(i+1,j+1,k,n,iside)+1) ;
1041 faces[nface][3] = GetEdgeVisibility(i,j,k,n 1044 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1)
1042 * ( GetNode(i ,j+1,k 1045 * ( GetNode(i ,j+1,k,n,iside)+1) ;
1043 } 1046 }
1044 } 1047 }
1045 } 1048 }
1046 } 1049 }
1047 1050