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