source: trunk/source/geometry/solids/specific/src/G4TwistTubsHypeSide.cc @ 1337

Last change on this file since 1337 was 1337, checked in by garnier, 14 years ago

tag geant4.9.4 beta 1 + modifs locales
File size: 33.0 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: G4TwistTubsHypeSide.cc,v 1.6 2007/05/18 07:39:56 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-01 $
29//
30//
31// --------------------------------------------------------------------
32// GEANT 4 class source file
33//
34//
35// G4TwistTubsHypeSide.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 "G4TwistTubsHypeSide.hh"
46#include "G4GeometryTolerance.hh"
47
48//=====================================================================
49//* constructors ------------------------------------------------------
50
51G4TwistTubsHypeSide::G4TwistTubsHypeSide(const G4String         &name,
52                                         const G4RotationMatrix &rot,
53                                         const G4ThreeVector    &tlate,
54                                         const G4int             handedness,
55                                         const G4double          kappa,
56                                         const G4double          tanstereo,
57                                         const G4double          r0,
58                                         const EAxis             axis0,
59                                         const EAxis             axis1,
60                                               G4double          axis0min,
61                                               G4double          axis1min,
62                                               G4double          axis0max,
63                                               G4double          axis1max )
64  : G4VTwistSurface(name, rot, tlate, handedness, axis0, axis1,
65               axis0min, axis1min, axis0max, axis1max),
66    fKappa(kappa), fTanStereo(tanstereo),
67    fTan2Stereo(tanstereo*tanstereo), fR0(r0), fR02(r0*r0)
68{
69   if (axis0 == kZAxis && axis1 == kPhi) {
70      G4Exception("G4TwistTubsHypeSide::G4TwistTubsHypeSide()", "InvalidSetup",
71                  FatalException, "Should swap axis0 and axis1!");
72   }
73   
74   fInside.gp.set(kInfinity, kInfinity, kInfinity);
75   fInside.inside = kOutside;
76   fIsValidNorm = false;
77   
78   SetCorners();
79   SetBoundaries();
80
81}
82
83G4TwistTubsHypeSide::G4TwistTubsHypeSide(const G4String      &name,
84                                         G4double         EndInnerRadius[2],
85                                         G4double         EndOuterRadius[2],
86                                         G4double         DPhi,
87                                         G4double         EndPhi[2],
88                                         G4double         EndZ[2], 
89                                         G4double         InnerRadius,
90                                         G4double         OuterRadius,
91                                         G4double         Kappa,
92                                         G4double         TanInnerStereo,
93                                         G4double         TanOuterStereo,
94                                         G4int            handedness)
95   : G4VTwistSurface(name)
96{
97
98   fHandedness = handedness;   // +z = +ve, -z = -ve
99   fAxis[0]    = kPhi;
100   fAxis[1]    = kZAxis;
101   fAxisMin[0] = kInfinity;         // we cannot fix boundary min of Phi,
102   fAxisMax[0] = kInfinity;         // because it depends on z.
103   fAxisMin[1] = EndZ[0];
104   fAxisMax[1] = EndZ[1];
105   fKappa      = Kappa;
106   fDPhi       = DPhi ;
107
108   if (handedness < 0) { // inner hyperbolic surface
109      fTanStereo  = TanInnerStereo;
110      fR0         = InnerRadius;
111   } else {              // outer hyperbolic surface
112      fTanStereo  = TanOuterStereo;
113      fR0         = OuterRadius;
114   }
115   fTan2Stereo = fTanStereo * fTanStereo;
116   fR02        = fR0 * fR0;
117   
118   fTrans.set(0, 0, 0);
119   fIsValidNorm = false;
120
121   fInside.gp.set(kInfinity, kInfinity, kInfinity);
122   fInside.inside = kOutside;
123   
124   SetCorners(EndInnerRadius, EndOuterRadius, DPhi, EndPhi, EndZ) ; 
125
126   SetBoundaries();
127}
128
129//=====================================================================
130//* Fake default constructor ------------------------------------------
131
132G4TwistTubsHypeSide::G4TwistTubsHypeSide( __void__& a )
133  : G4VTwistSurface(a)
134{
135}
136
137//=====================================================================
138//* destructor --------------------------------------------------------
139
140G4TwistTubsHypeSide::~G4TwistTubsHypeSide()
141{
142}
143
144//=====================================================================
145//* GetNormal ---------------------------------------------------------
146
147G4ThreeVector G4TwistTubsHypeSide::GetNormal(const G4ThreeVector &tmpxx, 
148                                                   G4bool isGlobal) 
149{
150   // GetNormal returns a normal vector at a surface (or very close
151   // to surface) point at tmpxx.
152   // If isGlobal=true, it returns the normal in global coordinate.
153   //
154   
155   G4ThreeVector xx;
156   if (isGlobal) {
157      xx = ComputeLocalPoint(tmpxx);
158      if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance) {
159         return ComputeGlobalDirection(fCurrentNormal.normal);
160      }
161   } else {
162      xx = tmpxx;
163      if (xx == fCurrentNormal.p) {
164         return fCurrentNormal.normal;
165      }
166   }
167   
168   fCurrentNormal.p = xx;
169
170   G4ThreeVector normal( xx.x(), xx.y(), -xx.z() * fTan2Stereo);
171   normal *= fHandedness;
172   normal = normal.unit();
173
174   if (isGlobal) {
175      fCurrentNormal.normal = ComputeLocalDirection(normal);
176   } else {
177      fCurrentNormal.normal = normal;
178   }
179   return fCurrentNormal.normal;
180}
181
182//=====================================================================
183//* Inside() ----------------------------------------------------------
184
185EInside G4TwistTubsHypeSide::Inside(const G4ThreeVector &gp) 
186{
187   // Inside returns
188   static const G4double halftol
189     = 0.5 * G4GeometryTolerance::GetInstance()->GetRadialTolerance();
190
191   if (fInside.gp == gp) {
192      return fInside.inside;
193   }
194   fInside.gp = gp;
195   
196   G4ThreeVector p = ComputeLocalPoint(gp);
197 
198
199   if (p.mag() < DBL_MIN) {
200      fInside.inside = kOutside;
201      return fInside.inside;
202   }
203   
204   G4double rhohype = GetRhoAtPZ(p);
205   G4double distanceToOut = fHandedness * (rhohype - p.getRho());
206                            // +ve : inside
207
208   if (distanceToOut < -halftol) {
209
210     fInside.inside = kOutside;
211
212   } else {
213
214      G4int areacode = GetAreaCode(p);
215      if (IsOutside(areacode)) {
216         fInside.inside = kOutside;
217      } else if (IsBoundary(areacode)) {
218         fInside.inside = kSurface;
219      } else if (IsInside(areacode)) {
220         if (distanceToOut <= halftol) {
221            fInside.inside = kSurface;
222         } else {
223            fInside.inside = kInside;
224         }
225      } else {
226         G4cout << "WARNING - G4TwistTubsHypeSide::Inside()" << G4endl
227                << "          Invalid option !" << G4endl
228                << "          name, areacode, distanceToOut = "
229                << GetName() << ", " << std::hex << areacode << std::dec << ", "
230                << distanceToOut << G4endl;
231      }
232   }
233   
234   return fInside.inside; 
235}
236
237//=====================================================================
238//* DistanceToSurface -------------------------------------------------
239
240G4int G4TwistTubsHypeSide::DistanceToSurface(const G4ThreeVector &gp,
241                                             const G4ThreeVector &gv,
242                                                   G4ThreeVector  gxx[],
243                                                   G4double       distance[],
244                                                   G4int          areacode[],
245                                                   G4bool         isvalid[],
246                                                   EValidate      validate)
247{
248   //
249   // Decide if and where a line intersects with a hyperbolic
250   // surface (of infinite extent)
251   //
252   // Arguments:
253   //     p       - (in) Point on trajectory
254   //     v       - (in) Vector along trajectory
255   //     r2      - (in) Square of radius at z = 0
256   //     tan2phi - (in) std::tan(stereo)**2
257   //     s       - (out) Up to two points of intersection, where the
258   //                     intersection point is p + s*v, and if there are
259   //                     two intersections, s[0] < s[1]. May be negative.
260   // Returns:
261   //     The number of intersections. If 0, the trajectory misses.
262   //
263   //
264   // Equation of a line:
265   //
266   //       x = x0 + s*tx      y = y0 + s*ty      z = z0 + s*tz
267   //
268   // Equation of a hyperbolic surface:
269   //
270   //       x**2 + y**2 = r**2 + (z*tanPhi)**2
271   //
272   // Solution is quadratic:
273   //
274   //  a*s**2 + b*s + c = 0
275   //
276   // where:
277   //
278   //  a = tx**2 + ty**2 - (tz*tanPhi)**2
279   //
280   //  b = 2*( x0*tx + y0*ty - z0*tz*tanPhi**2 )
281   //
282   //  c = x0**2 + y0**2 - r**2 - (z0*tanPhi)**2
283   //
284     
285   fCurStatWithV.ResetfDone(validate, &gp, &gv);
286
287   if (fCurStatWithV.IsDone()) {
288      G4int i;
289      for (i=0; i<fCurStatWithV.GetNXX(); i++) {
290         gxx[i] = fCurStatWithV.GetXX(i);
291         distance[i] = fCurStatWithV.GetDistance(i);
292         areacode[i] = fCurStatWithV.GetAreacode(i);
293         isvalid[i]  = fCurStatWithV.IsValid(i);
294      }
295      return fCurStatWithV.GetNXX();
296   } else {
297      // initialize
298      G4int i;
299      for (i=0; i<2; i++) {
300         distance[i] = kInfinity;
301         areacode[i] = sOutside;
302         isvalid[i]  = false;
303         gxx[i].set(kInfinity, kInfinity, kInfinity);
304      }
305   }
306   
307   G4ThreeVector p = ComputeLocalPoint(gp);
308   G4ThreeVector v = ComputeLocalDirection(gv);
309   G4ThreeVector xx[2]; 
310   
311   //
312   // special case!  p is on origin.
313   //
314
315   if (p.mag() == 0) {
316      // p is origin.
317      // unique solution of 2-dimension question in r-z plane
318      // Equations:
319      //    r^2 = fR02 + z^2*fTan2Stere0
320      //    r = beta*z
321      //        where
322      //        beta = vrho / vz
323      // Solution (z value of intersection point):
324      //    xxz = +- std::sqrt (fR02 / (beta^2 - fTan2Stereo))
325      //
326
327      G4double vz    = v.z();
328      G4double absvz = std::abs(vz);
329      G4double vrho  = v.getRho();       
330      G4double vslope = vrho/vz;
331      G4double vslope2 = vslope * vslope;
332      if (vrho == 0 || (vrho/absvz) <= (absvz*std::fabs(fTanStereo)/absvz)) {
333         // vz/vrho is bigger than slope of asymptonic line
334         distance[0] = kInfinity;
335         fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
336                                        isvalid[0], 0, validate, &gp, &gv);
337         return 0;
338      }
339       
340      if (vz) { 
341         G4double xxz  = std::sqrt(fR02 / (vslope2 - fTan2Stereo)) 
342                        * (vz / std::fabs(vz)) ;
343         G4double t = xxz / vz;
344         xx[0].set(t*v.x(), t*v.y(), xxz);
345      } else {
346         // p.z = 0 && v.z =0
347         xx[0].set(v.x()*fR0, v.y()*fR0, 0);  // v is a unit vector.
348      }
349      distance[0] = xx[0].mag();
350      gxx[0]      = ComputeGlobalPoint(xx[0]);
351
352      if (validate == kValidateWithTol) {
353         areacode[0] = GetAreaCode(xx[0]);
354         if (!IsOutside(areacode[0])) {
355            if (distance[0] >= 0) isvalid[0] = true;
356         }
357      } else if (validate == kValidateWithoutTol) {
358         areacode[0] = GetAreaCode(xx[0], false);
359         if (IsInside(areacode[0])) {
360            if (distance[0] >= 0) isvalid[0] = true;
361         }
362      } else { // kDontValidate                       
363         areacode[0] = sInside;
364            if (distance[0] >= 0) isvalid[0] = true;
365      }
366                 
367      fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
368                                        isvalid[0], 1, validate, &gp, &gv);
369      return 1;
370   }
371
372   //
373   // special case end.
374   //
375
376   G4double a = v.x()*v.x() + v.y()*v.y() - v.z()*v.z()*fTan2Stereo;
377   G4double b = 2.0 * ( p.x() * v.x() + p.y() * v.y() - p.z() * v.z() * fTan2Stereo );
378   G4double c = p.x()*p.x() + p.y()*p.y() - fR02 - p.z()*p.z()*fTan2Stereo;
379   G4double D = b*b - 4*a*c;          //discriminant
380   
381   if (std::fabs(a) < DBL_MIN) {
382      if (std::fabs(b) > DBL_MIN) {           // single solution
383
384         distance[0] = -c/b;
385         xx[0] = p + distance[0]*v;
386         gxx[0] = ComputeGlobalPoint(xx[0]);
387
388         if (validate == kValidateWithTol) {
389            areacode[0] = GetAreaCode(xx[0]);
390            if (!IsOutside(areacode[0])) {
391               if (distance[0] >= 0) isvalid[0] = true;
392            }
393         } else if (validate == kValidateWithoutTol) {
394            areacode[0] = GetAreaCode(xx[0], false);
395            if (IsInside(areacode[0])) {
396               if (distance[0] >= 0) isvalid[0] = true;
397            }
398         } else { // kDontValidate                       
399            areacode[0] = sInside;
400               if (distance[0] >= 0) isvalid[0] = true;
401         }
402                 
403         fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
404                                        isvalid[0], 1, validate, &gp, &gv);
405         return 1;
406         
407      } else {
408         // if a=b=0 and c != 0, p is origin and v is parallel to asymptotic line.
409         // if a=b=c=0, p is on surface and v is paralell to stereo wire.
410         // return distance = infinity.
411
412         fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
413                                        isvalid[0], 0, validate, &gp, &gv);
414
415         return 0;
416      }
417     
418   } else if (D > DBL_MIN) {         // double solutions
419     
420      D = std::sqrt(D);
421      G4double      factor = 0.5/a;
422      G4double      tmpdist[2] = {kInfinity, kInfinity};
423      G4ThreeVector tmpxx[2] ;
424      G4int         tmpareacode[2] = {sOutside, sOutside};
425      G4bool        tmpisvalid[2]  = {false, false};
426      G4int i;
427
428      for (i=0; i<2; i++) {
429         tmpdist[i] = factor*(-b - D);
430         D = -D;
431         tmpxx[i] = p + tmpdist[i]*v;
432       
433         if (validate == kValidateWithTol) {
434            tmpareacode[i] = GetAreaCode(tmpxx[i]);
435            if (!IsOutside(tmpareacode[i])) {
436               if (tmpdist[i] >= 0) tmpisvalid[i] = true;
437               continue;
438            }
439         } else if (validate == kValidateWithoutTol) {
440            tmpareacode[i] = GetAreaCode(tmpxx[i], false);
441            if (IsInside(tmpareacode[i])) {
442               if (tmpdist[i] >= 0) tmpisvalid[i] = true;
443               continue;
444            }
445         } else { // kDontValidate
446            tmpareacode[i] = sInside;
447               if (tmpdist[i] >= 0) tmpisvalid[i] = true;
448            continue;
449         }
450      }     
451
452      if (tmpdist[0] <= tmpdist[1]) {
453          distance[0] = tmpdist[0];
454          distance[1] = tmpdist[1];
455          xx[0]       = tmpxx[0];
456          xx[1]       = tmpxx[1];
457          gxx[0]      = ComputeGlobalPoint(tmpxx[0]);
458          gxx[1]      = ComputeGlobalPoint(tmpxx[1]);
459          areacode[0] = tmpareacode[0];
460          areacode[1] = tmpareacode[1];
461          isvalid[0]  = tmpisvalid[0];
462          isvalid[1]  = tmpisvalid[1];
463      } else {
464          distance[0] = tmpdist[1];
465          distance[1] = tmpdist[0];
466          xx[0]       = tmpxx[1];
467          xx[1]       = tmpxx[0];
468          gxx[0]      = ComputeGlobalPoint(tmpxx[1]);
469          gxx[1]      = ComputeGlobalPoint(tmpxx[0]);
470          areacode[0] = tmpareacode[1];
471          areacode[1] = tmpareacode[0];
472          isvalid[0]  = tmpisvalid[1];
473          isvalid[1]  = tmpisvalid[0];
474      }
475         
476      fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
477                                     isvalid[0], 2, validate, &gp, &gv);
478      fCurStatWithV.SetCurrentStatus(1, gxx[1], distance[1], areacode[1],
479                                     isvalid[1], 2, validate, &gp, &gv);
480      return 2;
481     
482   } else {
483      // if D<0, no solution
484      // if D=0, just grazing the surfaces, return kInfinity
485
486      fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
487                                     isvalid[0], 0, validate, &gp, &gv);
488      return 0;
489   }
490   G4Exception("G4TwistTubsHypeSide::DistanceToSurface(p,v)",
491               "InvalidCondition", FatalException, "Illegal operation !");
492   return 1;
493}
494
495   
496//=====================================================================
497//* DistanceToSurface -------------------------------------------------
498
499G4int G4TwistTubsHypeSide::DistanceToSurface(const G4ThreeVector &gp,
500                                                   G4ThreeVector  gxx[],
501                                                   G4double       distance[],
502                                                   G4int          areacode[])
503{
504    // Find the approximate distance of a point of a hyperbolic surface.
505    // The distance must be an underestimate.
506    // It will also be nice (although not necessary) that the estimate is
507    // always finite no matter how close the point is.
508    //
509    // We arranged G4Hype::ApproxDistOutside and G4Hype::ApproxDistInside
510    // for this function. See these discriptions.
511   
512   static const G4double halftol
513     = 0.5 * G4GeometryTolerance::GetInstance()->GetRadialTolerance();
514
515   fCurStat.ResetfDone(kDontValidate, &gp);
516
517   if (fCurStat.IsDone()) {
518      for (G4int i=0; i<fCurStat.GetNXX(); i++) {
519         gxx[i] = fCurStat.GetXX(i);
520         distance[i] = fCurStat.GetDistance(i);
521         areacode[i] = fCurStat.GetAreacode(i);
522      }
523      return fCurStat.GetNXX();
524   } else {
525      // initialize
526      for (G4int i=0; i<2; i++) {
527         distance[i] = kInfinity;
528         areacode[i] = sOutside;
529         gxx[i].set(kInfinity, kInfinity, kInfinity);
530      }
531   }
532   
533
534   G4ThreeVector p = ComputeLocalPoint(gp);
535   G4ThreeVector xx;
536
537   //
538   // special case!
539   // If p is on surface, return distance = 0 immediatery .
540   //
541   G4ThreeVector  lastgxx[2];
542   G4double       distfromlast[2];
543   for (G4int i=0; i<2; i++) {
544      lastgxx[i] = fCurStatWithV.GetXX(i);
545      distfromlast[i] = (gp - lastgxx[i]).mag();
546   }
547
548   if ((gp - lastgxx[0]).mag() < halftol || (gp - lastgxx[1]).mag() < halftol) {
549      // last winner, or last poststep point is on the surface.
550      xx = p;             
551      gxx[0] = gp;
552      distance[0] = 0;     
553
554      G4bool isvalid = true;
555      fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
556                                isvalid, 1, kDontValidate, &gp);
557
558      return 1;
559
560   }
561   //
562   // special case end
563   //
564       
565   G4double prho       = p.getRho();
566   G4double pz         = std::fabs(p.z());           // use symmetry
567   G4double r1         = std::sqrt(fR02 + pz * pz * fTan2Stereo);
568   
569   G4ThreeVector pabsz(p.x(), p.y(), pz);
570   
571   if (prho > r1 + halftol) {  // p is outside of Hyperbolic surface
572
573      // First point xx1
574      G4double t = r1 / prho;
575      G4ThreeVector xx1(t * pabsz.x(), t * pabsz.y() , pz);
576     
577      // Second point xx2
578      G4double z2 = (prho * fTanStereo + pz) / (1 + fTan2Stereo);
579      G4double r2 = std::sqrt(fR02 + z2 * z2 * fTan2Stereo);
580      t = r2 / prho;
581      G4ThreeVector xx2(t * pabsz.x(), t * pabsz.y() , z2);
582           
583      G4double len = (xx2 - xx1).mag();
584      if (len < DBL_MIN) {
585         // xx2 = xx1?? I guess we
586         // must have really bracketed the normal
587         distance[0] = (pabsz - xx1).mag();
588         xx = xx1;
589      } else {
590         distance[0] = DistanceToLine(pabsz, xx1, (xx2 - xx1) , xx);
591      }
592     
593   } else if (prho < r1 - halftol) { // p is inside of Hyperbolic surface.
594           
595      // First point xx1
596      G4double t;
597      G4ThreeVector xx1;
598      if (prho < DBL_MIN) {
599         xx1.set(r1, 0. , pz);
600      } else {
601         t = r1 / prho;
602         xx1.set(t * pabsz.x(), t * pabsz.y() , pz);
603      }
604     
605      // dr, dz is tangential vector of Hyparbolic surface at xx1
606      // dr = r, dz = z*tan2stereo
607      G4double dr  = pz * fTan2Stereo;
608      G4double dz  = r1;
609      G4double tanbeta   = dr / dz;
610      G4double pztanbeta = pz * tanbeta;
611     
612      // Second point xx2
613      // xx2 is intersection between x-axis and tangential vector
614      G4double r2 = r1 - pztanbeta;
615      G4ThreeVector xx2;
616      if (prho < DBL_MIN) {
617         xx2.set(r2, 0. , 0.);
618      } else {
619         t  = r2 / prho;
620         xx2.set(t * pabsz.x(), t * pabsz.y() , 0.);
621      }
622     
623      G4ThreeVector d = xx2 - xx1;
624      distance[0] = DistanceToLine(pabsz, xx1, d, xx);
625         
626   } else {  // p is on Hyperbolic surface.
627   
628      distance[0] = 0;
629      xx.set(p.x(), p.y(), pz);
630
631   }
632
633   if (p.z() < 0) {
634      G4ThreeVector tmpxx(xx.x(), xx.y(), -xx.z());
635      xx = tmpxx;
636   }
637       
638   gxx[0] = ComputeGlobalPoint(xx);
639   areacode[0]    = sInside;
640   G4bool isvalid = true;
641   fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
642                             isvalid, 1, kDontValidate, &gp);
643   return 1;
644}
645
646//=====================================================================
647//* GetAreaCode -------------------------------------------------------
648
649G4int G4TwistTubsHypeSide::GetAreaCode(const G4ThreeVector &xx, 
650                                             G4bool         withTol)
651{
652   static const G4double ctol = 0.5 * kCarTolerance;
653   G4int areacode = sInside;
654
655   if ((fAxis[0] == kPhi && fAxis[1] == kZAxis))  {
656      //G4int phiaxis = 0;
657      G4int zaxis   = 1;
658     
659      if (withTol) {
660
661         G4bool isoutside      = false;
662         G4int  phiareacode    = GetAreaCodeInPhi(xx);
663         G4bool isoutsideinphi = IsOutside(phiareacode);
664
665         // test boundary of phiaxis
666
667         if ((phiareacode & sAxisMin) == sAxisMin) {
668
669            areacode |= (sAxis0 & (sAxisPhi | sAxisMin)) | sBoundary;
670            if (isoutsideinphi) isoutside = true;
671
672         } else if ((phiareacode & sAxisMax)  == sAxisMax) {
673
674            areacode |= (sAxis0 & (sAxisPhi | sAxisMax)) | sBoundary;
675            if (isoutsideinphi) isoutside = true;
676
677         }
678
679         // test boundary of zaxis
680
681         if (xx.z() < fAxisMin[zaxis] + ctol) {
682
683            areacode |= (sAxis1 & (sAxisZ | sAxisMin));
684            if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
685            else                        areacode |= sBoundary;
686
687            if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true;
688
689         } else if (xx.z() > fAxisMax[zaxis] - ctol) {
690
691            areacode |= (sAxis1 & (sAxisZ | sAxisMax));
692            if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
693            else                        areacode |= sBoundary;
694
695            if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true;
696         }
697
698         // if isoutside = true, clear sInside bit.
699         // if not on boundary, add boundary information.
700
701         if (isoutside) {
702            G4int tmpareacode = areacode & (~sInside);
703            areacode = tmpareacode;
704         } else if ((areacode & sBoundary) != sBoundary) {
705            areacode |= (sAxis0 & sAxisPhi) | (sAxis1 & sAxisZ);
706         }
707
708         return areacode;
709     
710      } else {
711
712         G4int phiareacode = GetAreaCodeInPhi(xx, false);
713         
714         // test boundary of z-axis
715
716         if (xx.z() < fAxisMin[zaxis]) {
717
718            areacode |= (sAxis1 & (sAxisZ | sAxisMin)) | sBoundary;
719
720         } else if (xx.z() > fAxisMax[zaxis]) {
721
722            areacode |= (sAxis1 & (sAxisZ | sAxisMax)) | sBoundary;
723
724         }
725
726         // boundary of phi-axis
727
728         if (phiareacode == sAxisMin) {
729
730            areacode |= (sAxis0 & (sAxisPhi | sAxisMin));
731            if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
732            else                        areacode |= sBoundary; 
733             
734         } else if (phiareacode == sAxisMax) {
735
736            areacode |= (sAxis0 & (sAxisPhi | sAxisMax));
737            if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
738            else                        areacode |= sBoundary; 
739           
740         }
741
742         // if not on boundary, add boundary information.
743
744         if ((areacode & sBoundary) != sBoundary) {
745            areacode |= (sAxis0 & sAxisPhi) | (sAxis1 & sAxisZ);
746         }
747         return areacode;
748      }
749   } else {
750      G4cerr << "ERROR - G4TwistTubsHypeSide::GetAreaCode()" << G4endl
751             << "        fAxis[0] = " << fAxis[0] << G4endl
752             << "        fAxis[1] = " << fAxis[1] << G4endl;
753      G4Exception("G4TwistTubsHypeSide::GetAreaCode()",
754                  "NotImplemented", FatalException,
755                  "Feature NOT implemented !");
756   }
757   return areacode;
758}
759
760//=====================================================================
761//* GetAreaCodeInPhi --------------------------------------------------
762
763G4int G4TwistTubsHypeSide::GetAreaCodeInPhi(const G4ThreeVector &xx,
764                                                  G4bool withTol)
765{
766   
767   G4ThreeVector lowerlimit; // lower phi-boundary limit at z = xx.z()
768   G4ThreeVector upperlimit; // upper phi-boundary limit at z = xx.z()
769   lowerlimit = GetBoundaryAtPZ(sAxis0 & sAxisMin, xx);
770   upperlimit = GetBoundaryAtPZ(sAxis0 & sAxisMax, xx);
771
772   G4int  areacode  = sInside;
773   G4bool isoutside = false; 
774   
775   if (withTol) {
776         
777      if (AmIOnLeftSide(xx, lowerlimit) >= 0) {        // xx is on lowerlimit
778         areacode |= (sAxisMin | sBoundary);
779         if (AmIOnLeftSide(xx, lowerlimit) > 0) isoutside = true; 
780
781      } else if (AmIOnLeftSide(xx, upperlimit) <= 0) { // xx is on upperlimit
782         areacode |= (sAxisMax | sBoundary);
783         if (AmIOnLeftSide(xx, upperlimit) < 0) isoutside = true; 
784      }
785
786      // if isoutside = true, clear inside bit.
787
788      if (isoutside) {
789         G4int tmpareacode = areacode & (~sInside);
790         areacode = tmpareacode;
791      }
792
793
794   } else {
795   
796      if (AmIOnLeftSide(xx, lowerlimit, false) >= 0) {
797         areacode |= (sAxisMin | sBoundary);
798      } else if (AmIOnLeftSide(xx, upperlimit, false) <= 0) {
799         areacode |= (sAxisMax | sBoundary);
800      }
801   }
802
803   return areacode;
804   
805}
806
807//=====================================================================
808//* SetCorners(EndInnerRadius, EndOuterRadius,DPhi,EndPhi,EndZ) -------
809
810void G4TwistTubsHypeSide::SetCorners(
811                                     G4double         EndInnerRadius[2],
812                                     G4double         EndOuterRadius[2],
813                                     G4double         DPhi,
814                                     G4double         endPhi[2],
815                                     G4double         endZ[2] 
816                                     )
817{
818   // Set Corner points in local coodinate.
819
820   if (fAxis[0] == kPhi && fAxis[1] == kZAxis) {
821
822      G4int i;
823      G4double endRad[2];
824      G4double halfdphi = 0.5*DPhi;
825     
826      for (i=0; i<2; i++) { // i=0,1 : -ve z, +ve z
827         endRad[i] = (fHandedness == 1 ? EndOuterRadius[i]
828                                      : EndInnerRadius[i]);
829      }
830
831      G4int zmin = 0 ;  // at -ve z
832      G4int zmax = 1 ;  // at +ve z
833
834      G4double x, y, z;
835     
836      // corner of Axis0min and Axis1min
837      x = endRad[zmin]*std::cos(endPhi[zmin] - halfdphi);
838      y = endRad[zmin]*std::sin(endPhi[zmin] - halfdphi);
839      z = endZ[zmin];
840      SetCorner(sC0Min1Min, x, y, z);
841     
842      // corner of Axis0max and Axis1min
843      x = endRad[zmin]*std::cos(endPhi[zmin] + halfdphi);
844      y = endRad[zmin]*std::sin(endPhi[zmin] + halfdphi);
845      z = endZ[zmin];
846      SetCorner(sC0Max1Min, x, y, z);
847     
848      // corner of Axis0max and Axis1max
849      x = endRad[zmax]*std::cos(endPhi[zmax] + halfdphi);
850      y = endRad[zmax]*std::sin(endPhi[zmax] + halfdphi);
851      z = endZ[zmax];
852      SetCorner(sC0Max1Max, x, y, z);
853     
854      // corner of Axis0min and Axis1max
855      x = endRad[zmax]*std::cos(endPhi[zmax] - halfdphi);
856      y = endRad[zmax]*std::sin(endPhi[zmax] - halfdphi);
857      z = endZ[zmax];
858      SetCorner(sC0Min1Max, x, y, z);
859
860   } else {
861      G4cerr << "ERROR - G4TwistTubsFlatSide::SetCorners()" << G4endl
862             << "        fAxis[0] = " << fAxis[0] << G4endl
863             << "        fAxis[1] = " << fAxis[1] << G4endl;
864      G4Exception("G4TwistTubsHypeSide::SetCorners()",
865                  "NotImplemented", FatalException,
866                  "Feature NOT implemented !");
867   }
868}
869
870
871//=====================================================================
872//* SetCorners() ------------------------------------------------------
873
874void G4TwistTubsHypeSide::SetCorners()
875{
876   G4Exception("G4TwistTubsHypeSide::SetCorners()",
877               "NotImplemented", FatalException,
878               "Method NOT implemented !");
879}
880
881//=====================================================================
882//* SetBoundaries() ---------------------------------------------------
883
884void G4TwistTubsHypeSide::SetBoundaries()
885{
886   // Set direction-unit vector of phi-boundary-lines in local coodinate.
887   // sAxis0 must be kPhi.
888   // This fanction set lower phi-boundary and upper phi-boundary.
889
890   if (fAxis[0] == kPhi && fAxis[1] == kZAxis) {
891
892      G4ThreeVector direction;
893      // sAxis0 & sAxisMin
894      direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
895      direction = direction.unit();
896      SetBoundary(sAxis0 & (sAxisPhi | sAxisMin), direction, 
897                   GetCorner(sC0Min1Min), sAxisZ);
898
899      // sAxis0 & sAxisMax
900      direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
901      direction = direction.unit();
902      SetBoundary(sAxis0 & (sAxisPhi | sAxisMax), direction, 
903                  GetCorner(sC0Max1Min), sAxisZ);
904
905      // sAxis1 & sAxisMin
906      direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
907      direction = direction.unit();
908      SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction, 
909                   GetCorner(sC0Min1Min), sAxisPhi);
910
911      // sAxis1 & sAxisMax
912      direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
913      direction = direction.unit();
914      SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction, 
915                  GetCorner(sC0Min1Max), sAxisPhi);
916   } else {
917      G4cerr << "ERROR - G4TwistTubsHypeSide::SetBoundaries()" << G4endl
918             << "        fAxis[0] = " << fAxis[0] << G4endl
919             << "        fAxis[1] = " << fAxis[1] << G4endl;
920      G4Exception("G4TwistTubsHypeSide::SetBoundaries()",
921                  "NotImplemented", FatalException,
922                  "Feature NOT implemented !");
923   }
924}
925
926//=====================================================================
927//* GetFacets() -------------------------------------------------------
928
929void G4TwistTubsHypeSide::GetFacets( G4int m, G4int n, G4double xyz[][3],
930                                     G4int faces[][4], G4int iside ) 
931{
932
933  G4double z ;     // the two parameters for the surface equation
934  G4double x,xmin,xmax ;
935
936  G4ThreeVector p ;  // a point on the surface, given by (z,u)
937
938  G4int nnode ;
939  G4int nface ;
940
941  // calculate the (n-1)*(m-1) vertices
942
943  G4int i,j ;
944
945  for ( i = 0 ; i<n ; i++ ) {
946
947    z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin[1])/(n-1) ;
948
949    for ( j = 0 ; j<m ; j++ )
950    {
951      nnode = GetNode(i,j,m,n,iside) ;
952
953      xmin = GetBoundaryMin(z) ; 
954      xmax = GetBoundaryMax(z) ;
955
956      if (fHandedness < 0) { // inner hyperbolic surface
957        x = xmin + j*(xmax-xmin)/(m-1) ;
958      } else {               // outer hyperbolic surface
959        x = xmax - j*(xmax-xmin)/(m-1) ;
960      }
961
962      p = SurfacePoint(x,z,true) ;  // surface point in global coord.system
963
964      xyz[nnode][0] = p.x() ;
965      xyz[nnode][1] = p.y() ;
966      xyz[nnode][2] = p.z() ;
967
968      if ( i<n-1 && j<m-1 ) {   // clock wise filling
969       
970        nface = GetFace(i,j,m,n,iside) ;
971
972        faces[nface][0] = GetEdgeVisibility(i,j,m,n,0,1) * ( GetNode(i  ,j  ,m,n,iside)+1) ; 
973        faces[nface][1] = GetEdgeVisibility(i,j,m,n,1,1) * ( GetNode(i+1,j  ,m,n,iside)+1) ;
974        faces[nface][2] = GetEdgeVisibility(i,j,m,n,2,1) * ( GetNode(i+1,j+1,m,n,iside)+1) ;
975        faces[nface][3] = GetEdgeVisibility(i,j,m,n,3,1) * ( GetNode(i  ,j+1,m,n,iside)+1) ;
976
977      }
978    }
979  }
980}
Note: See TracBrowser for help on using the repository browser.