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