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

Last change on this file since 1202 was 1058, checked in by garnier, 15 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.