Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/specific/src/G4TwistTrapFlatSide.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/G4TwistTrapFlatSide.cc (Version 11.3.0) and /geometry/solids/specific/src/G4TwistTrapFlatSide.cc (Version 9.0.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 // G4TwistTrapFlatSide implementation          << 
 27 //                                                 26 //
 28 // Author: 30-Aug-2002 - O.Link (Oliver.Link@c <<  27 // $Id: G4TwistTrapFlatSide.cc,v 1.6 2007/05/23 09:31:02 gcosmo Exp $
                                                   >>  28 // GEANT4 tag $Name: geant4-09-00-patch-01 $
                                                   >>  29 //
                                                   >>  30 // 
                                                   >>  31 // --------------------------------------------------------------------
                                                   >>  32 // GEANT 4 class source file
                                                   >>  33 //
                                                   >>  34 //
                                                   >>  35 // G4TwistTrapFlatSide.cc
                                                   >>  36 //
                                                   >>  37 // Author: 
                                                   >>  38 //   30-Aug-2002 - O.Link (Oliver.Link@cern.ch)
                                                   >>  39 //
 29 // -------------------------------------------     40 // --------------------------------------------------------------------
 30                                                    41 
 31 #include "G4TwistTrapFlatSide.hh"                  42 #include "G4TwistTrapFlatSide.hh"
 32                                                    43 
 33 //============================================     44 //=====================================================================
 34 //* constructors -----------------------------     45 //* constructors ------------------------------------------------------
 35                                                    46 
 36 G4TwistTrapFlatSide::G4TwistTrapFlatSide( cons <<  47 G4TwistTrapFlatSide::G4TwistTrapFlatSide( const G4String        &name,
 37                                                <<  48                               G4double      PhiTwist,
 38                                                <<  49                               G4double      pDx1,
 39                                                <<  50                               G4double      pDx2,
 40                                                <<  51                               G4double      pDy,
 41                                                <<  52                               G4double      pDz,
 42                                                <<  53                               G4double      pAlpha,
 43                                                <<  54                               G4double      pPhi,
 44                                                <<  55                               G4double      pTheta,
 45                                                <<  56                               G4int         handedness) 
                                                   >>  57 
 46   : G4VTwistSurface(name)                          58   : G4VTwistSurface(name)
 47 {                                                  59 {
 48    fHandedness = handedness;   // +z = +ve, -z     60    fHandedness = handedness;   // +z = +ve, -z = -ve
 49                                                    61 
 50    fDx1 = pDx1 ;                                   62    fDx1 = pDx1 ;
 51    fDx2 = pDx2 ;                                   63    fDx2 = pDx2 ;
 52    fDy = pDy ;                                     64    fDy = pDy ;
 53    fDz = pDz ;                                     65    fDz = pDz ;
 54    fAlpha = pAlpha ;                               66    fAlpha = pAlpha ;
 55    fTAlph = std::tan(fAlpha) ;                     67    fTAlph = std::tan(fAlpha) ;
 56    fPhi  = pPhi ;                                  68    fPhi  = pPhi ;
 57    fTheta = pTheta ;                               69    fTheta = pTheta ;
 58                                                    70 
 59    fdeltaX = 2 * fDz * std::tan(fTheta) * std:     71    fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi)  ;
 60      // dx in surface equation                     72      // dx in surface equation
 61    fdeltaY = 2 * fDz * std::tan(fTheta) * std:     73    fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi)  ;
 62      // dy in surface equation                     74      // dy in surface equation
 63                                                    75 
 64    fPhiTwist = PhiTwist ;                          76    fPhiTwist = PhiTwist ;
 65                                                    77 
 66    fCurrentNormal.normal.set( 0, 0, (fHandedne     78    fCurrentNormal.normal.set( 0, 0, (fHandedness < 0 ? -1 : 1)); 
 67          // Unit vector, in local coordinate s     79          // Unit vector, in local coordinate system
 68    fRot.rotateZ( fHandedness > 0                   80    fRot.rotateZ( fHandedness > 0 
 69                  ? 0.5 * fPhiTwist                 81                  ? 0.5 * fPhiTwist
 70                  : -0.5 * fPhiTwist );             82                  : -0.5 * fPhiTwist );
 71                                                    83 
 72    fTrans.set(                                     84    fTrans.set(
 73               fHandedness > 0 ? 0.5*fdeltaX :      85               fHandedness > 0 ? 0.5*fdeltaX : -0.5*fdeltaX , 
 74               fHandedness > 0 ? 0.5*fdeltaY :      86               fHandedness > 0 ? 0.5*fdeltaY : -0.5*fdeltaY ,
 75               fHandedness > 0 ? fDz : -fDz ) ;     87               fHandedness > 0 ? fDz : -fDz ) ;
 76                                                    88 
 77    fIsValidNorm = true;                            89    fIsValidNorm = true;
 78                                                    90 
 79                                                    91 
 80    fAxis[0] = kXAxis ;                             92    fAxis[0] = kXAxis ;
 81    fAxis[1] = kYAxis ;                             93    fAxis[1] = kYAxis ;
 82    fAxisMin[0] = kInfinity ;  // x-Axis cannot     94    fAxisMin[0] = kInfinity ;  // x-Axis cannot be fixed, because it 
 83    fAxisMax[0] =  kInfinity ; // depends on y      95    fAxisMax[0] =  kInfinity ; // depends on y
 84    fAxisMin[1] = -fDy ;  // y - axis               96    fAxisMin[1] = -fDy ;  // y - axis
 85    fAxisMax[1] =  fDy ;                            97    fAxisMax[1] =  fDy ;
 86                                                    98 
 87    SetCorners();                                   99    SetCorners();
 88    SetBoundaries();                               100    SetBoundaries();
 89 }                                                 101 }
 90                                                   102 
 91                                                   103 
 92 //============================================    104 //=====================================================================
 93 //* Fake default constructor -----------------    105 //* Fake default constructor ------------------------------------------
 94                                                   106 
 95 G4TwistTrapFlatSide::G4TwistTrapFlatSide( __vo    107 G4TwistTrapFlatSide::G4TwistTrapFlatSide( __void__& a )
 96   : G4VTwistSurface(a)                            108   : G4VTwistSurface(a)
 97 {                                                 109 {
 98 }                                                 110 }
 99                                                   111 
100                                                   112 
101 //============================================    113 //=====================================================================
102 //* destructor -------------------------------    114 //* destructor --------------------------------------------------------
103                                                   115 
104 G4TwistTrapFlatSide::~G4TwistTrapFlatSide() =  << 116 G4TwistTrapFlatSide::~G4TwistTrapFlatSide()
                                                   >> 117 {
                                                   >> 118 }
105                                                   119 
106 //============================================    120 //=====================================================================
107 //* GetNormal --------------------------------    121 //* GetNormal ---------------------------------------------------------
108                                                   122 
109 G4ThreeVector G4TwistTrapFlatSide::GetNormal(c << 123 G4ThreeVector G4TwistTrapFlatSide::GetNormal(const G4ThreeVector & /* xx */ , 
110                                              G    124                                              G4bool isGlobal)
111 {                                                 125 {
112    if (isGlobal)                               << 126    if (isGlobal) {
113    {                                           << 
114       return ComputeGlobalDirection(fCurrentNo    127       return ComputeGlobalDirection(fCurrentNormal.normal);
115    }                                           << 128    } else {
116    else                                        << 
117    {                                           << 
118       return fCurrentNormal.normal;               129       return fCurrentNormal.normal;
119    }                                              130    }
120 }                                                 131 }
121                                                   132 
122 //============================================    133 //=====================================================================
123 //* DistanceToSurface(p, v) ------------------    134 //* DistanceToSurface(p, v) -------------------------------------------
124                                                   135 
125 G4int G4TwistTrapFlatSide::DistanceToSurface(c << 136 G4int G4TwistTrapFlatSide::DistanceToSurface(const G4ThreeVector &gp,
126                                              c << 137                                        const G4ThreeVector &gv,
127                                                << 138                                              G4ThreeVector  gxx[],
128                                                << 139                                              G4double       distance[],
129                                                << 140                                              G4int          areacode[],
130                                                << 141                                              G4bool         isvalid[],
131                                                << 142                                              EValidate      validate) 
132 {                                                 143 {
133    fCurStatWithV.ResetfDone(validate, &gp, &gv    144    fCurStatWithV.ResetfDone(validate, &gp, &gv);
134                                                   145    
135    if (fCurStatWithV.IsDone())                 << 146    if (fCurStatWithV.IsDone()) {
136    {                                           << 147       G4int i;
137       for (G4int i=0; i<fCurStatWithV.GetNXX() << 148       for (i=0; i<fCurStatWithV.GetNXX(); i++) {
138       {                                        << 
139          gxx[i] = fCurStatWithV.GetXX(i);         149          gxx[i] = fCurStatWithV.GetXX(i);
140          distance[i] = fCurStatWithV.GetDistan    150          distance[i] = fCurStatWithV.GetDistance(i);
141          areacode[i] = fCurStatWithV.GetAreaco    151          areacode[i] = fCurStatWithV.GetAreacode(i);
142          isvalid[i]  = fCurStatWithV.IsValid(i    152          isvalid[i]  = fCurStatWithV.IsValid(i);
143       }                                           153       }
144       return fCurStatWithV.GetNXX();              154       return fCurStatWithV.GetNXX();
145    }                                           << 155    } else {
146    else   // initialize                        << 156       // initialize
147    {                                           << 157       G4int i;
148       for (auto i=0; i<2; ++i)                 << 158       for (i=0; i<2; i++) {
149       {                                        << 
150          distance[i] = kInfinity;                 159          distance[i] = kInfinity;
151          areacode[i] = sOutside;                  160          areacode[i] = sOutside;
152          isvalid[i]  = false;                     161          isvalid[i]  = false;
153          gxx[i].set(kInfinity, kInfinity, kInf    162          gxx[i].set(kInfinity, kInfinity, kInfinity);
154       }                                           163       }
155    }                                              164    }
156                                                   165 
157    G4ThreeVector p = ComputeLocalPoint(gp);       166    G4ThreeVector p = ComputeLocalPoint(gp);
158    G4ThreeVector v = ComputeLocalDirection(gv)    167    G4ThreeVector v = ComputeLocalDirection(gv);
159                                                   168 
160    //                                             169    //
161    // special case!                               170    // special case!
162    // if p is on surface, distance = 0.           171    // if p is on surface, distance = 0. 
163    //                                             172    //
164                                                   173 
165    if (std::fabs(p.z()) == 0.)     // if p is  << 174    if (std::fabs(p.z()) == 0.) {   // if p is on the plane
166    {                                           << 
167       distance[0] = 0;                            175       distance[0] = 0;
168       G4ThreeVector xx = p;                       176       G4ThreeVector xx = p;
169       gxx[0] = ComputeGlobalPoint(xx);            177       gxx[0] = ComputeGlobalPoint(xx);
170                                                   178       
171       if (validate == kValidateWithTol)        << 179       if (validate == kValidateWithTol) {
172       {                                        << 
173          areacode[0] = GetAreaCode(xx);           180          areacode[0] = GetAreaCode(xx);
174          if (!IsOutside(areacode[0]))          << 181          if (!IsOutside(areacode[0])) {
175          {                                     << 
176             isvalid[0] = true;                    182             isvalid[0] = true;
177          }                                        183          }
178       }                                        << 184       } else if (validate == kValidateWithoutTol) {
179       else if (validate == kValidateWithoutTol << 
180       {                                        << 
181          areacode[0] = GetAreaCode(xx, false);    185          areacode[0] = GetAreaCode(xx, false);
182          if (IsInside(areacode[0]))            << 186          if (IsInside(areacode[0])) {
183          {                                     << 
184             isvalid[0] = true;                    187             isvalid[0] = true;
185          }                                        188          }
186       }                                        << 189       } else { // kDontValidate
187       else   // kDontValidate                  << 
188       {                                        << 
189          areacode[0] = sInside;                   190          areacode[0] = sInside;
190          isvalid[0] = true;                       191          isvalid[0] = true;
191       }                                           192       }
                                                   >> 193 
192       return 1;                                   194       return 1;
193    }                                              195    }
194    //                                             196    //
195    // special case end                            197    // special case end
196    //                                             198    //
197                                                   199    
198    if (v.z() == 0) {                              200    if (v.z() == 0) { 
199                                                   201 
200       fCurStatWithV.SetCurrentStatus(0, gxx[0]    202       fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0], 
201                                      isvalid[0    203                                      isvalid[0], 0, validate, &gp, &gv);
202       return 0;                                   204       return 0;
203    }                                              205    }
204                                                   206    
205    distance[0] = - (p.z() / v.z());               207    distance[0] = - (p.z() / v.z());
206                                                   208       
207    G4ThreeVector xx = p + distance[0]*v;          209    G4ThreeVector xx = p + distance[0]*v;
208    gxx[0] = ComputeGlobalPoint(xx);               210    gxx[0] = ComputeGlobalPoint(xx);
209                                                   211 
210    if (validate == kValidateWithTol)           << 212    if (validate == kValidateWithTol) {
211    {                                           << 
212       areacode[0] = GetAreaCode(xx);              213       areacode[0] = GetAreaCode(xx);
213       if (!IsOutside(areacode[0]))             << 214       if (!IsOutside(areacode[0])) {
214       {                                        << 
215          if (distance[0] >= 0) isvalid[0] = tr    215          if (distance[0] >= 0) isvalid[0] = true;
216       }                                           216       }
217    }                                           << 217    } else if (validate == kValidateWithoutTol) {
218    else if (validate == kValidateWithoutTol)   << 
219    {                                           << 
220       areacode[0] = GetAreaCode(xx, false);       218       areacode[0] = GetAreaCode(xx, false);
221       if (IsInside(areacode[0]))               << 219       if (IsInside(areacode[0])) {
222       {                                        << 
223          if (distance[0] >= 0) isvalid[0] = tr    220          if (distance[0] >= 0) isvalid[0] = true;
224       }                                           221       }
225    }                                           << 222    } else { // kDontValidate
226    else   // kDontValidate                     << 
227    {                                           << 
228       areacode[0] = sInside;                      223       areacode[0] = sInside;
229          if (distance[0] >= 0) isvalid[0] = tr    224          if (distance[0] >= 0) isvalid[0] = true;
230    }                                              225    }
231                                                   226 
                                                   >> 227 
232    fCurStatWithV.SetCurrentStatus(0, gxx[0], d    228    fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
233                                   isvalid[0],     229                                   isvalid[0], 1, validate, &gp, &gv);
234                                                   230 
235 #ifdef G4TWISTDEBUG                               231 #ifdef G4TWISTDEBUG
236    G4cerr << "ERROR - G4TwistTrapFlatSide::Dis    232    G4cerr << "ERROR - G4TwistTrapFlatSide::DistanceToSurface(p,v)" << G4endl;
237    G4cerr << "        Name        : " << GetNa    233    G4cerr << "        Name        : " << GetName() << G4endl;
238    G4cerr << "        xx          : " << xx <<    234    G4cerr << "        xx          : " << xx << G4endl;
239    G4cerr << "        gxx[0]      : " << gxx[0    235    G4cerr << "        gxx[0]      : " << gxx[0] << G4endl;
240    G4cerr << "        dist[0]     : " << dista    236    G4cerr << "        dist[0]     : " << distance[0] << G4endl;
241    G4cerr << "        areacode[0] : " << areac    237    G4cerr << "        areacode[0] : " << areacode[0] << G4endl;
242    G4cerr << "        isvalid[0]  : " << isval    238    G4cerr << "        isvalid[0]  : " << isvalid[0]  << G4endl;
243 #endif                                            239 #endif
244    return 1;                                      240    return 1;
245 }                                                 241 }
246                                                   242 
247 //============================================    243 //=====================================================================
248 //* DistanceToSurface(p) ---------------------    244 //* DistanceToSurface(p) ----------------------------------------------
249                                                   245 
250 G4int G4TwistTrapFlatSide::DistanceToSurface(c << 246 G4int G4TwistTrapFlatSide::DistanceToSurface(const G4ThreeVector &gp,
251                                              G    247                                              G4ThreeVector  gxx[],
252                                              G    248                                              G4double       distance[],
253                                              G    249                                              G4int          areacode[])
254 {                                                 250 {
255    // Calculate distance to plane in local coo    251    // Calculate distance to plane in local coordinate,
256    // then return distance and global intersec    252    // then return distance and global intersection points.
257    //                                             253    //  
258                                                   254 
259    fCurStat.ResetfDone(kDontValidate, &gp);       255    fCurStat.ResetfDone(kDontValidate, &gp);
260                                                   256 
261    if (fCurStat.IsDone())                      << 257    if (fCurStat.IsDone()) {
262    {                                           << 258       G4int i;
263       for (G4int i=0; i<fCurStat.GetNXX(); ++i << 259       for (i=0; i<fCurStat.GetNXX(); i++) {
264       {                                        << 
265          gxx[i] = fCurStat.GetXX(i);              260          gxx[i] = fCurStat.GetXX(i);
266          distance[i] = fCurStat.GetDistance(i)    261          distance[i] = fCurStat.GetDistance(i);
267          areacode[i] = fCurStat.GetAreacode(i)    262          areacode[i] = fCurStat.GetAreacode(i);
268       }                                           263       }
269       return fCurStat.GetNXX();                   264       return fCurStat.GetNXX();
270    }                                           << 265    } else {
271    else  // initialize                         << 266       // initialize
272    {                                           << 267       G4int i;
273       for (auto i=0; i<2; ++i)                 << 268       for (i=0; i<2; i++) {
274       {                                        << 
275          distance[i] = kInfinity;                 269          distance[i] = kInfinity;
276          areacode[i] = sOutside;                  270          areacode[i] = sOutside;
277          gxx[i].set(kInfinity, kInfinity, kInf    271          gxx[i].set(kInfinity, kInfinity, kInfinity);
278       }                                           272       }
279    }                                              273    }
280                                                   274    
281    G4ThreeVector p = ComputeLocalPoint(gp);       275    G4ThreeVector p = ComputeLocalPoint(gp);
282    G4ThreeVector xx;                              276    G4ThreeVector xx;
283                                                   277 
284    // The plane is placed on origin with makin    278    // The plane is placed on origin with making its normal 
285    // parallel to z-axis.                         279    // parallel to z-axis. 
286    if (std::fabs(p.z()) <= 0.5 * kCarTolerance    280    if (std::fabs(p.z()) <= 0.5 * kCarTolerance)
287    {   // if p is on the plane, return 1          281    {   // if p is on the plane, return 1
288       distance[0] = 0;                            282       distance[0] = 0;
289       xx = p;                                     283       xx = p;
290    }                                           << 284    } else {
291    else                                        << 
292    {                                           << 
293       distance[0] = std::fabs(p.z());             285       distance[0] = std::fabs(p.z());
294       xx.set(p.x(), p.y(), 0);                    286       xx.set(p.x(), p.y(), 0);  
295    }                                              287    }
296                                                   288 
297    gxx[0] = ComputeGlobalPoint(xx);               289    gxx[0] = ComputeGlobalPoint(xx);
298    areacode[0] = sInside;                         290    areacode[0] = sInside;
299    G4bool isvalid = true;                         291    G4bool isvalid = true;
300    fCurStat.SetCurrentStatus(0, gxx[0], distan    292    fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
301                              isvalid, 1, kDont    293                              isvalid, 1, kDontValidate, &gp);
302    return 1;                                      294    return 1;
303                                                   295 
304 }                                                 296 }
305                                                   297 
306 //============================================ << 298 G4int G4TwistTrapFlatSide::GetAreaCode(const G4ThreeVector &xx, 
307 //* GetAreaCode() ---------------------------- << 299                                        G4bool withTol)
308                                                << 
309 G4int G4TwistTrapFlatSide::GetAreaCode(const G << 
310                                              G << 
311 {                                                 300 {
312                                                   301 
313   static const G4double ctol = 0.5 * kCarToler    302   static const G4double ctol = 0.5 * kCarTolerance;
314   G4int areacode = sInside;                       303   G4int areacode = sInside;
315                                                   304   
316   if (fAxis[0] == kXAxis && fAxis[1] == kYAxis << 305   if (fAxis[0] == kXAxis && fAxis[1] == kYAxis) {
317   {                                            << 306 
318     G4int yaxis = 1;                              307     G4int yaxis = 1;
319                                                   308     
320     G4double wmax = xAxisMax(xx.y(), fTAlph )     309     G4double wmax = xAxisMax(xx.y(), fTAlph ) ;
321     G4double wmin = -xAxisMax(xx.y(), -fTAlph     310     G4double wmin = -xAxisMax(xx.y(), -fTAlph ) ;
322                                                   311 
323     if (withTol)                               << 312     if (withTol) {
324     {                                          << 313       
325       G4bool isoutside   = false;                 314       G4bool isoutside   = false;
326                                                   315       
327       // test boundary of x-axis                  316       // test boundary of x-axis
328                                                   317       
329       if (xx.x() < wmin + ctol)                << 318       if (xx.x() < wmin + ctol) {
330       {                                        << 
331         areacode |= (sAxis0 & (sAxisX | sAxisM    319         areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary; 
332         if (xx.x() <= wmin - ctol) isoutside =    320         if (xx.x() <= wmin - ctol) isoutside = true;
333                                                   321         
334       }                                        << 322       } else if (xx.x() > wmax - ctol) {
335       else if (xx.x() > wmax - ctol)           << 
336       {                                        << 
337         areacode |= (sAxis0 & (sAxisX | sAxisM    323         areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
338         if (xx.x() >= wmax + ctol)  isoutside     324         if (xx.x() >= wmax + ctol)  isoutside = true;
339       }                                           325       }
340                                                   326       
341       // test boundary of y-axis                  327       // test boundary of y-axis
342                                                   328       
343       if (xx.y() < fAxisMin[yaxis] + ctol)     << 329       if (xx.y() < fAxisMin[yaxis] + ctol) {
344       {                                        << 
345         areacode |= (sAxis1 & (sAxisY | sAxisM    330         areacode |= (sAxis1 & (sAxisY | sAxisMin)); 
346                                                   331         
347         if   ((areacode & sBoundary) != 0) are << 332         if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
348         else                        areacode |    333         else                        areacode |= sBoundary;
349         if (xx.y() <= fAxisMin[yaxis] - ctol)     334         if (xx.y() <= fAxisMin[yaxis] - ctol) isoutside = true;
350                                                   335         
351       }                                        << 336       } else if (xx.y() > fAxisMax[yaxis] - ctol) {
352       else if (xx.y() > fAxisMax[yaxis] - ctol << 
353       {                                        << 
354         areacode |= (sAxis1 & (sAxisY | sAxisM    337         areacode |= (sAxis1 & (sAxisY | sAxisMax));
355                                                   338         
356         if   ((areacode & sBoundary) != 0) are << 339         if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
357         else                        areacode |    340         else                        areacode |= sBoundary; 
358         if (xx.y() >= fAxisMax[yaxis] + ctol)     341         if (xx.y() >= fAxisMax[yaxis] + ctol) isoutside = true;
359       }                                           342       }
360                                                   343       
361       // if isoutside = true, clear inside bit    344       // if isoutside = true, clear inside bit.             
362       // if not on boundary, add axis informat    345       // if not on boundary, add axis information.             
363                                                   346       
364       if (isoutside)                           << 347       if (isoutside) {
365       {                                        << 
366         G4int tmpareacode = areacode & (~sInsi    348         G4int tmpareacode = areacode & (~sInside);
367         areacode = tmpareacode;                   349         areacode = tmpareacode;
368       }                                        << 350       } else if ((areacode & sBoundary) != sBoundary) {
369       else if ((areacode & sBoundary) != sBoun << 
370       {                                        << 
371         areacode |= (sAxis0 & sAxisX) | (sAxis    351         areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisY);
372       }                                           352       }           
373     }                                          << 353       
374     else                                       << 354     } else {
375     {                                          << 355       
376       // boundary of x-axis                       356       // boundary of x-axis
377                                                << 357       
378       if (xx.x() < wmin )                      << 358       if (xx.x() < wmin ) {
379       {                                        << 
380         areacode |= (sAxis0 & (sAxisX | sAxisM    359         areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary;
381       }                                        << 360       } else if (xx.x() > wmax) {
382       else if (xx.x() > wmax)                  << 
383       {                                        << 
384         areacode |= (sAxis0 & (sAxisX | sAxisM    361         areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
385       }                                           362       }
386                                                   363       
387       // boundary of y-axis                       364       // boundary of y-axis
388                                                   365       
389       if (xx.y() < fAxisMin[yaxis])            << 366       if (xx.y() < fAxisMin[yaxis]) {
390       {                                        << 
391         areacode |= (sAxis1 & (sAxisY | sAxisM    367         areacode |= (sAxis1 & (sAxisY | sAxisMin));
392         if   ((areacode & sBoundary) != 0) are << 368         if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
393         else                        areacode |    369         else                        areacode |= sBoundary; 
394                                                   370         
395       }                                        << 371       } else if (xx.y() > fAxisMax[yaxis]) {
396       else if (xx.y() > fAxisMax[yaxis])       << 
397       {                                        << 
398         areacode |= (sAxis1 & (sAxisY | sAxisM    372         areacode |= (sAxis1 & (sAxisY | sAxisMax)) ;
399         if   ((areacode & sBoundary) != 0) are << 373         if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
400         else                        areacode |    374         else                        areacode |= sBoundary; 
401       }                                           375       }
402                                                   376       
403       if ((areacode & sBoundary) != sBoundary) << 377       if ((areacode & sBoundary) != sBoundary) {
404       {                                        << 
405         areacode |= (sAxis0 & sAxisX) | (sAxis    378         areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisY);
406       }                                           379       }           
407     }                                             380     }
408     return areacode;                              381     return areacode;
409   }                                            << 382   } else {
410   else                                         << 
411   {                                            << 
412     G4Exception("G4TwistTrapFlatSide::GetAreaC    383     G4Exception("G4TwistTrapFlatSide::GetAreaCode()",
413                 "GeomSolids0001", FatalExcepti << 384                 "NotImplemented", FatalException,
414                 "Feature NOT implemented !");     385                 "Feature NOT implemented !");
415   }                                               386   }
416                                                   387   
417   return areacode;                                388   return areacode;
418 }                                                 389 }
419                                                   390 
                                                   >> 391 
420 //============================================    392 //=====================================================================
421 //* SetCorners -------------------------------    393 //* SetCorners --------------------------------------------------------
422                                                   394 
423 void G4TwistTrapFlatSide::SetCorners()            395 void G4TwistTrapFlatSide::SetCorners()
424 {                                                 396 {
425    // Set Corner points in local coodinate.       397    // Set Corner points in local coodinate.
426                                                   398 
427    if (fAxis[0] == kXAxis && fAxis[1] == kYAxi << 399    if (fAxis[0] == kXAxis && fAxis[1] == kYAxis) {
428    {                                           << 400    
429      G4double x, y, z;                            401      G4double x, y, z;
430                                                   402       
431      // corner of Axis0min and Axis1min           403      // corner of Axis0min and Axis1min
432      x = -fDx1 + fDy * fTAlph ;                   404      x = -fDx1 + fDy * fTAlph ;
433      y = -fDy ;                                   405      y = -fDy ;
434      z = 0 ;                                      406      z = 0 ;
435      SetCorner(sC0Min1Min, x, y, z);              407      SetCorner(sC0Min1Min, x, y, z);
436                                                   408       
437      // corner of Axis0max and Axis1min           409      // corner of Axis0max and Axis1min
438      x = fDx1 + fDy * fTAlph ;                    410      x = fDx1 + fDy * fTAlph ;
439      y = -fDy ;                                   411      y = -fDy ;
440      z = 0 ;                                      412      z = 0 ;
441      SetCorner(sC0Max1Min, x, y, z);              413      SetCorner(sC0Max1Min, x, y, z);
442                                                   414      
443      // corner of Axis0max and Axis1max           415      // corner of Axis0max and Axis1max
444      x = fDx2 - fDy * fTAlph ;                    416      x = fDx2 - fDy * fTAlph ;
445      y = fDy ;                                    417      y = fDy ;
446      z = 0 ;                                      418      z = 0 ;
447      SetCorner(sC0Max1Max, x, y, z);              419      SetCorner(sC0Max1Max, x, y, z);
448                                                   420      
449      // corner of Axis0min and Axis1max           421      // corner of Axis0min and Axis1max
450      x = -fDx2 - fDy * fTAlph ;                   422      x = -fDx2 - fDy * fTAlph ;
451      y = fDy ;                                    423      y = fDy ;
452      z = 0 ;                                      424      z = 0 ;
453      SetCorner(sC0Min1Max, x, y, z);              425      SetCorner(sC0Min1Max, x, y, z);
454                                                   426      
455    }                                           << 427    } else {
456    else                                        << 428      G4cerr << "ERROR - G4TwistTrapFlatSide::SetCorners()" << G4endl
457    {                                           << 429             << "        fAxis[0] = " << fAxis[0] << G4endl
458      std::ostringstream message;               << 430             << "        fAxis[1] = " << fAxis[1] << G4endl;
459      message << "Feature NOT implemented !" << << 
460              << "        fAxis[0] = " << fAxis << 
461              << "        fAxis[1] = " << fAxis << 
462      G4Exception("G4TwistTrapFlatSide::SetCorn    431      G4Exception("G4TwistTrapFlatSide::SetCorners()",
463                  "GeomSolids0001", FatalExcept << 432                  "NotImplemented", FatalException,
                                                   >> 433                  "Feature NOT implemented !");
464    }                                              434    }
465 }                                                 435 }
466                                                   436 
467 //============================================    437 //=====================================================================
468 //* SetBoundaries() --------------------------    438 //* SetBoundaries() ---------------------------------------------------
469                                                   439 
470 void G4TwistTrapFlatSide::SetBoundaries()         440 void G4TwistTrapFlatSide::SetBoundaries()
471 {                                                 441 {
472    // Set direction-unit vector of phi-boundar    442    // Set direction-unit vector of phi-boundary-lines in local coodinate.
473    // Don't call the function twice.              443    // Don't call the function twice.
474                                                   444    
475   G4ThreeVector direction ;                       445   G4ThreeVector direction ;
476                                                   446 
477   if (fAxis[0] == kXAxis && fAxis[1] == kYAxis << 447   if (fAxis[0] == kXAxis && fAxis[1] == kYAxis) {
478   {                                            << 448    
479     // sAxis0 & sAxisMin                          449     // sAxis0 & sAxisMin
480     direction = - ( GetCorner(sC0Min1Max) - Ge    450     direction = - ( GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min) ) ;
481     direction = direction.unit();                 451     direction = direction.unit();
482     SetBoundary(sAxis0 & (sAxisX | sAxisMin),     452     SetBoundary(sAxis0 & (sAxisX | sAxisMin), direction, 
483                 GetCorner(sC0Min1Max), sAxisY)    453                 GetCorner(sC0Min1Max), sAxisY) ;
484                                                   454     
485     // sAxis0 & sAxisMax                          455     // sAxis0 & sAxisMax
486     direction = GetCorner(sC0Max1Max) - GetCor    456     direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min)  ; // inverse
487     direction = direction.unit();                 457     direction = direction.unit();
488     SetBoundary(sAxis0 & (sAxisX | sAxisMax),     458     SetBoundary(sAxis0 & (sAxisX | sAxisMax), direction, 
489                 GetCorner(sC0Max1Min), sAxisY)    459                 GetCorner(sC0Max1Min), sAxisY);
490                                                   460     
491     // sAxis1 & sAxisMin                          461     // sAxis1 & sAxisMin
492     direction = GetCorner(sC0Max1Min) - GetCor    462     direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
493     direction = direction.unit();                 463     direction = direction.unit();
494     SetBoundary(sAxis1 & (sAxisY | sAxisMin),     464     SetBoundary(sAxis1 & (sAxisY | sAxisMin), direction, 
495                 GetCorner(sC0Min1Min), sAxisX)    465                 GetCorner(sC0Min1Min), sAxisX);
496                                                   466     
497     // sAxis1 & sAxisMax                          467     // sAxis1 & sAxisMax
498     direction = - ( GetCorner(sC0Max1Max) - Ge    468     direction = - ( GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max) ) ;
499     direction = direction.unit();                 469     direction = direction.unit();
500     SetBoundary(sAxis1 & (sAxisY | sAxisMax),     470     SetBoundary(sAxis1 & (sAxisY | sAxisMax), direction, 
501                 GetCorner(sC0Max1Max), sAxisX)    471                 GetCorner(sC0Max1Max), sAxisX);
502                                                   472     
503   }                                            << 473   } else {
504   else                                         << 474     G4cerr << "ERROR - G4TwistTrapFlatSide::SetBoundaries()" << G4endl
505   {                                            << 475            << "        fAxis[0] = " << fAxis[0] << G4endl
506     std::ostringstream message;                << 476            << "        fAxis[1] = " << fAxis[1] << G4endl;
507     message << "Feature NOT implemented !" <<  << 
508             << "        fAxis[0] = " << fAxis[ << 
509             << "        fAxis[1] = " << fAxis[ << 
510     G4Exception("G4TwistTrapFlatSide::SetCorne    477     G4Exception("G4TwistTrapFlatSide::SetCorners()",
511                 "GeomSolids0001", FatalExcepti << 478                 "NotImplemented", FatalException,
512   }                                            << 479                 "Feature NOT implemented !");
                                                   >> 480    }
513 }                                                 481 }
514                                                   482 
515 //============================================    483 //=====================================================================
516 //* GetFacets() ------------------------------    484 //* GetFacets() -------------------------------------------------------
517                                                   485 
518 void G4TwistTrapFlatSide::GetFacets( G4int k,  << 486 void G4TwistTrapFlatSide::GetFacets( G4int m, G4int n, G4double xyz[][3],
519                                      G4int fac    487                                      G4int faces[][4], G4int iside ) 
520 {                                                 488 {
                                                   >> 489 
521   G4double x,y    ;     // the two parameters     490   G4double x,y    ;     // the two parameters for the surface equation
522   G4ThreeVector p ;  // a point on the surface    491   G4ThreeVector p ;  // a point on the surface, given by (z,u)
523                                                   492 
524   G4int nnode ;                                   493   G4int nnode ;
525   G4int nface ;                                   494   G4int nface ;
526                                                   495 
527   G4double xmin,xmax ;                            496   G4double xmin,xmax ;
528                                                   497 
529   // calculate the (n-1)*(k-1) vertices        << 498   // calculate the (n-1)*(m-1) vertices
                                                   >> 499 
                                                   >> 500   G4int i,j ;
                                                   >> 501 
                                                   >> 502   for ( i = 0 ; i<n ; i++ ) {
530                                                   503 
531   for ( G4int i = 0 ; i<n ; ++i )              << 
532   {                                            << 
533     y = -fDy + i*(2*fDy)/(n-1) ;                  504     y = -fDy + i*(2*fDy)/(n-1) ;
534                                                   505 
535     for ( G4int j = 0 ; j<k ; ++j )            << 506     for ( j = 0 ; j<m ; j++ ) {
536     {                                          << 507 
537       xmin = GetBoundaryMin(y) ;                  508       xmin = GetBoundaryMin(y) ;
538       xmax = GetBoundaryMax(y) ;                  509       xmax = GetBoundaryMax(y) ;
539       x = xmin + j*(xmax-xmin)/(k-1) ;         << 510       x = xmin + j*(xmax-xmin)/(m-1) ;
540                                                   511 
541       nnode = GetNode(i,j,k,n,iside) ;         << 512       nnode = GetNode(i,j,m,n,iside) ;
542       p = SurfacePoint(x,y,true) ;  // surface    513       p = SurfacePoint(x,y,true) ;  // surface point in global coordinate system
543                                                   514 
544       xyz[nnode][0] = p.x() ;                     515       xyz[nnode][0] = p.x() ;
545       xyz[nnode][1] = p.y() ;                     516       xyz[nnode][1] = p.y() ;
546       xyz[nnode][2] = p.z() ;                     517       xyz[nnode][2] = p.z() ;
547                                                   518 
548       if ( i<n-1 && j<k-1 )                    << 519       if ( i<n-1 && j<m-1 ) {   
549       {                                        << 520 
550         nface = GetFace(i,j,k,n,iside) ;       << 521         nface = GetFace(i,j,m,n,iside) ;
551                                                << 522 
552         if (fHandedness < 0)  // lower side    << 523         if (fHandedness < 0) {  // lower side 
553         {                                      << 524           faces[nface][0] = GetEdgeVisibility(i,j,m,n,0,1) * ( GetNode(i  ,j  ,m,n,iside)+1) ;  
554           faces[nface][0] = GetEdgeVisibility( << 525           faces[nface][1] = GetEdgeVisibility(i,j,m,n,1,1) * ( GetNode(i+1,j  ,m,n,iside)+1) ;
555                           * ( GetNode(i  ,j  , << 526           faces[nface][2] = GetEdgeVisibility(i,j,m,n,2,1) * ( GetNode(i+1,j+1,m,n,iside)+1) ;
556           faces[nface][1] = GetEdgeVisibility( << 527           faces[nface][3] = GetEdgeVisibility(i,j,m,n,3,1) * ( GetNode(i  ,j+1,m,n,iside)+1) ;
557                           * ( GetNode(i+1,j  , << 528         } else {                // upper side
558           faces[nface][2] = GetEdgeVisibility( << 529           faces[nface][0] = GetEdgeVisibility(i,j,m,n,0,-1) * ( GetNode(i  ,j  ,m,n,iside)+1) ;  
559                           * ( GetNode(i+1,j+1, << 530           faces[nface][1] = GetEdgeVisibility(i,j,m,n,1,-1) * ( GetNode(i  ,j+1,m,n,iside)+1) ;
560           faces[nface][3] = GetEdgeVisibility( << 531           faces[nface][2] = GetEdgeVisibility(i,j,m,n,2,-1) * ( GetNode(i+1,j+1,m,n,iside)+1) ;
561                           * ( GetNode(i  ,j+1, << 532           faces[nface][3] = GetEdgeVisibility(i,j,m,n,3,-1) * ( GetNode(i+1,j  ,m,n,iside)+1) ;
562         }                                      << 
563         else                  // upper side    << 
564         {                                      << 
565           faces[nface][0] = GetEdgeVisibility( << 
566                           * ( GetNode(i  ,j  , << 
567           faces[nface][1] = GetEdgeVisibility( << 
568                           * ( GetNode(i  ,j+1, << 
569           faces[nface][2] = GetEdgeVisibility( << 
570                           * ( GetNode(i+1,j+1, << 
571           faces[nface][3] = GetEdgeVisibility( << 
572                           * ( GetNode(i+1,j  , << 
573         }                                         533         }
                                                   >> 534 
574       }                                           535       }
575     }                                             536     }
576   }                                               537   }
577 }                                                 538 }
578                                                   539