| 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 // G4TwistTubsFlatSide implementation 26 // G4TwistTubsFlatSide implementation
27 // 27 //
28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepbu 28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp), created.
29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), 29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
30 // from original version in Jupi 30 // from original version in Jupiter-2.5.02 application.
31 // ------------------------------------------- 31 // --------------------------------------------------------------------
32 32
33 #include "G4TwistTubsFlatSide.hh" 33 #include "G4TwistTubsFlatSide.hh"
34 #include "G4GeometryTolerance.hh" 34 #include "G4GeometryTolerance.hh"
35 35
36 //============================================ 36 //=====================================================================
37 //* constructors ----------------------------- 37 //* constructors ------------------------------------------------------
38 38
39 G4TwistTubsFlatSide::G4TwistTubsFlatSide(const 39 G4TwistTubsFlatSide::G4TwistTubsFlatSide(const G4String& name,
40 const G4RotationM 40 const G4RotationMatrix& rot,
41 const G4ThreeVect 41 const G4ThreeVector& tlate,
42 const G4ThreeVect 42 const G4ThreeVector& n,
43 const EAxis 43 const EAxis axis0 ,
44 const EAxis 44 const EAxis axis1 ,
45 G4double 45 G4double axis0min,
46 G4double 46 G4double axis1min,
47 G4double 47 G4double axis0max,
48 G4double 48 G4double axis1max )
49 : G4VTwistSurface(name, rot, tlate, 0, axis0 49 : G4VTwistSurface(name, rot, tlate, 0, axis0, axis1,
50 axis0min, axis1min, axis0m 50 axis0min, axis1min, axis0max, axis1max)
51 { 51 {
52 if (axis0 == kPhi && axis1 == kRho) 52 if (axis0 == kPhi && axis1 == kRho)
53 { 53 {
54 G4Exception("G4TwistTubsFlatSide::G4Twis 54 G4Exception("G4TwistTubsFlatSide::G4TwistTubsFlatSide()",
55 "GeomSolids0002", FatalError 55 "GeomSolids0002", FatalErrorInArgument,
56 "Should swap axis0 and axis1 56 "Should swap axis0 and axis1!");
57 } 57 }
58 58
59 G4ThreeVector normal = rot.inverse()*n; 59 G4ThreeVector normal = rot.inverse()*n;
60 fCurrentNormal.normal = normal.unit(); // 60 fCurrentNormal.normal = normal.unit(); // in local coordinate system
61 fIsValidNorm = true; 61 fIsValidNorm = true;
62 62
63 SetCorners(); 63 SetCorners();
64 SetBoundaries(); 64 SetBoundaries();
65 65
66 fSurfaceArea = 1. ; // NOTE: not yet imple 66 fSurfaceArea = 1. ; // NOTE: not yet implemented!
67 } 67 }
68 68
69 G4TwistTubsFlatSide::G4TwistTubsFlatSide( cons 69 G4TwistTubsFlatSide::G4TwistTubsFlatSide( const G4String& name,
70 70 G4double EndInnerRadius[2],
71 71 G4double EndOuterRadius[2],
72 72 G4double DPhi,
73 73 G4double EndPhi[2],
74 74 G4double EndZ[2],
75 75 G4int handedness )
76 : G4VTwistSurface(name) 76 : G4VTwistSurface(name)
77 { 77 {
78 fHandedness = handedness; // +z = +ve, -z 78 fHandedness = handedness; // +z = +ve, -z = -ve
79 fAxis[0] = kRho; // in local coo 79 fAxis[0] = kRho; // in local coordinate system
80 fAxis[1] = kPhi; 80 fAxis[1] = kPhi;
81 G4int i = (handedness < 0 ? 0 : 1); 81 G4int i = (handedness < 0 ? 0 : 1);
82 fAxisMin[0] = EndInnerRadius[i]; // Inner- 82 fAxisMin[0] = EndInnerRadius[i]; // Inner-hype radius at z=0
83 fAxisMax[0] = EndOuterRadius[i]; // Outer- 83 fAxisMax[0] = EndOuterRadius[i]; // Outer-hype radius at z=0
84 fAxisMin[1] = -0.5*DPhi; 84 fAxisMin[1] = -0.5*DPhi;
85 fAxisMax[1] = -fAxisMin[1]; 85 fAxisMax[1] = -fAxisMin[1];
86 fCurrentNormal.normal.set(0, 0, (fHandednes 86 fCurrentNormal.normal.set(0, 0, (fHandedness < 0 ? -1 : 1));
87 // Unit vector, in local coordinate s 87 // Unit vector, in local coordinate system
88 fRot.rotateZ(EndPhi[i]); 88 fRot.rotateZ(EndPhi[i]);
89 fTrans.set(0, 0, EndZ[i]); 89 fTrans.set(0, 0, EndZ[i]);
90 fIsValidNorm = true; 90 fIsValidNorm = true;
91 91
92 SetCorners(); 92 SetCorners();
93 SetBoundaries(); 93 SetBoundaries();
94 94
95 fSurfaceArea = 0.5*DPhi * (EndOuterRadius[ 95 fSurfaceArea = 0.5*DPhi * (EndOuterRadius[i]*EndOuterRadius[i]
96 - EndInnerRadius[ 96 - EndInnerRadius[i]*EndInnerRadius[i] ) ;
97 } 97 }
98 98
99 //============================================ 99 //=====================================================================
100 //* Fake default constructor ----------------- 100 //* Fake default constructor ------------------------------------------
101 101
102 G4TwistTubsFlatSide::G4TwistTubsFlatSide( __vo 102 G4TwistTubsFlatSide::G4TwistTubsFlatSide( __void__& a )
103 : G4VTwistSurface(a) 103 : G4VTwistSurface(a)
104 { 104 {
105 } 105 }
106 106
107 //============================================ 107 //=====================================================================
108 //* destructor ------------------------------- 108 //* destructor --------------------------------------------------------
109 109
110 G4TwistTubsFlatSide::~G4TwistTubsFlatSide() = << 110 G4TwistTubsFlatSide::~G4TwistTubsFlatSide()
>> 111 {
>> 112 }
111 113
112 //============================================ 114 //=====================================================================
113 //* GetNormal -------------------------------- 115 //* GetNormal ---------------------------------------------------------
114 116
115 G4ThreeVector G4TwistTubsFlatSide::GetNormal(c 117 G4ThreeVector G4TwistTubsFlatSide::GetNormal(const G4ThreeVector& /* xx */ ,
116 118 G4bool isGlobal)
117 { 119 {
118 if (isGlobal) 120 if (isGlobal)
119 { 121 {
120 return ComputeGlobalDirection(fCurrentNo 122 return ComputeGlobalDirection(fCurrentNormal.normal);
121 } 123 }
122 else 124 else
123 { 125 {
124 return fCurrentNormal.normal; 126 return fCurrentNormal.normal;
125 } 127 }
126 } 128 }
127 129
128 //============================================ 130 //=====================================================================
129 //* DistanceToSurface(p, v) ------------------ 131 //* DistanceToSurface(p, v) -------------------------------------------
130 132
131 G4int G4TwistTubsFlatSide::DistanceToSurface(c 133 G4int G4TwistTubsFlatSide::DistanceToSurface(const G4ThreeVector& gp,
132 c 134 const G4ThreeVector& gv,
133 135 G4ThreeVector gxx[],
134 136 G4double distance[],
135 137 G4int areacode[],
136 138 G4bool isvalid[],
137 139 EValidate validate)
138 { 140 {
139 fCurStatWithV.ResetfDone(validate, &gp, &gv 141 fCurStatWithV.ResetfDone(validate, &gp, &gv);
140 142
141 if (fCurStatWithV.IsDone()) 143 if (fCurStatWithV.IsDone())
142 { 144 {
143 for (G4int i=0; i<fCurStatWithV.GetNXX() 145 for (G4int i=0; i<fCurStatWithV.GetNXX(); ++i)
144 { 146 {
145 gxx[i] = fCurStatWithV.GetXX(i); 147 gxx[i] = fCurStatWithV.GetXX(i);
146 distance[i] = fCurStatWithV.GetDistan 148 distance[i] = fCurStatWithV.GetDistance(i);
147 areacode[i] = fCurStatWithV.GetAreaco 149 areacode[i] = fCurStatWithV.GetAreacode(i);
148 isvalid[i] = fCurStatWithV.IsValid(i 150 isvalid[i] = fCurStatWithV.IsValid(i);
149 } 151 }
150 return fCurStatWithV.GetNXX(); 152 return fCurStatWithV.GetNXX();
151 } 153 }
152 else // initialize 154 else // initialize
153 { 155 {
154 for (G4int i=0; i<2; ++i) 156 for (G4int i=0; i<2; ++i)
155 { 157 {
156 distance[i] = kInfinity; 158 distance[i] = kInfinity;
157 areacode[i] = sOutside; 159 areacode[i] = sOutside;
158 isvalid[i] = false; 160 isvalid[i] = false;
159 gxx[i].set(kInfinity, kInfinity, kInf 161 gxx[i].set(kInfinity, kInfinity, kInfinity);
160 } 162 }
161 } 163 }
162 164
163 G4ThreeVector p = ComputeLocalPoint(gp); 165 G4ThreeVector p = ComputeLocalPoint(gp);
164 G4ThreeVector v = ComputeLocalDirection(gv) 166 G4ThreeVector v = ComputeLocalDirection(gv);
165 167
166 // 168 //
167 // special case! 169 // special case!
168 // if p is on surface, distance = 0. 170 // if p is on surface, distance = 0.
169 // 171 //
170 172
171 if (std::fabs(p.z()) == 0.) // if p is o 173 if (std::fabs(p.z()) == 0.) // if p is on the plane
172 { 174 {
173 distance[0] = 0; 175 distance[0] = 0;
174 G4ThreeVector xx = p; 176 G4ThreeVector xx = p;
175 gxx[0] = ComputeGlobalPoint(xx); 177 gxx[0] = ComputeGlobalPoint(xx);
176 178
177 if (validate == kValidateWithTol) 179 if (validate == kValidateWithTol)
178 { 180 {
179 areacode[0] = GetAreaCode(xx); 181 areacode[0] = GetAreaCode(xx);
180 if (!IsOutside(areacode[0])) 182 if (!IsOutside(areacode[0]))
181 { 183 {
182 isvalid[0] = true; 184 isvalid[0] = true;
183 } 185 }
184 } 186 }
185 else if (validate == kValidateWithoutTol 187 else if (validate == kValidateWithoutTol)
186 { 188 {
187 areacode[0] = GetAreaCode(xx, false); 189 areacode[0] = GetAreaCode(xx, false);
188 if (IsInside(areacode[0])) 190 if (IsInside(areacode[0]))
189 { 191 {
190 isvalid[0] = true; 192 isvalid[0] = true;
191 } 193 }
192 } 194 }
193 else // kDontValidate 195 else // kDontValidate
194 { 196 {
195 areacode[0] = sInside; 197 areacode[0] = sInside;
196 isvalid[0] = true; 198 isvalid[0] = true;
197 } 199 }
198 return 1; 200 return 1;
199 } 201 }
200 // 202 //
201 // special case end 203 // special case end
202 // 204 //
203 205
204 if (v.z() == 0) 206 if (v.z() == 0)
205 { 207 {
206 fCurStatWithV.SetCurrentStatus(0, gxx[0] 208 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
207 isvalid[0 209 isvalid[0], 0, validate, &gp, &gv);
208 return 0; 210 return 0;
209 } 211 }
210 212
211 distance[0] = - (p.z() / v.z()); 213 distance[0] = - (p.z() / v.z());
212 214
213 G4ThreeVector xx = p + distance[0]*v; 215 G4ThreeVector xx = p + distance[0]*v;
214 gxx[0] = ComputeGlobalPoint(xx); 216 gxx[0] = ComputeGlobalPoint(xx);
215 217
216 if (validate == kValidateWithTol) 218 if (validate == kValidateWithTol)
217 { 219 {
218 areacode[0] = GetAreaCode(xx); 220 areacode[0] = GetAreaCode(xx);
219 if (!IsOutside(areacode[0])) 221 if (!IsOutside(areacode[0]))
220 { 222 {
221 if (distance[0] >= 0) isvalid[0] = tr 223 if (distance[0] >= 0) isvalid[0] = true;
222 } 224 }
223 } 225 }
224 else if (validate == kValidateWithoutTol) 226 else if (validate == kValidateWithoutTol)
225 { 227 {
226 areacode[0] = GetAreaCode(xx, false); 228 areacode[0] = GetAreaCode(xx, false);
227 if (IsInside(areacode[0])) 229 if (IsInside(areacode[0]))
228 { 230 {
229 if (distance[0] >= 0) isvalid[0] = tr 231 if (distance[0] >= 0) isvalid[0] = true;
230 } 232 }
231 } 233 }
232 else // kDontValidate 234 else // kDontValidate
233 { 235 {
234 areacode[0] = sInside; 236 areacode[0] = sInside;
235 if (distance[0] >= 0) isvalid[0] = tr 237 if (distance[0] >= 0) isvalid[0] = true;
236 } 238 }
237 239
238 fCurStatWithV.SetCurrentStatus(0, gxx[0], d 240 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
239 isvalid[0], 241 isvalid[0], 1, validate, &gp, &gv);
240 242
241 #ifdef G4TWISTDEBUG 243 #ifdef G4TWISTDEBUG
242 G4cerr << "ERROR - G4TwistTubsFlatSide::Dis 244 G4cerr << "ERROR - G4TwistTubsFlatSide::DistanceToSurface(p,v)" << G4endl;
243 G4cerr << " Name : " << GetNa 245 G4cerr << " Name : " << GetName() << G4endl;
244 G4cerr << " xx : " << xx << 246 G4cerr << " xx : " << xx << G4endl;
245 G4cerr << " gxx[0] : " << gxx[0 247 G4cerr << " gxx[0] : " << gxx[0] << G4endl;
246 G4cerr << " dist[0] : " << dista 248 G4cerr << " dist[0] : " << distance[0] << G4endl;
247 G4cerr << " areacode[0] : " << areac 249 G4cerr << " areacode[0] : " << areacode[0] << G4endl;
248 G4cerr << " isvalid[0] : " << isval 250 G4cerr << " isvalid[0] : " << isvalid[0] << G4endl;
249 #endif 251 #endif
250 return 1; 252 return 1;
251 } 253 }
252 254
253 //============================================ 255 //=====================================================================
254 //* DistanceToSurface(p) --------------------- 256 //* DistanceToSurface(p) ----------------------------------------------
255 257
256 G4int G4TwistTubsFlatSide::DistanceToSurface(c 258 G4int G4TwistTubsFlatSide::DistanceToSurface(const G4ThreeVector& gp,
257 G 259 G4ThreeVector gxx[],
258 G 260 G4double distance[],
259 G 261 G4int areacode[])
260 { 262 {
261 // Calculate distance to plane in local coo 263 // Calculate distance to plane in local coordinate,
262 // then return distance and global intersec 264 // then return distance and global intersection points.
263 // 265 //
264 266
265 fCurStat.ResetfDone(kDontValidate, &gp); 267 fCurStat.ResetfDone(kDontValidate, &gp);
266 268
267 if (fCurStat.IsDone()) 269 if (fCurStat.IsDone())
268 { 270 {
269 for (G4int i=0; i<fCurStat.GetNXX(); ++i 271 for (G4int i=0; i<fCurStat.GetNXX(); ++i)
270 { 272 {
271 gxx[i] = fCurStat.GetXX(i); 273 gxx[i] = fCurStat.GetXX(i);
272 distance[i] = fCurStat.GetDistance(i) 274 distance[i] = fCurStat.GetDistance(i);
273 areacode[i] = fCurStat.GetAreacode(i) 275 areacode[i] = fCurStat.GetAreacode(i);
274 } 276 }
275 return fCurStat.GetNXX(); 277 return fCurStat.GetNXX();
276 } 278 }
277 else // initialize 279 else // initialize
278 { 280 {
279 for (auto i=0; i<2; ++i) 281 for (auto i=0; i<2; ++i)
280 { 282 {
281 distance[i] = kInfinity; 283 distance[i] = kInfinity;
282 areacode[i] = sOutside; 284 areacode[i] = sOutside;
283 gxx[i].set(kInfinity, kInfinity, kInf 285 gxx[i].set(kInfinity, kInfinity, kInfinity);
284 } 286 }
285 } 287 }
286 288
287 G4ThreeVector p = ComputeLocalPoint(gp); 289 G4ThreeVector p = ComputeLocalPoint(gp);
288 G4ThreeVector xx; 290 G4ThreeVector xx;
289 291
290 // The plane is placed on origin with makin 292 // The plane is placed on origin with making its normal
291 // parallel to z-axis. 293 // parallel to z-axis.
292 if (std::fabs(p.z()) <= 0.5 * kCarTolerance 294 if (std::fabs(p.z()) <= 0.5 * kCarTolerance) // if p is on plane, return 1
293 { 295 {
294 distance[0] = 0; 296 distance[0] = 0;
295 xx = p; 297 xx = p;
296 } 298 }
297 else 299 else
298 { 300 {
299 distance[0] = std::fabs(p.z()); 301 distance[0] = std::fabs(p.z());
300 xx.set(p.x(), p.y(), 0); 302 xx.set(p.x(), p.y(), 0);
301 } 303 }
302 304
303 gxx[0] = ComputeGlobalPoint(xx); 305 gxx[0] = ComputeGlobalPoint(xx);
304 areacode[0] = sInside; 306 areacode[0] = sInside;
305 G4bool isvalid = true; 307 G4bool isvalid = true;
306 fCurStat.SetCurrentStatus(0, gxx[0], distan 308 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
307 isvalid, 1, kDont 309 isvalid, 1, kDontValidate, &gp);
308 return 1; 310 return 1;
309 } 311 }
310 312
311 //============================================ 313 //=====================================================================
312 //* GetAreaCode ------------------------------ 314 //* GetAreaCode -------------------------------------------------------
313 315
314 G4int G4TwistTubsFlatSide::GetAreaCode(const G 316 G4int G4TwistTubsFlatSide::GetAreaCode(const G4ThreeVector &xx,
315 G 317 G4bool withTol)
316 { 318 {
317 const G4double rtol 319 const G4double rtol
318 = 0.5*G4GeometryTolerance::GetInstance()- 320 = 0.5*G4GeometryTolerance::GetInstance()->GetRadialTolerance();
319 321
320 G4int areacode = sInside; 322 G4int areacode = sInside;
321 323
322 if (fAxis[0] == kRho && fAxis[1] == kPhi) 324 if (fAxis[0] == kRho && fAxis[1] == kPhi)
323 { 325 {
324 G4int rhoaxis = 0; 326 G4int rhoaxis = 0;
325 327
326 G4ThreeVector dphimin; // direction of 328 G4ThreeVector dphimin; // direction of phi-minimum boundary
327 G4ThreeVector dphimax; // direction of 329 G4ThreeVector dphimax; // direction of phi-maximum boundary
328 dphimin = GetCorner(sC0Max1Min); 330 dphimin = GetCorner(sC0Max1Min);
329 dphimax = GetCorner(sC0Max1Max); 331 dphimax = GetCorner(sC0Max1Max);
330 332
331 if (withTol) 333 if (withTol)
332 { 334 {
333 G4bool isoutside = false; 335 G4bool isoutside = false;
334 336
335 // test boundary of rho-axis 337 // test boundary of rho-axis
336 338
337 if (xx.getRho() <= fAxisMin[rhoaxis] 339 if (xx.getRho() <= fAxisMin[rhoaxis] + rtol)
338 { 340 {
339 areacode |= (sAxis0 & (sAxisRho|sA 341 areacode |= (sAxis0 & (sAxisRho|sAxisMin)) | sBoundary; // rho-min
340 if (xx.getRho() < fAxisMin[rhoaxis 342 if (xx.getRho() < fAxisMin[rhoaxis] - rtol) isoutside = true;
341 343
342 } 344 }
343 else if (xx.getRho() >= fAxisMax[rhoa 345 else if (xx.getRho() >= fAxisMax[rhoaxis] - rtol)
344 { 346 {
345 areacode |= (sAxis0 & (sAxisRho|sA 347 areacode |= (sAxis0 & (sAxisRho|sAxisMax)) | sBoundary; // rho-max
346 if (xx.getRho() > fAxisMax[rhoaxis 348 if (xx.getRho() > fAxisMax[rhoaxis] + rtol) isoutside = true;
347 } 349 }
348 350
349 // test boundary of phi-axis 351 // test boundary of phi-axis
350 352
351 if (AmIOnLeftSide(xx, dphimin) >= 0) 353 if (AmIOnLeftSide(xx, dphimin) >= 0) // xx is on dphimin
352 { 354 {
353 areacode |= (sAxis1 & (sAxisPhi | 355 areacode |= (sAxis1 & (sAxisPhi | sAxisMin));
354 if ((areacode & sBoundary) != 0) << 356 if (areacode & sBoundary) areacode |= sCorner; // xx is on corner.
355 else areaco 357 else areacode |= sBoundary;
356 358
357 if (AmIOnLeftSide(xx, dphimin) > 0 359 if (AmIOnLeftSide(xx, dphimin) > 0) isoutside = true;
358 360
359 } 361 }
360 else if (AmIOnLeftSide(xx, dphimax) < 362 else if (AmIOnLeftSide(xx, dphimax) <= 0) // xx is on dphimax
361 { 363 {
362 areacode |= (sAxis1 & (sAxisPhi | 364 areacode |= (sAxis1 & (sAxisPhi | sAxisMax));
363 if ((areacode & sBoundary) != 0) << 365 if (areacode & sBoundary) areacode |= sCorner; // xx is on corner.
364 else areaco 366 else areacode |= sBoundary;
365 367
366 if (AmIOnLeftSide(xx, dphimax) < 0 368 if (AmIOnLeftSide(xx, dphimax) < 0) isoutside = true;
367 } 369 }
368 370
369 // if isoutside = true, clear inside 371 // if isoutside = true, clear inside bit.
370 // if not on boundary, add axis infor 372 // if not on boundary, add axis information.
371 373
372 if (isoutside) 374 if (isoutside)
373 { 375 {
374 G4int tmpareacode = areacode & (~s 376 G4int tmpareacode = areacode & (~sInside);
375 areacode = tmpareacode; 377 areacode = tmpareacode;
376 } 378 }
377 else if ((areacode & sBoundary) != sB 379 else if ((areacode & sBoundary) != sBoundary)
378 { 380 {
379 areacode |= (sAxis0 & sAxisRho) | 381 areacode |= (sAxis0 & sAxisRho) | (sAxis1 & sAxisPhi);
380 } 382 }
381 383
382 } 384 }
383 else 385 else
384 { 386 {
385 // out of boundary of rho-axis 387 // out of boundary of rho-axis
386 388
387 if (xx.getRho() < fAxisMin[rhoaxis]) 389 if (xx.getRho() < fAxisMin[rhoaxis])
388 { 390 {
389 areacode |= (sAxis0 & (sAxisRho | 391 areacode |= (sAxis0 & (sAxisRho | sAxisMin)) | sBoundary;
390 } 392 }
391 else if (xx.getRho() > fAxisMax[rhoax 393 else if (xx.getRho() > fAxisMax[rhoaxis])
392 { 394 {
393 areacode |= (sAxis0 & (sAxisRho | 395 areacode |= (sAxis0 & (sAxisRho | sAxisMax)) | sBoundary;
394 } 396 }
395 397
396 // out of boundary of phi-axis 398 // out of boundary of phi-axis
397 399
398 if (AmIOnLeftSide(xx, dphimin, false) 400 if (AmIOnLeftSide(xx, dphimin, false) >= 0) // xx is leftside or
399 { 401 {
400 areacode |= (sAxis1 & (sAxisPhi | s 402 areacode |= (sAxis1 & (sAxisPhi | sAxisMin)) ; // boundary of dphimin
401 if ((areacode & sBoundary) != 0) << 403 if (areacode & sBoundary) areacode |= sCorner; // xx is on corner
402 else areacod 404 else areacode |= sBoundary;
403 405
404 } 406 }
405 else if (AmIOnLeftSide(xx, dphimax, f 407 else if (AmIOnLeftSide(xx, dphimax, false) <= 0) // xx is rightside or
406 { 408 {
407 areacode |= (sAxis1 & (sAxisPhi | s 409 areacode |= (sAxis1 & (sAxisPhi | sAxisMax)); // boundary of dphimax
408 if ((areacode & sBoundary) != 0) << 410 if (areacode & sBoundary) areacode |= sCorner; // xx is on corner
409 else areacod 411 else areacode |= sBoundary;
410 412
411 } 413 }
412 414
413 if ((areacode & sBoundary) != sBounda 415 if ((areacode & sBoundary) != sBoundary) {
414 areacode |= (sAxis0 & sAxisRho) | 416 areacode |= (sAxis0 & sAxisRho) | (sAxis1 & sAxisPhi);
415 } 417 }
416 418
417 } 419 }
418 return areacode; 420 return areacode;
419 } 421 }
420 else 422 else
421 { 423 {
422 std::ostringstream message; 424 std::ostringstream message;
423 message << "Feature NOT implemented !" < 425 message << "Feature NOT implemented !" << G4endl
424 << " fAxis[0] = " << fAxi 426 << " fAxis[0] = " << fAxis[0] << G4endl
425 << " fAxis[1] = " << fAxi 427 << " fAxis[1] = " << fAxis[1];
426 G4Exception("G4TwistTubsFlatSide::GetAre 428 G4Exception("G4TwistTubsFlatSide::GetAreaCode()", "GeomSolids0001",
427 FatalException, message); 429 FatalException, message);
428 } 430 }
429 return areacode; 431 return areacode;
430 } 432 }
431 433
432 //============================================ 434 //=====================================================================
433 //* SetCorners ------------------------------- 435 //* SetCorners --------------------------------------------------------
434 436
435 void G4TwistTubsFlatSide::SetCorners() 437 void G4TwistTubsFlatSide::SetCorners()
436 { 438 {
437 // Set Corner points in local coodinate. 439 // Set Corner points in local coodinate.
438 440
439 if (fAxis[0] == kRho && fAxis[1] == kPhi) 441 if (fAxis[0] == kRho && fAxis[1] == kPhi)
440 { 442 {
441 G4int rhoaxis = 0; // kRho 443 G4int rhoaxis = 0; // kRho
442 G4int phiaxis = 1; // kPhi 444 G4int phiaxis = 1; // kPhi
443 445
444 G4double x, y, z; 446 G4double x, y, z;
445 // corner of Axis0min and Axis1min 447 // corner of Axis0min and Axis1min
446 x = fAxisMin[rhoaxis]*std::cos(fAxisM 448 x = fAxisMin[rhoaxis]*std::cos(fAxisMin[phiaxis]);
447 y = fAxisMin[rhoaxis]*std::sin(fAxisM 449 y = fAxisMin[rhoaxis]*std::sin(fAxisMin[phiaxis]);
448 z = 0; 450 z = 0;
449 SetCorner(sC0Min1Min, x, y, z); 451 SetCorner(sC0Min1Min, x, y, z);
450 // corner of Axis0max and Axis1min 452 // corner of Axis0max and Axis1min
451 x = fAxisMax[rhoaxis]*std::cos(fAxisM 453 x = fAxisMax[rhoaxis]*std::cos(fAxisMin[phiaxis]);
452 y = fAxisMax[rhoaxis]*std::sin(fAxisM 454 y = fAxisMax[rhoaxis]*std::sin(fAxisMin[phiaxis]);
453 z = 0; 455 z = 0;
454 SetCorner(sC0Max1Min, x, y, z); 456 SetCorner(sC0Max1Min, x, y, z);
455 // corner of Axis0max and Axis1max 457 // corner of Axis0max and Axis1max
456 x = fAxisMax[rhoaxis]*std::cos(fAxisM 458 x = fAxisMax[rhoaxis]*std::cos(fAxisMax[phiaxis]);
457 y = fAxisMax[rhoaxis]*std::sin(fAxisM 459 y = fAxisMax[rhoaxis]*std::sin(fAxisMax[phiaxis]);
458 z = 0; 460 z = 0;
459 SetCorner(sC0Max1Max, x, y, z); 461 SetCorner(sC0Max1Max, x, y, z);
460 // corner of Axis0min and Axis1max 462 // corner of Axis0min and Axis1max
461 x = fAxisMin[rhoaxis]*std::cos(fAxisM 463 x = fAxisMin[rhoaxis]*std::cos(fAxisMax[phiaxis]);
462 y = fAxisMin[rhoaxis]*std::sin(fAxisM 464 y = fAxisMin[rhoaxis]*std::sin(fAxisMax[phiaxis]);
463 z = 0; 465 z = 0;
464 SetCorner(sC0Min1Max, x, y, z); 466 SetCorner(sC0Min1Max, x, y, z);
465 467
466 } 468 }
467 else 469 else
468 { 470 {
469 std::ostringstream message; 471 std::ostringstream message;
470 message << "Feature NOT implemented !" < 472 message << "Feature NOT implemented !" << G4endl
471 << " fAxis[0] = " << fAxi 473 << " fAxis[0] = " << fAxis[0] << G4endl
472 << " fAxis[1] = " << fAxi 474 << " fAxis[1] = " << fAxis[1];
473 G4Exception("G4TwistTubsFlatSide::SetCor 475 G4Exception("G4TwistTubsFlatSide::SetCorners()", "GeomSolids0001",
474 FatalException, message); 476 FatalException, message);
475 } 477 }
476 } 478 }
477 479
478 //============================================ 480 //=====================================================================
479 //* SetBoundaries() -------------------------- 481 //* SetBoundaries() ---------------------------------------------------
480 482
481 void G4TwistTubsFlatSide::SetBoundaries() 483 void G4TwistTubsFlatSide::SetBoundaries()
482 { 484 {
483 // Set direction-unit vector of phi-boundar 485 // Set direction-unit vector of phi-boundary-lines in local coodinate.
484 // Don't call the function twice. 486 // Don't call the function twice.
485 487
486 if (fAxis[0] == kRho && fAxis[1] == kPhi) 488 if (fAxis[0] == kRho && fAxis[1] == kPhi)
487 { 489 {
488 G4ThreeVector direction; 490 G4ThreeVector direction;
489 // sAxis0 & sAxisMin 491 // sAxis0 & sAxisMin
490 direction = GetCorner(sC0Min1Max) - GetC 492 direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
491 direction = direction.unit(); 493 direction = direction.unit();
492 SetBoundary(sAxis0 & (sAxisPhi | sAxisMi 494 SetBoundary(sAxis0 & (sAxisPhi | sAxisMin), direction,
493 GetCorner(sC0Min1Min), sAxis 495 GetCorner(sC0Min1Min), sAxisPhi);
494 496
495 // sAxis0 & sAxisMax 497 // sAxis0 & sAxisMax
496 direction = GetCorner(sC0Max1Max) - GetC 498 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
497 direction = direction.unit(); 499 direction = direction.unit();
498 SetBoundary(sAxis0 & (sAxisPhi | sAxisMa 500 SetBoundary(sAxis0 & (sAxisPhi | sAxisMax), direction,
499 GetCorner(sC0Max1Min), sAxis 501 GetCorner(sC0Max1Min), sAxisPhi);
500 502
501 // sAxis1 & sAxisMin 503 // sAxis1 & sAxisMin
502 direction = GetCorner(sC0Max1Min) - GetC 504 direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
503 direction = direction.unit(); 505 direction = direction.unit();
504 SetBoundary(sAxis1 & (sAxisRho | sAxisMi 506 SetBoundary(sAxis1 & (sAxisRho | sAxisMin), direction,
505 GetCorner(sC0Min1Min), sAxis 507 GetCorner(sC0Min1Min), sAxisRho);
506 508
507 // sAxis1 & sAxisMax 509 // sAxis1 & sAxisMax
508 direction = GetCorner(sC0Max1Max) - GetC 510 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
509 direction = direction.unit(); 511 direction = direction.unit();
510 SetBoundary(sAxis1 & (sAxisRho | sAxisMa 512 SetBoundary(sAxis1 & (sAxisRho | sAxisMax), direction,
511 GetCorner(sC0Min1Max), sAxis 513 GetCorner(sC0Min1Max), sAxisPhi);
512 } 514 }
513 else 515 else
514 { 516 {
515 std::ostringstream message; 517 std::ostringstream message;
516 message << "Feature NOT implemented !" < 518 message << "Feature NOT implemented !" << G4endl
517 << " fAxis[0] = " << fAxi 519 << " fAxis[0] = " << fAxis[0] << G4endl
518 << " fAxis[1] = " << fAxi 520 << " fAxis[1] = " << fAxis[1];
519 G4Exception("G4TwistTubsFlatSide::SetBou 521 G4Exception("G4TwistTubsFlatSide::SetBoundaries()", "GeomSolids0001",
520 FatalException, message); 522 FatalException, message);
521 } 523 }
522 } 524 }
523 525
524 //============================================ 526 //=====================================================================
525 //* GetFacets() ------------------------------ 527 //* GetFacets() -------------------------------------------------------
526 528
527 void G4TwistTubsFlatSide::GetFacets( G4int k, 529 void G4TwistTubsFlatSide::GetFacets( G4int k, G4int n, G4double xyz[][3],
528 G4int fac 530 G4int faces[][4], G4int iside )
529 { 531 {
530 G4ThreeVector p ; 532 G4ThreeVector p ;
531 533
532 G4double rmin = fAxisMin[0] ; 534 G4double rmin = fAxisMin[0] ;
533 G4double rmax = fAxisMax[0] ; 535 G4double rmax = fAxisMax[0] ;
534 G4double phimin, phimax ; 536 G4double phimin, phimax ;
535 537
536 G4double r,phi ; 538 G4double r,phi ;
537 G4int nnode,nface ; 539 G4int nnode,nface ;
538 540
539 for ( G4int i = 0 ; i<n ; ++i ) 541 for ( G4int i = 0 ; i<n ; ++i )
540 { 542 {
541 r = rmin + i*(rmax-rmin)/(n-1) ; 543 r = rmin + i*(rmax-rmin)/(n-1) ;
542 544
543 phimin = GetBoundaryMin(r) ; 545 phimin = GetBoundaryMin(r) ;
544 phimax = GetBoundaryMax(r) ; 546 phimax = GetBoundaryMax(r) ;
545 547
546 for ( G4int j = 0 ; j<k ; ++j ) 548 for ( G4int j = 0 ; j<k ; ++j )
547 { 549 {
548 phi = phimin + j*(phimax-phimin)/(k-1) ; 550 phi = phimin + j*(phimax-phimin)/(k-1) ;
549 551
550 nnode = GetNode(i,j,k,n,iside) ; 552 nnode = GetNode(i,j,k,n,iside) ;
551 p = SurfacePoint(phi,r,true) ; // surfa 553 p = SurfacePoint(phi,r,true) ; // surface point in global coord.system
552 554
553 xyz[nnode][0] = p.x() ; 555 xyz[nnode][0] = p.x() ;
554 xyz[nnode][1] = p.y() ; 556 xyz[nnode][1] = p.y() ;
555 xyz[nnode][2] = p.z() ; 557 xyz[nnode][2] = p.z() ;
556 558
557 if ( i<n-1 && j<k-1 ) // conterclock 559 if ( i<n-1 && j<k-1 ) // conterclock wise filling
558 { 560 {
559 nface = GetFace(i,j,k,n,iside) ; 561 nface = GetFace(i,j,k,n,iside) ;
560 562
561 if (fHandedness < 0) // lower side 563 if (fHandedness < 0) // lower side
562 { 564 {
563 faces[nface][0] = GetEdgeVisibility( 565 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,-1)
564 * ( GetNode(i ,j , 566 * ( GetNode(i ,j ,k,n,iside)+1) ;
565 faces[nface][1] = GetEdgeVisibility( 567 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,-1)
566 * ( GetNode(i ,j+1, 568 * ( GetNode(i ,j+1,k,n,iside)+1) ;
567 faces[nface][2] = GetEdgeVisibility( 569 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,-1)
568 * ( GetNode(i+1,j+1, 570 * ( GetNode(i+1,j+1,k,n,iside)+1) ;
569 faces[nface][3] = GetEdgeVisibility( 571 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,-1)
570 * ( GetNode(i+1,j , 572 * ( GetNode(i+1,j ,k,n,iside)+1) ;
571 } 573 }
572 else // upper side 574 else // upper side
573 { 575 {
574 faces[nface][0] = GetEdgeVisibility( 576 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1)
575 * ( GetNode(i ,j , 577 * ( GetNode(i ,j ,k,n,iside)+1) ;
576 faces[nface][1] = GetEdgeVisibility( 578 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1)
577 * ( GetNode(i+1,j , 579 * ( GetNode(i+1,j ,k,n,iside)+1) ;
578 faces[nface][2] = GetEdgeVisibility( 580 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1)
579 * ( GetNode(i+1,j+1, 581 * ( GetNode(i+1,j+1,k,n,iside)+1) ;
580 faces[nface][3] = GetEdgeVisibility( 582 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1)
581 * ( GetNode(i ,j+1, 583 * ( GetNode(i ,j+1,k,n,iside)+1) ;
582 } 584 }
583 } 585 }
584 } 586 }
585 } 587 }
586 } 588 }
587 589