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