Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/specific/src/G4TwistTubsFlatSide.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /geometry/solids/specific/src/G4TwistTubsFlatSide.cc (Version 11.3.0) and /geometry/solids/specific/src/G4TwistTubsFlatSide.cc (Version 10.4.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // G4TwistTubsFlatSide implementation          << 
 27 //                                                 26 //
 28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepbu <<  27 // $Id: G4TwistTubsFlatSide.cc 72937 2013-08-14 13:20:38Z gcosmo $
 29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), <<  28 //
 30 //               from original version in Jupi <<  29 // 
                                                   >>  30 // --------------------------------------------------------------------
                                                   >>  31 // GEANT 4 class source file
                                                   >>  32 //
                                                   >>  33 //
                                                   >>  34 // G4TwistTubsFlatSide.cc
                                                   >>  35 //
                                                   >>  36 // Author: 
                                                   >>  37 //   01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp)
                                                   >>  38 //
                                                   >>  39 // History:
                                                   >>  40 //   13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
                                                   >>  41 //                 from original version in Jupiter-2.5.02 application.
 31 // -------------------------------------------     42 // --------------------------------------------------------------------
 32                                                    43 
 33 #include "G4TwistTubsFlatSide.hh"                  44 #include "G4TwistTubsFlatSide.hh"
 34 #include "G4GeometryTolerance.hh"                  45 #include "G4GeometryTolerance.hh"
 35                                                    46 
 36 //============================================     47 //=====================================================================
 37 //* constructors -----------------------------     48 //* constructors ------------------------------------------------------
 38                                                    49 
 39 G4TwistTubsFlatSide::G4TwistTubsFlatSide(const <<  50 G4TwistTubsFlatSide::G4TwistTubsFlatSide(const G4String         &name,
 40                              const G4RotationM <<  51                              const G4RotationMatrix &rot,
 41                              const G4ThreeVect <<  52                              const G4ThreeVector    &tlate,
 42                              const G4ThreeVect <<  53                              const G4ThreeVector    &n,
 43                              const EAxis           54                              const EAxis             axis0 ,
 44                              const EAxis           55                              const EAxis             axis1 ,
 45                                    G4double        56                                    G4double          axis0min,
 46                                    G4double        57                                    G4double          axis1min,
 47                                    G4double        58                                    G4double          axis0max,
 48                                    G4double        59                                    G4double          axis1max )
 49   : G4VTwistSurface(name, rot, tlate, 0, axis0     60   : G4VTwistSurface(name, rot, tlate, 0, axis0, axis1,
 50                     axis0min, axis1min, axis0m <<  61                axis0min, axis1min, axis0max, axis1max)
 51 {                                                  62 {   
 52    if (axis0 == kPhi && axis1 == kRho)         <<  63    if (axis0 == kPhi && axis1 == kRho) {
 53    {                                           << 
 54       G4Exception("G4TwistTubsFlatSide::G4Twis     64       G4Exception("G4TwistTubsFlatSide::G4TwistTubsFlatSide()",
 55                   "GeomSolids0002", FatalError     65                   "GeomSolids0002", FatalErrorInArgument,
 56                   "Should swap axis0 and axis1     66                   "Should swap axis0 and axis1!");
 57    }                                               67    }
 58                                                    68    
 59    G4ThreeVector normal = rot.inverse()*n;         69    G4ThreeVector normal = rot.inverse()*n;
 60    fCurrentNormal.normal = normal.unit();   //     70    fCurrentNormal.normal = normal.unit();   // in local coordinate system
 61    fIsValidNorm = true;                            71    fIsValidNorm = true;
 62                                                    72 
 63    SetCorners();                                   73    SetCorners();
 64    SetBoundaries();                                74    SetBoundaries();
 65                                                    75 
 66    fSurfaceArea = 1. ;  // NOTE: not yet imple <<  76    fSurfaceArea = 1 ;  // not yet implemented. This is NOT a problem for tracking
                                                   >>  77 
 67 }                                                  78 }
 68                                                    79 
 69 G4TwistTubsFlatSide::G4TwistTubsFlatSide( cons <<  80 
 70                                                <<  81 
 71                                                <<  82 G4TwistTubsFlatSide::G4TwistTubsFlatSide( const G4String        &name,
 72                                                <<  83                                     G4double         EndInnerRadius[2],
 73                                                <<  84                                     G4double         EndOuterRadius[2],
 74                                                <<  85                                     G4double         DPhi,
 75                                                <<  86                                     G4double         EndPhi[2],
                                                   >>  87                                     G4double         EndZ[2], 
                                                   >>  88                                     G4int            handedness ) 
 76   : G4VTwistSurface(name)                          89   : G4VTwistSurface(name)
 77 {                                                  90 {
 78    fHandedness = handedness;   // +z = +ve, -z     91    fHandedness = handedness;   // +z = +ve, -z = -ve
 79    fAxis[0]    = kRho;         // in local coo     92    fAxis[0]    = kRho;         // in local coordinate system
 80    fAxis[1]    = kPhi;                             93    fAxis[1]    = kPhi;
 81    G4int i     = (handedness < 0 ? 0 : 1);         94    G4int i     = (handedness < 0 ? 0 : 1);
 82    fAxisMin[0] = EndInnerRadius[i];  // Inner-     95    fAxisMin[0] = EndInnerRadius[i];  // Inner-hype radius at z=0
 83    fAxisMax[0] = EndOuterRadius[i];  // Outer-     96    fAxisMax[0] = EndOuterRadius[i];  // Outer-hype radius at z=0
 84    fAxisMin[1] = -0.5*DPhi;                        97    fAxisMin[1] = -0.5*DPhi;
 85    fAxisMax[1] = -fAxisMin[1];                     98    fAxisMax[1] = -fAxisMin[1];
 86    fCurrentNormal.normal.set(0, 0, (fHandednes     99    fCurrentNormal.normal.set(0, 0, (fHandedness < 0 ? -1 : 1)); 
 87          // Unit vector, in local coordinate s    100          // Unit vector, in local coordinate system
 88    fRot.rotateZ(EndPhi[i]);                       101    fRot.rotateZ(EndPhi[i]);
 89    fTrans.set(0, 0, EndZ[i]);                     102    fTrans.set(0, 0, EndZ[i]);
 90    fIsValidNorm = true;                           103    fIsValidNorm = true;
 91                                                   104 
 92    SetCorners();                                  105    SetCorners();
 93    SetBoundaries();                               106    SetBoundaries();
 94                                                   107 
 95    fSurfaceArea =  0.5*DPhi * (EndOuterRadius[    108    fSurfaceArea =  0.5*DPhi * (EndOuterRadius[i]*EndOuterRadius[i]  
 96                              - EndInnerRadius[ << 109                                - EndInnerRadius[i]*EndInnerRadius[i] ) ; 
                                                   >> 110 
 97 }                                                 111 }
 98                                                   112 
                                                   >> 113 
 99 //============================================    114 //=====================================================================
100 //* Fake default constructor -----------------    115 //* Fake default constructor ------------------------------------------
101                                                   116 
102 G4TwistTubsFlatSide::G4TwistTubsFlatSide( __vo    117 G4TwistTubsFlatSide::G4TwistTubsFlatSide( __void__& a )
103   : G4VTwistSurface(a)                         << 118   : G4VTwistSurface(a), fSurfaceArea(0.)
104 {                                                 119 {
105 }                                                 120 }
106                                                   121 
                                                   >> 122 
107 //============================================    123 //=====================================================================
108 //* destructor -------------------------------    124 //* destructor --------------------------------------------------------
109                                                   125 
110 G4TwistTubsFlatSide::~G4TwistTubsFlatSide() =  << 126 G4TwistTubsFlatSide::~G4TwistTubsFlatSide()
                                                   >> 127 {
                                                   >> 128 }
111                                                   129 
112 //============================================    130 //=====================================================================
113 //* GetNormal --------------------------------    131 //* GetNormal ---------------------------------------------------------
114                                                   132 
115 G4ThreeVector G4TwistTubsFlatSide::GetNormal(c << 133 G4ThreeVector G4TwistTubsFlatSide::GetNormal(const G4ThreeVector & /* xx */ , 
116                                                << 134                                              G4bool isGlobal)
117 {                                                 135 {
118    if (isGlobal)                               << 136    if (isGlobal) {
119    {                                           << 
120       return ComputeGlobalDirection(fCurrentNo    137       return ComputeGlobalDirection(fCurrentNormal.normal);
121    }                                           << 138    } else {
122    else                                        << 
123    {                                           << 
124       return fCurrentNormal.normal;               139       return fCurrentNormal.normal;
125    }                                              140    }
126 }                                                 141 }
127                                                   142 
128 //============================================    143 //=====================================================================
129 //* DistanceToSurface(p, v) ------------------    144 //* DistanceToSurface(p, v) -------------------------------------------
130                                                   145 
131 G4int G4TwistTubsFlatSide::DistanceToSurface(c << 146 G4int G4TwistTubsFlatSide::DistanceToSurface(const G4ThreeVector &gp,
132                                              c << 147                                        const G4ThreeVector &gv,
133                                                << 148                                              G4ThreeVector  gxx[],
134                                                << 149                                              G4double       distance[],
135                                                << 150                                              G4int          areacode[],
136                                                << 151                                              G4bool         isvalid[],
137                                                << 152                                              EValidate      validate) 
138 {                                                 153 {
139    fCurStatWithV.ResetfDone(validate, &gp, &gv    154    fCurStatWithV.ResetfDone(validate, &gp, &gv);
140                                                   155    
141    if (fCurStatWithV.IsDone())                 << 156    if (fCurStatWithV.IsDone()) {
142    {                                           << 157       G4int i;
143       for (G4int i=0; i<fCurStatWithV.GetNXX() << 158       for (i=0; i<fCurStatWithV.GetNXX(); i++) {
144       {                                        << 
145          gxx[i] = fCurStatWithV.GetXX(i);         159          gxx[i] = fCurStatWithV.GetXX(i);
146          distance[i] = fCurStatWithV.GetDistan    160          distance[i] = fCurStatWithV.GetDistance(i);
147          areacode[i] = fCurStatWithV.GetAreaco    161          areacode[i] = fCurStatWithV.GetAreacode(i);
148          isvalid[i]  = fCurStatWithV.IsValid(i    162          isvalid[i]  = fCurStatWithV.IsValid(i);
149       }                                           163       }
150       return fCurStatWithV.GetNXX();              164       return fCurStatWithV.GetNXX();
151    }                                           << 165    } else {
152    else  // initialize                         << 166       // initialize
153    {                                           << 167       G4int i;
154       for (G4int i=0; i<2; ++i)                << 168       for (i=0; i<2; i++) {
155       {                                        << 
156          distance[i] = kInfinity;                 169          distance[i] = kInfinity;
157          areacode[i] = sOutside;                  170          areacode[i] = sOutside;
158          isvalid[i]  = false;                     171          isvalid[i]  = false;
159          gxx[i].set(kInfinity, kInfinity, kInf    172          gxx[i].set(kInfinity, kInfinity, kInfinity);
160       }                                           173       }
161    }                                              174    }
162                                                   175 
163    G4ThreeVector p = ComputeLocalPoint(gp);       176    G4ThreeVector p = ComputeLocalPoint(gp);
164    G4ThreeVector v = ComputeLocalDirection(gv)    177    G4ThreeVector v = ComputeLocalDirection(gv);
165                                                   178 
166    //                                             179    //
167    // special case!                               180    // special case!
168    // if p is on surface, distance = 0.           181    // if p is on surface, distance = 0. 
169    //                                             182    //
170                                                   183 
171    if (std::fabs(p.z()) == 0.)    // if p is o << 184    if (std::fabs(p.z()) == 0.) {   // if p is on the plane
172    {                                           << 
173       distance[0] = 0;                            185       distance[0] = 0;
174       G4ThreeVector xx = p;                       186       G4ThreeVector xx = p;
175       gxx[0] = ComputeGlobalPoint(xx);            187       gxx[0] = ComputeGlobalPoint(xx);
176                                                   188       
177       if (validate == kValidateWithTol)        << 189       if (validate == kValidateWithTol) {
178       {                                        << 
179          areacode[0] = GetAreaCode(xx);           190          areacode[0] = GetAreaCode(xx);
180          if (!IsOutside(areacode[0]))          << 191          if (!IsOutside(areacode[0])) {
181          {                                     << 
182             isvalid[0] = true;                    192             isvalid[0] = true;
183          }                                        193          }
184       }                                        << 194       } else if (validate == kValidateWithoutTol) {
185       else if (validate == kValidateWithoutTol << 
186       {                                        << 
187          areacode[0] = GetAreaCode(xx, false);    195          areacode[0] = GetAreaCode(xx, false);
188          if (IsInside(areacode[0]))            << 196          if (IsInside(areacode[0])) {
189          {                                     << 
190             isvalid[0] = true;                    197             isvalid[0] = true;
191          }                                        198          }
192       }                                        << 199       } else { // kDontValidate
193       else  // kDontValidate                   << 
194       {                                        << 
195          areacode[0] = sInside;                   200          areacode[0] = sInside;
196          isvalid[0] = true;                       201          isvalid[0] = true;
197       }                                           202       }
                                                   >> 203 
198       return 1;                                   204       return 1;
199    }                                              205    }
200    //                                             206    //
201    // special case end                            207    // special case end
202    //                                             208    //
203                                                   209    
204    if (v.z() == 0)                             << 210    if (v.z() == 0) { 
205    {                                           << 211 
206       fCurStatWithV.SetCurrentStatus(0, gxx[0]    212       fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0], 
207                                      isvalid[0    213                                      isvalid[0], 0, validate, &gp, &gv);
208       return 0;                                   214       return 0;
209    }                                              215    }
210                                                   216    
211    distance[0] = - (p.z() / v.z());               217    distance[0] = - (p.z() / v.z());
212                                                   218       
213    G4ThreeVector xx = p + distance[0]*v;          219    G4ThreeVector xx = p + distance[0]*v;
214    gxx[0] = ComputeGlobalPoint(xx);               220    gxx[0] = ComputeGlobalPoint(xx);
215                                                   221 
216    if (validate == kValidateWithTol)           << 222    if (validate == kValidateWithTol) {
217    {                                           << 
218       areacode[0] = GetAreaCode(xx);              223       areacode[0] = GetAreaCode(xx);
219       if (!IsOutside(areacode[0]))             << 224       if (!IsOutside(areacode[0])) {
220       {                                        << 
221          if (distance[0] >= 0) isvalid[0] = tr    225          if (distance[0] >= 0) isvalid[0] = true;
222       }                                           226       }
223    }                                           << 227    } else if (validate == kValidateWithoutTol) {
224    else if (validate == kValidateWithoutTol)   << 
225    {                                           << 
226       areacode[0] = GetAreaCode(xx, false);       228       areacode[0] = GetAreaCode(xx, false);
227       if (IsInside(areacode[0]))               << 229       if (IsInside(areacode[0])) {
228       {                                        << 
229          if (distance[0] >= 0) isvalid[0] = tr    230          if (distance[0] >= 0) isvalid[0] = true;
230       }                                           231       }
231    }                                           << 232    } else { // kDontValidate
232    else  // kDontValidate                      << 
233    {                                           << 
234       areacode[0] = sInside;                      233       areacode[0] = sInside;
235          if (distance[0] >= 0) isvalid[0] = tr    234          if (distance[0] >= 0) isvalid[0] = true;
236    }                                              235    }
237                                                   236 
238    fCurStatWithV.SetCurrentStatus(0, gxx[0], d    237    fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
239                                   isvalid[0],     238                                   isvalid[0], 1, validate, &gp, &gv);
240                                                   239 
241 #ifdef G4TWISTDEBUG                               240 #ifdef G4TWISTDEBUG
242    G4cerr << "ERROR - G4TwistTubsFlatSide::Dis    241    G4cerr << "ERROR - G4TwistTubsFlatSide::DistanceToSurface(p,v)" << G4endl;
243    G4cerr << "        Name        : " << GetNa    242    G4cerr << "        Name        : " << GetName() << G4endl;
244    G4cerr << "        xx          : " << xx <<    243    G4cerr << "        xx          : " << xx << G4endl;
245    G4cerr << "        gxx[0]      : " << gxx[0    244    G4cerr << "        gxx[0]      : " << gxx[0] << G4endl;
246    G4cerr << "        dist[0]     : " << dista    245    G4cerr << "        dist[0]     : " << distance[0] << G4endl;
247    G4cerr << "        areacode[0] : " << areac    246    G4cerr << "        areacode[0] : " << areacode[0] << G4endl;
248    G4cerr << "        isvalid[0]  : " << isval    247    G4cerr << "        isvalid[0]  : " << isvalid[0]  << G4endl;
249 #endif                                            248 #endif
250    return 1;                                      249    return 1;
251 }                                                 250 }
252                                                   251 
253 //============================================    252 //=====================================================================
254 //* DistanceToSurface(p) ---------------------    253 //* DistanceToSurface(p) ----------------------------------------------
255                                                   254 
256 G4int G4TwistTubsFlatSide::DistanceToSurface(c << 255 G4int G4TwistTubsFlatSide::DistanceToSurface(const G4ThreeVector &gp,
257                                              G    256                                              G4ThreeVector  gxx[],
258                                              G    257                                              G4double       distance[],
259                                              G    258                                              G4int          areacode[])
260 {                                                 259 {
261    // Calculate distance to plane in local coo    260    // Calculate distance to plane in local coordinate,
262    // then return distance and global intersec    261    // then return distance and global intersection points.
263    //                                             262    //  
264                                                   263 
265    fCurStat.ResetfDone(kDontValidate, &gp);       264    fCurStat.ResetfDone(kDontValidate, &gp);
266                                                   265 
267    if (fCurStat.IsDone())                      << 266    if (fCurStat.IsDone()) {
268    {                                           << 267       G4int i;
269       for (G4int i=0; i<fCurStat.GetNXX(); ++i << 268       for (i=0; i<fCurStat.GetNXX(); i++) {
270       {                                        << 
271          gxx[i] = fCurStat.GetXX(i);              269          gxx[i] = fCurStat.GetXX(i);
272          distance[i] = fCurStat.GetDistance(i)    270          distance[i] = fCurStat.GetDistance(i);
273          areacode[i] = fCurStat.GetAreacode(i)    271          areacode[i] = fCurStat.GetAreacode(i);
274       }                                           272       }
275       return fCurStat.GetNXX();                   273       return fCurStat.GetNXX();
276    }                                           << 274    } else {
277    else   // initialize                        << 275       // initialize
278    {                                           << 276       G4int i;
279       for (auto i=0; i<2; ++i)                 << 277       for (i=0; i<2; i++) {
280       {                                        << 
281          distance[i] = kInfinity;                 278          distance[i] = kInfinity;
282          areacode[i] = sOutside;                  279          areacode[i] = sOutside;
283          gxx[i].set(kInfinity, kInfinity, kInf    280          gxx[i].set(kInfinity, kInfinity, kInfinity);
284       }                                           281       }
285    }                                              282    }
286                                                   283    
287    G4ThreeVector p = ComputeLocalPoint(gp);       284    G4ThreeVector p = ComputeLocalPoint(gp);
288    G4ThreeVector xx;                              285    G4ThreeVector xx;
289                                                   286 
290    // The plane is placed on origin with makin    287    // The plane is placed on origin with making its normal 
291    // parallel to z-axis.                         288    // parallel to z-axis. 
292    if (std::fabs(p.z()) <= 0.5 * kCarTolerance << 289    if (std::fabs(p.z()) <= 0.5 * kCarTolerance) {   // if p is on the plane, return 1
293    {                                           << 290       distance[0] = 0;
294      distance[0] = 0;                          << 
295       xx = p;                                     291       xx = p;
296    }                                           << 292    } else {
297    else                                        << 
298    {                                           << 
299       distance[0] = std::fabs(p.z());             293       distance[0] = std::fabs(p.z());
300       xx.set(p.x(), p.y(), 0);                    294       xx.set(p.x(), p.y(), 0);  
301    }                                              295    }
302                                                   296 
303    gxx[0] = ComputeGlobalPoint(xx);               297    gxx[0] = ComputeGlobalPoint(xx);
304    areacode[0] = sInside;                         298    areacode[0] = sInside;
305    G4bool isvalid = true;                         299    G4bool isvalid = true;
306    fCurStat.SetCurrentStatus(0, gxx[0], distan    300    fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
307                              isvalid, 1, kDont    301                              isvalid, 1, kDontValidate, &gp);
308    return 1;                                      302    return 1;
                                                   >> 303 
309 }                                                 304 }
310                                                   305 
311 //============================================    306 //=====================================================================
312 //* GetAreaCode ------------------------------    307 //* GetAreaCode -------------------------------------------------------
313                                                   308 
314 G4int G4TwistTubsFlatSide::GetAreaCode(const G    309 G4int G4TwistTubsFlatSide::GetAreaCode(const G4ThreeVector &xx, 
315                                              G << 310                                        G4bool withTol)
316 {                                                 311 {
317    const G4double rtol                            312    const G4double rtol
318      = 0.5*G4GeometryTolerance::GetInstance()-    313      = 0.5*G4GeometryTolerance::GetInstance()->GetRadialTolerance();
319                                                   314    
320    G4int areacode = sInside;                      315    G4int areacode = sInside;
321                                                   316 
322    if (fAxis[0] == kRho && fAxis[1] == kPhi)   << 317    if (fAxis[0] == kRho && fAxis[1] == kPhi) {
323    {                                           << 
324       G4int rhoaxis = 0;                          318       G4int rhoaxis = 0;
                                                   >> 319      // G4int phiaxis = 0;
325                                                   320 
326       G4ThreeVector dphimin;   // direction of    321       G4ThreeVector dphimin;   // direction of phi-minimum boundary
327       G4ThreeVector dphimax;   // direction of    322       G4ThreeVector dphimax;   // direction of phi-maximum boundary
328       dphimin = GetCorner(sC0Max1Min);            323       dphimin = GetCorner(sC0Max1Min);
329       dphimax = GetCorner(sC0Max1Max);            324       dphimax = GetCorner(sC0Max1Max);   
330                                                   325       
331       if (withTol)                             << 326       if (withTol) {
332       {                                        << 327 
333          G4bool isoutside = false;                328          G4bool isoutside = false;
334                                                   329 
335          // test boundary of rho-axis             330          // test boundary of rho-axis
336                                                   331 
337          if (xx.getRho() <= fAxisMin[rhoaxis]  << 332          if (xx.getRho() <= fAxisMin[rhoaxis] + rtol) {
338          {                                     << 333 
339             areacode |= (sAxis0 & (sAxisRho|sA << 334             areacode |= (sAxis0 & (sAxisRho | sAxisMin)) | sBoundary; // rho-min
340             if (xx.getRho() < fAxisMin[rhoaxis    335             if (xx.getRho() < fAxisMin[rhoaxis] - rtol) isoutside = true; 
341                                                   336             
342          }                                     << 337          } else if (xx.getRho() >= fAxisMax[rhoaxis] - rtol) {
343          else if (xx.getRho() >= fAxisMax[rhoa << 338 
344          {                                     << 339             areacode |= (sAxis0 & (sAxisRho | sAxisMax)) | sBoundary; // rho-max
345             areacode |= (sAxis0 & (sAxisRho|sA << 340             if (xx.getRho() > fAxisMax[rhoaxis] + rtol) isoutside = true; 
346             if (xx.getRho() > fAxisMax[rhoaxis << 341 
347          }                                        342          }         
348                                                   343          
349          // test boundary of phi-axis             344          // test boundary of phi-axis
350                                                   345 
351          if (AmIOnLeftSide(xx, dphimin) >= 0)  << 346          if (AmIOnLeftSide(xx, dphimin) >= 0) {           // xx is on dphimin
352          {                                     << 347 
353             areacode |= (sAxis1 & (sAxisPhi |     348             areacode |= (sAxis1 & (sAxisPhi | sAxisMin)); 
354             if   ((areacode & sBoundary) != 0) << 349             if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
355             else                        areaco    350             else                        areacode |= sBoundary;
356                                                   351 
357             if (AmIOnLeftSide(xx, dphimin) > 0    352             if (AmIOnLeftSide(xx, dphimin) > 0) isoutside = true; 
358                                                   353 
359          }                                     << 354          } else if (AmIOnLeftSide(xx, dphimax) <= 0) {    // xx is on dphimax
360          else if (AmIOnLeftSide(xx, dphimax) < << 355 
361          {                                     << 
362             areacode |= (sAxis1 & (sAxisPhi |     356             areacode |= (sAxis1 & (sAxisPhi | sAxisMax)); 
363             if   ((areacode & sBoundary) != 0) << 357             if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
364             else                        areaco    358             else                        areacode |= sBoundary;
365                                                   359 
366             if (AmIOnLeftSide(xx, dphimax) < 0 << 360             if (AmIOnLeftSide(xx, dphimax) < 0) isoutside = true; 
                                                   >> 361 
367          }                                        362          }
368                                                   363          
369          // if isoutside = true, clear inside     364          // if isoutside = true, clear inside bit.
370          // if not on boundary, add axis infor    365          // if not on boundary, add axis information. 
371                                                   366 
372          if (isoutside)                        << 367          if (isoutside) {
373          {                                     << 
374             G4int tmpareacode = areacode & (~s    368             G4int tmpareacode = areacode & (~sInside);
375             areacode = tmpareacode;               369             areacode = tmpareacode;
376          }                                     << 370          } else if ((areacode & sBoundary) != sBoundary) {
377          else if ((areacode & sBoundary) != sB << 
378          {                                     << 
379             areacode |= (sAxis0 & sAxisRho) |     371             areacode |= (sAxis0 & sAxisRho) | (sAxis1 & sAxisPhi);
380          }                                        372          }
381                                                   373 
382       }                                        << 374       } else {
383       else                                     << 375 
384       {                                        << 
385          // out of boundary of rho-axis           376          // out of boundary of rho-axis
386                                                   377 
387          if (xx.getRho() < fAxisMin[rhoaxis])  << 378          if (xx.getRho() < fAxisMin[rhoaxis]) {
388          {                                     << 
389             areacode |= (sAxis0 & (sAxisRho |     379             areacode |= (sAxis0 & (sAxisRho | sAxisMin)) | sBoundary;
390          }                                     << 380          } else if (xx.getRho() > fAxisMax[rhoaxis]) {
391          else if (xx.getRho() > fAxisMax[rhoax << 
392          {                                     << 
393             areacode |= (sAxis0 & (sAxisRho |     381             areacode |= (sAxis0 & (sAxisRho | sAxisMax)) | sBoundary;
394          }                                        382          }
395                                                   383          
396          // out of boundary of phi-axis           384          // out of boundary of phi-axis
397                                                   385 
398          if (AmIOnLeftSide(xx, dphimin, false) << 386          if (AmIOnLeftSide(xx, dphimin, false) >= 0) {       // xx is leftside or
399          {                                     << 387             areacode |= (sAxis1 & (sAxisPhi | sAxisMin)) ;   // boundary of dphimin
400            areacode |= (sAxis1 & (sAxisPhi | s << 388             if   (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
401            if   ((areacode & sBoundary) != 0)  << 389             else                        areacode |= sBoundary;
402            else                        areacod << 
403                                                   390 
404          }                                     << 391          } else if (AmIOnLeftSide(xx, dphimax, false) <= 0) { // xx is rightside or
405          else if (AmIOnLeftSide(xx, dphimax, f << 392             areacode |= (sAxis1 & (sAxisPhi | sAxisMax)) ;    // boundary of dphimax
406          {                                     << 393             if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
407            areacode |= (sAxis1 & (sAxisPhi | s << 394             else                        areacode |= sBoundary;
408            if   ((areacode & sBoundary) != 0)  << 
409            else                        areacod << 
410                                                   395            
411          }                                        396          }
412                                                   397 
413          if ((areacode & sBoundary) != sBounda    398          if ((areacode & sBoundary) != sBoundary) {
414             areacode |= (sAxis0 & sAxisRho) |     399             areacode |= (sAxis0 & sAxisRho) | (sAxis1 & sAxisPhi);
415          }                                        400          }
416                                                   401 
417       }                                           402       }
418       return areacode;                            403       return areacode;
419    }                                           << 404    } else {
420    else                                        << 405 
421    {                                           << 
422       std::ostringstream message;                 406       std::ostringstream message;
423       message << "Feature NOT implemented !" <    407       message << "Feature NOT implemented !" << G4endl
424               << "        fAxis[0] = " << fAxi    408               << "        fAxis[0] = " << fAxis[0] << G4endl
425               << "        fAxis[1] = " << fAxi    409               << "        fAxis[1] = " << fAxis[1];
426       G4Exception("G4TwistTubsFlatSide::GetAre    410       G4Exception("G4TwistTubsFlatSide::GetAreaCode()", "GeomSolids0001",
427                   FatalException, message);       411                   FatalException, message);
428    }                                              412    }
429    return areacode;                               413    return areacode;
430 }                                                 414 }
431                                                << 415                                       
                                                   >> 416                                            
432 //============================================    417 //=====================================================================
433 //* SetCorners -------------------------------    418 //* SetCorners --------------------------------------------------------
434                                                   419 
435 void G4TwistTubsFlatSide::SetCorners()            420 void G4TwistTubsFlatSide::SetCorners()
436 {                                                 421 {
437    // Set Corner points in local coodinate.       422    // Set Corner points in local coodinate.
438                                                   423    
439    if (fAxis[0] == kRho && fAxis[1] == kPhi)   << 424    if (fAxis[0] == kRho && fAxis[1] == kPhi) {
440    {                                           << 425        
441       G4int rhoaxis = 0;  // kRho                 426       G4int rhoaxis = 0;  // kRho
442       G4int phiaxis = 1;  // kPhi                 427       G4int phiaxis = 1;  // kPhi
443                                                   428       
444       G4double x, y, z;                           429       G4double x, y, z;
445       // corner of Axis0min and Axis1min          430       // corner of Axis0min and Axis1min
446          x = fAxisMin[rhoaxis]*std::cos(fAxisM    431          x = fAxisMin[rhoaxis]*std::cos(fAxisMin[phiaxis]);
447          y = fAxisMin[rhoaxis]*std::sin(fAxisM    432          y = fAxisMin[rhoaxis]*std::sin(fAxisMin[phiaxis]);
448          z = 0;                                   433          z = 0;
449          SetCorner(sC0Min1Min, x, y, z);          434          SetCorner(sC0Min1Min, x, y, z);
450       // corner of Axis0max and Axis1min          435       // corner of Axis0max and Axis1min
451          x = fAxisMax[rhoaxis]*std::cos(fAxisM    436          x = fAxisMax[rhoaxis]*std::cos(fAxisMin[phiaxis]);
452          y = fAxisMax[rhoaxis]*std::sin(fAxisM    437          y = fAxisMax[rhoaxis]*std::sin(fAxisMin[phiaxis]);
453          z = 0;                                   438          z = 0;
454          SetCorner(sC0Max1Min, x, y, z);          439          SetCorner(sC0Max1Min, x, y, z);
455       // corner of Axis0max and Axis1max          440       // corner of Axis0max and Axis1max
456          x = fAxisMax[rhoaxis]*std::cos(fAxisM    441          x = fAxisMax[rhoaxis]*std::cos(fAxisMax[phiaxis]);
457          y = fAxisMax[rhoaxis]*std::sin(fAxisM    442          y = fAxisMax[rhoaxis]*std::sin(fAxisMax[phiaxis]);
458          z = 0;                                   443          z = 0;
459          SetCorner(sC0Max1Max, x, y, z);          444          SetCorner(sC0Max1Max, x, y, z);
460       // corner of Axis0min and Axis1max          445       // corner of Axis0min and Axis1max
461          x = fAxisMin[rhoaxis]*std::cos(fAxisM    446          x = fAxisMin[rhoaxis]*std::cos(fAxisMax[phiaxis]);
462          y = fAxisMin[rhoaxis]*std::sin(fAxisM    447          y = fAxisMin[rhoaxis]*std::sin(fAxisMax[phiaxis]);
463          z = 0;                                   448          z = 0;
464          SetCorner(sC0Min1Max, x, y, z);          449          SetCorner(sC0Min1Max, x, y, z);
465                                                   450        
466    }                                           << 451    } else {
467    else                                        << 
468    {                                           << 
469       std::ostringstream message;                 452       std::ostringstream message;
470       message << "Feature NOT implemented !" <    453       message << "Feature NOT implemented !" << G4endl
471               << "        fAxis[0] = " << fAxi    454               << "        fAxis[0] = " << fAxis[0] << G4endl
472               << "        fAxis[1] = " << fAxi    455               << "        fAxis[1] = " << fAxis[1];
473       G4Exception("G4TwistTubsFlatSide::SetCor    456       G4Exception("G4TwistTubsFlatSide::SetCorners()", "GeomSolids0001",
474                   FatalException, message);       457                   FatalException, message);
475    }                                              458    }
476 }                                                 459 }
477                                                   460 
478 //============================================    461 //=====================================================================
479 //* SetBoundaries() --------------------------    462 //* SetBoundaries() ---------------------------------------------------
480                                                   463 
481 void G4TwistTubsFlatSide::SetBoundaries()         464 void G4TwistTubsFlatSide::SetBoundaries()
482 {                                                 465 {
483    // Set direction-unit vector of phi-boundar    466    // Set direction-unit vector of phi-boundary-lines in local coodinate.
484    // Don't call the function twice.              467    // Don't call the function twice.
485                                                   468    
486    if (fAxis[0] == kRho && fAxis[1] == kPhi)   << 469    if (fAxis[0] == kRho && fAxis[1] == kPhi) {
487    {                                           << 470    
488       G4ThreeVector direction;                    471       G4ThreeVector direction;
489       // sAxis0 & sAxisMin                        472       // sAxis0 & sAxisMin
490       direction = GetCorner(sC0Min1Max) - GetC    473       direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
491       direction = direction.unit();               474       direction = direction.unit();
492       SetBoundary(sAxis0 & (sAxisPhi | sAxisMi    475       SetBoundary(sAxis0 & (sAxisPhi | sAxisMin), direction,
493                   GetCorner(sC0Min1Min), sAxis    476                   GetCorner(sC0Min1Min), sAxisPhi);
494                                                   477                   
495       // sAxis0 & sAxisMax                        478       // sAxis0 & sAxisMax
496       direction = GetCorner(sC0Max1Max) - GetC    479       direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
497       direction = direction.unit();               480       direction = direction.unit();
498       SetBoundary(sAxis0 & (sAxisPhi | sAxisMa    481       SetBoundary(sAxis0 & (sAxisPhi | sAxisMax), direction,
499                   GetCorner(sC0Max1Min), sAxis    482                   GetCorner(sC0Max1Min), sAxisPhi);
500                                                   483 
501       // sAxis1 & sAxisMin                        484       // sAxis1 & sAxisMin
502       direction = GetCorner(sC0Max1Min) - GetC    485       direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
503       direction = direction.unit();               486       direction = direction.unit();
504       SetBoundary(sAxis1 & (sAxisRho | sAxisMi    487       SetBoundary(sAxis1 & (sAxisRho | sAxisMin), direction,
505                   GetCorner(sC0Min1Min), sAxis    488                   GetCorner(sC0Min1Min), sAxisRho);
506                                                   489       
507       // sAxis1 & sAxisMax                        490       // sAxis1 & sAxisMax
508       direction = GetCorner(sC0Max1Max) - GetC    491       direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
509       direction = direction.unit();               492       direction = direction.unit();
510       SetBoundary(sAxis1 & (sAxisRho | sAxisMa    493       SetBoundary(sAxis1 & (sAxisRho | sAxisMax), direction,
511                   GetCorner(sC0Min1Max), sAxis    494                   GetCorner(sC0Min1Max), sAxisPhi);
512    }                                           << 495    } else {
513    else                                        << 
514    {                                           << 
515       std::ostringstream message;                 496       std::ostringstream message;
516       message << "Feature NOT implemented !" <    497       message << "Feature NOT implemented !" << G4endl
517               << "        fAxis[0] = " << fAxi    498               << "        fAxis[0] = " << fAxis[0] << G4endl
518               << "        fAxis[1] = " << fAxi    499               << "        fAxis[1] = " << fAxis[1];
519       G4Exception("G4TwistTubsFlatSide::SetBou    500       G4Exception("G4TwistTubsFlatSide::SetBoundaries()", "GeomSolids0001",
520                   FatalException, message);       501                   FatalException, message);
521    }                                              502    }
522 }                                                 503 }
523                                                   504 
524 //============================================    505 //=====================================================================
525 //* GetFacets() ------------------------------    506 //* GetFacets() -------------------------------------------------------
526                                                   507 
527 void G4TwistTubsFlatSide::GetFacets( G4int k,     508 void G4TwistTubsFlatSide::GetFacets( G4int k, G4int n, G4double xyz[][3],
528                                      G4int fac    509                                      G4int faces[][4], G4int iside ) 
529 {                                                 510 {
                                                   >> 511 
530   G4ThreeVector p ;                               512   G4ThreeVector p ;
531                                                   513 
532   G4double rmin = fAxisMin[0] ;                   514   G4double rmin = fAxisMin[0] ;
533   G4double rmax = fAxisMax[0] ;                   515   G4double rmax = fAxisMax[0] ;
534   G4double phimin, phimax ;                       516   G4double phimin, phimax ;
535                                                   517 
536   G4double r,phi ;                                518   G4double r,phi ;
                                                   >> 519 
                                                   >> 520   G4int i,j ;
                                                   >> 521 
537   G4int nnode,nface ;                             522   G4int nnode,nface ;
538                                                   523 
539   for ( G4int i = 0 ; i<n ; ++i )              << 524   for ( i = 0 ; i<n ; i++ ) {
540   {                                            << 525 
541     r = rmin + i*(rmax-rmin)/(n-1) ;              526     r = rmin + i*(rmax-rmin)/(n-1) ;
542                                                   527 
543     phimin = GetBoundaryMin(r) ;                  528     phimin = GetBoundaryMin(r) ; 
544     phimax = GetBoundaryMax(r) ;                  529     phimax = GetBoundaryMax(r) ;
545                                                   530 
546     for ( G4int j = 0 ; j<k ; ++j )            << 531     for ( j = 0 ; j<k ; j++ )
547     {                                             532     {
548       phi = phimin + j*(phimax-phimin)/(k-1) ;    533       phi = phimin + j*(phimax-phimin)/(k-1) ;
549                                                   534 
550       nnode = GetNode(i,j,k,n,iside) ;            535       nnode = GetNode(i,j,k,n,iside) ;
551       p = SurfacePoint(phi,r,true) ;  // surfa    536       p = SurfacePoint(phi,r,true) ;  // surface point in global coord.system
552                                                   537 
553       xyz[nnode][0] = p.x() ;                     538       xyz[nnode][0] = p.x() ;
554       xyz[nnode][1] = p.y() ;                     539       xyz[nnode][1] = p.y() ;
555       xyz[nnode][2] = p.z() ;                     540       xyz[nnode][2] = p.z() ;
556                                                   541 
557       if ( i<n-1 && j<k-1 )    // conterclock  << 542       if ( i<n-1 && j<k-1 ) {   // conterclock wise filling
558       {                                        << 543         
559         nface = GetFace(i,j,k,n,iside) ;          544         nface = GetFace(i,j,k,n,iside) ;
560                                                   545 
561         if (fHandedness < 0)   // lower side   << 546         if (fHandedness < 0) {  // lower side 
562         {                                      << 547           faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,-1) * ( GetNode(i  ,j  ,k,n,iside)+1) ;  
563           faces[nface][0] = GetEdgeVisibility( << 548           faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,-1) * ( GetNode(i  ,j+1,k,n,iside)+1) ;
564                           * ( GetNode(i  ,j  , << 549           faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,-1) * ( GetNode(i+1,j+1,k,n,iside)+1) ;
565           faces[nface][1] = GetEdgeVisibility( << 550           faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,-1) * ( GetNode(i+1,j  ,k,n,iside)+1) ;
566                           * ( GetNode(i  ,j+1, << 551         } else {                // upper side
567           faces[nface][2] = GetEdgeVisibility( << 552           faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1) * ( GetNode(i  ,j  ,k,n,iside)+1) ;  
568                           * ( GetNode(i+1,j+1, << 553           faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1) * ( GetNode(i+1,j  ,k,n,iside)+1) ;
569           faces[nface][3] = GetEdgeVisibility( << 554           faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1) * ( GetNode(i+1,j+1,k,n,iside)+1) ;
570                           * ( GetNode(i+1,j  , << 555           faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1) * ( GetNode(i  ,j+1,k,n,iside)+1) ;
571         }                                      << 556 
572         else                   // upper side   << 
573         {                                      << 
574           faces[nface][0] = GetEdgeVisibility( << 
575                           * ( GetNode(i  ,j  , << 
576           faces[nface][1] = GetEdgeVisibility( << 
577                           * ( GetNode(i+1,j  , << 
578           faces[nface][2] = GetEdgeVisibility( << 
579                           * ( GetNode(i+1,j+1, << 
580           faces[nface][3] = GetEdgeVisibility( << 
581                           * ( GetNode(i  ,j+1, << 
582         }                                         557         }
                                                   >> 558 
                                                   >> 559 
                                                   >> 560 
583       }                                           561       }
584     }                                             562     }
585   }                                               563   }
586 }                                                 564 }
587                                                   565