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