Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/CSG/src/G4Para.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/CSG/src/G4Para.cc (Version 11.3.0) and /geometry/solids/CSG/src/G4Para.cc (Version 9.4)


  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 //
                                                   >>  27 // $Id: G4Para.cc,v 1.43 2010/10/19 15:42:10 gcosmo Exp $
                                                   >>  28 // GEANT4 tag $Name: geant4-09-04 $
                                                   >>  29 //
                                                   >>  30 // class G4Para
                                                   >>  31 //
 26 // Implementation for G4Para class                 32 // Implementation for G4Para class
 27 //                                                 33 //
 28 // 21.03.95 P.Kent: Modified for `tolerant' ge <<  34 // History:
                                                   >>  35 //
                                                   >>  36 // 23.10.05 V.Grichine: bug fixed in DistanceToOut(p,v,...) for the v.x()<0 case 
                                                   >>  37 // 28.04.05 V.Grichine: new SurfaceNormal according to J. Apostolakis proposal 
                                                   >>  38 // 30.11.04 V.Grichine: modifications in SurfaceNormal for edges/vertices and
                                                   >>  39 //                      in constructor with vertices
                                                   >>  40 // 14.02.02 V.Grichine: bug fixed in Inside according to proposal of D.Wright
                                                   >>  41 // 18.11.99 V.Grichine: kUndef was added to ESide
 29 // 31.10.96 V.Grichine: Modifications accordin     42 // 31.10.96 V.Grichine: Modifications according G4Box/Tubs before to commit
 30 // 28.04.05 V.Grichine: new SurfaceNormal acco <<  43 // 21.03.95 P.Kent: Modified for `tolerant' geom
 31 // 29.05.17 E.Tcherniaev: complete revision, s <<  44 //
 32 //////////////////////////////////////////////     45 ////////////////////////////////////////////////////////////////////////////
 33                                                    46 
 34 #include "G4Para.hh"                               47 #include "G4Para.hh"
 35                                                    48 
 36 #if !defined(G4GEOM_USE_UPARA)                 << 
 37                                                << 
 38 #include "G4VoxelLimits.hh"                        49 #include "G4VoxelLimits.hh"
 39 #include "G4AffineTransform.hh"                    50 #include "G4AffineTransform.hh"
 40 #include "G4BoundingEnvelope.hh"               << 
 41 #include "Randomize.hh"                            51 #include "Randomize.hh"
 42                                                    52 
 43 #include "G4VPVParameterisation.hh"                53 #include "G4VPVParameterisation.hh"
 44                                                    54 
 45 #include "G4VGraphicsScene.hh"                     55 #include "G4VGraphicsScene.hh"
                                                   >>  56 #include "G4Polyhedron.hh"
                                                   >>  57 #include "G4NURBS.hh"
                                                   >>  58 #include "G4NURBSbox.hh"
 46                                                    59 
 47 using namespace CLHEP;                             60 using namespace CLHEP;
 48                                                    61 
 49 ////////////////////////////////////////////// <<  62 // Private enum: Not for external use 
                                                   >>  63     
                                                   >>  64 enum ESide {kUndef,kPX,kMX,kPY,kMY,kPZ,kMZ};
                                                   >>  65 
                                                   >>  66 // used internally for normal routine
                                                   >>  67 
                                                   >>  68 enum ENSide {kNZ,kNX,kNY};
                                                   >>  69 
                                                   >>  70 /////////////////////////////////////////////////////////////////////
                                                   >>  71 //
                                                   >>  72 // Constructor - check and set half-widths
                                                   >>  73 
                                                   >>  74 void G4Para::SetAllParameters( G4double pDx, G4double pDy, G4double pDz, 
                                                   >>  75                                G4double pAlpha, G4double pTheta, G4double pPhi )
                                                   >>  76 {
                                                   >>  77   if ( pDx > 0 && pDy > 0 && pDz > 0 )
                                                   >>  78   {
                                                   >>  79     fDx         = pDx;
                                                   >>  80     fDy         = pDy;
                                                   >>  81     fDz         = pDz;
                                                   >>  82     fTalpha     = std::tan(pAlpha);
                                                   >>  83     fTthetaCphi = std::tan(pTheta)*std::cos(pPhi);
                                                   >>  84     fTthetaSphi = std::tan(pTheta)*std::sin(pPhi);
                                                   >>  85   }
                                                   >>  86   else
                                                   >>  87   {
                                                   >>  88     G4cerr << "ERROR - G4Para()::SetAllParameters(): " << GetName() << G4endl
                                                   >>  89            << "        Invalid dimensions ! - "
                                                   >>  90            << pDx << ", " << pDy << ", " << pDz << G4endl;
                                                   >>  91     G4Exception("G4Para::SetAllParameters()", "InvalidSetup",
                                                   >>  92                 FatalException, "Invalid Length Parameters.");
                                                   >>  93   }
                                                   >>  94   fCubicVolume = 0.;
                                                   >>  95   fSurfaceArea = 0.;
                                                   >>  96   fpPolyhedron = 0;
                                                   >>  97 }
                                                   >>  98 
                                                   >>  99 ///////////////////////////////////////////////////////////////////////////
 50 //                                                100 //
 51 //  Constructor - set & check half widths      << 
 52                                                   101 
 53 G4Para::G4Para(const G4String& pName,             102 G4Para::G4Para(const G4String& pName,
 54                      G4double pDx, G4double pD    103                      G4double pDx, G4double pDy, G4double pDz,
 55                      G4double pAlpha, G4double    104                      G4double pAlpha, G4double pTheta, G4double pPhi)
 56   : G4CSGSolid(pName), halfCarTolerance(0.5*kC << 105   : G4CSGSolid(pName)
 57 {                                                 106 {
 58   SetAllParameters(pDx, pDy, pDz, pAlpha, pThe << 107   if ((pDx<=0) || (pDy<=0) || (pDz<=0))
 59   fRebuildPolyhedron = false;  // default valu << 108   {
                                                   >> 109     G4cerr << "ERROR - G4Para()::G4Para(): " << GetName() << G4endl
                                                   >> 110            << "        Invalid dimensions ! - "
                                                   >> 111            << pDx << ", " << pDy << ", " << pDz << G4endl;
                                                   >> 112     G4Exception("G4Para::G4Para()", "InvalidSetup",
                                                   >> 113                 FatalException, "Invalid Length Parameters.");
                                                   >> 114   }
                                                   >> 115   SetAllParameters( pDx, pDy, pDz, pAlpha, pTheta, pPhi);
 60 }                                                 116 }
 61                                                   117 
 62 ////////////////////////////////////////////// << 118 ////////////////////////////////////////////////////////////////////////
 63 //                                                119 //
 64 // Constructor - design of trapezoid based on  << 120 // Constructor - Design of trapezoid based on 8 G4ThreeVector parameters, 
                                                   >> 121 // which are its vertices. Checking of planarity with preparation of 
                                                   >> 122 // fPlanes[] and than calculation of other members
 65                                                   123 
 66 G4Para::G4Para( const G4String& pName,            124 G4Para::G4Para( const G4String& pName,
 67                 const G4ThreeVector pt[8] )       125                 const G4ThreeVector pt[8] )
 68   : G4CSGSolid(pName), halfCarTolerance(0.5*kC << 126   : G4CSGSolid(pName)
 69 {                                                 127 {
 70   // Find dimensions and trigonometric values  << 128   if (!( pt[0].z()<0 && pt[0].z()==pt[1].z() && pt[0].z()==pt[2].z() &&
 71   //                                           << 129          pt[0].z()==pt[3].z() && pt[4].z()>0 && pt[4].z()==pt[5].z() &&
 72   fDx = (pt[3].x() - pt[2].x())*0.5;           << 130          pt[4].z()==pt[6].z() && pt[4].z()==pt[7].z()           &&
 73   fDy = (pt[2].y() - pt[1].y())*0.5;           << 131         (pt[0].z()+pt[4].z())==0                                &&
 74   fDz = pt[7].z();                             << 132          pt[0].y()==pt[1].y() && pt[2].y()==pt[3].y()           &&
 75   CheckParameters(); // check dimensions       << 133          pt[4].y()==pt[5].y() && pt[6].y()==pt[7].y()           &&
 76                                                << 134        ( pt[0].y() + pt[2].y() + pt[4].y() + pt[6].y() ) == 0   && 
 77   fTalpha = (pt[2].x() + pt[3].x() - pt[1].x() << 135        ( pt[0].x() + pt[1].x() + pt[4].x() + pt[5].x() ) == 0) )
 78   fTthetaCphi = (pt[4].x() + fDy*fTalpha + fDx << 136   {
 79   fTthetaSphi = (pt[4].y() + fDy)/fDz;         << 137     G4cerr << "ERROR - G4Para()::G4Para(): " << GetName() << G4endl
 80   MakePlanes();                                << 138            << "        Invalid dimensions !" << G4endl;
 81                                                << 139     G4Exception("G4Para::G4Para()", "InvalidSetup",
 82   // Recompute vertices                        << 140                 FatalException, "Invalid vertice coordinates.");
 83   //                                           << 141   }    
 84   G4ThreeVector v[8];                          << 142   fDx = ((pt[3]).x()-(pt[2]).x())*0.5;
 85   G4double DyTalpha = fDy*fTalpha;             << 143   fDy = ((pt[2]).y()-(pt[1]).y())*0.5;
 86   G4double DzTthetaSphi = fDz*fTthetaSphi;     << 144   fDz = (pt[7]).z();
 87   G4double DzTthetaCphi = fDz*fTthetaCphi;     << 145   fTalpha = ((pt[2]).x()+(pt[3]).x()-(pt[1]).x()-(pt[0]).x())*0.25/fDy ;
 88   v[0].set(-DzTthetaCphi-DyTalpha-fDx, -DzTthe << 146   fTthetaCphi = ((pt[4]).x()+fDy*fTalpha+fDx)/fDz ;
 89   v[1].set(-DzTthetaCphi-DyTalpha+fDx, -DzTthe << 147   fTthetaSphi = ((pt[4]).y()+fDy)/fDz ;
 90   v[2].set(-DzTthetaCphi+DyTalpha-fDx, -DzTthe << 148   fCubicVolume = 0.;
 91   v[3].set(-DzTthetaCphi+DyTalpha+fDx, -DzTthe << 149   fSurfaceArea = 0.;
 92   v[4].set( DzTthetaCphi-DyTalpha-fDx,  DzTthe << 150   fpPolyhedron = 0;
 93   v[5].set( DzTthetaCphi-DyTalpha+fDx,  DzTthe << 
 94   v[6].set( DzTthetaCphi+DyTalpha-fDx,  DzTthe << 
 95   v[7].set( DzTthetaCphi+DyTalpha+fDx,  DzTthe << 
 96                                                << 
 97   // Compare with original vertices            << 
 98   //                                           << 
 99   for (G4int i=0; i<8; ++i)                    << 
100   {                                            << 
101     G4double delx = std::abs(pt[i].x() - v[i]. << 
102     G4double dely = std::abs(pt[i].y() - v[i]. << 
103     G4double delz = std::abs(pt[i].z() - v[i]. << 
104     G4double discrepancy = std::max(std::max(d << 
105     if (discrepancy > 0.1*kCarTolerance)       << 
106     {                                          << 
107       std::ostringstream message;              << 
108       G4long oldprc = message.precision(16);   << 
109       message << "Invalid vertice coordinates  << 
110               << "\nVertix #" << i << ", discr << 
111               << "\n  original   : " << pt[i]  << 
112               << "\n  recomputed : " << v[i];  << 
113       G4cout.precision(oldprc);                << 
114       G4Exception("G4Para::G4Para()", "GeomSol << 
115                   FatalException, message);    << 
116                                                << 
117     }                                          << 
118   }                                            << 
119 }                                                 151 }
120                                                   152 
121 ////////////////////////////////////////////// << 153 ///////////////////////////////////////////////////////////////////////
122 //                                                154 //
123 // Fake default constructor - sets only member    155 // Fake default constructor - sets only member data and allocates memory
124 //                            for usage restri << 156 //                            for usage restricted to object persistency.
125                                                << 157 //
126 G4Para::G4Para( __void__& a )                     158 G4Para::G4Para( __void__& a )
127   : G4CSGSolid(a), halfCarTolerance(0.5*kCarTo << 159   : G4CSGSolid(a), fDx(0.), fDy(0.), fDz(0.),
                                                   >> 160     fTalpha(0.), fTthetaCphi(0.), fTthetaSphi(0.)
128 {                                                 161 {
129   SetAllParameters(1., 1., 1., 0., 0., 0.);    << 
130   fRebuildPolyhedron = false; // default value << 
131 }                                                 162 }
132                                                   163 
133 //////////////////////////////////////////////    164 //////////////////////////////////////////////////////////////////////////
134 //                                                165 //
135 // Destructor                                  << 
136                                                   166 
137 G4Para::~G4Para() = default;                   << 167 G4Para::~G4Para()
                                                   >> 168 {
                                                   >> 169 }
138                                                   170 
139 //////////////////////////////////////////////    171 //////////////////////////////////////////////////////////////////////////
140 //                                                172 //
141 // Copy constructor                               173 // Copy constructor
142                                                   174 
143 G4Para::G4Para(const G4Para& rhs)                 175 G4Para::G4Para(const G4Para& rhs)
144   : G4CSGSolid(rhs), halfCarTolerance(rhs.half << 176   : G4CSGSolid(rhs), fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz),
145     fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz),  << 177     fTalpha(rhs.fTalpha), fTthetaCphi(rhs.fTthetaCphi),
146     fTthetaCphi(rhs.fTthetaCphi),fTthetaSphi(r << 178     fTthetaSphi(rhs.fTthetaSphi)
147 {                                                 179 {
148   for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs << 
149 }                                                 180 }
150                                                   181 
151 //////////////////////////////////////////////    182 //////////////////////////////////////////////////////////////////////////
152 //                                                183 //
153 // Assignment operator                            184 // Assignment operator
154                                                   185 
155 G4Para& G4Para::operator = (const G4Para& rhs) << 186 G4Para& G4Para::operator = (const G4Para& rhs) 
156 {                                                 187 {
157    // Check assignment to self                    188    // Check assignment to self
158    //                                             189    //
159    if (this == &rhs)  { return *this; }           190    if (this == &rhs)  { return *this; }
160                                                   191 
161    // Copy base class data                        192    // Copy base class data
162    //                                             193    //
163    G4CSGSolid::operator=(rhs);                    194    G4CSGSolid::operator=(rhs);
164                                                   195 
165    // Copy data                                   196    // Copy data
166    //                                             197    //
167    halfCarTolerance = rhs.halfCarTolerance;    << 198    fDx = rhs.fDx; fDy = rhs.fDy; fDz = rhs.fDz;
168    fDx = rhs.fDx;                              << 199    fTalpha = rhs.fTalpha; fTthetaCphi = rhs.fTthetaCphi;
169    fDy = rhs.fDy;                              << 
170    fDz = rhs.fDz;                              << 
171    fTalpha = rhs.fTalpha;                      << 
172    fTthetaCphi = rhs.fTthetaCphi;              << 
173    fTthetaSphi = rhs.fTthetaSphi;                 200    fTthetaSphi = rhs.fTthetaSphi;
174    for (G4int i=0; i<4; ++i) { fPlanes[i] = rh << 
175                                                   201 
176    return *this;                                  202    return *this;
177 }                                                 203 }
178                                                   204 
179 //////////////////////////////////////////////    205 //////////////////////////////////////////////////////////////////////////
180 //                                                206 //
181 // Set all parameters, as for constructor - se << 
182                                                << 
183 void G4Para::SetAllParameters(G4double pDx, G4 << 
184                               G4double pAlpha, << 
185 {                                              << 
186   // Reset data of the base class              << 
187   fCubicVolume = 0;                            << 
188   fSurfaceArea = 0;                            << 
189   fRebuildPolyhedron = true;                   << 
190                                                << 
191   // Set parameters                            << 
192   fDx = pDx;                                   << 
193   fDy = pDy;                                   << 
194   fDz = pDz;                                   << 
195   fTalpha = std::tan(pAlpha);                  << 
196   fTthetaCphi = std::tan(pTheta)*std::cos(pPhi << 
197   fTthetaSphi = std::tan(pTheta)*std::sin(pPhi << 
198                                                << 
199   CheckParameters();                           << 
200   MakePlanes();                                << 
201 }                                              << 
202                                                << 
203 ////////////////////////////////////////////// << 
204 //                                             << 
205 // Check dimensions                            << 
206                                                << 
207 void G4Para::CheckParameters()                 << 
208 {                                              << 
209   if (fDx < 2*kCarTolerance ||                 << 
210       fDy < 2*kCarTolerance ||                 << 
211       fDz < 2*kCarTolerance)                   << 
212   {                                            << 
213     std::ostringstream message;                << 
214     message << "Invalid (too small or negative << 
215             << GetName()                       << 
216             << "\n  X - " << fDx               << 
217             << "\n  Y - " << fDy               << 
218             << "\n  Z - " << fDz;              << 
219     G4Exception("G4Para::CheckParameters()", " << 
220                 FatalException, message);      << 
221   }                                            << 
222 }                                              << 
223                                                << 
224 ////////////////////////////////////////////// << 
225 //                                             << 
226 // Set side planes                             << 
227                                                << 
228 void G4Para::MakePlanes()                      << 
229 {                                              << 
230   G4ThreeVector vx(1, 0, 0);                   << 
231   G4ThreeVector vy(fTalpha, 1, 0);             << 
232   G4ThreeVector vz(fTthetaCphi, fTthetaSphi, 1 << 
233                                                << 
234   // Set -Y & +Y planes                        << 
235   //                                           << 
236   G4ThreeVector ynorm = (vx.cross(vz)).unit(); << 
237                                                << 
238   fPlanes[0].a = 0.;                           << 
239   fPlanes[0].b = ynorm.y();                    << 
240   fPlanes[0].c = ynorm.z();                    << 
241   fPlanes[0].d = fPlanes[0].b*fDy; // point (0 << 
242                                                << 
243   fPlanes[1].a =  0.;                          << 
244   fPlanes[1].b = -fPlanes[0].b;                << 
245   fPlanes[1].c = -fPlanes[0].c;                << 
246   fPlanes[1].d =  fPlanes[0].d;                << 
247                                                << 
248   // Set -X & +X planes                        << 
249   //                                           << 
250   G4ThreeVector xnorm = (vz.cross(vy)).unit(); << 
251                                                << 
252   fPlanes[2].a = xnorm.x();                    << 
253   fPlanes[2].b = xnorm.y();                    << 
254   fPlanes[2].c = xnorm.z();                    << 
255   fPlanes[2].d = fPlanes[2].a*fDx; // point (f << 
256                                                << 
257   fPlanes[3].a = -fPlanes[2].a;                << 
258   fPlanes[3].b = -fPlanes[2].b;                << 
259   fPlanes[3].c = -fPlanes[2].c;                << 
260   fPlanes[3].d =  fPlanes[2].d;                << 
261 }                                              << 
262                                                << 
263 ////////////////////////////////////////////// << 
264 //                                             << 
265 // Get volume                                  << 
266                                                << 
267 G4double G4Para::GetCubicVolume()              << 
268 {                                              << 
269   // It is like G4Box, since para transformati << 
270   if (fCubicVolume == 0)                       << 
271   {                                            << 
272     fCubicVolume = 8*fDx*fDy*fDz;              << 
273   }                                            << 
274   return fCubicVolume;                         << 
275 }                                              << 
276                                                << 
277 ////////////////////////////////////////////// << 
278 //                                             << 
279 // Get surface area                            << 
280                                                << 
281 G4double G4Para::GetSurfaceArea()              << 
282 {                                              << 
283   if(fSurfaceArea == 0)                        << 
284   {                                            << 
285     G4ThreeVector vx(fDx, 0, 0);               << 
286     G4ThreeVector vy(fDy*fTalpha, fDy, 0);     << 
287     G4ThreeVector vz(fDz*fTthetaCphi, fDz*fTth << 
288                                                << 
289     G4double sxy = fDx*fDy; // (vx.cross(vy)). << 
290     G4double sxz = (vx.cross(vz)).mag();       << 
291     G4double syz = (vy.cross(vz)).mag();       << 
292                                                << 
293     fSurfaceArea = 8*(sxy+sxz+syz);            << 
294   }                                            << 
295   return fSurfaceArea;                         << 
296 }                                              << 
297                                                << 
298 ////////////////////////////////////////////// << 
299 //                                             << 
300 // Dispatch to parameterisation for replicatio    207 // Dispatch to parameterisation for replication mechanism dimension
301 // computation & modification                  << 208 // computation & modification.
302                                                   209 
303 void G4Para::ComputeDimensions(      G4VPVPara    210 void G4Para::ComputeDimensions(      G4VPVParameterisation* p,
304                                 const G4int n,    211                                 const G4int n,
305                                 const G4VPhysi    212                                 const G4VPhysicalVolume* pRep )
306 {                                                 213 {
307   p->ComputeDimensions(*this,n,pRep);             214   p->ComputeDimensions(*this,n,pRep);
308 }                                                 215 }
309                                                   216 
310 ////////////////////////////////////////////// << 
311 //                                             << 
312 // Get bounding box                            << 
313                                                << 
314 void G4Para::BoundingLimits(G4ThreeVector& pMi << 
315 {                                              << 
316   G4double dz = GetZHalfLength();              << 
317   G4double dx = GetXHalfLength();              << 
318   G4double dy = GetYHalfLength();              << 
319                                                << 
320   G4double x0 = dz*fTthetaCphi;                << 
321   G4double x1 = dy*GetTanAlpha();              << 
322   G4double xmin =                              << 
323     std::min(                                  << 
324     std::min(                                  << 
325     std::min(-x0-x1-dx,-x0+x1-dx),x0-x1-dx),x0 << 
326   G4double xmax =                              << 
327     std::max(                                  << 
328     std::max(                                  << 
329     std::max(-x0-x1+dx,-x0+x1+dx),x0-x1+dx),x0 << 
330                                                << 
331   G4double y0 = dz*fTthetaSphi;                << 
332   G4double ymin = std::min(-y0-dy,y0-dy);      << 
333   G4double ymax = std::max(-y0+dy,y0+dy);      << 
334                                                << 
335   pMin.set(xmin,ymin,-dz);                     << 
336   pMax.set(xmax,ymax, dz);                     << 
337                                                << 
338   // Check correctness of the bounding box     << 
339   //                                           << 
340   if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 
341   {                                            << 
342     std::ostringstream message;                << 
343     message << "Bad bounding box (min >= max)  << 
344             << GetName() << " !"               << 
345             << "\npMin = " << pMin             << 
346             << "\npMax = " << pMax;            << 
347     G4Exception("G4Para::BoundingLimits()", "G << 
348                 JustWarning, message);         << 
349     DumpInfo();                                << 
350   }                                            << 
351 }                                              << 
352                                                   217 
353 ////////////////////////////////////////////// << 218 //////////////////////////////////////////////////////////////
354 //                                                219 //
355 // Calculate extent under transform and specif    220 // Calculate extent under transform and specified limit
356                                                   221 
357 G4bool G4Para::CalculateExtent( const EAxis pA    222 G4bool G4Para::CalculateExtent( const EAxis pAxis,
358                                 const G4VoxelL    223                                 const G4VoxelLimits& pVoxelLimit,
359                                 const G4Affine    224                                 const G4AffineTransform& pTransform,
360                                      G4double&    225                                      G4double& pMin, G4double& pMax ) const
361 {                                                 226 {
362   G4ThreeVector bmin, bmax;                    << 227   G4bool flag;
363   G4bool exist;                                << 
364                                                   228 
365   // Check bounding box (bbox)                 << 229   if (!pTransform.IsRotated())
366   //                                           << 230   {  
367   BoundingLimits(bmin,bmax);                   << 231     // Special case handling for unrotated trapezoids
368   G4BoundingEnvelope bbox(bmin,bmax);          << 232     // Compute z/x/y/ mins and maxs respecting limits, with early returns
369 #ifdef G4BBOX_EXTENT                           << 233     // if outside limits. Then switch() on pAxis
370   return bbox.CalculateExtent(pAxis,pVoxelLimi << 234 
371 #endif                                         << 235     G4int i ; 
372   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 236     G4double xoffset,xMin,xMax;
373   {                                            << 237     G4double yoffset,yMin,yMax;
374     return exist = pMin < pMax;                << 238     G4double zoffset,zMin,zMax;
375   }                                            << 239     G4double temp[8] ;       // some points for intersection with zMin/zMax
                                                   >> 240     
                                                   >> 241     xoffset=pTransform.NetTranslation().x();      
                                                   >> 242     yoffset=pTransform.NetTranslation().y();
                                                   >> 243     zoffset=pTransform.NetTranslation().z();
                                                   >> 244  
                                                   >> 245     G4ThreeVector pt[8];   // vertices after translation
                                                   >> 246     pt[0]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy*fTalpha-fDx,
                                                   >> 247                         yoffset-fDz*fTthetaSphi-fDy,zoffset-fDz);
                                                   >> 248     pt[1]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy*fTalpha+fDx,
                                                   >> 249                         yoffset-fDz*fTthetaSphi-fDy,zoffset-fDz);
                                                   >> 250     pt[2]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy*fTalpha-fDx,
                                                   >> 251                         yoffset-fDz*fTthetaSphi+fDy,zoffset-fDz);
                                                   >> 252     pt[3]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy*fTalpha+fDx,
                                                   >> 253                         yoffset-fDz*fTthetaSphi+fDy,zoffset-fDz);
                                                   >> 254     pt[4]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy*fTalpha-fDx,
                                                   >> 255                         yoffset+fDz*fTthetaSphi-fDy,zoffset+fDz);
                                                   >> 256     pt[5]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy*fTalpha+fDx,
                                                   >> 257                         yoffset+fDz*fTthetaSphi-fDy,zoffset+fDz);
                                                   >> 258     pt[6]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy*fTalpha-fDx,
                                                   >> 259                         yoffset+fDz*fTthetaSphi+fDy,zoffset+fDz);
                                                   >> 260     pt[7]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy*fTalpha+fDx,
                                                   >> 261                         yoffset+fDz*fTthetaSphi+fDy,zoffset+fDz);
                                                   >> 262     zMin=zoffset-fDz;
                                                   >> 263     zMax=zoffset+fDz;
                                                   >> 264     if ( pVoxelLimit.IsZLimited() )
                                                   >> 265     {
                                                   >> 266       if   ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
                                                   >> 267           || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
                                                   >> 268       {
                                                   >> 269         return false;
                                                   >> 270       }
                                                   >> 271       else
                                                   >> 272       {
                                                   >> 273         if (zMin<pVoxelLimit.GetMinZExtent())
                                                   >> 274         {
                                                   >> 275           zMin=pVoxelLimit.GetMinZExtent();
                                                   >> 276         }
                                                   >> 277         if (zMax>pVoxelLimit.GetMaxZExtent())
                                                   >> 278         {
                                                   >> 279           zMax=pVoxelLimit.GetMaxZExtent();
                                                   >> 280         }
                                                   >> 281       }
                                                   >> 282     }
376                                                   283 
377   // Set bounding envelope (benv) and calculat << 284     temp[0] = pt[0].y()+(pt[4].y()-pt[0].y())
378   //                                           << 285                        *(zMin-pt[0].z())/(pt[4].z()-pt[0].z()) ;
379   G4double dz = GetZHalfLength();              << 286     temp[1] = pt[0].y()+(pt[4].y()-pt[0].y())
380   G4double dx = GetXHalfLength();              << 287                        *(zMax-pt[0].z())/(pt[4].z()-pt[0].z()) ;
381   G4double dy = GetYHalfLength();              << 288     temp[2] = pt[2].y()+(pt[6].y()-pt[2].y())
                                                   >> 289                        *(zMin-pt[2].z())/(pt[6].z()-pt[2].z()) ;
                                                   >> 290     temp[3] = pt[2].y()+(pt[6].y()-pt[2].y())
                                                   >> 291                        *(zMax-pt[2].z())/(pt[6].z()-pt[2].z()) ;        
                                                   >> 292     yMax = yoffset - std::fabs(fDz*fTthetaSphi) - fDy - fDy ;
                                                   >> 293     yMin = -yMax ;
                                                   >> 294     for(i=0;i<4;i++)
                                                   >> 295     {
                                                   >> 296       if(temp[i] > yMax) yMax = temp[i] ;
                                                   >> 297       if(temp[i] < yMin) yMin = temp[i] ;
                                                   >> 298     }
                                                   >> 299       
                                                   >> 300     if (pVoxelLimit.IsYLimited())
                                                   >> 301     {
                                                   >> 302       if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
                                                   >> 303         || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
                                                   >> 304       {
                                                   >> 305         return false;
                                                   >> 306       }
                                                   >> 307       else
                                                   >> 308       {
                                                   >> 309         if (yMin<pVoxelLimit.GetMinYExtent())
                                                   >> 310         {
                                                   >> 311           yMin=pVoxelLimit.GetMinYExtent();
                                                   >> 312         }
                                                   >> 313         if (yMax>pVoxelLimit.GetMaxYExtent())
                                                   >> 314         {
                                                   >> 315           yMax=pVoxelLimit.GetMaxYExtent();
                                                   >> 316         }
                                                   >> 317       }
                                                   >> 318     }
382                                                   319 
383   G4double x0 = dz*fTthetaCphi;                << 320     temp[0] = pt[0].x()+(pt[4].x()-pt[0].x())
384   G4double x1 = dy*GetTanAlpha();              << 321                        *(zMin-pt[0].z())/(pt[4].z()-pt[0].z()) ;
385   G4double y0 = dz*fTthetaSphi;                << 322     temp[1] = pt[0].x()+(pt[4].x()-pt[0].x())
                                                   >> 323                        *(zMax-pt[0].z())/(pt[4].z()-pt[0].z()) ;
                                                   >> 324     temp[2] = pt[2].x()+(pt[6].x()-pt[2].x())
                                                   >> 325                        *(zMin-pt[2].z())/(pt[6].z()-pt[2].z()) ;
                                                   >> 326     temp[3] = pt[2].x()+(pt[6].x()-pt[2].x())
                                                   >> 327                        *(zMax-pt[2].z())/(pt[6].z()-pt[2].z()) ;
                                                   >> 328     temp[4] = pt[3].x()+(pt[7].x()-pt[3].x())
                                                   >> 329                        *(zMin-pt[3].z())/(pt[7].z()-pt[3].z()) ;
                                                   >> 330     temp[5] = pt[3].x()+(pt[7].x()-pt[3].x())
                                                   >> 331                        *(zMax-pt[3].z())/(pt[7].z()-pt[3].z()) ;
                                                   >> 332     temp[6] = pt[1].x()+(pt[5].x()-pt[1].x())
                                                   >> 333                        *(zMin-pt[1].z())/(pt[5].z()-pt[1].z()) ;
                                                   >> 334     temp[7] = pt[1].x()+(pt[5].x()-pt[1].x())
                                                   >> 335                        *(zMax-pt[1].z())/(pt[5].z()-pt[1].z()) ;
                                                   >> 336 
                                                   >> 337     xMax = xoffset - std::fabs(fDz*fTthetaCphi) - fDx - fDx -fDx - fDx;
                                                   >> 338     xMin = -xMax ;
                                                   >> 339     for(i=0;i<8;i++)
                                                   >> 340     {
                                                   >> 341       if(temp[i] > xMax) xMax = temp[i] ;
                                                   >> 342       if(temp[i] < xMin) xMin = temp[i] ;
                                                   >> 343     }
                                                   >> 344       // xMax/Min = f(yMax/Min) ?
                                                   >> 345     if (pVoxelLimit.IsXLimited())
                                                   >> 346     {
                                                   >> 347       if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
                                                   >> 348         || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
                                                   >> 349       {
                                                   >> 350         return false;
                                                   >> 351       }
                                                   >> 352       else
                                                   >> 353       {
                                                   >> 354         if (xMin<pVoxelLimit.GetMinXExtent())
                                                   >> 355         {
                                                   >> 356           xMin=pVoxelLimit.GetMinXExtent();
                                                   >> 357         }
                                                   >> 358         if (xMax>pVoxelLimit.GetMaxXExtent())
                                                   >> 359         {
                                                   >> 360           xMax=pVoxelLimit.GetMaxXExtent();
                                                   >> 361         }
                                                   >> 362       }
                                                   >> 363     }
                                                   >> 364 
                                                   >> 365     switch (pAxis)
                                                   >> 366     {
                                                   >> 367       case kXAxis:
                                                   >> 368         pMin=xMin;
                                                   >> 369         pMax=xMax;
                                                   >> 370         break;
                                                   >> 371       case kYAxis:
                                                   >> 372         pMin=yMin;
                                                   >> 373         pMax=yMax;
                                                   >> 374         break;
                                                   >> 375       case kZAxis:
                                                   >> 376         pMin=zMin;
                                                   >> 377         pMax=zMax;
                                                   >> 378         break;
                                                   >> 379       default:
                                                   >> 380         break;
                                                   >> 381     }
386                                                   382 
387   G4ThreeVectorList baseA(4), baseB(4);        << 383     pMin-=kCarTolerance;
388   baseA[0].set(-x0-x1-dx,-y0-dy,-dz);          << 384     pMax+=kCarTolerance;
389   baseA[1].set(-x0-x1+dx,-y0-dy,-dz);          << 385     flag = true;
390   baseA[2].set(-x0+x1+dx,-y0+dy,-dz);          << 386   }
391   baseA[3].set(-x0+x1-dx,-y0+dy,-dz);          << 387   else
                                                   >> 388   {
                                                   >> 389     // General rotated case - create and clip mesh to boundaries
392                                                   390 
393   baseB[0].set(+x0-x1-dx, y0-dy, dz);          << 391     G4bool existsAfterClip=false;
394   baseB[1].set(+x0-x1+dx, y0-dy, dz);          << 392     G4ThreeVectorList *vertices;
395   baseB[2].set(+x0+x1+dx, y0+dy, dz);          << 
396   baseB[3].set(+x0+x1-dx, y0+dy, dz);          << 
397                                                   393 
398   std::vector<const G4ThreeVectorList *> polyg << 394     pMin=+kInfinity;
399   polygons[0] = &baseA;                        << 395     pMax=-kInfinity;
400   polygons[1] = &baseB;                        << 
401                                                   396 
402   G4BoundingEnvelope benv(bmin,bmax,polygons); << 397     // Calculate rotated vertex coordinates
403   exist = benv.CalculateExtent(pAxis,pVoxelLim << 398 
404   return exist;                                << 399     vertices=CreateRotatedVertices(pTransform);
                                                   >> 400     ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
                                                   >> 401     ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax);
                                                   >> 402     ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
                                                   >> 403       
                                                   >> 404     if (pMin!=kInfinity||pMax!=-kInfinity)
                                                   >> 405     {
                                                   >> 406       existsAfterClip=true;
                                                   >> 407         
                                                   >> 408       // Add 2*tolerance to avoid precision troubles
                                                   >> 409       //
                                                   >> 410       pMin-=kCarTolerance;
                                                   >> 411       pMax+=kCarTolerance;
                                                   >> 412     }
                                                   >> 413     else
                                                   >> 414     {
                                                   >> 415       // Check for case where completely enveloping clipping volume
                                                   >> 416       // If point inside then we are confident that the solid completely
                                                   >> 417       // envelopes the clipping volume. Hence set min/max extents according
                                                   >> 418       // to clipping volume extents along the specified axis.
                                                   >> 419        
                                                   >> 420       G4ThreeVector clipCentre(
                                                   >> 421         (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
                                                   >> 422         (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
                                                   >> 423         (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
                                                   >> 424         
                                                   >> 425       if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
                                                   >> 426       {
                                                   >> 427         existsAfterClip=true;
                                                   >> 428         pMin=pVoxelLimit.GetMinExtent(pAxis);
                                                   >> 429         pMax=pVoxelLimit.GetMaxExtent(pAxis);
                                                   >> 430       }
                                                   >> 431     }
                                                   >> 432     delete vertices ;          //  'new' in the function called
                                                   >> 433     flag = existsAfterClip ;
                                                   >> 434   }
                                                   >> 435   return flag;
405 }                                                 436 }
406                                                   437 
407 ////////////////////////////////////////////// << 438 /////////////////////////////////////////////////////////////////////////////
408 //                                             << 
409 // Determine where is point p, inside/on_surfa << 
410 //                                                439 //
                                                   >> 440 // Check in p is inside/on surface/outside solid
411                                                   441 
412 EInside G4Para::Inside( const G4ThreeVector& p    442 EInside G4Para::Inside( const G4ThreeVector& p ) const
413 {                                                 443 {
414   G4double xx = fPlanes[2].a*p.x()+fPlanes[2]. << 444   G4double xt, yt, yt1;
415   G4double dx = std::abs(xx) + fPlanes[2].d;   << 445   EInside  in = kOutside;
                                                   >> 446 
                                                   >> 447   yt1 = p.y() - fTthetaSphi*p.z();
                                                   >> 448   yt  = std::fabs(yt1) ;
416                                                   449 
417   G4double yy = fPlanes[0].b*p.y()+fPlanes[0]. << 450   // xt = std::fabs( p.x() - fTthetaCphi*p.z() - fTalpha*yt );
418   G4double dy = std::abs(yy) + fPlanes[0].d;   << 
419   G4double dxy = std::max(dx,dy);              << 
420                                                   451 
421   G4double dz = std::abs(p.z())-fDz;           << 452   xt = std::fabs( p.x() - fTthetaCphi*p.z() - fTalpha*yt1 );
422   G4double dist = std::max(dxy,dz);            << 
423                                                   453 
424   if (dist > halfCarTolerance) return kOutside << 454   if ( std::fabs( p.z() ) <= fDz - kCarTolerance*0.5)
425   return (dist > -halfCarTolerance) ? kSurface << 455   {
                                                   >> 456     if (yt <= fDy - kCarTolerance*0.5)
                                                   >> 457     {
                                                   >> 458       if      ( xt <= fDx - kCarTolerance*0.5 ) in = kInside;
                                                   >> 459       else if ( xt <= fDx + kCarTolerance*0.5 ) in = kSurface;
                                                   >> 460     }
                                                   >> 461     else if ( yt <= fDy + kCarTolerance*0.5)
                                                   >> 462     {
                                                   >> 463       if ( xt <= fDx + kCarTolerance*0.5 ) in = kSurface;  
                                                   >> 464     }
                                                   >> 465   }
                                                   >> 466   else  if ( std::fabs(p.z()) <= fDz + kCarTolerance*0.5 )
                                                   >> 467   {
                                                   >> 468     if ( yt <= fDy + kCarTolerance*0.5)
                                                   >> 469     {
                                                   >> 470       if ( xt <= fDx + kCarTolerance*0.5 ) in = kSurface;  
                                                   >> 471     }
                                                   >> 472   }
                                                   >> 473   return in;
426 }                                                 474 }
427                                                   475 
428 ////////////////////////////////////////////// << 476 ///////////////////////////////////////////////////////////////////////////
429 //                                                477 //
430 // Determine side where point is, and return c << 478 // Calculate side nearest to p, and return normal
                                                   >> 479 // If 2+ sides equidistant, first side's normal returned (arbitrarily)
431                                                   480 
432 G4ThreeVector G4Para::SurfaceNormal( const G4T    481 G4ThreeVector G4Para::SurfaceNormal( const G4ThreeVector& p ) const
433 {                                                 482 {
434   G4int nsurf = 0; // number of surfaces where << 483   G4ThreeVector norm, sumnorm(0.,0.,0.);
                                                   >> 484   G4int noSurfaces = 0; 
                                                   >> 485   G4double distx,disty,distz;
                                                   >> 486   G4double newpx,newpy,xshift;
                                                   >> 487   G4double calpha,salpha;      // Sin/Cos(alpha) - needed to recalc G4Parameter
                                                   >> 488   G4double tntheta,cosntheta;  // tan and cos of normal's theta component
                                                   >> 489   G4double ycomp;
                                                   >> 490   G4double delta = 0.5*kCarTolerance;
                                                   >> 491 
                                                   >> 492   newpx  = p.x()-fTthetaCphi*p.z();
                                                   >> 493   newpy  = p.y()-fTthetaSphi*p.z();
                                                   >> 494 
                                                   >> 495   calpha = 1/std::sqrt(1+fTalpha*fTalpha);
                                                   >> 496   if (fTalpha) {salpha = -calpha/fTalpha;} // NOTE: using MINUS std::sin(alpha)
                                                   >> 497   else         {salpha = 0.;}
                                                   >> 498   
                                                   >> 499   //  xshift = newpx*calpha+newpy*salpha;
                                                   >> 500   xshift = newpx - newpy*fTalpha;
435                                                   501 
436   // Check Z faces                             << 502   //  distx  = std::fabs(std::fabs(xshift)-fDx*calpha);
437   //                                           << 503   distx  = std::fabs(std::fabs(xshift)-fDx);
438   G4double nz = 0;                             << 504   disty  = std::fabs(std::fabs(newpy)-fDy);
439   G4double dz = std::abs(p.z()) - fDz;         << 505   distz  = std::fabs(std::fabs(p.z())-fDz);
440   if (std::abs(dz) <= halfCarTolerance)        << 
441   {                                            << 
442     nz = (p.z() < 0) ? -1 : 1;                 << 
443     ++nsurf;                                   << 
444   }                                            << 
445                                                   506 
446   // Check Y faces                             << 507   tntheta   = fTthetaCphi*calpha + fTthetaSphi*salpha;
447   //                                           << 508   cosntheta = 1/std::sqrt(1+tntheta*tntheta);
448   G4double ny = 0;                             << 509   ycomp     = 1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
449   G4double yy = fPlanes[0].b*p.y()+fPlanes[0]. << 510 
450   if (std::abs(fPlanes[0].d + yy) <= halfCarTo << 511   G4ThreeVector nX  = G4ThreeVector( calpha*cosntheta,
451   {                                            << 512                                      salpha*cosntheta,
452     ny  = fPlanes[0].b;                        << 513                                     -tntheta*cosntheta);
453     nz += fPlanes[0].c;                        << 514   G4ThreeVector nY  = G4ThreeVector( 0, ycomp,-fTthetaSphi*ycomp);
454     ++nsurf;                                   << 515   G4ThreeVector nZ  = G4ThreeVector( 0, 0,  1.0);
455   }                                            << 516 
456   else if (std::abs(fPlanes[1].d - yy) <= half << 517   if (distx <= delta)      
457   {                                               518   {
458     ny  = fPlanes[1].b;                        << 519     noSurfaces ++;
459     nz += fPlanes[1].c;                        << 520     if ( xshift >= 0.) {sumnorm += nX;}
460     ++nsurf;                                   << 521     else               {sumnorm -= nX;}
461   }                                               522   }
462                                                << 523   if (disty <= delta)
463   // Check X faces                             << 
464   //                                           << 
465   G4double nx = 0;                             << 
466   G4double xx = fPlanes[2].a*p.x()+fPlanes[2]. << 
467   if (std::abs(fPlanes[2].d + xx) <= halfCarTo << 
468   {                                               524   {
469     nx  = fPlanes[2].a;                        << 525     noSurfaces ++;
470     ny += fPlanes[2].b;                        << 526     if ( newpy >= 0.)  {sumnorm += nY;}
471     nz += fPlanes[2].c;                        << 527     else               {sumnorm -= nY;}
472     ++nsurf;                                   << 
473   }                                               528   }
474   else if (std::abs(fPlanes[3].d - xx) <= half << 529   if (distz <= delta)  
475   {                                               530   {
476     nx  = fPlanes[3].a;                        << 531     noSurfaces ++;
477     ny += fPlanes[3].b;                        << 532     if ( p.z() >= 0.)  {sumnorm += nZ;}
478     nz += fPlanes[3].c;                        << 533     else               {sumnorm -= nZ;}
479     ++nsurf;                                   << 
480   }                                               534   }
481                                                << 535   if ( noSurfaces == 0 )
482   // Return normal                             << 
483   //                                           << 
484   if (nsurf == 1)      return {nx,ny,nz};      << 
485   else if (nsurf != 0) return G4ThreeVector(nx << 
486   else                                         << 
487   {                                               536   {
488     // Point is not on the surface             << 
489     //                                         << 
490 #ifdef G4CSGDEBUG                                 537 #ifdef G4CSGDEBUG
491     std::ostringstream message;                << 538     G4Exception("G4Para::SurfaceNormal(p)", "Notification", JustWarning, 
492     G4int oldprc = message.precision(16);      << 539                 "Point p is not on surface !?" );
493     message << "Point p is not on surface (!?) << 540 #endif 
494             << GetName() << G4endl;            << 541      norm = ApproxSurfaceNormal(p);
495     message << "Position:\n";                  << 
496     message << "   p.x() = " << p.x()/mm << "  << 
497     message << "   p.y() = " << p.y()/mm << "  << 
498     message << "   p.z() = " << p.z()/mm << "  << 
499     G4cout.precision(oldprc) ;                 << 
500     G4Exception("G4Para::SurfaceNormal(p)", "G << 
501                 JustWarning, message );        << 
502     DumpInfo();                                << 
503 #endif                                         << 
504     return ApproxSurfaceNormal(p);             << 
505   }                                               542   }
                                                   >> 543   else if ( noSurfaces == 1 ) {norm = sumnorm;}
                                                   >> 544   else                        {norm = sumnorm.unit();}
                                                   >> 545 
                                                   >> 546   return norm;
506 }                                                 547 }
507                                                   548 
508 ////////////////////////////////////////////// << 549 
                                                   >> 550 ////////////////////////////////////////////////////////////////////////
509 //                                                551 //
510 // Algorithm for SurfaceNormal() following the    552 // Algorithm for SurfaceNormal() following the original specification
511 // for points not on the surface                  553 // for points not on the surface
512                                                   554 
513 G4ThreeVector G4Para::ApproxSurfaceNormal( con    555 G4ThreeVector G4Para::ApproxSurfaceNormal( const G4ThreeVector& p ) const
514 {                                                 556 {
515   G4double dist = -DBL_MAX;                    << 557   ENSide  side;
516   G4int iside = 0;                             << 558   G4ThreeVector norm;
517   for (G4int i=0; i<4; ++i)                    << 559   G4double distx,disty,distz;
                                                   >> 560   G4double newpx,newpy,xshift;
                                                   >> 561   G4double calpha,salpha;  // Sin/Cos(alpha) - needed to recalc G4Parameter 
                                                   >> 562   G4double tntheta,cosntheta;  // tan and cos of normal's theta component
                                                   >> 563   G4double ycomp;
                                                   >> 564 
                                                   >> 565   newpx=p.x()-fTthetaCphi*p.z();
                                                   >> 566   newpy=p.y()-fTthetaSphi*p.z();
                                                   >> 567 
                                                   >> 568   calpha=1/std::sqrt(1+fTalpha*fTalpha);
                                                   >> 569   if (fTalpha)
518   {                                               570   {
519     G4double d = fPlanes[i].a*p.x() +          << 571     salpha=-calpha/fTalpha;  // NOTE: actually use MINUS std::sin(alpha)
520                  fPlanes[i].b*p.y() +          << 572   }
521                  fPlanes[i].c*p.z() + fPlanes[ << 573   else
522     if (d > dist) { dist = d; iside = i; }     << 574   {
                                                   >> 575     salpha=0;
523   }                                               576   }
524                                                   577 
525   G4double distz = std::abs(p.z()) - fDz;      << 578   xshift=newpx*calpha+newpy*salpha;
526   if (dist > distz)                            << 579 
527     return { fPlanes[iside].a, fPlanes[iside]. << 580   distx=std::fabs(std::fabs(xshift)-fDx*calpha);
                                                   >> 581   disty=std::fabs(std::fabs(newpy)-fDy);
                                                   >> 582   distz=std::fabs(std::fabs(p.z())-fDz);
                                                   >> 583     
                                                   >> 584   if (distx<disty)
                                                   >> 585   {
                                                   >> 586     if (distx<distz) {side=kNX;}
                                                   >> 587     else {side=kNZ;}
                                                   >> 588   }
528   else                                            589   else
529     return { 0, 0, (G4double)((p.z() < 0) ? -1 << 590   {
                                                   >> 591     if (disty<distz) {side=kNY;}
                                                   >> 592     else {side=kNZ;}
                                                   >> 593   }
                                                   >> 594 
                                                   >> 595   switch (side)
                                                   >> 596   {
                                                   >> 597     case kNX:
                                                   >> 598       tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
                                                   >> 599       if (xshift<0)
                                                   >> 600       {
                                                   >> 601         cosntheta=-1/std::sqrt(1+tntheta*tntheta);
                                                   >> 602       }
                                                   >> 603       else
                                                   >> 604       {
                                                   >> 605         cosntheta=1/std::sqrt(1+tntheta*tntheta);
                                                   >> 606       }
                                                   >> 607       norm=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
                                                   >> 608       break;
                                                   >> 609     case kNY:
                                                   >> 610       if (newpy<0)
                                                   >> 611       {
                                                   >> 612         ycomp=-1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
                                                   >> 613       }
                                                   >> 614       else
                                                   >> 615       {
                                                   >> 616         ycomp=1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
                                                   >> 617       }
                                                   >> 618       norm=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
                                                   >> 619       break;
                                                   >> 620     case kNZ:           // Closest to Z
                                                   >> 621       if (p.z()>=0)
                                                   >> 622       {
                                                   >> 623         norm=G4ThreeVector(0,0,1);
                                                   >> 624       }
                                                   >> 625       else
                                                   >> 626       {
                                                   >> 627         norm=G4ThreeVector(0,0,-1);
                                                   >> 628       }
                                                   >> 629       break;
                                                   >> 630   }
                                                   >> 631   return norm;
530 }                                                 632 }
531                                                   633 
532 ////////////////////////////////////////////// << 634 //////////////////////////////////////////////////////////////////////////////
533 //                                                635 //
534 // Calculate distance to shape from outside       636 // Calculate distance to shape from outside
535 //  - return kInfinity if no intersection      << 637 // - return kInfinity if no intersection
536                                                << 638 //
537 G4double G4Para::DistanceToIn(const G4ThreeVec << 639 // ALGORITHM:
538                               const G4ThreeVec << 640 // For each component, calculate pair of minimum and maximum intersection
539 {                                              << 641 // values for which the particle is in the extent of the shape
540   // Z intersections                           << 642 // - The smallest (MAX minimum) allowed distance of the pairs is intersect
                                                   >> 643 // - Z plane intersectin uses tolerance
                                                   >> 644 // - XZ YZ planes use logic & *SLIGHTLY INCORRECT* tolerance
                                                   >> 645 //   (this saves at least 1 sqrt, 1 multiply and 1 divide... in applicable
                                                   >> 646 //    cases)
                                                   >> 647 // - Note: XZ and YZ planes each divide space into four regions,
                                                   >> 648 //   characterised by ss1 ss2
                                                   >> 649 
                                                   >> 650 G4double G4Para::DistanceToIn( const G4ThreeVector& p,
                                                   >> 651                                const G4ThreeVector& v ) const
                                                   >> 652 {
                                                   >> 653   G4double snxt;    // snxt = default return value
                                                   >> 654   G4double smin,smax;
                                                   >> 655   G4double tmin,tmax;
                                                   >> 656   G4double yt,vy,xt,vx;
                                                   >> 657   G4double max;
541   //                                              658   //
542   if ((std::abs(p.z()) - fDz) >= -halfCarToler << 659   // Z Intersection range
543     return kInfinity;                          << 
544   G4double invz = (-v.z() == 0) ? DBL_MAX : -1 << 
545   G4double dz = (invz < 0) ? fDz : -fDz;       << 
546   G4double tzmin = (p.z() + dz)*invz;          << 
547   G4double tzmax = (p.z() - dz)*invz;          << 
548                                                << 
549   // Y intersections                           << 
550   //                                              660   //
551   G4double tmin0 = tzmin, tmax0 = tzmax;       << 661   if (v.z()>0)
552   G4double cos0 = fPlanes[0].b*v.y() + fPlanes << 662   {
553   G4double disy = fPlanes[0].b*p.y() + fPlanes << 663     max=fDz-p.z();
554   G4double dis0 = fPlanes[0].d + disy;         << 664     if (max>kCarTolerance*0.5)
555   if (dis0 >= -halfCarTolerance)               << 665     {
                                                   >> 666       smax=max/v.z();
                                                   >> 667       smin=(-fDz-p.z())/v.z();
                                                   >> 668     }
                                                   >> 669     else
                                                   >> 670     {
                                                   >> 671       return snxt=kInfinity;
                                                   >> 672     }
                                                   >> 673   }
                                                   >> 674   else if (v.z()<0)
556   {                                               675   {
557     if (cos0 >= 0) return kInfinity;           << 676     max=-fDz-p.z();
558     G4double tmp  = -dis0/cos0;                << 677     if (max<-kCarTolerance*0.5)
559     if (tmin0 < tmp) tmin0 = tmp;              << 678     {
                                                   >> 679       smax=max/v.z();
                                                   >> 680       smin=(fDz-p.z())/v.z();
                                                   >> 681     }
                                                   >> 682     else
                                                   >> 683     {
                                                   >> 684       return snxt=kInfinity;
                                                   >> 685     }
560   }                                               686   }
561   else if (cos0 > 0)                           << 687   else
562   {                                               688   {
563     G4double tmp  = -dis0/cos0;                << 689     if (std::fabs(p.z())<=fDz) // Inside
564     if (tmax0 > tmp) tmax0 = tmp;              << 690     {
                                                   >> 691       smin=0;
                                                   >> 692       smax=kInfinity;
                                                   >> 693     }
                                                   >> 694     else
                                                   >> 695     {
                                                   >> 696       return snxt=kInfinity;
                                                   >> 697     }
565   }                                               698   }
                                                   >> 699     
                                                   >> 700   //
                                                   >> 701   // Y G4Parallel planes intersection
                                                   >> 702   //
                                                   >> 703 
                                                   >> 704   yt=p.y()-fTthetaSphi*p.z();
                                                   >> 705   vy=v.y()-fTthetaSphi*v.z();
566                                                   706 
567   G4double tmin1 = tmin0, tmax1 = tmax0;       << 707   if (vy>0)
568   G4double cos1 = -cos0;                       << 
569   G4double dis1 = fPlanes[1].d - disy;         << 
570   if (dis1 >= -halfCarTolerance)               << 
571   {                                               708   {
572     if (cos1 >= 0) return kInfinity;           << 709     max=fDy-yt;
573     G4double tmp  = -dis1/cos1;                << 710     if (max>kCarTolerance*0.5)
574     if (tmin1 < tmp) tmin1 = tmp;              << 711     {
                                                   >> 712       tmax=max/vy;
                                                   >> 713       tmin=(-fDy-yt)/vy;
                                                   >> 714     }
                                                   >> 715     else
                                                   >> 716     {
                                                   >> 717       return snxt=kInfinity;
                                                   >> 718     }
                                                   >> 719   }
                                                   >> 720   else if (vy<0)
                                                   >> 721   {
                                                   >> 722     max=-fDy-yt;
                                                   >> 723     if (max<-kCarTolerance*0.5)
                                                   >> 724     {
                                                   >> 725       tmax=max/vy;
                                                   >> 726       tmin=(fDy-yt)/vy;
                                                   >> 727     }
                                                   >> 728     else
                                                   >> 729     {
                                                   >> 730       return snxt=kInfinity;
                                                   >> 731     }
575   }                                               732   }
576   else if (cos1 > 0)                           << 733   else
577   {                                               734   {
578     G4double tmp  = -dis1/cos1;                << 735     if (std::fabs(yt)<=fDy)
579     if (tmax1 > tmp) tmax1 = tmp;              << 736     {
                                                   >> 737       tmin=0;
                                                   >> 738       tmax=kInfinity;
                                                   >> 739     }
                                                   >> 740     else
                                                   >> 741     {
                                                   >> 742       return snxt=kInfinity;
                                                   >> 743     }
580   }                                               744   }
581                                                   745 
582   // X intersections                           << 746   // Re-Calc valid intersection range
583   //                                              747   //
584   G4double tmin2 = tmin1, tmax2 = tmax1;       << 748   if (tmin>smin) smin=tmin;
585   G4double cos2 = fPlanes[2].a*v.x() + fPlanes << 749   if (tmax<smax) smax=tmax;
586   G4double disx = fPlanes[2].a*p.x() + fPlanes << 750   if (smax<=smin)
587   G4double dis2 = fPlanes[2].d + disx;         << 
588   if (dis2 >= -halfCarTolerance)               << 
589   {                                               751   {
590     if (cos2 >= 0) return kInfinity;           << 752     return snxt=kInfinity;
591     G4double tmp  = -dis2/cos2;                << 
592     if (tmin2 < tmp) tmin2 = tmp;              << 
593   }                                               753   }
594   else if (cos2 > 0)                           << 754   else
595   {                                               755   {
596     G4double tmp  = -dis2/cos2;                << 756     //
597     if (tmax2 > tmp) tmax2 = tmp;              << 757     // X G4Parallel planes intersection
                                                   >> 758     //
                                                   >> 759     xt=p.x()-fTthetaCphi*p.z()-fTalpha*yt;
                                                   >> 760     vx=v.x()-fTthetaCphi*v.z()-fTalpha*vy;
                                                   >> 761     if (vx>0)
                                                   >> 762     {
                                                   >> 763       max=fDx-xt;
                                                   >> 764       if (max>kCarTolerance*0.5)
                                                   >> 765       {
                                                   >> 766         tmax=max/vx;
                                                   >> 767         tmin=(-fDx-xt)/vx;
                                                   >> 768       }
                                                   >> 769       else
                                                   >> 770       {
                                                   >> 771         return snxt=kInfinity;
                                                   >> 772       }
                                                   >> 773     }
                                                   >> 774     else if (vx<0)
                                                   >> 775     {
                                                   >> 776       max=-fDx-xt;
                                                   >> 777       if (max<-kCarTolerance*0.5)
                                                   >> 778       {
                                                   >> 779         tmax=max/vx;
                                                   >> 780         tmin=(fDx-xt)/vx;
                                                   >> 781       }
                                                   >> 782       else
                                                   >> 783       {
                                                   >> 784         return snxt=kInfinity;
                                                   >> 785       }
                                                   >> 786     }
                                                   >> 787     else
                                                   >> 788     {
                                                   >> 789       if (std::fabs(xt)<=fDx)
                                                   >> 790       {
                                                   >> 791         tmin=0;
                                                   >> 792         tmax=kInfinity;
                                                   >> 793       }
                                                   >> 794       else
                                                   >> 795       {
                                                   >> 796         return snxt=kInfinity;
                                                   >> 797       }
                                                   >> 798     }
                                                   >> 799     if (tmin>smin) smin=tmin;
                                                   >> 800     if (tmax<smax) smax=tmax;
598   }                                               801   }
599                                                   802 
600   G4double tmin3 = tmin2, tmax3 = tmax2;       << 803   if (smax>0&&smin<smax)
601   G4double cos3 = -cos2;                       << 
602   G4double dis3 = fPlanes[3].d - disx;         << 
603   if (dis3 >= -halfCarTolerance)               << 
604   {                                               804   {
605     if (cos3 >= 0) return kInfinity;           << 805     if (smin>0)
606     G4double tmp  = -dis3/cos3;                << 806     {
607     if (tmin3 < tmp) tmin3 = tmp;              << 807       snxt=smin;
                                                   >> 808     }
                                                   >> 809     else
                                                   >> 810     {
                                                   >> 811       snxt=0;
                                                   >> 812     }
608   }                                               813   }
609   else if (cos3 > 0)                           << 814   else
610   {                                               815   {
611     G4double tmp  = -dis3/cos3;                << 816     snxt=kInfinity;
612     if (tmax3 > tmp) tmax3 = tmp;              << 
613   }                                               817   }
614                                                << 818   return snxt;
615   // Find distance                             << 
616   //                                           << 
617   G4double tmin = tmin3, tmax = tmax3;         << 
618   if (tmax <= tmin + halfCarTolerance) return  << 
619   return (tmin < halfCarTolerance ) ? 0. : tmi << 
620 }                                                 819 }
621                                                   820 
622 ////////////////////////////////////////////// << 821 ////////////////////////////////////////////////////////////////////////////
623 //                                                822 //
624 // Calculate exact shortest distance to any bo    823 // Calculate exact shortest distance to any boundary from outside
625 // - returns 0 is point inside                 << 824 // - Returns 0 is point inside
626                                                   825 
627 G4double G4Para::DistanceToIn( const G4ThreeVe    826 G4double G4Para::DistanceToIn( const G4ThreeVector& p ) const
628 {                                                 827 {
629   G4double xx = fPlanes[2].a*p.x()+fPlanes[2]. << 828   G4double safe=0.0;
630   G4double dx = std::abs(xx) + fPlanes[2].d;   << 829   G4double distz1,distz2,disty1,disty2,distx1,distx2;
                                                   >> 830   G4double trany,cosy,tranx,cosx;
631                                                   831 
632   G4double yy = fPlanes[0].b*p.y()+fPlanes[0]. << 832   // Z planes
633   G4double dy = std::abs(yy) + fPlanes[0].d;   << 833   //
634   G4double dxy = std::max(dx,dy);              << 834   distz1=p.z()-fDz;
                                                   >> 835   distz2=-fDz-p.z();
                                                   >> 836   if (distz1>distz2)
                                                   >> 837   {
                                                   >> 838     safe=distz1;
                                                   >> 839   }
                                                   >> 840   else
                                                   >> 841   {
                                                   >> 842     safe=distz2;
                                                   >> 843   }
635                                                   844 
636   G4double dz = std::abs(p.z())-fDz;           << 845   trany=p.y()-fTthetaSphi*p.z(); // Transformed y into `box' system
637   G4double dist = std::max(dxy,dz);            << 846 
                                                   >> 847   // Transformed x into `box' system
                                                   >> 848   //
                                                   >> 849   cosy=1.0/std::sqrt(1.0+fTthetaSphi*fTthetaSphi);
                                                   >> 850   disty1=(trany-fDy)*cosy;
                                                   >> 851   disty2=(-fDy-trany)*cosy;
                                                   >> 852     
                                                   >> 853   if (disty1>safe) safe=disty1;
                                                   >> 854   if (disty2>safe) safe=disty2;
638                                                   855 
639   return (dist > 0) ? dist : 0.;               << 856   tranx=p.x()-fTthetaCphi*p.z()-fTalpha*trany;
                                                   >> 857   cosx=1.0/std::sqrt(1.0+fTalpha*fTalpha+fTthetaCphi*fTthetaCphi);
                                                   >> 858   distx1=(tranx-fDx)*cosx;
                                                   >> 859   distx2=(-fDx-tranx)*cosx;
                                                   >> 860     
                                                   >> 861   if (distx1>safe) safe=distx1;
                                                   >> 862   if (distx2>safe) safe=distx2;
                                                   >> 863     
                                                   >> 864   if (safe<0) safe=0;
                                                   >> 865   return safe;  
640 }                                                 866 }
641                                                   867 
642 //////////////////////////////////////////////    868 //////////////////////////////////////////////////////////////////////////
643 //                                                869 //
644 // Calculate distance to surface of shape from << 870 // Calculate distance to surface of shape from inside
645 // find normal at exit point                   << 871 // Calculate distance to x/y/z planes - smallest is exiting distance
646 // - when leaving the surface, return 0        << 
647                                                   872 
648 G4double G4Para::DistanceToOut(const G4ThreeVe    873 G4double G4Para::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
649                                const G4bool ca    874                                const G4bool calcNorm,
650                                      G4bool* v << 875                                G4bool *validNorm, G4ThreeVector *n) const
651 {                                                 876 {
652   // Z intersections                           << 877   ESide side = kUndef;
                                                   >> 878   G4double snxt;    // snxt = return value
                                                   >> 879   G4double max,tmax;
                                                   >> 880   G4double yt,vy,xt,vx;
                                                   >> 881 
                                                   >> 882   G4double ycomp,calpha,salpha,tntheta,cosntheta;
                                                   >> 883 
                                                   >> 884   //
                                                   >> 885   // Z Intersections
653   //                                              886   //
654   if ((std::abs(p.z()) - fDz) >= -halfCarToler << 887 
                                                   >> 888   if (v.z()>0)
655   {                                               889   {
656     if (calcNorm)                              << 890     max=fDz-p.z();
                                                   >> 891     if (max>kCarTolerance*0.5)
657     {                                             892     {
658       *validNorm = true;                       << 893       snxt=max/v.z();
659       n->set(0, 0, (p.z() < 0) ? -1 : 1);      << 894       side=kPZ;
                                                   >> 895     }
                                                   >> 896     else
                                                   >> 897     {
                                                   >> 898       if (calcNorm)
                                                   >> 899       {
                                                   >> 900         *validNorm=true;
                                                   >> 901         *n=G4ThreeVector(0,0,1);
                                                   >> 902       }
                                                   >> 903       return snxt=0;
660     }                                             904     }
661     return 0.;                                 << 
662   }                                               905   }
663   G4double vz = v.z();                         << 906   else if (v.z()<0)
664   G4double tmax = (vz == 0) ? DBL_MAX : (std:: << 
665   G4int iside = (vz < 0) ? -4 : -2; // little  << 
666                                                << 
667   // Y intersections                           << 
668   //                                           << 
669   G4double cos0 = fPlanes[0].b*v.y() + fPlanes << 
670   if (cos0 > 0)                                << 
671   {                                               907   {
672     G4double dis0 = fPlanes[0].b*p.y() + fPlan << 908     max=-fDz-p.z();
673     if (dis0 >= -halfCarTolerance)             << 909     if (max<-kCarTolerance*0.5)
                                                   >> 910     {
                                                   >> 911       snxt=max/v.z();
                                                   >> 912       side=kMZ;
                                                   >> 913     }
                                                   >> 914     else
674     {                                             915     {
675       if (calcNorm)                               916       if (calcNorm)
676       {                                           917       {
677         *validNorm = true;                     << 918         *validNorm=true;
678         n->set(0, fPlanes[0].b, fPlanes[0].c); << 919         *n=G4ThreeVector(0,0,-1);
679       }                                           920       }
680       return 0.;                               << 921       return snxt=0;
681     }                                             922     }
682     G4double tmp = -dis0/cos0;                 << 
683     if (tmax > tmp) { tmax = tmp; iside = 0; } << 
684   }                                               923   }
                                                   >> 924   else
                                                   >> 925   {
                                                   >> 926     snxt=kInfinity;
                                                   >> 927   }
                                                   >> 928     
                                                   >> 929   //
                                                   >> 930   // Y plane intersection
                                                   >> 931   //
                                                   >> 932 
                                                   >> 933   yt=p.y()-fTthetaSphi*p.z();
                                                   >> 934   vy=v.y()-fTthetaSphi*v.z();
685                                                   935 
686   G4double cos1 = -cos0;                       << 936   if (vy>0)
687   if (cos1 > 0)                                << 
688   {                                               937   {
689     G4double dis1 = fPlanes[1].b*p.y() + fPlan << 938     max=fDy-yt;
690     if (dis1 >= -halfCarTolerance)             << 939     if (max>kCarTolerance*0.5)
                                                   >> 940     {
                                                   >> 941       tmax=max/vy;
                                                   >> 942       if (tmax<snxt)
                                                   >> 943       {
                                                   >> 944         snxt=tmax;
                                                   >> 945         side=kPY;
                                                   >> 946       }
                                                   >> 947     }
                                                   >> 948     else
                                                   >> 949     {
                                                   >> 950       if (calcNorm)
                                                   >> 951       {      
                                                   >> 952         *validNorm=true; // Leaving via plus Y
                                                   >> 953         ycomp=1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
                                                   >> 954         *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
                                                   >> 955       }
                                                   >> 956       return snxt=0;
                                                   >> 957     }
                                                   >> 958   }
                                                   >> 959   else if (vy<0)
                                                   >> 960   {
                                                   >> 961     max=-fDy-yt;
                                                   >> 962     if (max<-kCarTolerance*0.5)
                                                   >> 963     {
                                                   >> 964       tmax=max/vy;
                                                   >> 965       if (tmax<snxt)
                                                   >> 966       {
                                                   >> 967         snxt=tmax;
                                                   >> 968         side=kMY;
                                                   >> 969       }
                                                   >> 970     }
                                                   >> 971     else
691     {                                             972     {
692       if (calcNorm)                               973       if (calcNorm)
693       {                                           974       {
694         *validNorm = true;                     << 975         *validNorm=true; // Leaving via minus Y
695         n->set(0, fPlanes[1].b, fPlanes[1].c); << 976         ycomp=-1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
                                                   >> 977         *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
696       }                                           978       }
697       return 0.;                               << 979       return snxt=0;
698     }                                             980     }
699     G4double tmp = -dis1/cos1;                 << 
700     if (tmax > tmp) { tmax = tmp; iside = 1; } << 
701   }                                               981   }
702                                                   982 
703   // X intersections                           << 
704   //                                              983   //
705   G4double cos2 = fPlanes[2].a*v.x() + fPlanes << 984   // X plane intersection
706   if (cos2 > 0)                                << 985   //
                                                   >> 986 
                                                   >> 987   xt=p.x()-fTthetaCphi*p.z()-fTalpha*yt;
                                                   >> 988   vx=v.x()-fTthetaCphi*v.z()-fTalpha*vy;
                                                   >> 989   if (vx>0)
707   {                                               990   {
708     G4double dis2 = fPlanes[2].a*p.x()+fPlanes << 991     max=fDx-xt;
709     if (dis2 >= -halfCarTolerance)             << 992     if (max>kCarTolerance*0.5)
                                                   >> 993     {
                                                   >> 994       tmax=max/vx;
                                                   >> 995       if (tmax<snxt)
                                                   >> 996       {
                                                   >> 997         snxt=tmax;
                                                   >> 998         side=kPX;
                                                   >> 999       }
                                                   >> 1000     }
                                                   >> 1001     else
710     {                                             1002     {
711       if (calcNorm)                               1003       if (calcNorm)
712       {                                           1004       {
713          *validNorm = true;                    << 1005         *validNorm=true; // Leaving via plus X
714          n->set(fPlanes[2].a, fPlanes[2].b, fP << 1006         calpha=1/std::sqrt(1+fTalpha*fTalpha);
                                                   >> 1007         if (fTalpha)
                                                   >> 1008         {
                                                   >> 1009           salpha=-calpha/fTalpha;  // NOTE: actually use MINUS std::sin(alpha)
                                                   >> 1010         }
                                                   >> 1011         else
                                                   >> 1012         {
                                                   >> 1013           salpha=0;
                                                   >> 1014         }
                                                   >> 1015         tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
                                                   >> 1016         cosntheta=1/std::sqrt(1+tntheta*tntheta);
                                                   >> 1017         *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
715       }                                           1018       }
716       return 0.;                               << 1019       return snxt=0;
717     }                                             1020     }
718     G4double tmp = -dis2/cos2;                 << 
719     if (tmax > tmp) { tmax = tmp; iside = 2; } << 
720   }                                               1021   }
721                                                << 1022   else if (vx<0)
722   G4double cos3 = -cos2;                       << 
723   if (cos3 > 0)                                << 
724   {                                               1023   {
725     G4double dis3 = fPlanes[3].a*p.x()+fPlanes << 1024     max=-fDx-xt;
726     if (dis3 >= -halfCarTolerance)             << 1025     if (max<-kCarTolerance*0.5)
                                                   >> 1026     {
                                                   >> 1027       tmax=max/vx;
                                                   >> 1028       if (tmax<snxt)
                                                   >> 1029       {
                                                   >> 1030         snxt=tmax;
                                                   >> 1031         side=kMX;
                                                   >> 1032       }
                                                   >> 1033     }
                                                   >> 1034     else
727     {                                             1035     {
728       if (calcNorm)                               1036       if (calcNorm)
729       {                                           1037       {
730          *validNorm = true;                    << 1038         *validNorm=true; // Leaving via minus X
731          n->set(fPlanes[3].a, fPlanes[3].b, fP << 1039         calpha=1/std::sqrt(1+fTalpha*fTalpha);
                                                   >> 1040         if (fTalpha)
                                                   >> 1041         {
                                                   >> 1042           salpha=-calpha/fTalpha;  // NOTE: actually use MINUS std::sin(alpha)
                                                   >> 1043         }
                                                   >> 1044         else
                                                   >> 1045         {
                                                   >> 1046           salpha=0;
                                                   >> 1047         }
                                                   >> 1048         tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
                                                   >> 1049         cosntheta=-1/std::sqrt(1+tntheta*tntheta);
                                                   >> 1050         *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
732       }                                           1051       }
733       return 0.;                               << 1052       return snxt=0;
734     }                                             1053     }
735     G4double tmp = -dis3/cos3;                 << 
736     if (tmax > tmp) { tmax = tmp; iside = 3; } << 
737   }                                               1054   }
738                                                   1055 
739   // Set normal, if required, and return dista << 1056   if (calcNorm)
740   //                                           << 
741   if (calcNorm)                                << 
742   {                                               1057   {
743     *validNorm = true;                         << 1058     *validNorm=true;
744     if (iside < 0)                             << 1059     switch (side)
745       n->set(0, 0, iside + 3); // (-4+3)=-1, ( << 1060     {
746     else                                       << 1061       case kMZ:
747       n->set(fPlanes[iside].a, fPlanes[iside]. << 1062         *n=G4ThreeVector(0,0,-1);
                                                   >> 1063         break;
                                                   >> 1064       case kPZ:
                                                   >> 1065         *n=G4ThreeVector(0,0,1);
                                                   >> 1066         break;
                                                   >> 1067       case kMY:
                                                   >> 1068         ycomp=-1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
                                                   >> 1069         *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
                                                   >> 1070         break;        
                                                   >> 1071       case kPY:
                                                   >> 1072         ycomp=1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
                                                   >> 1073         *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
                                                   >> 1074         break;        
                                                   >> 1075       case kMX:
                                                   >> 1076         calpha=1/std::sqrt(1+fTalpha*fTalpha);
                                                   >> 1077         if (fTalpha)
                                                   >> 1078         {
                                                   >> 1079           salpha=-calpha/fTalpha;  // NOTE: actually use MINUS std::sin(alpha)
                                                   >> 1080         }
                                                   >> 1081         else
                                                   >> 1082         {
                                                   >> 1083           salpha=0;
                                                   >> 1084         }
                                                   >> 1085         tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
                                                   >> 1086         cosntheta=-1/std::sqrt(1+tntheta*tntheta);
                                                   >> 1087         *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
                                                   >> 1088         break;
                                                   >> 1089       case kPX:
                                                   >> 1090         calpha=1/std::sqrt(1+fTalpha*fTalpha);
                                                   >> 1091         if (fTalpha)
                                                   >> 1092         {
                                                   >> 1093           salpha=-calpha/fTalpha;  // NOTE: actually use MINUS std::sin(alpha)
                                                   >> 1094         }
                                                   >> 1095         else
                                                   >> 1096         {
                                                   >> 1097           salpha=0;
                                                   >> 1098         }
                                                   >> 1099         tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
                                                   >> 1100         cosntheta=1/std::sqrt(1+tntheta*tntheta);
                                                   >> 1101         *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
                                                   >> 1102         break;
                                                   >> 1103       default:
                                                   >> 1104         DumpInfo();
                                                   >> 1105         G4Exception("G4Para::DistanceToOut(p,v,..)","Notification",JustWarning,
                                                   >> 1106                     "Undefined side for valid surface normal to solid.");
                                                   >> 1107         break;
                                                   >> 1108     }
748   }                                               1109   }
749   return tmax;                                 << 1110   return snxt;
750 }                                                 1111 }
751                                                   1112 
752 ////////////////////////////////////////////// << 1113 /////////////////////////////////////////////////////////////////////////////
753 //                                                1114 //
754 // Calculate exact shortest distance to any bo    1115 // Calculate exact shortest distance to any boundary from inside
755 // - returns 0 is point outside                << 1116 // - Returns 0 is point outside
756                                                   1117 
757 G4double G4Para::DistanceToOut( const G4ThreeV    1118 G4double G4Para::DistanceToOut( const G4ThreeVector& p ) const
758 {                                                 1119 {
                                                   >> 1120   G4double safe=0.0;
                                                   >> 1121   G4double distz1,distz2,disty1,disty2,distx1,distx2;
                                                   >> 1122   G4double trany,cosy,tranx,cosx;
                                                   >> 1123 
759 #ifdef G4CSGDEBUG                                 1124 #ifdef G4CSGDEBUG
760   if( Inside(p) == kOutside )                     1125   if( Inside(p) == kOutside )
761   {                                               1126   {
762     std::ostringstream message;                << 1127      G4int oldprc = G4cout.precision(16) ;
763     G4int oldprc = message.precision(16);      << 1128      G4cout << G4endl ;
764     message << "Point p is outside (!?) of sol << 1129      DumpInfo();
765     message << "Position:\n";                  << 1130      G4cout << "Position:"  << G4endl << G4endl ;
766     message << "   p.x() = " << p.x()/mm << "  << 1131      G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
767     message << "   p.y() = " << p.y()/mm << "  << 1132      G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
768     message << "   p.z() = " << p.z()/mm << "  << 1133      G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
769     G4cout.precision(oldprc) ;                 << 1134      G4cout.precision(oldprc) ;
770     G4Exception("G4Para::DistanceToOut(p)", "G << 1135      G4Exception("G4Para::DistanceToOut(p)", "Notification",
771                 JustWarning, message );        << 1136                  JustWarning, "Point p is outside !?" );
772     DumpInfo();                                << 1137   }
773     }                                          << 
774 #endif                                            1138 #endif
775   G4double xx = fPlanes[2].a*p.x()+fPlanes[2]. << 
776   G4double dx = std::abs(xx) + fPlanes[2].d;   << 
777                                                   1139 
778   G4double yy = fPlanes[0].b*p.y()+fPlanes[0]. << 1140   // Z planes
779   G4double dy = std::abs(yy) + fPlanes[0].d;   << 1141   //
780   G4double dxy = std::max(dx,dy);              << 1142   distz1=fDz-p.z();
                                                   >> 1143   distz2=fDz+p.z();
                                                   >> 1144   if (distz1<distz2)
                                                   >> 1145   {
                                                   >> 1146     safe=distz1;
                                                   >> 1147   }
                                                   >> 1148   else
                                                   >> 1149   {
                                                   >> 1150     safe=distz2;
                                                   >> 1151   }
                                                   >> 1152 
                                                   >> 1153   trany=p.y()-fTthetaSphi*p.z(); // Transformed y into `box' system
781                                                   1154 
782   G4double dz = std::abs(p.z())-fDz;           << 1155   // Transformed x into `box' system
783   G4double dist = std::max(dxy,dz);            << 1156   //
                                                   >> 1157   cosy=1.0/std::sqrt(1.0+fTthetaSphi*fTthetaSphi);
                                                   >> 1158   disty1=(fDy-trany)*cosy;
                                                   >> 1159   disty2=(fDy+trany)*cosy;
                                                   >> 1160     
                                                   >> 1161   if (disty1<safe) safe=disty1;
                                                   >> 1162   if (disty2<safe) safe=disty2;
784                                                   1163 
785   return (dist < 0) ? -dist : 0.;              << 1164   tranx=p.x()-fTthetaCphi*p.z()-fTalpha*trany;
                                                   >> 1165   cosx=1.0/std::sqrt(1.0+fTalpha*fTalpha+fTthetaCphi*fTthetaCphi);
                                                   >> 1166   distx1=(fDx-tranx)*cosx;
                                                   >> 1167   distx2=(fDx+tranx)*cosx;
                                                   >> 1168     
                                                   >> 1169   if (distx1<safe) safe=distx1;
                                                   >> 1170   if (distx2<safe) safe=distx2;
                                                   >> 1171     
                                                   >> 1172   if (safe<0) safe=0;
                                                   >> 1173   return safe;  
786 }                                                 1174 }
787                                                   1175 
788 ////////////////////////////////////////////// << 1176 ////////////////////////////////////////////////////////////////////////////////
789 //                                                1177 //
790 // GetEntityType                               << 1178 // Create a List containing the transformed vertices
791                                                << 1179 // Ordering [0-3] -fDz cross section
792 G4GeometryType G4Para::GetEntityType() const   << 1180 //          [4-7] +fDz cross section such that [0] is below [4],
793 {                                              << 1181 //                                             [1] below [5] etc.
794   return {"G4Para"};                           << 1182 // Note:
                                                   >> 1183 //  Caller has deletion resposibility
                                                   >> 1184 
                                                   >> 1185 G4ThreeVectorList*
                                                   >> 1186 G4Para::CreateRotatedVertices( const G4AffineTransform& pTransform ) const
                                                   >> 1187 {
                                                   >> 1188   G4ThreeVectorList *vertices;
                                                   >> 1189   vertices=new G4ThreeVectorList();
                                                   >> 1190   if (vertices)
                                                   >> 1191   {
                                                   >> 1192     vertices->reserve(8);
                                                   >> 1193     G4ThreeVector vertex0(-fDz*fTthetaCphi-fDy*fTalpha-fDx,
                                                   >> 1194                           -fDz*fTthetaSphi-fDy, -fDz);
                                                   >> 1195     G4ThreeVector vertex1(-fDz*fTthetaCphi-fDy*fTalpha+fDx,
                                                   >> 1196                           -fDz*fTthetaSphi-fDy, -fDz);
                                                   >> 1197     G4ThreeVector vertex2(-fDz*fTthetaCphi+fDy*fTalpha-fDx,
                                                   >> 1198                           -fDz*fTthetaSphi+fDy, -fDz);
                                                   >> 1199     G4ThreeVector vertex3(-fDz*fTthetaCphi+fDy*fTalpha+fDx,
                                                   >> 1200                           -fDz*fTthetaSphi+fDy, -fDz);
                                                   >> 1201     G4ThreeVector vertex4(+fDz*fTthetaCphi-fDy*fTalpha-fDx,
                                                   >> 1202                           +fDz*fTthetaSphi-fDy, +fDz);
                                                   >> 1203     G4ThreeVector vertex5(+fDz*fTthetaCphi-fDy*fTalpha+fDx,
                                                   >> 1204                           +fDz*fTthetaSphi-fDy, +fDz);
                                                   >> 1205     G4ThreeVector vertex6(+fDz*fTthetaCphi+fDy*fTalpha-fDx,
                                                   >> 1206                           +fDz*fTthetaSphi+fDy, +fDz);
                                                   >> 1207     G4ThreeVector vertex7(+fDz*fTthetaCphi+fDy*fTalpha+fDx,
                                                   >> 1208                           +fDz*fTthetaSphi+fDy, +fDz);
                                                   >> 1209 
                                                   >> 1210     vertices->push_back(pTransform.TransformPoint(vertex0));
                                                   >> 1211     vertices->push_back(pTransform.TransformPoint(vertex1));
                                                   >> 1212     vertices->push_back(pTransform.TransformPoint(vertex2));
                                                   >> 1213     vertices->push_back(pTransform.TransformPoint(vertex3));
                                                   >> 1214     vertices->push_back(pTransform.TransformPoint(vertex4));
                                                   >> 1215     vertices->push_back(pTransform.TransformPoint(vertex5));
                                                   >> 1216     vertices->push_back(pTransform.TransformPoint(vertex6));
                                                   >> 1217     vertices->push_back(pTransform.TransformPoint(vertex7));
                                                   >> 1218   }
                                                   >> 1219   else
                                                   >> 1220   {
                                                   >> 1221     DumpInfo();
                                                   >> 1222     G4Exception("G4Para::CreateRotatedVertices()",
                                                   >> 1223                 "FatalError", FatalException,
                                                   >> 1224                 "Error in allocation of vertices. Out of memory !");
                                                   >> 1225   }
                                                   >> 1226   return vertices;
795 }                                                 1227 }
796                                                   1228 
797 //////////////////////////////////////////////    1229 //////////////////////////////////////////////////////////////////////////
798 //                                                1230 //
799 // IsFaceted                                   << 1231 // GetEntityType
800                                                   1232 
801 G4bool G4Para::IsFaceted() const               << 1233 G4GeometryType G4Para::GetEntityType() const
802 {                                                 1234 {
803   return true;                                 << 1235   return G4String("G4Para");
804 }                                                 1236 }
805                                                   1237 
806 //////////////////////////////////////////////    1238 //////////////////////////////////////////////////////////////////////////
807 //                                                1239 //
808 // Make a clone of the object                     1240 // Make a clone of the object
809 //                                                1241 //
810 G4VSolid* G4Para::Clone() const                   1242 G4VSolid* G4Para::Clone() const
811 {                                                 1243 {
812   return new G4Para(*this);                       1244   return new G4Para(*this);
813 }                                                 1245 }
814                                                   1246 
815 //////////////////////////////////////////////    1247 //////////////////////////////////////////////////////////////////////////
816 //                                                1248 //
817 // Stream object contents to an output stream     1249 // Stream object contents to an output stream
818                                                   1250 
819 std::ostream& G4Para::StreamInfo( std::ostream    1251 std::ostream& G4Para::StreamInfo( std::ostream& os ) const
820 {                                                 1252 {
821   G4double alpha = std::atan(fTalpha);         << 
822   G4double theta = std::atan(std::sqrt(fTtheta << 
823                                        fTtheta << 
824   G4double phi   = std::atan2(fTthetaSphi,fTth << 
825                                                << 
826   G4long oldprc = os.precision(16);            << 
827   os << "-------------------------------------    1253   os << "-----------------------------------------------------------\n"
828      << "    *** Dump for solid - " << GetName    1254      << "    *** Dump for solid - " << GetName() << " ***\n"
829      << "    =================================    1255      << "    ===================================================\n"
830      << " Solid type: G4Para\n"                   1256      << " Solid type: G4Para\n"
831      << " Parameters:\n"                       << 1257      << " Parameters: \n"
832      << "    half length X: " << fDx/mm << " m << 1258      << "    half length X: " << fDx/mm << " mm \n"
833      << "    half length Y: " << fDy/mm << " m << 1259      << "    half length Y: " << fDy/mm << " mm \n"
834      << "    half length Z: " << fDz/mm << " m << 1260      << "    half length Z: " << fDz/mm << " mm \n"
835      << "    alpha: " << alpha/degree << "degr << 1261      << "    std::tan(alpha)         : " << fTalpha/degree << " degrees \n"
836      << "    theta: " << theta/degree << "degr << 1262      << "    std::tan(theta)*std::cos(phi): " << fTthetaCphi/degree
837      << "    phi: " << phi/degree << "degrees\ << 1263      << " degrees \n"
                                                   >> 1264      << "    std::tan(theta)*std::sin(phi): " << fTthetaSphi/degree
                                                   >> 1265      << " degrees \n"
838      << "-------------------------------------    1266      << "-----------------------------------------------------------\n";
839   os.precision(oldprc);                        << 
840                                                   1267 
841   return os;                                      1268   return os;
842 }                                                 1269 }
843                                                   1270 
844 ////////////////////////////////////////////// << 1271 //////////////////////////////////////////////////////////////////////////////
                                                   >> 1272 //
                                                   >> 1273 // GetPointOnPlane
                                                   >> 1274 // Auxiliary method for Get Point on Surface
845 //                                                1275 //
846 // Return a point randomly and uniformly selec << 
847                                                   1276 
848 G4ThreeVector G4Para::GetPointOnSurface() cons << 1277 G4ThreeVector G4Para::GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1, 
                                                   >> 1278                                       G4ThreeVector p2, G4ThreeVector p3, 
                                                   >> 1279                                       G4double& area) const
849 {                                                 1280 {
850   G4double DyTalpha = fDy*fTalpha;             << 1281   G4double lambda1, lambda2, chose, aOne, aTwo;
851   G4double DzTthetaSphi = fDz*fTthetaSphi;     << 1282   G4ThreeVector t, u, v, w, Area, normal;
852   G4double DzTthetaCphi = fDz*fTthetaCphi;     << 1283   
853                                                << 1284   t = p1 - p0;
854   // Set vertices                              << 1285   u = p2 - p1;
855   //                                           << 1286   v = p3 - p2;
856   G4ThreeVector pt[8];                         << 1287   w = p0 - p3;
857   pt[0].set(-DzTthetaCphi-DyTalpha-fDx, -DzTth << 1288 
858   pt[1].set(-DzTthetaCphi-DyTalpha+fDx, -DzTth << 1289   Area = G4ThreeVector(w.y()*v.z() - w.z()*v.y(),
859   pt[2].set(-DzTthetaCphi+DyTalpha-fDx, -DzTth << 1290                        w.z()*v.x() - w.x()*v.z(),
860   pt[3].set(-DzTthetaCphi+DyTalpha+fDx, -DzTth << 1291                        w.x()*v.y() - w.y()*v.x());
861   pt[4].set( DzTthetaCphi-DyTalpha-fDx,  DzTth << 
862   pt[5].set( DzTthetaCphi-DyTalpha+fDx,  DzTth << 
863   pt[6].set( DzTthetaCphi+DyTalpha-fDx,  DzTth << 
864   pt[7].set( DzTthetaCphi+DyTalpha+fDx,  DzTth << 
865                                                << 
866   // Set areas (-Z, -Y, +Y, -X, +X, +Z)        << 
867   //                                           << 
868   G4ThreeVector vx(fDx, 0, 0);                 << 
869   G4ThreeVector vy(DyTalpha, fDy, 0);          << 
870   G4ThreeVector vz(DzTthetaCphi, DzTthetaSphi, << 
871                                                << 
872   G4double sxy = fDx*fDy; // (vx.cross(vy)).ma << 
873   G4double sxz = (vx.cross(vz)).mag();         << 
874   G4double syz = (vy.cross(vz)).mag();         << 
875                                                   1292   
876   G4double sface[6] = { sxy, syz, syz, sxz, sx << 1293   aOne = 0.5*Area.mag();
877   for (G4int i=1; i<6; ++i) { sface[i] += sfac << 1294   
                                                   >> 1295   Area = G4ThreeVector(t.y()*u.z() - t.z()*u.y(),
                                                   >> 1296                        t.z()*u.x() - t.x()*u.z(),
                                                   >> 1297                        t.x()*u.y() - t.y()*u.x());
                                                   >> 1298   
                                                   >> 1299   aTwo = 0.5*Area.mag();
                                                   >> 1300   
                                                   >> 1301   area = aOne + aTwo;
                                                   >> 1302   
                                                   >> 1303   chose = RandFlat::shoot(0.,aOne+aTwo);
878                                                   1304 
879   // Select face                               << 1305   if( (chose>=0.) && (chose < aOne) )
880   //                                           << 1306   {
881   G4double select = sface[5]*G4UniformRand();  << 1307     lambda1 = RandFlat::shoot(0.,1.);
882   G4int k = 5;                                 << 1308     lambda2 = RandFlat::shoot(0.,lambda1);
883   if (select <= sface[4]) k = 4;               << 1309     return (p2+lambda1*v+lambda2*w);    
884   if (select <= sface[3]) k = 3;               << 1310   }
885   if (select <= sface[2]) k = 2;               << 1311 
886   if (select <= sface[1]) k = 1;               << 1312   // else
887   if (select <= sface[0]) k = 0;               << 1313 
888                                                << 1314   lambda1 = RandFlat::shoot(0.,1.);
889   // Generate point                            << 1315   lambda2 = RandFlat::shoot(0.,lambda1);
890   //                                           << 1316   return (p0+lambda1*t+lambda2*u);    
891   G4int ip[6][3] = {{0,1,2}, {0,4,1}, {2,3,6}, << 
892   G4double u = G4UniformRand();                << 
893   G4double v = G4UniformRand();                << 
894   return (1.-u-v)*pt[ip[k][0]] + u*pt[ip[k][1] << 
895 }                                                 1317 }
896                                                   1318 
897 ////////////////////////////////////////////// << 1319 /////////////////////////////////////////////////////////////////////////
                                                   >> 1320 //
                                                   >> 1321 // GetPointOnSurface
                                                   >> 1322 //
                                                   >> 1323 // Return a point (G4ThreeVector) randomly and uniformly
                                                   >> 1324 // selected on the solid surface
                                                   >> 1325 
                                                   >> 1326 G4ThreeVector G4Para::GetPointOnSurface() const
                                                   >> 1327 {
                                                   >> 1328   G4ThreeVector One, Two, Three, Four, Five, Six;
                                                   >> 1329   G4ThreeVector pt[8] ;
                                                   >> 1330   G4double chose, aOne, aTwo, aThree, aFour, aFive, aSix;
                                                   >> 1331 
                                                   >> 1332   pt[0] = G4ThreeVector(-fDz*fTthetaCphi-fDy*fTalpha-fDx,
                                                   >> 1333                         -fDz*fTthetaSphi-fDy, -fDz);
                                                   >> 1334   pt[1] = G4ThreeVector(-fDz*fTthetaCphi-fDy*fTalpha+fDx,
                                                   >> 1335                         -fDz*fTthetaSphi-fDy, -fDz);
                                                   >> 1336   pt[2] = G4ThreeVector(-fDz*fTthetaCphi+fDy*fTalpha-fDx,
                                                   >> 1337                         -fDz*fTthetaSphi+fDy, -fDz);
                                                   >> 1338   pt[3] = G4ThreeVector(-fDz*fTthetaCphi+fDy*fTalpha+fDx,
                                                   >> 1339                         -fDz*fTthetaSphi+fDy, -fDz);
                                                   >> 1340   pt[4] = G4ThreeVector(+fDz*fTthetaCphi-fDy*fTalpha-fDx,
                                                   >> 1341                         +fDz*fTthetaSphi-fDy, +fDz);
                                                   >> 1342   pt[5] = G4ThreeVector(+fDz*fTthetaCphi-fDy*fTalpha+fDx,
                                                   >> 1343                         +fDz*fTthetaSphi-fDy, +fDz);
                                                   >> 1344   pt[6] = G4ThreeVector(+fDz*fTthetaCphi+fDy*fTalpha-fDx,
                                                   >> 1345                         +fDz*fTthetaSphi+fDy, +fDz);
                                                   >> 1346   pt[7] = G4ThreeVector(+fDz*fTthetaCphi+fDy*fTalpha+fDx,
                                                   >> 1347                         +fDz*fTthetaSphi+fDy, +fDz);
                                                   >> 1348 
                                                   >> 1349   // make sure we provide the points in a clockwise fashion
                                                   >> 1350 
                                                   >> 1351   One   = GetPointOnPlane(pt[0],pt[1],pt[3],pt[2], aOne);
                                                   >> 1352   Two   = GetPointOnPlane(pt[4],pt[5],pt[7],pt[6], aTwo);
                                                   >> 1353   Three = GetPointOnPlane(pt[6],pt[7],pt[3],pt[2], aThree);
                                                   >> 1354   Four  = GetPointOnPlane(pt[4],pt[5],pt[1],pt[0], aFour); 
                                                   >> 1355   Five  = GetPointOnPlane(pt[0],pt[2],pt[6],pt[4], aFive);
                                                   >> 1356   Six   = GetPointOnPlane(pt[1],pt[3],pt[7],pt[5], aSix);
                                                   >> 1357 
                                                   >> 1358   chose = RandFlat::shoot(0.,aOne+aTwo+aThree+aFour+aFive+aSix);
                                                   >> 1359   
                                                   >> 1360   if( (chose>=0.) && (chose<aOne) )                    
                                                   >> 1361     { return One; }
                                                   >> 1362   else if(chose>=aOne && chose<aOne+aTwo)  
                                                   >> 1363     { return Two; }
                                                   >> 1364   else if(chose>=aOne+aTwo && chose<aOne+aTwo+aThree)
                                                   >> 1365     { return Three; }
                                                   >> 1366   else if(chose>=aOne+aTwo+aThree && chose<aOne+aTwo+aThree+aFour)
                                                   >> 1367     { return Four; }
                                                   >> 1368   else if(chose>=aOne+aTwo+aThree+aFour && chose<aOne+aTwo+aThree+aFour+aFive)
                                                   >> 1369     { return Five; }
                                                   >> 1370   return Six;
                                                   >> 1371 }
                                                   >> 1372 
                                                   >> 1373 ////////////////////////////////////////////////////////////////////////////
898 //                                                1374 //
899 // Methods for visualisation                      1375 // Methods for visualisation
900                                                   1376 
901 void G4Para::DescribeYourselfTo ( G4VGraphicsS    1377 void G4Para::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
902 {                                                 1378 {
903   scene.AddSolid (*this);                         1379   scene.AddSolid (*this);
904 }                                                 1380 }
905                                                   1381 
906 G4Polyhedron* G4Para::CreatePolyhedron () cons    1382 G4Polyhedron* G4Para::CreatePolyhedron () const
907 {                                                 1383 {
908   G4double phi = std::atan2(fTthetaSphi, fTthe    1384   G4double phi = std::atan2(fTthetaSphi, fTthetaCphi);
909   G4double alpha = std::atan(fTalpha);            1385   G4double alpha = std::atan(fTalpha);
910   G4double theta = std::atan(std::sqrt(fTtheta << 1386   G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi
911                                        fTtheta << 1387                             +fTthetaSphi*fTthetaSphi));
912                                                   1388     
913   return new G4PolyhedronPara(fDx, fDy, fDz, a    1389   return new G4PolyhedronPara(fDx, fDy, fDz, alpha, theta, phi);
914 }                                                 1390 }
915 #endif                                         << 1391 
                                                   >> 1392 G4NURBS* G4Para::CreateNURBS () const
                                                   >> 1393 {
                                                   >> 1394   // return new G4NURBSbox (fDx, fDy, fDz);
                                                   >> 1395   return 0 ;
                                                   >> 1396 }
916                                                   1397