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