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

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

file release beta
File size: 44.7 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: G4VTwistSurface.cc,v 1.9 2007/05/31 13:52:48 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-02-ref-02 $
29//
30//
31// --------------------------------------------------------------------
32// GEANT 4 class source file
33//
34//
35// G4VTwistSurface.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 <iomanip>
46
47#include "G4VTwistSurface.hh"
48#include "G4GeometryTolerance.hh"
49
50const G4int G4VTwistSurface::sOutside = 0x00000000;
51const G4int G4VTwistSurface::sInside = 0x10000000;
52const G4int G4VTwistSurface::sBoundary = 0x20000000;
53const G4int G4VTwistSurface::sCorner = 0x40000000;
54const G4int G4VTwistSurface::sC0Min1Min = 0x40000101;
55const G4int G4VTwistSurface::sC0Max1Min = 0x40000201;
56const G4int G4VTwistSurface::sC0Max1Max = 0x40000202;
57const G4int G4VTwistSurface::sC0Min1Max = 0x40000102;
58const G4int G4VTwistSurface::sAxisMin = 0x00000101;
59const G4int G4VTwistSurface::sAxisMax = 0x00000202;
60const G4int G4VTwistSurface::sAxisX = 0x00000404;
61const G4int G4VTwistSurface::sAxisY = 0x00000808;
62const G4int G4VTwistSurface::sAxisZ = 0x00000C0C;
63const G4int G4VTwistSurface::sAxisRho = 0x00001010;
64const G4int G4VTwistSurface::sAxisPhi = 0x00001414;
65
66// mask
67const G4int G4VTwistSurface::sAxis0 = 0x0000FF00;
68const G4int G4VTwistSurface::sAxis1 = 0x000000FF;
69const G4int G4VTwistSurface::sSizeMask = 0x00000303;
70const G4int G4VTwistSurface::sAxisMask = 0x0000FCFC;
71const G4int G4VTwistSurface::sAreaMask = 0XF0000000;
72
73//=====================================================================
74//* constructors ------------------------------------------------------
75
76G4VTwistSurface::G4VTwistSurface(const G4String &name)
77 : fName(name)
78{
79
80 fAxis[0] = kUndefined;
81 fAxis[1] = kUndefined;
82 fAxisMin[0] = kInfinity;
83 fAxisMin[1] = kInfinity;
84 fAxisMax[0] = kInfinity;
85 fAxisMax[1] = kInfinity;
86 fHandedness = 1;
87
88 G4int i;
89 for (i=0; i<4; i++) {
90 fCorners[i].set(kInfinity, kInfinity, kInfinity);
91 fNeighbours[i] = 0;
92 }
93
94 fCurrentNormal.p.set(kInfinity, kInfinity, kInfinity);
95
96 fAmIOnLeftSide.me.set(kInfinity, kInfinity, kInfinity);
97 fAmIOnLeftSide.vec.set(kInfinity, kInfinity, kInfinity);
98 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
99}
100
101G4VTwistSurface::G4VTwistSurface(const G4String &name,
102 const G4RotationMatrix &rot,
103 const G4ThreeVector &tlate,
104 G4int handedness,
105 const EAxis axis0 ,
106 const EAxis axis1 ,
107 G4double axis0min,
108 G4double axis1min,
109 G4double axis0max,
110 G4double axis1max )
111 : fName(name)
112{
113 fAxis[0] = axis0;
114 fAxis[1] = axis1;
115 fAxisMin[0] = axis0min;
116 fAxisMin[1] = axis1min;
117 fAxisMax[0] = axis0max;
118 fAxisMax[1] = axis1max;
119 fHandedness = handedness;
120 fRot = rot;
121 fTrans = tlate;
122
123 G4int i;
124 for (i=0; i<4; i++) {
125 fCorners[i].set(kInfinity, kInfinity, kInfinity);
126 fNeighbours[i] = 0;
127 }
128
129 fCurrentNormal.p.set(kInfinity, kInfinity, kInfinity);
130
131 fAmIOnLeftSide.me.set(kInfinity, kInfinity, kInfinity);
132 fAmIOnLeftSide.vec.set(kInfinity, kInfinity, kInfinity);
133 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
134}
135
136//=====================================================================
137//* Fake default constructor ------------------------------------------
138
139G4VTwistSurface::G4VTwistSurface( __void__& )
140 : fName("")
141{
142}
143
144//=====================================================================
145//* destructor --------------------------------------------------------
146
147G4VTwistSurface::~G4VTwistSurface()
148{
149}
150
151//=====================================================================
152//* AmIOnLeftSide -----------------------------------------------------
153
154G4int G4VTwistSurface::AmIOnLeftSide(const G4ThreeVector &me,
155 const G4ThreeVector &vec,
156 G4bool withtol)
157{
158 // AmIOnLeftSide returns phi-location of "me"
159 // (phi relation between me and vec projected on z=0 plane).
160 // If "me" is on -ve-phi-side of "vec", it returns 1.
161 // On the other hand, if "me" is on +ve-phi-side of "vec",
162 // it returns -1.
163 // (The return value represents z-coordinate of normal vector
164 // of me.cross(vec).)
165 // If me is on boundary of vec, return 0.
166
167 static const G4double kAngTolerance
168 = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
169
170 static G4RotationMatrix unitrot; // unit matrix
171 static const G4RotationMatrix rottol = unitrot.rotateZ(0.5*kAngTolerance);
172 static const G4RotationMatrix invrottol = unitrot.rotateZ(-1.*kAngTolerance);
173
174 if (fAmIOnLeftSide.me == me
175 && fAmIOnLeftSide.vec == vec
176 && fAmIOnLeftSide.withTol == withtol) {
177 return fAmIOnLeftSide.amIOnLeftSide;
178 }
179
180 fAmIOnLeftSide.me = me;
181 fAmIOnLeftSide.vec = vec;
182 fAmIOnLeftSide.withTol = withtol;
183
184 G4ThreeVector met = (G4ThreeVector(me.x(), me.y(), 0.)).unit();
185 G4ThreeVector vect = (G4ThreeVector(vec.x(), vec.y(), 0.)).unit();
186
187 G4ThreeVector ivect = invrottol * vect;
188 G4ThreeVector rvect = rottol * vect;
189
190 G4double metcrossvect = met.x() * vect.y() - met.y() * vect.x();
191
192 if (withtol) {
193 if (met.x() * ivect.y() - met.y() * ivect.x() > 0 &&
194 metcrossvect >= 0) {
195 fAmIOnLeftSide.amIOnLeftSide = 1;
196 } else if (met.x() * rvect.y() - met.y() * rvect.x() < 0 &&
197 metcrossvect <= 0) {
198 fAmIOnLeftSide.amIOnLeftSide = -1;
199 } else {
200 fAmIOnLeftSide.amIOnLeftSide = 0;
201 }
202 } else {
203 if (metcrossvect > 0) {
204 fAmIOnLeftSide.amIOnLeftSide = 1;
205 } else if (metcrossvect < 0 ) {
206 fAmIOnLeftSide.amIOnLeftSide = -1;
207 } else {
208 fAmIOnLeftSide.amIOnLeftSide = 0;
209 }
210 }
211
212#ifdef G4TWISTDEBUG
213 G4cout << " === G4VTwistSurface::AmIOnLeftSide() =============="
214 << G4endl;
215 G4cout << " Name , returncode : " << fName << " "
216 << fAmIOnLeftSide.amIOnLeftSide << G4endl;
217 G4cout << " me, vec : " << std::setprecision(14) << me
218 << " " << vec << G4endl;
219 G4cout << " met, vect : " << met << " " << vect << G4endl;
220 G4cout << " ivec, rvec : " << ivect << " " << rvect << G4endl;
221 G4cout << " met x vect : " << metcrossvect << G4endl;
222 G4cout << " met x ivec : " << met.cross(ivect) << G4endl;
223 G4cout << " met x rvec : " << met.cross(rvect) << G4endl;
224 G4cout << " =============================================="
225 << G4endl;
226#endif
227
228 return fAmIOnLeftSide.amIOnLeftSide;
229}
230
231//=====================================================================
232//* DistanceToBoundary ------------------------------------------------
233
234G4double G4VTwistSurface::DistanceToBoundary(G4int areacode,
235 G4ThreeVector &xx,
236 const G4ThreeVector &p)
237{
238 // DistanceToBoundary
239 //
240 // return distance to nearest boundary from arbitrary point p
241 // in local coodinate.
242 // Argument areacode must be one of them:
243 // sAxis0 & sAxisMin, sAxis0 & sAxisMax,
244 // sAxis1 & sAxisMin, sAxis1 & sAxisMax.
245 //
246
247 G4ThreeVector d; // direction vector of the boundary
248 G4ThreeVector x0; // reference point of the boundary
249 G4double dist = kInfinity;
250 G4int boundarytype;
251
252 if (IsAxis0(areacode) && IsAxis1(areacode)) {
253 G4cerr << "ERROR - G4VTwistSurface::DistanceToBoundary()" << G4endl
254 << " Point is in the corner area. This function returns"
255 << G4endl
256 << " a direction vector of a boundary line." << G4endl
257 << " areacode = " << areacode << G4endl;
258 G4Exception("G4VTwistSurface::DistanceToBoundary()", "InvalidSetup",
259 FatalException, "Point is in the corner area.");
260 } else if (IsAxis0(areacode) || IsAxis1(areacode)) {
261 GetBoundaryParameters(areacode, d, x0, boundarytype);
262 if (boundarytype == sAxisPhi) {
263 G4double t = x0.getRho() / p.getRho();
264 xx.set(t*p.x(), t*p.y(), x0.z());
265 dist = (xx - p).mag();
266 } else {
267 // linear boundary
268 // sAxisX, sAxisY, sAxisZ, sAxisRho
269 dist = DistanceToLine(p, x0, d, xx);
270 }
271 } else {
272 G4cerr << "ERROR - G4VTwistSurface::DistanceToBoundary()" << G4endl
273 << " areacode = " << areacode << G4endl;
274 G4Exception("G4VTwistSurface::DistanceToBoundary()", "InvalidSetup",
275 FatalException, "Bad areacode of boundary.");
276 }
277 return dist;
278}
279
280//=====================================================================
281//* DistanceToIn ------------------------------------------------------
282
283G4double G4VTwistSurface::DistanceToIn(const G4ThreeVector &gp,
284 const G4ThreeVector &gv,
285 G4ThreeVector &gxxbest)
286{
287#ifdef G4TWISTDEBUG
288 G4cout << " ~~~~~ G4VTwistSurface::DistanceToIn(p,v) - Start ~~~~~" << G4endl;
289 G4cout << " Name : " << fName << G4endl;
290 G4cout << " gp : " << gp << G4endl;
291 G4cout << " gv : " << gv << G4endl;
292 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
293#endif
294
295 G4ThreeVector gxx[G4VSURFACENXX];
296 G4double distance[G4VSURFACENXX] ;
297 G4int areacode[G4VSURFACENXX] ;
298 G4bool isvalid[G4VSURFACENXX] ;
299
300 for (G4int i = 0 ; i<G4VSURFACENXX ; i++ ) {
301 distance[i] = kInfinity ;
302 areacode[i] = sOutside ;
303 isvalid[i] = false ;
304 }
305
306 G4double bestdistance = kInfinity;
307 G4int besti = -1;
308 G4ThreeVector bestgxx(kInfinity, kInfinity, kInfinity);
309
310 G4int nxx = DistanceToSurface(gp, gv, gxx, distance, areacode,
311 isvalid, kValidateWithTol);
312
313 for (G4int i=0; i< nxx; i++) {
314
315 // skip this intersection if:
316 // - invalid intersection
317 // - particle goes outword the surface
318
319 if (!isvalid[i]) {
320 // xx[i] is sOutside or distance[i] < 0
321 continue;
322 }
323
324 G4ThreeVector normal = GetNormal(gxx[i], true);
325
326 if ((normal * gv) >= 0) {
327
328#ifdef G4TWISTDEBUG
329 G4cout << " G4VTwistSurface::DistanceToIn(p,v): "
330 << "particle goes outword the surface." << G4endl;
331#endif
332 continue;
333 }
334
335 //
336 // accept this intersection if the intersection is inside.
337 //
338
339 if (IsInside(areacode[i])) {
340 if (distance[i] < bestdistance) {
341 bestdistance = distance[i];
342 bestgxx = gxx[i];
343 besti = i;
344
345#ifdef G4TWISTDEBUG
346 G4cout << " G4VTwistSurface::DistanceToIn(p,v): "
347 << " areacode sInside name, distance = "
348 << fName << " "<< bestdistance << G4endl;
349#endif
350 }
351
352 //
353 // else, the intersection is on boundary or corner.
354 //
355
356 } else {
357
358 G4VTwistSurface *neighbours[2];
359 G4bool isaccepted[2] = {false, false};
360 G4int nneighbours = GetNeighbours(areacode[i], neighbours);
361
362 for (G4int j=0; j< nneighbours; j++) {
363 // if on corner, nneighbours = 2.
364 // if on boundary, nneighbours = 1.
365
366 G4ThreeVector tmpgxx[G4VSURFACENXX];
367 G4double tmpdist[G4VSURFACENXX] ;
368 G4int tmpareacode[G4VSURFACENXX] ;
369 G4bool tmpisvalid[G4VSURFACENXX] ;
370
371 for (G4int l = 0 ; l<G4VSURFACENXX ; l++ ) {
372 tmpdist[l] = kInfinity ;
373 tmpareacode[l] = sOutside ;
374 tmpisvalid[l] = false ;
375 }
376
377 G4int tmpnxx = neighbours[j]->DistanceToSurface(
378 gp, gv, tmpgxx, tmpdist,
379 tmpareacode, tmpisvalid,
380 kValidateWithTol);
381 G4ThreeVector neighbournormal;
382
383 for (G4int k=0; k< tmpnxx; k++) {
384
385 //
386 // if tmpxx[k] is valid && sInside, the final winner must
387 // be neighbour surface. return kInfinity.
388 // else , choose tmpxx on same boundary of xx, then check normal
389 //
390
391 if (IsInside(tmpareacode[k])) {
392
393#ifdef G4TWISTDEBUG
394 G4cout << " G4VTwistSurface:DistanceToIn(p,v): "
395 << " intersection "<< tmpgxx[k] << G4endl
396 << " is inside of neighbour surface of " << fName
397 << " . returning kInfinity." << G4endl;
398 G4cout << "~~~~~ G4VTwistSurface::DistanceToIn(p,v) - return ~~~~"
399 << G4endl;
400 G4cout << " No intersections " << G4endl;
401 G4cout << " Name : " << fName << G4endl;
402 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
403 << G4endl;
404#endif
405 if (tmpisvalid[k]) return kInfinity;
406 continue;
407
408 //
409 // if tmpxx[k] is valid && sInside, the final winner must
410 // be neighbour surface. return .
411 //
412
413 } else if (IsSameBoundary(this,areacode[i],
414 neighbours[j], tmpareacode[k])) {
415 // tmpxx[k] is same boundary (or corner) of xx.
416
417 neighbournormal = neighbours[j]->GetNormal(tmpgxx[k], true);
418 if (neighbournormal * gv < 0) isaccepted[j] = true;
419 }
420 }
421
422 // if nneighbours = 1, chabge isaccepted[1] before
423 // exiting neighboursurface loop.
424
425 if (nneighbours == 1) isaccepted[1] = true;
426
427 } // neighboursurface loop end
428
429 // now, we can accept xx intersection
430
431 if (isaccepted[0] == true && isaccepted[1] == true) {
432 if (distance[i] < bestdistance) {
433 bestdistance = distance[i];
434 gxxbest = gxx[i];
435 besti = i;
436#ifdef G4TWISTDEBUG
437 G4cout << " G4VTwistSurface::DistanceToIn(p,v): "
438 << " areacode sBoundary & sBoundary distance = "
439 << fName << " " << distance[i] << G4endl;
440#endif
441 }
442 }
443
444 } // else end
445 } // intersection loop end
446
447 gxxbest = bestgxx;
448
449#ifdef G4TWISTDEBUG
450 if (besti < 0) {
451 G4cout << "~~~~~ G4VTwistSurface::DistanceToIn(p,v) - return ~~~~" << G4endl;
452 G4cout << " No intersections " << G4endl;
453 G4cout << " Name : " << fName << G4endl;
454 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
455 } else {
456 G4cout << "~~~~~ G4VTwistSurface::DistanceToIn(p,v) : return ~~~~" << G4endl;
457 G4cout << " Name, i : " << fName << " , " << besti << G4endl;
458 G4cout << " gxx[i] : " << gxxbest << G4endl;
459 G4cout << " bestdist : " << bestdistance << G4endl;
460 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
461 }
462
463#endif
464
465 return bestdistance;
466}
467
468//=====================================================================
469//* DistanceToOut(p, v) -----------------------------------------------
470
471G4double G4VTwistSurface::DistanceToOut(const G4ThreeVector &gp,
472 const G4ThreeVector &gv,
473 G4ThreeVector &gxxbest)
474{
475#ifdef G4TWISTDEBUG
476 G4cout << "~~~~~ G4VTwistSurface::DistanceToOut(p,v) - Start ~~~~" << G4endl;
477 G4cout << " Name : " << fName << G4endl;
478 G4cout << " gp : " << gp << G4endl;
479 G4cout << " gv : " << gv << G4endl;
480 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
481#endif
482
483 G4ThreeVector gxx[G4VSURFACENXX];
484 G4double distance[G4VSURFACENXX] ;
485 G4int areacode[G4VSURFACENXX] ;
486 G4bool isvalid[G4VSURFACENXX] ;
487 G4int i;
488
489 for ( i = 0 ; i<G4VSURFACENXX ; i++ )
490 {
491 distance[i] = kInfinity ;
492 areacode[i] = sOutside ;
493 isvalid[i] = false ;
494 }
495
496 G4int nxx;
497 G4double bestdistance = kInfinity;
498 G4int besti = -1;
499
500 nxx = DistanceToSurface(gp, gv, gxx, distance, areacode,
501 isvalid, kValidateWithTol);
502
503 for (i=0; i<nxx; i++) {
504 if (!(isvalid[i])) {
505 continue;
506 }
507
508 G4ThreeVector normal = GetNormal(gxx[i], true);
509 if (normal * gv <= 0) {
510 // particle goes toword inside of solid, return kInfinity
511#ifdef G4TWISTDEBUG
512 G4cout << " G4VTwistSurface::DistanceToOut(p,v): normal*gv < 0, normal "
513 << fName << " " << normal
514 << G4endl;
515#endif
516 } else {
517 // gxx[i] is accepted.
518 if (distance[i] < bestdistance) {
519 bestdistance = distance[i];
520 gxxbest = gxx[i];
521 besti = i;
522 }
523 }
524 }
525
526#ifdef G4TWISTDEBUG
527 if (besti < 0) {
528 G4cout << "~~~~~ G4VTwistSurface::DistanceToOut(p,v) - return ~~~" << G4endl;
529 G4cout << " No intersections " << G4endl;
530 G4cout << " Name : " << fName << G4endl;
531 G4cout << " bestdist : " << bestdistance << G4endl;
532 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
533 } else {
534 G4cout << "~~~~~ G4VTwistSurface::DistanceToOut(p,v) : return ~~~" << G4endl;
535 G4cout << " Name, i : " << fName << " , " << i << G4endl;
536 G4cout << " gxx[i] : " << gxxbest << G4endl;
537 G4cout << " bestdist : " << bestdistance << G4endl;
538 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
539 }
540#endif
541
542 return bestdistance;
543}
544
545//=====================================================================
546//* DistanceTo(p) -----------------------------------------------------
547
548G4double G4VTwistSurface::DistanceTo(const G4ThreeVector &gp,
549 G4ThreeVector &gxxbest)
550{
551#ifdef G4TWISTDEBUG
552 G4cout << "~~~~~ G4VTwistSurface::DistanceTo(p) - Start ~~~~~~~~~" << G4endl;
553 G4cout << " Name : " << fName << G4endl;
554 G4cout << " gp : " << gp << G4endl;
555 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
556#endif
557
558
559 G4ThreeVector gxx[G4VSURFACENXX];
560 G4double distance[G4VSURFACENXX] ;
561 G4int areacode[G4VSURFACENXX] ;
562
563 for (G4int i = 0 ; i<G4VSURFACENXX ; i++ ) {
564 distance[i] = kInfinity ;
565 areacode[i] = sOutside ;
566 }
567
568 G4int nxx;
569
570 nxx = DistanceToSurface(gp, gxx, distance, areacode);
571 gxxbest = gxx[0];
572
573#ifdef G4TWISTDEBUG
574 G4cout << "~~~~~ G4VTwistSurface::DistanceTo(p) - return ~~~~~~~~" << G4endl;
575 G4cout << " Name : " << fName << G4endl;
576 G4cout << " gxx : " << gxxbest << G4endl;
577 G4cout << " bestdist : " << distance[0] << G4endl;
578 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
579#endif
580
581 return distance[0];
582}
583
584//=====================================================================
585//* IsSameBoundary ----------------------------------------------------
586
587G4bool G4VTwistSurface::IsSameBoundary(G4VTwistSurface *surface1, G4int areacode1,
588 G4VTwistSurface *surface2, G4int areacode2 ) const
589{
590 //
591 // IsSameBoundary
592 //
593 // checking tool whether two boundaries on different surfaces are same or not.
594 //
595
596 G4bool testbitmode = true;
597 G4bool iscorner[2] = {IsCorner(areacode1, testbitmode),
598 IsCorner(areacode2, testbitmode)};
599
600 if (iscorner[0] && iscorner[1]) {
601 // on corner
602 G4ThreeVector corner1 =
603 surface1->ComputeGlobalPoint(surface1->GetCorner(areacode1));
604 G4ThreeVector corner2 =
605 surface2->ComputeGlobalPoint(surface2->GetCorner(areacode2));
606
607 if ((corner1 - corner2).mag() < kCarTolerance) {
608 return true;
609 } else {
610 return false;
611 }
612
613 } else if ((IsBoundary(areacode1, testbitmode) && (!iscorner[0])) &&
614 (IsBoundary(areacode2, testbitmode) && (!iscorner[1]))) {
615 // on boundary
616 G4ThreeVector d1, d2, ld1, ld2;
617 G4ThreeVector x01, x02, lx01, lx02;
618 G4int type1, type2;
619 surface1->GetBoundaryParameters(areacode1, ld1, lx01, type1);
620 surface2->GetBoundaryParameters(areacode2, ld2, lx02, type2);
621
622 x01 = surface1->ComputeGlobalPoint(lx01);
623 x02 = surface2->ComputeGlobalPoint(lx02);
624 d1 = surface1->ComputeGlobalDirection(ld1);
625 d2 = surface2->ComputeGlobalDirection(ld2);
626
627 if ((x01 - x02).mag() < kCarTolerance &&
628 (d1 - d2).mag() < kCarTolerance) {
629 return true;
630 } else {
631 return false;
632 }
633
634 } else {
635 return false;
636 }
637}
638
639//=====================================================================
640//* GetBoundaryParameters ---------------------------------------------
641
642void G4VTwistSurface::GetBoundaryParameters(const G4int &areacode,
643 G4ThreeVector &d,
644 G4ThreeVector &x0,
645 G4int &boundarytype) const
646{
647 // areacode must be one of them:
648 // sAxis0 & sAxisMin, sAxis0 & sAxisMax,
649 // sAxis1 & sAxisMin, sAxis1 & sAxisMax.
650
651 G4int i;
652 for (i=0; i<4; i++) {
653 if (fBoundaries[i].GetBoundaryParameters(areacode, d, x0,
654 boundarytype)) {
655 return;
656 }
657 }
658
659 G4cerr << "ERROR - G4VTwistSurface::GetBoundaryParameters()" << G4endl
660 << " Boundary at areacode " << std::hex << areacode
661 << std::dec << G4endl
662 << " is not be registered." << G4endl;
663 G4Exception("G4VTwistSurface::GetBoundaryParameters()", "InvalidSetup",
664 FatalException, "Not registered boundary.");
665}
666
667//=====================================================================
668//* GetBoundaryAtPZ ---------------------------------------------------
669
670G4ThreeVector G4VTwistSurface::GetBoundaryAtPZ(G4int areacode,
671 const G4ThreeVector &p) const
672{
673 // areacode must be one of them:
674 // sAxis0 & sAxisMin, sAxis0 & sAxisMax,
675 // sAxis1 & sAxisMin, sAxis1 & sAxisMax.
676
677 if (areacode & sAxis0 && areacode & sAxis1) {
678 G4cerr << "ERROR - G4VTwistSurface::GetBoundaryAtPZ()" << G4endl
679 << " Point is in the corner area. This function returns"
680 << G4endl
681 << " a direction vector of a boundary line." << G4endl
682 << " areacode = " << areacode << G4endl;
683 G4Exception("G4VTwistSurface::GetBoundaryAtPZ()", "InvalidCondition",
684 FatalException, "Point is in the corner area.");
685 }
686
687 G4ThreeVector d;
688 G4ThreeVector x0;
689 G4int boundarytype;
690 G4bool found = false;
691
692 for (G4int i=0; i<4; i++) {
693 if (fBoundaries[i].GetBoundaryParameters(areacode, d, x0,
694 boundarytype)){
695 found = true;
696 continue;
697 }
698 }
699
700 if (!found) {
701 G4cerr << "ERROR - G4VTwistSurface::GetBoundaryAtPZ()" << G4endl
702 << " Boundary at areacode " << areacode << G4endl
703 << " is not be registered." << G4endl;
704 G4Exception("G4VTwistSurface::GetBoundaryAtPZ()", "InvalidSetup",
705 FatalException, "Not registered boundary.");
706 }
707
708 if (((boundarytype & sAxisPhi) == sAxisPhi) ||
709 ((boundarytype & sAxisRho) == sAxisRho)) {
710 G4cerr << "ERROR - G4VTwistSurface::GetBoundaryAtPZ()" << G4endl
711 << " Boundary at areacode " << areacode << G4endl
712 << " is not a z-depended line." << G4endl;
713 G4Exception("G4VTwistSurface::GetBoundaryAtPZ()", "InvalidSetup",
714 FatalException, "Not a z-depended line boundary.");
715 }
716 return ((p.z() - x0.z()) / d.z()) * d + x0;
717}
718
719//=====================================================================
720//* SetCorner ---------------------------------------------------------
721
722void G4VTwistSurface::SetCorner(G4int areacode, G4double x, G4double y, G4double z)
723{
724 if ((areacode & sCorner) != sCorner){
725 G4cerr << "ERROR - G4VTwistSurface::SetCorner()" << G4endl
726 << " areacode " << areacode << G4endl;
727 G4Exception("G4VTwistSurface::SetCorner()", "InvalidSetup",
728 FatalException, "Area code must represents corner.");
729 }
730
731 if ((areacode & sC0Min1Min) == sC0Min1Min) {
732 fCorners[0].set(x, y, z);
733 } else if ((areacode & sC0Max1Min) == sC0Max1Min) {
734 fCorners[1].set(x, y, z);
735 } else if ((areacode & sC0Max1Max) == sC0Max1Max) {
736 fCorners[2].set(x, y, z);
737 } else if ((areacode & sC0Min1Max) == sC0Min1Max) {
738 fCorners[3].set(x, y, z);
739 }
740}
741
742//=====================================================================
743//* SetBoundaryAxis ---------------------------------------------------
744
745void G4VTwistSurface::GetBoundaryAxis(G4int areacode, EAxis axis[]) const
746{
747 if ((areacode & sBoundary) != sBoundary) {
748 G4Exception("G4VTwistSurface::GetBoundaryAxis()", "InvalidCondition",
749 FatalException, "Not located on a boundary!");
750 }
751 G4int i;
752 for (i=0; i<2; i++) {
753
754 G4int whichaxis = 0 ;
755 if (i == 0) {
756 whichaxis = sAxis0;
757 } else if (i == 1) {
758 whichaxis = sAxis1;
759 }
760
761 // extracted axiscode of whichaxis
762 G4int axiscode = whichaxis & sAxisMask & areacode ;
763 if (axiscode) {
764 if (axiscode == (whichaxis & sAxisX)) {
765 axis[i] = kXAxis;
766 } else if (axiscode == (whichaxis & sAxisY)) {
767 axis[i] = kYAxis;
768 } else if (axiscode == (whichaxis & sAxisZ)) {
769 axis[i] = kZAxis;
770 } else if (axiscode == (whichaxis & sAxisRho)) {
771 axis[i] = kRho;
772 } else if (axiscode == (whichaxis & sAxisPhi)) {
773 axis[i] = kPhi;
774 } else {
775 G4cerr << "ERROR - G4VTwistSurface::GetBoundaryAxis()" << G4endl
776 << " areacode " << areacode << G4endl;
777 G4Exception("G4VTwistSurface::GetBoundaryAxis()", "InvalidSetup",
778 FatalException, "Not supported areacode.");
779 }
780 }
781 }
782}
783
784//=====================================================================
785//* SetBoundaryLimit --------------------------------------------------
786
787void G4VTwistSurface::GetBoundaryLimit(G4int areacode, G4double limit[]) const
788{
789 if (areacode & sCorner) {
790 if (areacode & sC0Min1Max) {
791 limit[0] = fAxisMin[0];
792 limit[1] = fAxisMin[1];
793 } else if (areacode & sC0Max1Min) {
794 limit[0] = fAxisMax[0];
795 limit[1] = fAxisMin[1];
796 } else if (areacode & sC0Max1Max) {
797 limit[0] = fAxisMax[0];
798 limit[1] = fAxisMax[1];
799 } else if (areacode & sC0Min1Max) {
800 limit[0] = fAxisMin[0];
801 limit[1] = fAxisMax[1];
802 }
803 } else if (areacode & sBoundary) {
804 if (areacode & (sAxis0 | sAxisMin)) {
805 limit[0] = fAxisMin[0];
806 } else if (areacode & (sAxis1 | sAxisMin)) {
807 limit[0] = fAxisMin[1];
808 } else if (areacode & (sAxis0 | sAxisMax)) {
809 limit[0] = fAxisMax[0];
810 } else if (areacode & (sAxis1 | sAxisMax)) {
811 limit[0] = fAxisMax[1];
812 }
813 } else {
814 G4cerr << "WARNING - G4VTwistSurface::GetBoundaryAxis()" << G4endl
815 << " areacode " << areacode << G4endl;
816 G4Exception("G4VTwistSurface::GetBoundaryLimit()", "InvalidCondition",
817 JustWarning, "Not located on a boundary!");
818 }
819}
820
821//=====================================================================
822//* SetBoundary -------------------------------------------------------
823
824void G4VTwistSurface::SetBoundary(const G4int &axiscode,
825 const G4ThreeVector &direction,
826 const G4ThreeVector &x0,
827 const G4int &boundarytype)
828{
829 G4int code = (~sAxisMask) & axiscode;
830 if ((code == (sAxis0 & sAxisMin)) ||
831 (code == (sAxis0 & sAxisMax)) ||
832 (code == (sAxis1 & sAxisMin)) ||
833 (code == (sAxis1 & sAxisMax))) {
834
835 G4int i;
836 G4bool done = false;
837 for (i=0; i<4; i++) {
838 if (fBoundaries[i].IsEmpty()) {
839 fBoundaries[i].SetFields(axiscode, direction,
840 x0, boundarytype);
841 done = true;
842 break;
843 }
844 }
845
846 if (!done) {
847 G4Exception("G4VTwistSurface::SetBoundary()", "InvalidCondition",
848 FatalException, "Number of boundary exceeding 4!");
849 }
850 } else {
851 G4cerr << "ERROR - G4VTwistSurface::SetBoundary()" << G4endl
852 << " invalid axiscode. axiscode = "
853 << std::hex << axiscode << std::dec << G4endl;
854 G4Exception("G4VTwistSurface::SetBoundary()", "InvalidCondition",
855 FatalException, "Invalid axis-code.");
856 }
857}
858
859//=====================================================================
860//* GetFace -----------------------------------------------------------
861
862G4int G4VTwistSurface::GetFace( G4int i, G4int j, G4int m,
863 G4int n, G4int iside )
864{
865 // this is the face mapping function
866 // (i,j) -> face number
867
868 if ( iside == 0 ) {
869 return i * ( m - 1 ) + j ;
870 }
871
872 else if ( iside == 1 ) {
873 return (m-1)*(m-1) + i*(m-1) + j ;
874 }
875
876 else if ( iside == 2 ) {
877 return 2*(m-1)*(m-1) + i*(m-1) + j ;
878 }
879
880 else if ( iside == 3 ) {
881 return 2*(m-1)*(m-1) + (n-1)*(m-1) + i*(m-1) + j ;
882 }
883
884 else if ( iside == 4 ) {
885 return 2*(m-1)*(m-1) + 2*(n-1)*(m-1) + i*(m-1) + j ;
886 }
887
888 else if ( iside == 5 ) {
889 return 2*(m-1)*(m-1) + 3*(n-1)*(m-1) + i*(m-1) + j ;
890 }
891
892 else {
893
894 G4cerr << "ERROR - G4VTwistSurface::GetFace(): "
895 << GetName() << G4endl
896 << " Not correct side number! - " << G4endl
897 << "iside is " << iside << " but should be "
898 << "0,1,2,3,4 or 5" << "." << G4endl ;
899 G4Exception("G4TwistSurface::G4GetFace()", "InvalidSetup",
900 FatalException, "Not correct side number.");
901
902
903 }
904
905 return -1 ; // wrong face
906}
907
908//=====================================================================
909//* GetNode -----------------------------------------------------------
910
911G4int G4VTwistSurface::GetNode( G4int i, G4int j, G4int m,
912 G4int n, G4int iside )
913{
914 // this is the node mapping function
915 // (i,j) -> node number
916 // Depends on the side iside and the used meshing of the surface
917
918 if ( iside == 0 ) {
919 // lower endcap is mxm squared.
920 // n = m
921 return i * m + j ;
922 }
923
924 if ( iside == 1 ) {
925 // upper endcap is mxm squared. Shift by m*m
926 // n = m
927 return m*m + i*m + j ;
928 }
929
930 else if ( iside == 2 ) {
931 // front side.
932 if ( i == 0 ) { return j ; }
933 else if ( i == n-1 ) { return m*m + j ; }
934 else { return 2*m*m + 4*(i-1)*(m-1) + j ; }
935 }
936
937 else if ( iside == 3 ) {
938 // right side
939 if ( i == 0 ) { return (j+1)*m - 1 ; }
940 else if ( i == n-1 ) { return m*m + (j+1)*m - 1 ; }
941 else { return 2*m*m + 4*(i-1)*(m-1) + (m-1) + j ; }
942 }
943 else if ( iside == 4 ) {
944 // back side.
945 if ( i == 0 ) { return m*m - 1 - j ; } // reversed order
946 else if ( i == n-1 ) { return 2*m*m - 1 - j ; } // reversed order
947 else { return 2*m*m + 4*(i-1)*(m-1) + 2*(m-1) + j ; // normal order
948 }
949 }
950 else if ( iside == 5 ) {
951 // left side
952 if ( i == 0 ) { return m*m - (j+1)*m ; } // reversed order
953 else if ( i == n-1) { return 2*m*m - (j+1)*m ; } // reverded order
954 else {
955 if ( j == m-1 ) { return 2*m*m + 4*(i-1)*(m-1) ; } // special case
956 else { return 2*m*m + 4*(i-1)*(m-1) + 3*(m-1) + j ; } // normal order
957 }
958 }
959
960 else {
961
962 G4cerr << "ERROR - G4VTwistSurface::GetNode(): "
963 << GetName() << G4endl
964 << " Not correct side number! - " << G4endl
965 << "iside is " << iside << " but should be "
966 << "0,1,2,3,4 or 5" << "." << G4endl ;
967 G4Exception("G4TwistSurface::G4GetNode()", "InvalidSetup",
968 FatalException, "Not correct side number.");
969 }
970 return -1 ; // wrong node
971}
972
973//=====================================================================
974//* GetEdgeVisiblility ------------------------------------------------
975
976G4int G4VTwistSurface::GetEdgeVisibility( G4int i, G4int j, G4int m, G4int n, G4int number, G4int orientation)
977{
978
979 // clockwise filling -> positive orientation
980 // counter clockwise filling -> negative orientation
981
982 //
983 // d C c
984 // +------+
985 // | |
986 // | |
987 // | |
988 // D | |B
989 // | |
990 // | |
991 // | |
992 // +------+
993 // a A b
994 //
995 // a = +--+ A = ---+
996 // b = --++ B = --+-
997 // c = -++- C = -+--
998 // d = ++-- D = +---
999
1000
1001 // check first invisible faces
1002
1003 if ( ( i>0 && i<n-2 ) && ( j>0 && j<m-2 ) ) {
1004 return -1 ; // always invisible, signs: ----
1005 }
1006
1007 // change first the vertex number (depends on the orientation)
1008 // 0,1,2,3 -> 3,2,1,0
1009 if ( orientation < 0 ) { number = ( 3 - number ) ; }
1010
1011 // check true edges
1012 if ( ( j>=1 && j<=m-3 ) ) {
1013
1014 if ( i == 0 ) { // signs (A): ---+
1015 return ( number == 3 ) ? 1 : -1 ;
1016 }
1017
1018 else if ( i == n-2 ) { // signs (C): -+--
1019 return ( number == 1 ) ? 1 : -1 ;
1020 }
1021
1022 else {
1023 G4cerr << "ERROR - G4VTwistSurface::GetEdgeVisibliity(): "
1024 << GetName() << G4endl
1025 << " Not correct face number! - " << G4endl ;
1026 G4Exception("G4TwistSurface::G4GetEdgeVisibility()", "InvalidSetup",
1027 FatalException, "Not correct face number.");
1028 }
1029 }
1030
1031 if ( ( i>=1 && i<=n-3 ) ) {
1032
1033 if ( j == 0 ) { // signs (D): +---
1034 return ( number == 0 ) ? 1 : -1 ;
1035 }
1036
1037 else if ( j == m-2 ) { // signs (B): --+-
1038 return ( number == 2 ) ? 1 : -1 ;
1039 }
1040
1041 else {
1042 G4cerr << "ERROR - G4VTwistSurface::GetEdgeVisibliity(): "
1043 << GetName() << G4endl
1044 << " Not correct face number! - " << G4endl ;
1045 G4Exception("G4TwistSurface::G4GetEdgeVisibility()", "InvalidSetup",
1046 FatalException, "Not correct face number.");
1047 }
1048 }
1049
1050 // now the corners
1051 if ( i == 0 && j == 0 ) { // signs (a) : +--+
1052 return ( number == 0 || number == 3 ) ? 1 : -1 ;
1053 }
1054 else if ( i == 0 && j == m-2 ) { // signs (b) : --++
1055 return ( number == 2 || number == 3 ) ? 1 : -1 ;
1056 }
1057 else if ( i == n-2 && j == m-2 ) { // signs (c) : -++-
1058 return ( number == 1 || number == 2 ) ? 1 : -1 ;
1059 }
1060 else if ( i == n-2 && j == j ) { // signs (d) : ++--
1061 return ( number == 0 || number == 1 ) ? 1 : -1 ;
1062 }
1063 else {
1064 G4cerr << "ERROR - G4VTwistSurface::GetEdgeVisibliity(): "
1065 << GetName() << G4endl
1066 << " Not correct face number! - " << G4endl ;
1067 G4Exception("G4TwistSurface::G4GetEdgeVisibility()", "InvalidSetup",
1068 FatalException, "Not correct face number.");
1069 }
1070
1071 G4cerr << "ERROR - G4VTwistSurface::GetEdgeVisibliity(): "
1072 << GetName() << G4endl
1073 << " Not correct face number! - " << G4endl ;
1074 G4Exception("G4TwistSurface::G4GetEdgeVisibility()", "InvalidSetup",
1075 FatalException, "Not correct face number.");
1076
1077 return 0 ;
1078
1079}
1080
1081
1082//=====================================================================
1083//* DebugPrint --------------------------------------------------------
1084
1085void G4VTwistSurface::DebugPrint() const
1086{
1087 G4ThreeVector A = fRot * GetCorner(sC0Min1Min) + fTrans;
1088 G4ThreeVector B = fRot * GetCorner(sC0Max1Min) + fTrans;
1089 G4ThreeVector C = fRot * GetCorner(sC0Max1Max) + fTrans;
1090 G4ThreeVector D = fRot * GetCorner(sC0Min1Max) + fTrans;
1091
1092 G4cout << "/* G4VTwistSurface::DebugPrint():-------------------------------"
1093 << G4endl;
1094 G4cout << "/* Name = " << fName << G4endl;
1095 G4cout << "/* Axis = " << std::hex << fAxis[0] << " "
1096 << std::hex << fAxis[1]
1097 << " (0,1,2,3,5 = kXAxis,kYAxis,kZAxis,kRho,kPhi)"
1098 << std::dec << G4endl;
1099 G4cout << "/* BoundaryLimit(in local) fAxis0(min, max) = ("<<fAxisMin[0]
1100 << ", " << fAxisMax[0] << ")" << G4endl;
1101 G4cout << "/* BoundaryLimit(in local) fAxis1(min, max) = ("<<fAxisMin[1]
1102 << ", " << fAxisMax[1] << ")" << G4endl;
1103 G4cout << "/* Cornar point sC0Min1Min = " << A << G4endl;
1104 G4cout << "/* Cornar point sC0Max1Min = " << B << G4endl;
1105 G4cout << "/* Cornar point sC0Max1Max = " << C << G4endl;
1106 G4cout << "/* Cornar point sC0Min1Max = " << D << G4endl;
1107 G4cout << "/*---------------------------------------------------------"
1108 << G4endl;
1109}
1110
1111//=====================================================================
1112// G4VTwistSurface::CurrentStatus class
1113//=====================================================================
1114
1115//=====================================================================
1116//* CurrentStatus::CurrentStatus --------------------------------------
1117
1118G4VTwistSurface::CurrentStatus::CurrentStatus()
1119{
1120 for (size_t i=0; i<G4VSURFACENXX; i++)
1121 {
1122 fDistance[i] = kInfinity;
1123 fAreacode[i] = sOutside;
1124 fIsValid[i] = false;
1125 fXX[i].set(kInfinity, kInfinity, kInfinity);
1126 }
1127 fNXX = 0;
1128 fLastp.set(kInfinity, kInfinity, kInfinity);
1129 fLastv.set(kInfinity, kInfinity, kInfinity);
1130 fLastValidate = kUninitialized;
1131 fDone = false;
1132}
1133
1134//=====================================================================
1135//* CurrentStatus::~CurrentStatus -------------------------------------
1136
1137G4VTwistSurface::CurrentStatus::~CurrentStatus()
1138{
1139}
1140
1141//=====================================================================
1142//* CurrentStatus::SetCurrentStatus -----------------------------------
1143
1144void
1145G4VTwistSurface::CurrentStatus::SetCurrentStatus(G4int i,
1146 G4ThreeVector &xx,
1147 G4double &dist,
1148 G4int &areacode,
1149 G4bool &isvalid,
1150 G4int nxx,
1151 EValidate validate,
1152 const G4ThreeVector *p,
1153 const G4ThreeVector *v)
1154{
1155 fDistance[i] = dist;
1156 fAreacode[i] = areacode;
1157 fIsValid[i] = isvalid;
1158 fXX[i] = xx;
1159 fNXX = nxx;
1160 fLastValidate = validate;
1161 if (p)
1162 {
1163 fLastp = *p;
1164 }
1165 else
1166 {
1167 G4Exception("G4VTwistSurface::CurrentStatus::CurrentStatus()",
1168 "InvalidCondition", FatalException,
1169 "SetCurrentStatus: p = 0!");
1170 }
1171 if (v)
1172 {
1173 fLastv = *v;
1174 }
1175 else
1176 {
1177 fLastv.set(kInfinity, kInfinity, kInfinity);
1178 }
1179 fDone = true;
1180}
1181
1182//=====================================================================
1183//* CurrentStatus::ResetfDone -----------------------------------------
1184
1185void
1186G4VTwistSurface::CurrentStatus::ResetfDone(EValidate validate,
1187 const G4ThreeVector *p,
1188 const G4ThreeVector *v)
1189
1190{
1191 if (validate == fLastValidate && p && *p == fLastp)
1192 {
1193 if (!v || (v && *v == fLastv)) return;
1194 }
1195 G4ThreeVector xx(kInfinity, kInfinity, kInfinity);
1196 for (size_t i=0; i<G4VSURFACENXX; i++)
1197 {
1198 fDistance[i] = kInfinity;
1199 fAreacode[i] = sOutside;
1200 fIsValid[i] = false;
1201 fXX[i] = xx; // bug in old code ( was fXX[i] = xx[i] )
1202 }
1203 fNXX = 0;
1204 fLastp.set(kInfinity, kInfinity, kInfinity);
1205 fLastv.set(kInfinity, kInfinity, kInfinity);
1206 fLastValidate = kUninitialized;
1207 fDone = false;
1208}
1209
1210//=====================================================================
1211//* CurrentStatus::DebugPrint -----------------------------------------
1212
1213void
1214G4VTwistSurface::CurrentStatus::DebugPrint() const
1215{
1216 G4cout << "CurrentStatus::Dist0,1= " << fDistance[0]
1217 << " " << fDistance[1] << " areacode = " << fAreacode[0]
1218 << " " << fAreacode[1] << G4endl;
1219}
1220
1221//=====================================================================
1222// G4VTwistSurface::Boundary class
1223//=====================================================================
1224
1225//=====================================================================
1226//* Boundary::Boundary ------------------------------------------------
1227
1228G4VTwistSurface::Boundary::Boundary()
1229 : fBoundaryAcode(-1), fBoundaryType(0)
1230{
1231}
1232
1233//=====================================================================
1234//* Boundary::~Boundary -----------------------------------------------
1235
1236G4VTwistSurface::Boundary::~Boundary()
1237{
1238}
1239
1240//=====================================================================
1241//* Boundary::SetFields -----------------------------------------------
1242
1243void
1244G4VTwistSurface::Boundary::SetFields(const G4int &areacode,
1245 const G4ThreeVector &d,
1246 const G4ThreeVector &x0,
1247 const G4int &boundarytype)
1248{
1249 fBoundaryAcode = areacode;
1250 fBoundaryDirection = d;
1251 fBoundaryX0 = x0;
1252 fBoundaryType = boundarytype;
1253}
1254
1255//=====================================================================
1256//* Boundary::IsEmpty -------------------------------------------------
1257
1258G4bool G4VTwistSurface::Boundary::IsEmpty() const
1259{
1260 if (fBoundaryAcode == -1) return true;
1261 return false;
1262}
1263
1264//=====================================================================
1265//* Boundary::GetBoundaryParameters -----------------------------------
1266
1267G4bool
1268G4VTwistSurface::Boundary::GetBoundaryParameters(const G4int &areacode,
1269 G4ThreeVector &d,
1270 G4ThreeVector &x0,
1271 G4int &boundarytype) const
1272{
1273 // areacode must be one of them:
1274 // sAxis0 & sAxisMin, sAxis0 & sAxisMax,
1275 // sAxis1 & sAxisMin, sAxis1 & sAxisMax
1276 if ((areacode & sAxis0) && (areacode & sAxis1))
1277 {
1278 G4cerr << "ERROR - G4VTwistSurface::Boundary::GetBoundaryParameters()"
1279 << G4endl
1280 << " Located in the corner area. This function"
1281 << " returns a direction vector of a boundary line."
1282 << " areacode = " << areacode << G4endl;
1283 G4Exception("G4VTwistSurface::Boundary::GetBoundaryParameters()",
1284 "InvalidCondition", FatalException,
1285 "Located in the corner area.");
1286 }
1287 if ((areacode & sSizeMask) != (fBoundaryAcode & sSizeMask))
1288 {
1289 return false;
1290 }
1291 d = fBoundaryDirection;
1292 x0 = fBoundaryX0;
1293 boundarytype = fBoundaryType;
1294 return true;
1295}
Note: See TracBrowser for help on using the repository browser.