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