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