Ignore:
Timestamp:
Sep 10, 2008, 5:40:37 PM (16 years ago)
Author:
garnier
Message:

geant4.8.2 beta

Location:
trunk/source/geometry/solids/specific/src
Files:
36 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/geometry/solids/specific/src/G4ClippablePolygon.cc

    r831 r850  
    2626//
    2727// $Id: G4ClippablePolygon.cc,v 1.12 2007/05/11 13:54:28 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/specific/src/G4Ellipsoid.cc

    r831 r850  
    2525//
    2626// $Id: G4Ellipsoid.cc,v 1.14 2007/05/18 07:39:56 gcosmo Exp $
    27 // GEANT4 tag $Name: $
     27// GEANT4 tag $Name: HEAD $
    2828//
    2929// class G4Ellipsoid

  • trunk/source/geometry/solids/specific/src/G4EllipticalCone.cc

    r831 r850  
    2424// ********************************************************************
    2525//
    26 // $Id: G4EllipticalCone.cc,v 1.15.2.1 2008/04/25 09:11:34 gcosmo Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4EllipticalCone.cc,v 1.16 2008/04/25 08:45:26 gcosmo Exp $
     27// GEANT4 tag $Name: HEAD $
    2828//
    2929// Implementation of G4EllipticalCone class

  • trunk/source/geometry/solids/specific/src/G4EllipticalTube.cc

    r831 r850  
    2626//
    2727// $Id: G4EllipticalTube.cc,v 1.27 2006/10/20 13:45:21 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/specific/src/G4EnclosingCylinder.cc

    r831 r850  
    2626//
    2727// $Id: G4EnclosingCylinder.cc,v 1.10 2007/05/11 13:54:29 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/specific/src/G4ExtrudedSolid.cc

    r831 r850  
    2525//
    2626//
    27 // $Id: G4ExtrudedSolid.cc,v 1.11.2.1 2008/04/23 08:10:24 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     27// $Id: G4ExtrudedSolid.cc,v 1.17 2008/08/12 08:54:57 gcosmo Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//
    … …  
    6262  // General constructor
    6363
     64  G4String errorDescription = "InvalidSetup in \"";
     65  errorDescription += pName;
     66  errorDescription += "\"";
     67
    6468  // First check input parameters
    6569
    66   if ( fNv < 3 ) {
    67     G4Exception(
    68       "G4ExtrudedSolid::G4ExtrudedSolid()", "InvalidSetup",
    69       FatalException, "Number of polygon vertices < 3");
     70  if ( fNv < 3 )
     71  {
     72    G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", errorDescription,
     73                FatalException, "Number of polygon vertices < 3");
    7074  }
    7175     
    72   if ( fNz < 2 ) {
    73     G4Exception(
    74       "G4ExtrudedSolid::G4ExtrudedSolid()", "InvalidSetup",
    75       FatalException, "Number of z-sides < 2");
     76  if ( fNz < 2 )
     77  {
     78    G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", errorDescription,
     79                FatalException, "Number of z-sides < 2");
    7680  }
    7781     
    … …  
    8084    if ( zsections[i].fZ > zsections[i+1].fZ )
    8185    {
    82       G4Exception(
    83         "G4ExtrudedSolid::G4ExtrudedSolid()", "InvalidSetup",
     86      G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", errorDescription,
    8487        FatalException,
    8588        "Z-sections have to be ordered by z value (z0 < z1 < z2 ...)");
    8689    }
    87     if ( std::fabs( zsections[i+1].fZ - zsections[i].fZ ) < kCarTolerance )
     90    if ( std::fabs( zsections[i+1].fZ - zsections[i].fZ ) < kCarTolerance * 0.5 )
    8891    {
    89       G4Exception(
    90         "G4ExtrudedSolid::G4ExtrudedSolid()", "InvalidSetup",
     92      G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", errorDescription,
    9193        FatalException,
    9294        "Z-sections with the same z position are not supported.");
    … …  
    106108  if (!result)
    107109  {   
    108     G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "InvalidSetup",
     110    G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", errorDescription,
    109111                FatalException, "Making facets failed.");
    110112  }
    … …  
    134136  // Special constructor for solid with 2 z-sections
    135137
     138  G4String errorDescription = "InvalidSetup in \"";
     139  errorDescription += pName;
     140  errorDescription += "\"";
     141
    136142  // First check input parameters
    137143  //
    138144  if ( fNv < 3 )
    139145  {
    140     G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "InvalidSetup",
     146    G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", errorDescription,
    141147                FatalException, "Number of polygon vertices < 3");
    142148  }
    … …  
    154160  if (!result)
    155161  {   
    156     G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "InvalidSetup",
     162    G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", errorDescription,
    157163                FatalException, "Making facets failed.");
    158164  }
    … …  
    265271  if ( l1.x() == l2.x() )
    266272  {
    267     return std::fabs(p.x() - l1.x()) < kCarTolerance;
     273    return std::fabs(p.x() - l1.x()) < kCarTolerance * 0.5;
    268274  }
    269275
    270276  return std::fabs (p.y() - l1.y() - ((l2.y() - l1.y())/(l2.x() - l1.x()))
    271                                     *(p.x() - l1.x())) < kCarTolerance;
     277                                    *(p.x() - l1.x())) < kCarTolerance * 0.5;
    272278 }
    273279
    … …  
    280286  // l1 and l2
    281287
    282   if ( p.x() < std::min(l1.x(), l2.x()) - kCarTolerance ||
    283        p.x() > std::max(l1.x(), l2.x()) + kCarTolerance ||
    284        p.y() < std::min(l1.y(), l2.y()) - kCarTolerance||
    285        p.y() > std::max(l1.y(), l2.y()) + kCarTolerance )
     288  if ( p.x() < std::min(l1.x(), l2.x()) - kCarTolerance * 0.5 ||
     289       p.x() > std::max(l1.x(), l2.x()) + kCarTolerance * 0.5 ||
     290       p.y() < std::min(l1.y(), l2.y()) - kCarTolerance * 0.5 ||
     291       p.y() > std::max(l1.y(), l2.y()) + kCarTolerance * 0.5 )
    286292  {
    287293    return false;
    … …  
    309315                                      G4TwoVector c, G4TwoVector p) const
    310316{
    311   // Return true if p is inside of triangle abc, else returns false
     317  // Return true if p is inside of triangle abc or on its edges,
     318  // else returns false
    312319
    313320  // Check extent first
    … …  
    318325       ( p.y() > a.y() && p.y() > b.y() && p.y() > c.y() ) ) return false;
    319326 
    320   return   IsSameSide(p, a, b, c)
    321         && IsSameSide(p, b, a, c)
    322         && IsSameSide(p, c, a, b);
     327  G4bool inside
     328    = IsSameSide(p, a, b, c)
     329      && IsSameSide(p, b, a, c)
     330      && IsSameSide(p, c, a, b);
     331
     332  G4bool onEdge
     333    = IsSameLineSegment(p, a, b)
     334      || IsSameLineSegment(p, b, c)
     335      || IsSameLineSegment(p, c, a);
     336     
     337  return inside || onEdge;   
    323338}     
    324339
    … …  
    333348  G4TwoVector t2 = pb - po;
    334349 
    335   G4double result
    336     = (atan2(t1.y(), t1.x()) - atan2(t2.y(), t2.x()));
     350  G4double result = (std::atan2(t1.y(), t1.x()) - std::atan2(t2.y(), t2.x()));
    337351
    338352  if ( result < 0 ) result += 2*pi;
    … …  
    439453    //
    440454    G4double angle = GetAngle(c2->first, c3->first, c1->first);
    441     //G4cout << angle << G4endl;
    442     if ( angle > pi ) {
     455
     456    if ( angle > pi )
     457    {
    443458      // G4cout << "Skipping concave vertex " << c2->second << G4endl;
    444459
    … …  
    449464      ++c3;
    450465      if ( c3 == verticesToBeDone.end() ) { c3 = verticesToBeDone.begin(); }
     466
     467      // G4cout << "Looking at triangle : "
     468      //        << c1->second << "  " << c2->second
     469      //        << "  " << c3->second << G4endl; 
     470
    451471    }
    452472
    … …  
    623643  // Check first if outside extent
    624644  //
    625   if ( p.x() < GetMinXExtent() - kCarTolerance ||
    626        p.x() > GetMaxXExtent() + kCarTolerance ||
    627        p.y() < GetMinYExtent() - kCarTolerance ||
    628        p.y() > GetMaxYExtent() + kCarTolerance ||
    629        p.z() < GetMinZExtent() - kCarTolerance ||
    630        p.z() > GetMaxZExtent() + kCarTolerance )
     645  if ( p.x() < GetMinXExtent() - kCarTolerance * 0.5 ||
     646       p.x() > GetMaxXExtent() + kCarTolerance * 0.5 ||
     647       p.y() < GetMinYExtent() - kCarTolerance * 0.5 ||
     648       p.y() > GetMaxYExtent() + kCarTolerance * 0.5 ||
     649       p.z() < GetMinZExtent() - kCarTolerance * 0.5 ||
     650       p.z() > GetMaxZExtent() + kCarTolerance * 0.5 )
    631651  {
    632652    // G4cout << "G4ExtrudedSolid::Outside extent: " << p << G4endl;
    … …  
    660680    if ( IsPointInside(fPolygon[(*it)[0]], fPolygon[(*it)[1]],
    661681                       fPolygon[(*it)[2]], pscaled) )  { inside = true; }
    662                        
    663     if ( IsSameLineSegment(pscaled, fPolygon[(*it)[0]], fPolygon[(*it)[1]]) ||
    664          IsSameLineSegment(pscaled, fPolygon[(*it)[1]], fPolygon[(*it)[2]]) ||
    665          IsSameLineSegment(pscaled, fPolygon[(*it)[2]], fPolygon[(*it)[0]]) )
    666                                                        { inside = true; }
    667682    ++it;
    668683  } while ( (inside == false) && (it != fTriangles.end()) );
    … …  
    672687    // Check if on surface of z sides
    673688    //
    674     if ( std::fabs( p.z() - fZSections[0].fZ ) < kCarTolerance ||
    675          std::fabs( p.z() - fZSections[fNz-1].fZ ) < kCarTolerance )
     689    if ( std::fabs( p.z() - fZSections[0].fZ ) < kCarTolerance * 0.5 ||
     690         std::fabs( p.z() - fZSections[fNz-1].fZ ) < kCarTolerance * 0.5 )
    676691    {
    677692      // G4cout << "G4ExtrudedSolid::Inside return Surface (on z side)"

  • trunk/source/geometry/solids/specific/src/G4Hype.cc

    r831 r850  
    2525//
    2626//
    27 // $Id: G4Hype.cc,v 1.25.8.1 2008/04/23 08:10:24 gcosmo Exp $
     27// $Id: G4Hype.cc,v 1.27 2008/04/14 08:49:28 gcosmo Exp $
    2828// $Original: G4Hype.cc,v 1.0 1998/06/09 16:57:50 safai Exp $
    29 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     29// GEANT4 tag $Name: HEAD $
    3030//
    3131//

  • trunk/source/geometry/solids/specific/src/G4IntersectingCone.cc

    r831 r850  
    2525//
    2626//
    27 // $Id: G4IntersectingCone.cc,v 1.8.8.2 2008/04/28 09:03:16 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     27// $Id: G4IntersectingCone.cc,v 1.12 2008/04/28 08:59:47 gcosmo Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/specific/src/G4Paraboloid.cc

    r831 r850  
    2424// ********************************************************************
    2525//
    26 // $Id: G4Paraboloid.cc,v 1.5.2.1 2008/04/23 08:10:24 gcosmo Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: G4Paraboloid.cc,v 1.8 2008/07/17 07:33:00 gcosmo Exp $
     27// GEANT4 tag $Name: HEAD $
    2828//
    2929// class G4Paraboloid
    3030//
    3131// Implementation for G4Paraboloid class
    32 //
    33 // History:
    3432//
    3533// Author : Lukas Lindroos (CERN), July 2007
    … …  
    8280  }
    8381
    84 // r1^2 = k1 * (-dz) + k2
    85 // r2^2 = k1 * ( dz) + k2
    86 // => r1^2 + r2^2 = k2 + k2 => k2 = (r2^2 + r1^2) / 2
    87 // and r2^2 - r1^2 = k1 * dz - k1 * (-dz) => k1 = (r2^2 - r1^2) / 2 / dz
     82  // r1^2 = k1 * (-dz) + k2
     83  // r2^2 = k1 * ( dz) + k2
     84  // => r1^2 + r2^2 = k2 + k2 => k2 = (r2^2 + r1^2) / 2
     85  // and r2^2 - r1^2 = k1 * dz - k1 * (-dz) => k1 = (r2^2 - r1^2) / 2 / dz
    8886
    8987  k1 = (r2 * r2 - r1 * r1) / 2 / dz;
    … …  
    283281           rhoSurfTimesTol2  = (k1 * p.z() + k2) * sqr(kCarTolerance),
    284282           A = rho2 - ((k1 *p.z() + k2) + 0.25 * kCarTolerance * kCarTolerance);
    285 
     283 
    286284  if(A < 0 && sqr(A) > rhoSurfTimesTol2)
    287285  {
    288286    // Actually checking rho < radius of paraboloid at z = p.z().
    289287    // We're either inside or in lower/upper cutoff area.
    290 
     288   
    291289    if(std::fabs(p.z()) > dz - 0.5 * kCarTolerance)
    292290    {
    … …  
    432430      }
    433431    }
    434     else // Direction away, no posibility of intersection
    435       { return kInfinity; }
     432    else  // Direction away, no possibility of intersection
     433    {
     434      return kInfinity;
     435    }
    436436  }
    437437  else if(r1 && p.z() < tolh - dz)
    … …  
    455455      }
    456456    }
    457     else// Direction away, no posibility of intersection
    458       { return kInfinity; }
     457    else  // Direction away, no possibility of intersection
     458    {
     459      return kInfinity;
     460    }
    459461  }
    460462
    … …  
    577579///////////////////////////////////////////////////////////////////////////////
    578580//
    579 // Calculate distance to surface of shape from `inside'
     581// Calculate distance to surface of shape from 'inside'
    580582
    581583G4double G4Paraboloid::DistanceToOut(const G4ThreeVector& p,
    … …  
    598600  // to include all z) x^2 + y^2 = k1 * z + k2 => .. =>
    599601  // => s = (A +- std::sqrt(A^2 + B)) / vRho2
    600   // where
     602  // where:
     603  //
    601604  G4double A = k1 / 2 * v.z() - p.x() * v.x() - p.y() * v.y();
    602   // and
     605  //
     606  // and:
     607  //
    603608  G4double B = (-rho2 + paraRho2) * vRho2;
    604609
    … …  
    676681    else if( ((A <= 0) && (B >= sqr(A) * (sqr(vRho2) - 1))) || (A >= 0))
    677682    {
    678       intersection = (A + std::sqrt(B + sqr(A))) / vRho2;
     683      // intersection = (A + std::sqrt(B + sqr(A))) / vRho2;
     684      // The above calculation has a precision problem:
     685      // known problem of solving quadratic equation with small A 
     686
     687      A = A/vRho2;
     688      B = (k1 * p.z() + k2 - rho2)/vRho2;
     689      intersection = B/(-A + std::sqrt(B + sqr(A)));
    679690      if(calcNorm)
    680691      {
    … …  
    700711  {
    701712    // If this is true we're somewhere in the border.
    702 
     713   
    703714    G4ThreeVector normal = G4ThreeVector (p.x(), p.y(), -k1/2);
    704715
    … …  
    706717    {
    707718      // We're in the lower or upper edge
     719      //
    708720      if( ((v.z() > 0) && (p.z() > 0)) || ((v.z() < 0) && (p.z() < 0)) )
    709       // If we're headig out of the object that is treated here
    710       {
     721      {             // If we're heading out of the object that is treated here
    711722        if(calcNorm)
    712723        {
    … …  
    742753      }
    743754    }
    744     else if(normal.dot(v) >= 0)
    745     {
    746       if(calcNorm)
    747       {
    748         *validNorm = true;
    749         *n = normal.unit();
    750       }
    751       return 0;
    752     }
     755    //
     756    // Problem in the Logic :: Following condition for point on upper surface
     757    //                         and Vz<0  will return 0 (Problem #1015), but
     758    //                         it has to return intersection with parabolic
     759    //                         surface or with lower plane surface (z = -dz)
     760    // The logic has to be  :: If not found intersection until now,
     761    // do not exit but continue to search for possible intersection.
     762    // Only for point situated on both borders (Z and parabolic)
     763    // this condition has to be taken into account and done later
     764    //
     765    //
     766    // else if(normal.dot(v) >= 0)
     767    // {
     768    //   if(calcNorm)
     769    //   {
     770    //     *validNorm = true;
     771    //     *n = normal.unit();
     772    //   }
     773    //   return 0;
     774    // }
    753775
    754776    if(v.z() > 0)
    … …  
    780802      }
    781803    }
    782     if(r1 && v.z() < 0)
     804    if( v.z() < 0)
    783805    {
    784806      // Check for collision with lower edge.
    … …  
    809831    }
    810832
    811     if(vRho2 != 0)
    812       { intersection = (A + std::sqrt(B + sqr(A))) / vRho2; }
     833    // Note: comparison with zero below would not be correct !
     834    //
     835    if(std::fabs(vRho2) > tol2) // precision error in the calculation of
     836    {                           // intersection = (A+std::sqrt(B+sqr(A)))/vRho2
     837      A = A/vRho2;
     838      B = (k1 * p.z() + k2 - rho2);
     839      if(std::fabs(B)>kCarTolerance)
     840      {
     841        B = (B)/vRho2;
     842        intersection = B/(-A + std::sqrt(B + sqr(A)));
     843      }
     844      else                      // Point is On both borders: Z and parabolic
     845      {                         // solution depends on normal.dot(v) sign
     846        if(normal.dot(v) >= 0)
     847        {
     848          if(calcNorm)
     849          {
     850            *validNorm = true;
     851            *n = normal.unit();
     852          }
     853          return 0;
     854        }
     855        intersection = 2.*A;
     856      }
     857    }
    813858    else
    814       { intersection = ((rho2 - k2) / k1 - p.z()) / v.z(); }
     859    {
     860      intersection = ((rho2 - k2) / k1 - p.z()) / v.z();
     861    }
    815862
    816863    if(calcNorm)

  • trunk/source/geometry/solids/specific/src/G4PolyPhiFace.cc

    r831 r850  
    2525//
    2626//
    27 // $Id: G4PolyPhiFace.cc,v 1.13 2007/07/19 12:57:14 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: G4PolyPhiFace.cc,v 1.15 2008/05/15 11:41:59 gcosmo Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//
    … …  
    4747#include "G4GeometryTolerance.hh"
    4848
     49#include "Randomize.hh"
     50#include "G4TwoVector.hh"
     51
    4952//
    5053// Constructor
    … …  
    6366{
    6467  kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
     68  fSurfaceArea = 0.; 
    6569
    6670  numEdges = rz->NumVertices();
    … …  
    103107  //
    104108  corners = new G4PolyPhiFaceVertex[numEdges];
    105 
    106109  //
    107110  // Fill them
    … …  
    110113 
    111114  G4PolyPhiFaceVertex *corn = corners;
     115  G4PolyPhiFaceVertex *helper=corners;
     116
    112117  iterRZ.Begin();
    113118  do
    … …  
    117122    corn->x = corn->r*radial.x();
    118123    corn->y = corn->r*radial.y();
     124
     125    // Add pointer on prev corner
     126    //
     127    if( corn == corners )
     128      { corn->prev = corners+numEdges-1;}
     129    else
     130      { corn->prev = helper; }
     131
     132    // Add pointer on next corner
     133    //
     134    if( corn < corners+numEdges-1 )
     135      { corn->next = corn+1;}
     136    else
     137      { corn->next = corners; }
     138
     139    helper = corn;
    119140  } while( ++corn, iterRZ.Next() );
    120141
    … …  
    322343  normal    = source.normal;
    323344  radial    = source.radial;
    324   surface    = source.surface;
     345  surface   = source.surface;
    325346  rMin    = source.rMin;
    326347  rMax    = source.rMax;
    … …  
    330351
    331352  kCarTolerance = source.kCarTolerance;
     353  fSurfaceArea = source.fSurfaceArea;
    332354
    333355  //
    … …  
    895917  return answer;
    896918}
     919
     920//
     921// Calculation of Surface Area of a Triangle
     922// In the same time Random Point in Triangle is given
     923//
     924G4double G4PolyPhiFace::SurfaceTriangle( G4ThreeVector p1,
     925                                         G4ThreeVector p2,
     926                                         G4ThreeVector p3,
     927                                         G4ThreeVector *p4 )
     928{
     929  G4ThreeVector v, w;
     930 
     931  v = p3 - p1;
     932  w = p1 - p2;
     933  G4double lambda1 = G4UniformRand();
     934  G4double lambda2 = lambda1*G4UniformRand();
     935 
     936  *p4=p2 + lambda1*w + lambda2*v;
     937  return 0.5*(v.cross(w)).mag();
     938}
     939
     940//
     941// Compute surface area
     942//
     943G4double G4PolyPhiFace::SurfaceArea()
     944{
     945  if ( fSurfaceArea==0. ) { Triangulate(); }
     946  return fSurfaceArea;
     947}
     948
     949//
     950// Return random point on face
     951//
     952G4ThreeVector G4PolyPhiFace::GetPointOnFace()
     953{
     954  Triangulate();
     955  return surface_point;
     956}
     957
     958//
     959// Auxiliary Functions used for Finding the PointOnFace using Triangulation
     960//
     961
     962//
     963// Calculation of 2*Area of Triangle with Sign
     964//
     965G4double G4PolyPhiFace::Area2( G4TwoVector a,
     966                               G4TwoVector b,
     967                               G4TwoVector c )
     968{
     969  return ((b.x()-a.x())*(c.y()-a.y())-
     970          (c.x()-a.x())*(b.y()-a.y()));
     971}
     972
     973//
     974// Boolean function for sign of Surface
     975//
     976G4bool G4PolyPhiFace::Left( G4TwoVector a,
     977                            G4TwoVector b,
     978                            G4TwoVector c )
     979{
     980  return Area2(a,b,c)>0;
     981}
     982
     983//
     984// Boolean function for sign of Surface
     985//
     986G4bool G4PolyPhiFace::LeftOn( G4TwoVector a,
     987                              G4TwoVector b,
     988                              G4TwoVector c )
     989{
     990  return Area2(a,b,c)>=0;
     991}
     992
     993//
     994// Boolean function for sign of Surface
     995//
     996G4bool G4PolyPhiFace::Collinear( G4TwoVector a,
     997                                 G4TwoVector b,
     998                                 G4TwoVector c )
     999{
     1000  return Area2(a,b,c)==0;
     1001}
     1002
     1003//
     1004// Boolean function for finding "Proper" Intersection
     1005// That means Intersection of two lines segments (a,b) and (c,d)
     1006//
     1007G4bool G4PolyPhiFace::IntersectProp( G4TwoVector a,
     1008                                     G4TwoVector b,
     1009                                     G4TwoVector c, G4TwoVector d )
     1010{
     1011  if( Collinear(a,b,c) || Collinear(a,b,d)||
     1012      Collinear(c,d,a) || Collinear(c,d,b) )  { return false; }
     1013
     1014  G4bool Positive;
     1015  Positive = !(Left(a,b,c))^!(Left(a,b,d));
     1016  return Positive && (!Left(c,d,a)^!Left(c,d,b));
     1017}
     1018
     1019//
     1020// Boolean function for determining if Point c is between a and b
     1021// For the tree points(a,b,c) on the same line
     1022//
     1023G4bool G4PolyPhiFace::Between( G4TwoVector a, G4TwoVector b, G4TwoVector c )
     1024{
     1025  if( !Collinear(a,b,c) ) { return false; }
     1026
     1027  if(a.x()!=b.x())
     1028  {
     1029    return ((a.x()<=c.x())&&(c.x()<=b.x()))||
     1030           ((a.x()>=c.x())&&(c.x()>=b.x()));
     1031  }
     1032  else
     1033  {
     1034    return ((a.y()<=c.y())&&(c.y()<=b.y()))||
     1035           ((a.y()>=c.y())&&(c.y()>=b.y()));
     1036  }
     1037}
     1038
     1039//
     1040// Boolean function for finding Intersection "Proper" or not
     1041// Between two line segments (a,b) and (c,d)
     1042//
     1043G4bool G4PolyPhiFace::Intersect( G4TwoVector a,
     1044                                 G4TwoVector b,
     1045                                 G4TwoVector c, G4TwoVector d )
     1046{
     1047 if( IntersectProp(a,b,c,d) )
     1048   { return true; }
     1049 else if( Between(a,b,c)||
     1050          Between(a,b,d)||
     1051          Between(c,d,a)||
     1052          Between(c,d,b) )
     1053   { return true; }
     1054 else
     1055   { return false; }
     1056}
     1057
     1058//
     1059// Boolean Diagonalie help to determine
     1060// if diagonal s of segment (a,b) is convex or reflex
     1061//
     1062G4bool G4PolyPhiFace::Diagonalie( G4PolyPhiFaceVertex *a,
     1063                                  G4PolyPhiFaceVertex *b )
     1064{
     1065  G4PolyPhiFaceVertex   *corner = triangles;
     1066  G4PolyPhiFaceVertex   *corner_next=triangles;
     1067
     1068  // For each Edge (corner,corner_next)
     1069  do
     1070  {
     1071    corner_next=corner->next;
     1072
     1073    // Skip edges incident to a of b
     1074    //
     1075    if( (corner!=a)&&(corner_next!=a)
     1076      &&(corner!=b)&&(corner_next!=b) )
     1077    {
     1078       G4TwoVector rz1,rz2,rz3,rz4;
     1079       rz1 = G4TwoVector(a->r,a->z);
     1080       rz2 = G4TwoVector(b->r,b->z);
     1081       rz3 = G4TwoVector(corner->r,corner->z);
     1082       rz4 = G4TwoVector(corner_next->r,corner_next->z);
     1083       if( Intersect(rz1,rz2,rz3,rz4) ) { return false; }
     1084    }
     1085    corner=corner->next;   
     1086   
     1087  } while( corner != triangles );
     1088
     1089  return true;
     1090}
     1091
     1092//
     1093// Boolean function that determine if b is Inside Cone (a0,a,a1)
     1094// being a the center of the Cone
     1095//
     1096G4bool G4PolyPhiFace::InCone( G4PolyPhiFaceVertex *a, G4PolyPhiFaceVertex *b )
     1097{
     1098  // a0,a and a1 are consecutive vertices
     1099  //
     1100  G4PolyPhiFaceVertex *a0,*a1;
     1101  a1=a->next;
     1102  a0=a->prev;
     1103
     1104  G4TwoVector arz,arz0,arz1,brz;
     1105  arz=G4TwoVector(a->r,a->z);arz0=G4TwoVector(a0->r,a0->z);
     1106  arz1=G4TwoVector(a1->r,a1->z);brz=G4TwoVector(b->r,b->z);
     1107   
     1108 
     1109  if(LeftOn(arz,arz1,arz0))  // If a is convex vertex
     1110  {
     1111    return Left(arz,brz,arz0)&&Left(brz,arz,arz1);
     1112  }
     1113  else                       // Else a is reflex
     1114  {
     1115    return !( LeftOn(arz,brz,arz1)&&LeftOn(brz,arz,arz0));
     1116  }
     1117}
     1118
     1119//
     1120// Boolean function finding if Diagonal is possible
     1121// inside Polycone or PolyHedra
     1122//
     1123G4bool G4PolyPhiFace::Diagonal( G4PolyPhiFaceVertex *a, G4PolyPhiFaceVertex *b )
     1124{
     1125  return InCone(a,b) && InCone(b,a) && Diagonalie(a,b);
     1126}
     1127
     1128//
     1129// Initialisation for Triangulisation by ear tips
     1130// For details see "Computational Geometry in C" by Joseph O'Rourke
     1131//
     1132void G4PolyPhiFace::EarInit()
     1133{
     1134  G4PolyPhiFaceVertex   *corner = triangles;
     1135  G4PolyPhiFaceVertex *c_prev,*c_next;
     1136 
     1137  do
     1138  {
     1139     // We need to determine three consecutive vertices
     1140     //
     1141     c_next=corner->next;
     1142     c_prev=corner->prev;
     1143
     1144     // Calculation of ears
     1145     //
     1146     corner->ear=Diagonal(c_prev,c_next);   
     1147     corner=corner->next;
     1148
     1149  } while( corner!=triangles );
     1150}
     1151
     1152//
     1153// Triangulisation by ear tips for Polycone or Polyhedra
     1154// For details see "Computational Geometry in C" by Joseph O'Rourke
     1155//
     1156void G4PolyPhiFace::Triangulate()
     1157{
     1158  // The copy of Polycone is made and this copy is reordered in order to
     1159  // have a list of triangles. This list is used for GetPointOnFace().
     1160
     1161  G4PolyPhiFaceVertex *tri_help = new G4PolyPhiFaceVertex[numEdges];
     1162  triangles = tri_help;
     1163  G4PolyPhiFaceVertex *triang = triangles;
     1164
     1165  std::vector<G4double> areas;
     1166  std::vector<G4ThreeVector> points;
     1167  G4double area=0.;
     1168  G4PolyPhiFaceVertex *v0,*v1,*v2,*v3,*v4;
     1169  v2=triangles;
     1170
     1171  // Make copy for prev/next for triang=corners
     1172  //
     1173  G4PolyPhiFaceVertex *helper = corners;
     1174  G4PolyPhiFaceVertex *helper2 = corners;
     1175  do
     1176  {
     1177    triang->r = helper->r;
     1178    triang->z = helper->z;
     1179    triang->x = helper->x;
     1180    triang->y= helper->y;
     1181
     1182    // add pointer on prev corner
     1183    //
     1184    if( helper==corners )
     1185      { triang->prev=triangles+numEdges-1; }
     1186    else
     1187      { triang->prev=helper2; }
     1188
     1189    // add pointer on next corner
     1190    //
     1191    if( helper<corners+numEdges-1 )
     1192      { triang->next=triang+1; }
     1193    else
     1194      { triang->next=triangles; }
     1195    helper2=triang;
     1196    helper=helper->next;
     1197    triang=triang->next;
     1198
     1199  } while( helper!=corners );
     1200
     1201  EarInit();
     1202 
     1203  G4int n=numEdges;
     1204  G4int i=0;
     1205  G4ThreeVector p1,p2,p3,p4;
     1206  const G4int max_n_loops=numEdges*10000; // protection against infinite loop
     1207
     1208  // Each step of outer loop removes one ear
     1209  //
     1210  while(n>3)  // Inner loop searches for one ear
     1211  {
     1212    v2=triangles;
     1213    do
     1214    {
     1215      if(v2->ear)  // Ear found. Fill variables
     1216      {
     1217        // (v1,v3) is diagonal
     1218        //
     1219        v3=v2->next; v4=v3->next;
     1220        v1=v2->prev; v0=v1->prev;
     1221       
     1222        // Calculate areas and points
     1223
     1224        p1=G4ThreeVector((v2)->x,(v2)->y,(v2)->z);
     1225        p2=G4ThreeVector((v1)->x,(v1)->y,(v1)->z);
     1226        p3=G4ThreeVector((v3)->x,(v3)->y,(v3)->z);
     1227
     1228        G4double result1 = SurfaceTriangle(p1,p2,p3,&p4 );
     1229        points.push_back(p4);
     1230        areas.push_back(result1);
     1231        area=area+result1;
     1232
     1233        // Update earity of diagonal endpoints
     1234        //
     1235        v1->ear=Diagonal(v0,v3);
     1236        v3->ear=Diagonal(v1,v4);
     1237
     1238        // Cut off the ear v2
     1239        // Has to be done for a copy and not for real PolyPhiFace
     1240        //
     1241        v1->next=v3;
     1242        v3->prev=v1;
     1243        triangles=v3; // In case the head was v2
     1244        n--;
     1245 
     1246        break; // out of inner loop
     1247      }        // end if ear found
     1248
     1249      v2=v2->next;
     1250   
     1251    } while( v2!=triangles );
     1252
     1253    i++;
     1254    if(i>=max_n_loops)
     1255    {
     1256      G4Exception( "G4PolyPhiFace::Triangulation()",
     1257                   "Bad_Definition_of_Solid", FatalException,
     1258                   "Maximum number of steps is reached for triangulation!" );
     1259    }
     1260  }   // end outer while loop
     1261
     1262  if(v2->next)
     1263  {
     1264     // add last triangle
     1265     //
     1266     v2=v2->next;
     1267     p1=G4ThreeVector((v2)->x,(v2)->y,(v2)->z);
     1268     p2=G4ThreeVector((v2->next)->x,(v2->next)->y,(v2->next)->z);
     1269     p3=G4ThreeVector((v2->prev)->x,(v2->prev)->y,(v2->prev)->z);
     1270     G4double result1 = SurfaceTriangle(p1,p2,p3,&p4 );
     1271     points.push_back(p4);
     1272     areas.push_back(result1);
     1273     area=area+result1;
     1274  }
     1275 
     1276  // Surface Area is stored
     1277  //
     1278  fSurfaceArea = area;
     1279 
     1280  // Second Step: choose randomly one surface
     1281  //
     1282  G4double chose = area*G4UniformRand();
     1283   
     1284  // Third Step: Get a point on choosen surface
     1285  //
     1286  G4double Achose1, Achose2;
     1287  Achose1=0; Achose2=0.;
     1288  i=0;
     1289  do
     1290  {
     1291    Achose2+=areas[i];
     1292    if(chose>=Achose1 && chose<Achose2)
     1293    {
     1294      G4ThreeVector point;
     1295       point=points[i] ;
     1296       surface_point=point;
     1297       break;     
     1298    }
     1299    i++; Achose1=Achose2;
     1300  } while( i<numEdges-2 );
     1301 
     1302  delete [] tri_help;
     1303}

  • trunk/source/geometry/solids/specific/src/G4Polycone.cc

    r831 r850  
    2525//
    2626//
    27 // $Id: G4Polycone.cc,v 1.39 2007/10/02 09:50:46 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: G4Polycone.cc,v 1.43 2008/05/15 13:45:15 gcosmo Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//
    … …  
    781781//
    782782G4ThreeVector G4Polycone::GetPointOnSurface() const
    783 {
    784   G4double Area=0,totArea=0,Achose1=0,Achose2=0,phi,cosphi,sinphi,rRand;
    785   G4int i=0;
    786   G4int numPlanes = original_parameters->Num_z_planes;
    787  
    788   phi = RandFlat::shoot(startPhi,endPhi);
    789   cosphi = std::cos(phi);
    790   sinphi = std::sin(phi);
    791 
    792   rRand = RandFlat::shoot(original_parameters->Rmin[0],
    793                           original_parameters->Rmax[0]);
    794  
    795   std::vector<G4double> areas;       // (numPlanes+1);
    796   std::vector<G4ThreeVector> points; // (numPlanes-1);
    797  
    798   areas.push_back(pi*(sqr(original_parameters->Rmax[0])
    799                      -sqr(original_parameters->Rmin[0])));
    800 
    801   for(i=0; i<numPlanes-1; i++)
    802   {
    803     Area = (original_parameters->Rmin[i]+original_parameters->Rmin[i+1])*
    804       std::sqrt(sqr(original_parameters->Rmin[i]
    805                    -original_parameters->Rmin[i+1])+
    806                 sqr(original_parameters->Z_values[i+1]
    807                    -original_parameters->Z_values[i]));
     783{
     784  if (!genericPcon)  // Polycone by faces
     785  {
     786    G4double Area=0,totArea=0,Achose1=0,Achose2=0,phi,cosphi,sinphi,rRand;
     787    G4int i=0;
     788    G4int numPlanes = original_parameters->Num_z_planes;
     789 
     790    phi = RandFlat::shoot(startPhi,endPhi);
     791    cosphi = std::cos(phi);
     792    sinphi = std::sin(phi);
     793
     794    rRand = RandFlat::shoot(original_parameters->Rmin[0],
     795                            original_parameters->Rmax[0]);
     796 
     797    std::vector<G4double> areas;       // (numPlanes+1);
     798    std::vector<G4ThreeVector> points; // (numPlanes-1);
     799 
     800    areas.push_back(pi*(sqr(original_parameters->Rmax[0])
     801                       -sqr(original_parameters->Rmin[0])));
     802
     803    for(i=0; i<numPlanes-1; i++)
     804    {
     805      Area = (original_parameters->Rmin[i]+original_parameters->Rmin[i+1])
     806           * std::sqrt(sqr(original_parameters->Rmin[i]
     807                      -original_parameters->Rmin[i+1])+
     808                       sqr(original_parameters->Z_values[i+1]
     809                      -original_parameters->Z_values[i]));
     810
     811      Area += (original_parameters->Rmax[i]+original_parameters->Rmax[i+1])
     812            * std::sqrt(sqr(original_parameters->Rmax[i]
     813                       -original_parameters->Rmax[i+1])+
     814                        sqr(original_parameters->Z_values[i+1]
     815                       -original_parameters->Z_values[i]));
     816
     817      Area *= 0.5*(endPhi-startPhi);
    808818   
    809     Area += (original_parameters->Rmax[i]+original_parameters->Rmax[i+1])*
    810       std::sqrt(sqr(original_parameters->Rmax[i]
    811                    -original_parameters->Rmax[i+1])+
    812                 sqr(original_parameters->Z_values[i+1]
    813                    -original_parameters->Z_values[i]));
    814 
    815     Area *= 0.5*(endPhi-startPhi);
    816    
    817     if(startPhi==0.&& endPhi == twopi)
    818     {
    819       Area += std::fabs(original_parameters->Z_values[i+1]
    820                        -original_parameters->Z_values[i])*
    821                        (original_parameters->Rmax[i]
    822                        +original_parameters->Rmax[i+1]
    823                        -original_parameters->Rmin[i]
    824                        -original_parameters->Rmin[i+1]);
    825     }
    826     areas.push_back(Area);
    827     totArea += Area;
    828   }
    829  
    830   areas.push_back(pi*(sqr(original_parameters->Rmax[numPlanes-1])-
    831                       sqr(original_parameters->Rmin[numPlanes-1])));
    832  
    833   totArea += (areas[0]+areas[numPlanes]);
    834   G4double chose = RandFlat::shoot(0.,totArea);
    835 
    836   if( (chose>=0.) && (chose<areas[0]) )
    837   {
    838     return G4ThreeVector(rRand*cosphi, rRand*sinphi,
    839                          original_parameters->Z_values[0]);
    840   }
    841  
    842   for (i=0; i<numPlanes-1; i++)
    843   {
    844     Achose1 += areas[i];
    845     Achose2 = (Achose1+areas[i+1]);
    846     if(chose>=Achose1 && chose<Achose2)
    847       {// G4cout<<"will return Point On Cut"<<G4endl;
    848       return GetPointOnCut(original_parameters->Rmin[i],
    849                            original_parameters->Rmax[i],
    850                            original_parameters->Rmin[i+1],
    851                            original_parameters->Rmax[i+1],
    852                            original_parameters->Z_values[i],
    853                            original_parameters->Z_values[i+1], Area);
    854     }
    855   }
    856 
    857   rRand = RandFlat::shoot(original_parameters->Rmin[numPlanes-1],
    858                           original_parameters->Rmax[numPlanes-1]);
    859  
    860   return G4ThreeVector(rRand*cosphi,rRand*sinphi,
    861                        original_parameters->Z_values[numPlanes-1]); 
    862 }
    863 
     819      if(startPhi==0.&& endPhi == twopi)
     820      {
     821        Area += std::fabs(original_parameters->Z_values[i+1]
     822                         -original_parameters->Z_values[i])*
     823                         (original_parameters->Rmax[i]
     824                         +original_parameters->Rmax[i+1]
     825                         -original_parameters->Rmin[i]
     826                         -original_parameters->Rmin[i+1]);
     827      }
     828      areas.push_back(Area);
     829      totArea += Area;
     830    }
     831 
     832    areas.push_back(pi*(sqr(original_parameters->Rmax[numPlanes-1])-
     833                        sqr(original_parameters->Rmin[numPlanes-1])));
     834 
     835    totArea += (areas[0]+areas[numPlanes]);
     836    G4double chose = RandFlat::shoot(0.,totArea);
     837
     838    if( (chose>=0.) && (chose<areas[0]) )
     839    {
     840      return G4ThreeVector(rRand*cosphi, rRand*sinphi,
     841                           original_parameters->Z_values[0]);
     842    }
     843 
     844    for (i=0; i<numPlanes-1; i++)
     845    {
     846      Achose1 += areas[i];
     847      Achose2 = (Achose1+areas[i+1]);
     848      if(chose>=Achose1 && chose<Achose2)
     849      {
     850        return GetPointOnCut(original_parameters->Rmin[i],
     851                             original_parameters->Rmax[i],
     852                             original_parameters->Rmin[i+1],
     853                             original_parameters->Rmax[i+1],
     854                             original_parameters->Z_values[i],
     855                             original_parameters->Z_values[i+1], Area);
     856      }
     857    }
     858
     859    rRand = RandFlat::shoot(original_parameters->Rmin[numPlanes-1],
     860                            original_parameters->Rmax[numPlanes-1]);
     861 
     862    return G4ThreeVector(rRand*cosphi,rRand*sinphi,
     863                         original_parameters->Z_values[numPlanes-1]); 
     864
     865  }
     866  else  // Generic Polycone
     867  {
     868    return GetPointOnSurfaceGeneric(); 
     869  }
     870}
    864871
    865872//

  • trunk/source/geometry/solids/specific/src/G4PolyconeSide.cc

    r831 r850  
    2525//
    2626//
    27 // $Id: G4PolyconeSide.cc,v 1.17 2007/08/13 10:33:04 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: G4PolyconeSide.cc,v 1.19 2008/05/15 11:41:59 gcosmo Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//
    … …  
    4646#include "G4SolidExtentList.hh"
    4747#include "G4GeometryTolerance.hh"
     48
     49#include "Randomize.hh"
    4850
    4951//
    … …  
    6466{
    6567  kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
     68  fSurfaceArea = 0.0;
    6669
    6770  //
    … …  
    211214
    212215  kCarTolerance = source.kCarTolerance;
     216  fSurfaceArea = source.fSurfaceArea;
    213217 
    214218  cone    = new G4IntersectingCone( *source.cone );
    … …  
    10441048  y = 0.5*( y1+s1*ty1 + y2+s2*ty2 );
    10451049}
     1050
     1051//
     1052// Calculate surface area for GetPointOnSurface()
     1053//
     1054G4double G4PolyconeSide::SurfaceArea()
     1055{
     1056  if(fSurfaceArea==0)
     1057  {
     1058    fSurfaceArea = (r[0]+r[1])* std::sqrt(sqr(r[0]-r[1])+sqr(z[0]-z[1]));
     1059    fSurfaceArea *= 0.5*(deltaPhi);
     1060  } 
     1061  return fSurfaceArea;
     1062}
     1063
     1064//
     1065// GetPointOnFace
     1066//
     1067G4ThreeVector G4PolyconeSide::GetPointOnFace()
     1068{
     1069  G4double x,y,zz;
     1070  G4double rr,phi,dz,dr;
     1071  dr=r[1]-r[0];dz=z[1]-z[0];
     1072  phi=startPhi+deltaPhi*G4UniformRand();
     1073  rr=r[0]+dr*G4UniformRand();
     1074 
     1075  x=rr*std::cos(phi);
     1076  y=rr*std::sin(phi);
     1077
     1078  // PolyconeSide has a Ring Form
     1079  //
     1080  if (dz==0.)
     1081  {
     1082    zz=z[0];
     1083  }
     1084  else
     1085  {
     1086    if(dr==0.)  // PolyconeSide has a Tube Form
     1087    {
     1088      zz = z[0]+dz*G4UniformRand();
     1089    }
     1090    else
     1091    {
     1092      zz = z[0]+(rr-r[0])*dz/dr;
     1093    }
     1094  }
     1095
     1096  return G4ThreeVector(x,y,zz);
     1097}

  • trunk/source/geometry/solids/specific/src/G4Polyhedra.cc

    r831 r850  
    2525//
    2626//
    27 // $Id: G4Polyhedra.cc,v 1.36.8.1 2008/04/23 08:10:24 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     27// $Id: G4Polyhedra.cc,v 1.42 2008/05/15 13:45:15 gcosmo Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//
    … …  
    674674G4ThreeVector G4Polyhedra::GetPointOnSurface() const
    675675{
    676   G4int j, numPlanes = original_parameters->Num_z_planes, Flag=0;
    677   G4double chose, totArea=0., Achose1, Achose2,
    678            rad1, rad2, sinphi1, sinphi2, cosphi1, cosphi2;
    679   G4double a, b, l2, rang,
    680            totalPhi,ksi,
    681            area, aTop=0., aBottom=0.,zVal=0.;
    682   G4ThreeVector p0, p1, p2, p3;
    683   std::vector<G4double> aVector1;
    684   std::vector<G4double> aVector2;
    685   std::vector<G4double> aVector3;
    686 
    687    totalPhi= (phiIsOpen) ? (endPhi-startPhi) : twopi;
    688    ksi = totalPhi/numSide;
    689    G4double cosksi = std::cos(ksi/2.);
    690 
    691   // below we generate the areas relevant to our solid
    692   //
    693   for(j=0; j<numPlanes-1; j++)
    694   {
    695     a = original_parameters->Rmax[j+1];
    696     b = original_parameters->Rmax[j];
    697     l2 = sqr(original_parameters->Z_values[j]
    698             -original_parameters->Z_values[j+1]) + sqr(b-a);
    699     area = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
    700     aVector1.push_back(area);
    701   }
    702  
    703   for(j=0; j<numPlanes-1; j++)
    704   {
    705     a = original_parameters->Rmin[j+1];//*cosksi;
    706     b = original_parameters->Rmin[j];//*cosksi;
    707     l2 = sqr(original_parameters->Z_values[j]
    708             -original_parameters->Z_values[j+1]) + sqr(b-a);
    709     area = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
    710     aVector2.push_back(area);
    711   }
    712  
    713   for(j=0; j<numPlanes-1; j++)
    714   {
    715     if(phiIsOpen == true)
    716     {
    717       aVector3.push_back(0.5*(original_parameters->Rmax[j]
    718                              -original_parameters->Rmin[j]
    719                              +original_parameters->Rmax[j+1]
    720                              -original_parameters->Rmin[j+1])
    721       *std::fabs(original_parameters->Z_values[j+1]
    722                 -original_parameters->Z_values[j]));
    723     }
    724     else { aVector3.push_back(0.); }
    725   }
    726  
    727   for(j=0; j<numPlanes-1; j++)
    728   {
    729     totArea += numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
    730   }
    731  
    732   // must include top and bottom areas
    733   if(original_parameters->Rmax[numPlanes-1] != 0.)
    734   {
    735     a = original_parameters->Rmax[numPlanes-1];
    736     b = original_parameters->Rmin[numPlanes-1];
    737     l2 = sqr(a-b);
    738     aTop = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
    739   }
    740 
    741   if(original_parameters->Rmax[0] != 0.)
    742   {
    743     a = original_parameters->Rmax[0];
    744     b = original_parameters->Rmin[0];
    745     l2 = sqr(a-b);
    746     aBottom = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
    747   }
    748 
    749   Achose1 = 0.;
    750   Achose2 = numSide*(aVector1[0]+aVector2[0])+2.*aVector3[0];
    751 
    752   chose = RandFlat::shoot(0.,totArea+aTop+aBottom);
    753   if( (chose >= 0.) && (chose < aTop + aBottom) )
    754   { 
     676  if( !genericPgon )  // Polyhedra by faces
     677  {
     678    G4int j, numPlanes = original_parameters->Num_z_planes, Flag=0;
     679    G4double chose, totArea=0., Achose1, Achose2,
     680             rad1, rad2, sinphi1, sinphi2, cosphi1, cosphi2;
     681    G4double a, b, l2, rang, totalPhi, ksi,
     682             area, aTop=0., aBottom=0., zVal=0.;
     683
     684    G4ThreeVector p0, p1, p2, p3;
     685    std::vector<G4double> aVector1;
     686    std::vector<G4double> aVector2;
     687    std::vector<G4double> aVector3;
     688
     689    totalPhi= (phiIsOpen) ? (endPhi-startPhi) : twopi;
     690    ksi = totalPhi/numSide;
     691    G4double cosksi = std::cos(ksi/2.);
     692
     693    // Below we generate the areas relevant to our solid
     694    //
     695    for(j=0; j<numPlanes-1; j++)
     696    {
     697      a = original_parameters->Rmax[j+1];
     698      b = original_parameters->Rmax[j];
     699      l2 = sqr(original_parameters->Z_values[j]
     700              -original_parameters->Z_values[j+1]) + sqr(b-a);
     701      area = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
     702      aVector1.push_back(area);
     703    }
     704
     705    for(j=0; j<numPlanes-1; j++)
     706    {
     707      a = original_parameters->Rmin[j+1];//*cosksi;
     708      b = original_parameters->Rmin[j];//*cosksi;
     709      l2 = sqr(original_parameters->Z_values[j]
     710              -original_parameters->Z_values[j+1]) + sqr(b-a);
     711      area = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
     712      aVector2.push_back(area);
     713    }
     714 
     715    for(j=0; j<numPlanes-1; j++)
     716    {
     717      if(phiIsOpen == true)
     718      {
     719        aVector3.push_back(0.5*(original_parameters->Rmax[j]
     720                               -original_parameters->Rmin[j]
     721                               +original_parameters->Rmax[j+1]
     722                               -original_parameters->Rmin[j+1])
     723        *std::fabs(original_parameters->Z_values[j+1]
     724                  -original_parameters->Z_values[j]));
     725      }
     726      else { aVector3.push_back(0.); }
     727    }
     728 
     729    for(j=0; j<numPlanes-1; j++)
     730    {
     731      totArea += numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
     732    }
     733 
     734    // Must include top and bottom areas
     735    //
     736    if(original_parameters->Rmax[numPlanes-1] != 0.)
     737    {
     738      a = original_parameters->Rmax[numPlanes-1];
     739      b = original_parameters->Rmin[numPlanes-1];
     740      l2 = sqr(a-b);
     741      aTop = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
     742    }
     743
     744    if(original_parameters->Rmax[0] != 0.)
     745    {
     746      a = original_parameters->Rmax[0];
     747      b = original_parameters->Rmin[0];
     748      l2 = sqr(a-b);
     749      aBottom = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
     750    }
     751
     752    Achose1 = 0.;
     753    Achose2 = numSide*(aVector1[0]+aVector2[0])+2.*aVector3[0];
     754
     755    chose = RandFlat::shoot(0.,totArea+aTop+aBottom);
     756    if( (chose >= 0.) && (chose < aTop + aBottom) )
     757    {
     758      chose = RandFlat::shoot(startPhi,startPhi+totalPhi);
     759      rang = std::floor((chose-startPhi)/ksi-0.01);
     760      if(rang<0) { rang=0; }
     761      rang = std::fabs(rang); 
     762      sinphi1 = std::sin(startPhi+rang*ksi);
     763      sinphi2 = std::sin(startPhi+(rang+1)*ksi);
     764      cosphi1 = std::cos(startPhi+rang*ksi);
     765      cosphi2 = std::cos(startPhi+(rang+1)*ksi);
     766      chose = RandFlat::shoot(0., aTop + aBottom);
     767      if(chose>=0. && chose<aTop)
     768      {
     769        rad1 = original_parameters->Rmin[numPlanes-1];
     770        rad2 = original_parameters->Rmax[numPlanes-1];
     771        zVal = original_parameters->Z_values[numPlanes-1];
     772      }
     773      else
     774      {
     775        rad1 = original_parameters->Rmin[0];
     776        rad2 = original_parameters->Rmax[0];
     777        zVal = original_parameters->Z_values[0];
     778      }
     779      p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
     780      p1 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
     781      p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
     782      p3 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
     783      return GetPointOnPlane(p0,p1,p2,p3);
     784    }
     785    else
     786    {
     787      for (j=0; j<numPlanes-1; j++)
     788      {
     789        if( ((chose >= Achose1) && (chose < Achose2)) || (j == numPlanes-2) )
     790        {
     791          Flag = j; break;
     792        }
     793        Achose1 += numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
     794        Achose2 = Achose1 + numSide*(aVector1[j+1]+aVector2[j+1])
     795                          + 2.*aVector3[j+1];
     796      }
     797    }
     798
     799    // At this point we have chosen a subsection
     800    // between to adjacent plane cuts...
     801
     802    j = Flag;
    755803   
    756     chose = RandFlat::shoot(startPhi,startPhi+totalPhi);
    757     rang = std::floor((chose-startPhi)/ksi-0.01);
    758     if(rang<0)rang=0;
    759     rang = std::fabs(rang); 
    760     sinphi1 = std::sin(startPhi+rang*ksi);
    761     sinphi2 = std::sin(startPhi+(rang+1)*ksi);
    762     cosphi1 = std::cos(startPhi+rang*ksi);
    763     cosphi2 = std::cos(startPhi+(rang+1)*ksi);
    764     chose = RandFlat::shoot(0., aTop + aBottom);
    765     if(chose>=0. && chose<aTop)
    766     {
    767       rad1 = original_parameters->Rmin[numPlanes-1];
    768       rad2 = original_parameters->Rmax[numPlanes-1];
    769       zVal = original_parameters->Z_values[numPlanes-1];
    770     }
    771     else
    772     {
    773       rad1 = original_parameters->Rmin[0];
    774       rad2 = original_parameters->Rmax[0];
    775       zVal = original_parameters->Z_values[0];
    776     }
    777     p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
    778     p1 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
    779     p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
    780     p3 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
    781     return GetPointOnPlane(p0,p1,p2,p3);
    782   }
    783   else
    784   {
    785     for (j=0; j<numPlanes-1; j++)
    786     {
    787       if( ((chose >= Achose1) && (chose < Achose2)) || (j == numPlanes-2) )
    788       {
    789         Flag = j; break;
    790       }
    791       Achose1 += numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
    792       Achose2 = Achose1 + numSide*(aVector1[j+1]+aVector2[j+1])
    793                         + 2.*aVector3[j+1];
    794     }
    795   }
    796 
    797   // at this point we have chosen a subsection
    798   // between to adjacent plane cuts...
    799 
    800   j = Flag;
     804    totArea = numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
     805    chose = RandFlat::shoot(0.,totArea);
     806 
     807    if( (chose>=0.) && (chose<numSide*aVector1[j]) )
     808    {
     809      chose = RandFlat::shoot(startPhi,startPhi+totalPhi);
     810      rang = std::floor((chose-startPhi)/ksi-0.01);
     811      if(rang<0) { rang=0; }
     812      rang = std::fabs(rang);
     813      rad1 = original_parameters->Rmax[j];
     814      rad2 = original_parameters->Rmax[j+1];
     815      sinphi1 = std::sin(startPhi+rang*ksi);
     816      sinphi2 = std::sin(startPhi+(rang+1)*ksi);
     817      cosphi1 = std::cos(startPhi+rang*ksi);
     818      cosphi2 = std::cos(startPhi+(rang+1)*ksi);
     819      zVal = original_parameters->Z_values[j];
    801820   
    802   totArea = numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
    803   chose = RandFlat::shoot(0.,totArea);
    804  
    805   if( (chose>=0.) && (chose<numSide*aVector1[j]) )
    806   {
    807     chose = RandFlat::shoot(startPhi,startPhi+totalPhi);
    808     rang = std::floor((chose-startPhi)/ksi-0.01);
    809     if(rang<0)rang=0;
    810     rang = std::fabs(rang);
    811     rad1 = original_parameters->Rmax[j];
    812     rad2 = original_parameters->Rmax[j+1];
    813     sinphi1 = std::sin(startPhi+rang*ksi);
    814     sinphi2 = std::sin(startPhi+(rang+1)*ksi);
    815     cosphi1 = std::cos(startPhi+rang*ksi);
    816     cosphi2 = std::cos(startPhi+(rang+1)*ksi);
    817     zVal = original_parameters->Z_values[j];
     821      p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
     822      p1 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
     823
     824      zVal = original_parameters->Z_values[j+1];
     825
     826      p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
     827      p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
     828      return GetPointOnPlane(p0,p1,p2,p3);
     829    }
     830    else if ( (chose >= numSide*aVector1[j])
     831           && (chose <= numSide*(aVector1[j]+aVector2[j])) )
     832    {
     833      chose = RandFlat::shoot(startPhi,startPhi+totalPhi);
     834      rang = std::floor((chose-startPhi)/ksi-0.01);
     835      if(rang<0) { rang=0; }
     836      rang = std::fabs(rang);
     837      rad1 = original_parameters->Rmin[j];
     838      rad2 = original_parameters->Rmin[j+1];
     839      sinphi1 = std::sin(startPhi+rang*ksi);
     840      sinphi2 = std::sin(startPhi+(rang+1)*ksi);
     841      cosphi1 = std::cos(startPhi+rang*ksi);
     842      cosphi2 = std::cos(startPhi+(rang+1)*ksi);
     843      zVal = original_parameters->Z_values[j];
     844
     845      p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
     846      p1 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
     847
     848      zVal = original_parameters->Z_values[j+1];
    818849   
    819     p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
    820     p1 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
    821 
    822     zVal = original_parameters->Z_values[j+1];
    823 
    824     p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
    825     p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
     850      p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
     851      p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
     852      return GetPointOnPlane(p0,p1,p2,p3);
     853    }
     854
     855    chose = RandFlat::shoot(0.,2.2);
     856    if( (chose>=0.) && (chose < 1.) )
     857    {
     858      rang = startPhi;
     859    }
     860    else
     861    {
     862      rang = endPhi;
     863    }
     864
     865    cosphi1 = std::cos(rang); rad1 = original_parameters->Rmin[j];
     866    sinphi1 = std::sin(rang); rad2 = original_parameters->Rmax[j];
     867
     868    p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,
     869                       original_parameters->Z_values[j]);
     870    p1 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,
     871                       original_parameters->Z_values[j]);
     872
     873    rad1 = original_parameters->Rmax[j+1];
     874    rad2 = original_parameters->Rmin[j+1];
     875
     876    p2 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,
     877                       original_parameters->Z_values[j+1]);
     878    p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,
     879                       original_parameters->Z_values[j+1]);
    826880    return GetPointOnPlane(p0,p1,p2,p3);
    827881  }
    828   else if ( (chose >= numSide*aVector1[j])
    829          && (chose <= numSide*(aVector1[j]+aVector2[j])) )
    830   {
    831    
    832     chose = RandFlat::shoot(startPhi,startPhi+totalPhi);
    833     rang = std::floor((chose-startPhi)/ksi-0.01);
    834     if(rang<0)rang=0;
    835     rang = std::fabs(rang);
    836     rad1 = original_parameters->Rmin[j];
    837     rad2 = original_parameters->Rmin[j+1];
    838     sinphi1 = std::sin(startPhi+rang*ksi);
    839     sinphi2 = std::sin(startPhi+(rang+1)*ksi);
    840     cosphi1 = std::cos(startPhi+rang*ksi);
    841     cosphi2 = std::cos(startPhi+(rang+1)*ksi);
    842     zVal = original_parameters->Z_values[j];
    843    
    844     p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
    845     p1 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
    846 
    847     zVal = original_parameters->Z_values[j+1];
    848    
    849     p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
    850     p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
    851     return GetPointOnPlane(p0,p1,p2,p3);
    852   }
    853 
    854   chose = RandFlat::shoot(0.,2.2);
    855   if( (chose>=0.) && (chose < 1.) )
    856   {
    857     rang = startPhi;
    858   }
    859   else
    860   {
    861     rang = endPhi;
    862   }
    863 
    864   cosphi1 = std::cos(rang); rad1 = original_parameters->Rmin[j];
    865   sinphi1 = std::sin(rang); rad2 = original_parameters->Rmax[j];
    866    
    867   p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,
    868                      original_parameters->Z_values[j]);
    869   p1 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,
    870                      original_parameters->Z_values[j]);
    871    
    872   rad1 = original_parameters->Rmax[j+1];
    873   rad2 = original_parameters->Rmin[j+1];
    874      
    875   p2 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,
    876                      original_parameters->Z_values[j+1]);
    877   p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,
    878                      original_parameters->Z_values[j+1]);
    879   return GetPointOnPlane(p0,p1,p2,p3);
    880 }
    881 
     882  else  // Generic polyhedra
     883  {
     884    return GetPointOnSurfaceGeneric();
     885  }
     886}
    882887
    883888//

  • trunk/source/geometry/solids/specific/src/G4PolyhedraSide.cc

    r831 r850  
    2525//
    2626//
    27 // $Id: G4PolyhedraSide.cc,v 1.13 2007/05/31 13:52:48 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: G4PolyhedraSide.cc,v 1.15 2008/05/15 11:41:59 gcosmo Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//
    … …  
    3535// G4PolyhedraSide.cc
    3636//
    37 // Implemenation of the face representing one segmented side of a Polyhedra
     37// Implementation of the face representing one segmented side of a Polyhedra
    3838//
    3939// --------------------------------------------------------------------
    … …  
    4545#include "G4SolidExtentList.hh"
    4646#include "G4GeometryTolerance.hh"
     47
     48#include "Randomize.hh"
    4749
    4850//
    … …  
    6466
    6567  kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
    66 
     68  fSurfaceArea=0.;
    6769  //
    6870  // Record values
    … …  
    359361
    360362  kCarTolerance = source.kCarTolerance;
    361  
     363  fSurfaceArea = source.fSurfaceArea;
     364
    362365  cone = new G4IntersectingCone( *source.cone );
    363366
    … …  
    11401143  return std::sqrt( distFaceNorm*distFaceNorm + distOut2 );
    11411144}
     1145
     1146
     1147//
     1148// Calculation of surface area of a triangle.
     1149// At the same time a random point in the triangle is given
     1150//
     1151G4double G4PolyhedraSide::SurfaceTriangle( G4ThreeVector p1,
     1152                                           G4ThreeVector p2,
     1153                                           G4ThreeVector p3,
     1154                                           G4ThreeVector *p4 )
     1155{
     1156  G4ThreeVector v, w;
     1157 
     1158  v = p3 - p1;
     1159  w = p1 - p2;
     1160  G4double lambda1 = G4UniformRand();
     1161  G4double lambda2 = lambda1*G4UniformRand();
     1162 
     1163  *p4=p2 + lambda1*w + lambda2*v;
     1164  return 0.5*(v.cross(w)).mag();
     1165}
     1166
     1167
     1168//
     1169// GetPointOnPlane
     1170//
     1171// Auxiliary method for GetPointOnSurface()
     1172//
     1173G4ThreeVector
     1174G4PolyhedraSide::GetPointOnPlane( G4ThreeVector p0, G4ThreeVector p1,
     1175                                  G4ThreeVector p2, G4ThreeVector p3,
     1176                                  G4double *Area )
     1177{
     1178  G4double chose,aOne,aTwo;
     1179  G4ThreeVector point1,point2;
     1180  aOne = SurfaceTriangle(p0,p1,p2,&point1);
     1181  aTwo = SurfaceTriangle(p2,p3,p0,&point2);
     1182  *Area= aOne+aTwo;
     1183
     1184  chose = G4UniformRand()*(aOne+aTwo);
     1185  if( (chose>=0.) && (chose < aOne) )
     1186  {
     1187   return (point1);   
     1188  }
     1189  return (point2);
     1190}
     1191
     1192
     1193//
     1194// SurfaceArea()
     1195//
     1196G4double G4PolyhedraSide::SurfaceArea()
     1197{
     1198  if( fSurfaceArea==0. )
     1199  {
     1200    // Define the variables
     1201    //
     1202    G4double area,areas;
     1203    G4ThreeVector point1;
     1204    G4ThreeVector v1,v2,v3,v4;
     1205    G4PolyhedraSideVec *vec = vecs;
     1206    areas=0.;
     1207
     1208    // Do a loop on all SideEdge
     1209    //
     1210    do
     1211    {
     1212      // Define 4points for a Plane or Triangle
     1213      //
     1214      G4ThreeVector v1=vec->edges[0]->corner[0];
     1215      G4ThreeVector v2=vec->edges[0]->corner[1];
     1216      G4ThreeVector v3=vec->edges[1]->corner[1];
     1217      G4ThreeVector v4=vec->edges[1]->corner[0];
     1218      point1=GetPointOnPlane(v1,v2,v3,v4,&area);
     1219      areas+=area;
     1220    } while( ++vec < vecs + numSide);
     1221
     1222    fSurfaceArea=areas;
     1223  }
     1224  return fSurfaceArea;
     1225}
     1226
     1227
     1228//
     1229// GetPointOnFace()
     1230//
     1231G4ThreeVector G4PolyhedraSide::GetPointOnFace()
     1232{
     1233  // Define the variables
     1234  //
     1235  std::vector<G4double>areas;
     1236  std::vector<G4ThreeVector>points;
     1237  G4double area=0;
     1238  G4double result1;
     1239  G4ThreeVector point1;
     1240  G4ThreeVector v1,v2,v3,v4;
     1241  G4PolyhedraSideVec *vec = vecs;
     1242
     1243  // Do a loop on all SideEdge
     1244  //
     1245  do
     1246  {
     1247    // Define 4points for a Plane or Triangle
     1248    //
     1249    G4ThreeVector v1=vec->edges[0]->corner[0];
     1250    G4ThreeVector v2=vec->edges[0]->corner[1];
     1251    G4ThreeVector v3=vec->edges[1]->corner[1];
     1252    G4ThreeVector v4=vec->edges[1]->corner[0];
     1253    point1=GetPointOnPlane(v1,v2,v3,v4,&result1);
     1254    points.push_back(point1);
     1255    areas.push_back(result1);
     1256    area+=result1;
     1257  } while( ++vec < vecs+numSide );
     1258
     1259  // Choose randomly one of the surfaces and point on it
     1260  //
     1261  G4double chose = area*G4UniformRand();
     1262  G4double Achose1,Achose2;
     1263  Achose1=0;Achose2=0.;
     1264  G4int i=0;
     1265  do
     1266  {
     1267    Achose2+=areas[i];
     1268    if(chose>=Achose1 && chose<Achose2)
     1269    {
     1270      point1=points[i] ; break;     
     1271    }
     1272    i++; Achose1=Achose2;
     1273  } while( i<numSide );
     1274 
     1275  return point1;
     1276}

  • trunk/source/geometry/solids/specific/src/G4QuadrangularFacet.cc

    r831 r850  
    2626//
    2727// $Id: G4QuadrangularFacet.cc,v 1.6 2007/08/23 14:49:23 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  • trunk/source/geometry/solids/specific/src/G4ReduciblePolygon.cc

    r831 r850  
    2626//
    2727// $Id: G4ReduciblePolygon.cc,v 1.11 2006/06/29 18:48:53 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/specific/src/G4SolidExtentList.cc

    r831 r850  
    2626//
    2727// $Id: G4SolidExtentList.cc,v 1.5 2007/05/11 13:54:29 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/specific/src/G4TessellatedGeometryAlgorithms.cc

    r831 r850  
    2525//
    2626// $Id: G4TessellatedGeometryAlgorithms.cc,v 1.5 2007/12/12 16:51:12 gcosmo Exp $
    27 // GEANT4 tag $Name: $
     27// GEANT4 tag $Name: HEAD $
    2828//
    2929// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  • trunk/source/geometry/solids/specific/src/G4TessellatedSolid.cc

    r831 r850  
    2525// ********************************************************************
    2626//
    27 // $Id: G4TessellatedSolid.cc,v 1.14.2.1 2008/04/23 08:10:24 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     27// $Id: G4TessellatedSolid.cc,v 1.18 2008/03/13 11:58:28 gcosmo Exp $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  • trunk/source/geometry/solids/specific/src/G4Tet.cc

    r831 r850  
    2929//
    3030// $Id: G4Tet.cc,v 1.11 2006/11/13 08:58:03 gcosmo Exp $
    31 // GEANT4 tag $Name: $
     31// GEANT4 tag $Name: HEAD $
    3232//
    3333// class G4Tet

  • trunk/source/geometry/solids/specific/src/G4TriangularFacet.cc

    r831 r850  
    2626//
    2727// $Id: G4TriangularFacet.cc,v 1.10 2007/12/10 16:30:35 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  • trunk/source/geometry/solids/specific/src/G4TwistBoxSide.cc

    r831 r850  
    2626//
    2727// $Id: G4TwistBoxSide.cc,v 1.6 2007/05/23 09:31:02 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/specific/src/G4TwistTrapAlphaSide.cc

    r831 r850  
    2626//
    2727// $Id: G4TwistTrapAlphaSide.cc,v 1.8 2007/05/23 13:26:06 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/specific/src/G4TwistTrapFlatSide.cc

    r831 r850  
    2626//
    2727// $Id: G4TwistTrapFlatSide.cc,v 1.6 2007/05/23 09:31:02 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/specific/src/G4TwistTrapParallelSide.cc

    r831 r850  
    2626//
    2727// $Id: G4TwistTrapParallelSide.cc,v
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/specific/src/G4TwistTubsFlatSide.cc

    r831 r850  
    2626//
    2727// $Id: G4TwistTubsFlatSide.cc,v 1.7 2007/05/23 09:31:02 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/specific/src/G4TwistTubsHypeSide.cc

    r831 r850  
    2626//
    2727// $Id: G4TwistTubsHypeSide.cc,v 1.6 2007/05/18 07:39:56 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/specific/src/G4TwistTubsSide.cc

    r831 r850  
    2626//
    2727// $Id: G4TwistTubsSide.cc,v 1.5 2006/06/29 18:49:18 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/specific/src/G4TwistedBox.cc

    r831 r850  
    2626//
    2727// $Id: G4TwistedBox.cc,v 1.12 2006/06/29 18:49:20 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/specific/src/G4TwistedTrap.cc

    r831 r850  
    2626//
    2727// $Id: G4TwistedTrap.cc,v 1.14 2006/06/29 18:49:23 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/specific/src/G4TwistedTrd.cc

    r831 r850  
    2626//
    2727// $Id: G4TwistedTrd.cc,v 1.7 2006/06/29 18:49:25 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/specific/src/G4TwistedTubs.cc

    r831 r850  
    2626//
    2727// $Id: G4TwistedTubs.cc,v 1.24 2007/05/18 07:39:56 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/specific/src/G4VCSGfaceted.cc

    r831 r850  
    3030// and all its terms.
    3131//
    32 // $Id: G4VCSGfaceted.cc,v 1.20 2006/10/20 14:21:36 gcosmo Exp $
    33 // GEANT4 tag $Name: $
     32// $Id: G4VCSGfaceted.cc,v 1.25 2008/05/22 10:22:52 gcosmo Exp $
     33// GEANT4 tag $Name: HEAD $
    3434//
    3535//
    … …  
    5151#include "G4VoxelLimits.hh"
    5252#include "G4AffineTransform.hh"
     53
     54#include "Randomize.hh"
    5355
    5456#include "G4Polyhedron.hh"   
    … …  
    105107const G4VCSGfaceted &G4VCSGfaceted::operator=( const G4VCSGfaceted &source )
    106108{
    107   if (&source == this) return *this;
     109  if (&source == this) { return *this; }
    108110 
    109111  DeleteStuff();
    … …  
    122124{
    123125  numFace = source.numFace;
    124   if (numFace == 0) return;    // odd, but permissable?
     126  if (numFace == 0) { return; }    // odd, but permissable?
    125127 
    126128  faces = new G4VCSGface*[numFace];
    … …  
    128130  G4VCSGface **face = faces,
    129131       **sourceFace = source.faces;
    130   do {
     132  do
     133  {
    131134    *face = (*sourceFace)->Clone();
    132135  } while( ++sourceFace, ++face < faces+numFace );
    … …  
    146149  {
    147150    G4VCSGface **face = faces;
    148     do {
     151    do
     152    {
    149153      delete *face;
    150154    } while( ++face < faces + numFace );
    … …  
    170174  //
    171175  G4VCSGface **face = faces;
    172   do {
     176  do
     177  {
    173178    (*face)->CalculateExtent( axis, voxelLimit, transform, extentList );
    174179  } while( ++face < faces + numFace );
    … …  
    194199  G4VCSGface **face = faces;
    195200  G4double best = kInfinity;
    196   do {
     201  do
     202  {
    197203    G4double distance;
    198204    EInside result = (*face)->Inside( p, kCarTolerance/2, &distance );
    199     if (result == kSurface) return kSurface;
     205    if (result == kSurface) { return kSurface; }
    200206    if (distance < best)
    201207    {
    … …  
    217223  G4VCSGface **face = faces;
    218224  G4double best = kInfinity;
    219   do {
     225  do
     226  {
    220227    G4double distance;
    221228    G4ThreeVector normal = (*face)->Normal( p, &distance );
    … …  
    241248  G4VCSGface *bestFace=0;
    242249  G4VCSGface **face = faces;
    243   do {
     250  do
     251  {
    244252    G4double   faceDistance,
    245253               faceDistFromSurface;
    … …  
    258266        distFromSurface = faceDistFromSurface;
    259267        bestFace = *face;
    260         if (distFromSurface <= 0) return 0;
     268        if (distFromSurface <= 0) { return 0; }
    261269      }
    262270    }
    … …  
    265273  if (distance < kInfinity && distFromSurface<kCarTolerance/2)
    266274  {
    267     if (bestFace->Distance(p,false) < kCarTolerance/2) distance = 0;
     275    if (bestFace->Distance(p,false) < kCarTolerance/2)  { distance = 0; }
    268276  }
    269277
    … …  
    297305 
    298306  G4VCSGface **face = faces;
    299   do {
     307  do
     308  {
    300309    G4double  faceDistance,
    301310              faceDistFromSurface;
    … …  
    309318      // Intersecting face
    310319      //
    311       if ( (distance < kInfinity) || (!faceAllBehind) ) allBehind = false;
     320      if ( (distance < kInfinity) || (!faceAllBehind) )  { allBehind = false; }
    312321      if (faceDistance < distance)
    313322      {
    … …  
    316325        normal = faceNormal;
    317326        bestFace = *face;
    318         if (distFromSurface <= 0) break;
     327        if (distFromSurface <= 0)  { break; }
    319328      }
    320329    }
    … …  
    324333  {
    325334    if (distFromSurface <= 0)
     335    {
    326336      distance = 0;
     337    }
    327338    else if (distFromSurface<kCarTolerance/2)
    328339    {
    329       if (bestFace->Distance(p,true) < kCarTolerance/2) distance = 0;
     340      if (bestFace->Distance(p,true) < kCarTolerance/2)  { distance = 0; }
    330341    }
    331342
    … …  
    338349  else
    339350  {
    340     if (Inside(p) == kSurface) distance = 0;
    341     if (calcNorm) *validNorm = false;
     351    if (Inside(p) == kSurface)  { distance = 0; }
     352    if (calcNorm)  { *validNorm = false; }
    342353  }
    343354
    … …  
    365376  G4VCSGface **face = faces;
    366377  G4double best = kInfinity;
    367   do {
     378  do
     379  {
    368380    G4double distance = (*face)->Distance( p, outgoing );
    369     if (distance < best) best = distance;
     381    if (distance < best)  { best = distance; }
    370382  } while( ++face < faces + numFace );
    371383
    … …  
    400412
    401413  G4VCSGface **face = faces;
    402   do {   
     414  do
     415  {   
    403416    const G4ThreeVector **axis = axes+5 ;
    404417    G4double *answer = answers+5;
    405     do {
     418    do
     419    {
    406420      G4double testFace = (*face)->Extent( **axis );
    407       if (testFace > *answer) *answer = testFace;
     421      if (testFace > *answer)  { *answer = testFace; }
    408422    }
    409423    while( --axis, --answer >= answers );
    … …  
    412426 
    413427    return G4VisExtent( -answers[0], answers[1],
    414           -answers[2], answers[3],
    415           -answers[4], answers[5]  );
     428                        -answers[2], answers[3],
     429                        -answers[4], answers[5]  );
    416430}
    417431
    … …  
    524538G4double G4VCSGfaceted::GetCubicVolume()
    525539{
    526   if(fCubicVolume != 0.) ;
    527   else   fCubicVolume = EstimateCubicVolume(fStatistics,fCubVolEpsilon);
     540  if(fCubicVolume != 0.) {;}
     541  else   { fCubicVolume = EstimateCubicVolume(fStatistics,fCubVolEpsilon); }
    528542  return fCubicVolume;
    529543}
    … …  
    535549G4double G4VCSGfaceted::GetSurfaceArea()
    536550{
    537   if(fSurfaceArea != 0.) ;
    538   else   fSurfaceArea = EstimateCubicVolume(fStatistics,fAreaAccuracy);
     551  if(fSurfaceArea != 0.) {;}
     552  else   { fSurfaceArea = EstimateCubicVolume(fStatistics,fAreaAccuracy); }
    539553  return fSurfaceArea;
    540554}
    … …  
    549563      fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
    550564      fpPolyhedron->GetNumberOfRotationSteps())
    551     {
    552       delete fpPolyhedron;
    553       fpPolyhedron = CreatePolyhedron();
    554     }
     565  {
     566    delete fpPolyhedron;
     567    fpPolyhedron = CreatePolyhedron();
     568  }
    555569  return fpPolyhedron;
    556570}
     571
     572
     573//
     574// GetPointOnSurfaceGeneric proportional to Areas of faces
     575// in case of GenericPolycone or GenericPolyhedra
     576//
     577G4ThreeVector G4VCSGfaceted::GetPointOnSurfaceGeneric( ) const
     578{
     579  // Preparing variables
     580  //
     581  G4ThreeVector answer=G4ThreeVector(0.,0.,0.);
     582  G4VCSGface **face = faces;
     583  G4double area = 0;
     584  G4int i;
     585  std::vector<G4double> areas;
     586
     587  // First step: calculate surface areas
     588  //
     589  do
     590  {
     591    G4double result = (*face)->SurfaceArea( );
     592    areas.push_back(result);
     593    area=area+result;
     594  } while( ++face < faces + numFace );
     595
     596  // Second Step: choose randomly one surface
     597  //
     598  G4VCSGface **face1 = faces;
     599  G4double chose = area*G4UniformRand();
     600  G4double Achose1, Achose2;
     601  Achose1=0; Achose2=0.;
     602  i=0;
     603
     604  do
     605  {
     606    Achose2+=areas[i];
     607    if(chose>=Achose1 && chose<Achose2)
     608    {
     609      G4ThreeVector point;
     610      point= (*face1)->GetPointOnFace();
     611      return point;
     612    }
     613    i++;
     614    Achose1=Achose2;
     615  } while( ++face1 < faces + numFace );
     616
     617  return answer;
     618}

  • trunk/source/geometry/solids/specific/src/G4VFacet.cc

    r831 r850  
    2626//
    2727// $Id: G4VFacet.cc,v 1.6 2007/08/23 14:45:03 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  • trunk/source/geometry/solids/specific/src/G4VTwistSurface.cc

    r831 r850  
    2626//
    2727// $Id: G4VTwistSurface.cc,v 1.9 2007/05/31 13:52:48 gcosmo Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: HEAD $
    2929//
    3030//

  • trunk/source/geometry/solids/specific/src/G4VTwistedFaceted.cc

    r831 r850  
    2525//
    2626// $Id: G4VTwistedFaceted.cc,v 1.18 2007/05/25 09:42:34 gcosmo Exp $
    27 // GEANT4 tag $Name: $
     27// GEANT4 tag $Name: HEAD $
    2828//
    2929//
Note: See TracChangeset for help on using the changeset viewer.