source: trunk/source/geometry/solids/specific/src/G4TwistTubsFlatSide.cc@ 1202

Last change on this file since 1202 was 1058, checked in by garnier, 17 years ago

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