Changeset 1315 for trunk/source/geometry


Ignore:
Timestamp:
Jun 18, 2010, 11:42:07 AM (15 years ago)
Author:
garnier
Message:

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08

Location:
trunk/source/geometry
Files:
31 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/geometry/magneticfield/History

    r1231 r1315  
    1 $Id: History,v 1.147 2009/11/12 18:45:02 japost Exp $
     1$Id: History,v 1.146 2009/11/12 15:05:49 japost Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1818     ----------------------------------------------------------
    1919
    20 Nov 12th,  2009 J.Apostolakis - field-V09-02-09
    21 -----------------------------
    22 - G4MagIntegratorDriver:  Activate Peter Gumplinger's check on integration
    23                           error for spin.
    24 
    2520Nov 12th,  2009 J.Apostolakis - field-V09-02-08
    26 -----------------------------  (fix only)
     21-----------------------------
    2722- G4Nystrom:  Corrected interface method getField: array now has explicit dimension[4]
    2823                (Problem found by gcc 4.3 - it checked indices used in inline method! )

  • trunk/source/geometry/magneticfield/src/G4MagIntegratorDriver.cc

    r1231 r1315  
    2525//
    2626//
    27 // $Id: G4MagIntegratorDriver.cc,v 1.56 2009/11/12 18:41:03 japost Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: G4MagIntegratorDriver.cc,v 1.53 2009/11/05 22:31:43 japost Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030//
    … …  
    208208
    209209#ifdef G4DEBUG_FIELD
    210     G4double xSubStepStart= x;
    211210    for (i=0;i<nvar;i++)  { ySubStepStart[i] = y[i]; }
    212211    yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
    … …  
    228227      lastStepSucceeded= (hdid == h);   
    229228#ifdef G4DEBUG_FIELD
    230       if (dbg>2) {
    231         PrintStatus( ySubStepStart, xSubStepStart, y, x, h,  nstp); // Only
     229      if (dbg>2)
     230      {
     231        PrintStatus( ySubStepStart, xSubStart, y, x, h,  nstp); // Only
    232232      }
    233233#endif
    … …  
    291291    {
    292292      if( nstp==nStpPr )  { G4cout << "***** Many steps ****" << G4endl; }
    293       G4cout << "MagIntDrv: " ;
    294293      G4cout << "hdid="  << std::setw(12) << hdid  << " "
    295              << "hnext=" << std::setw(12) << hnext << " "
    296              << "hstep=" << std::setw(12) << hstep << " (requested) "
    297              << G4endl;
     294             << "hnext=" << std::setw(12) << hnext << " " << G4endl;
    298295      PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp);
    299296    }
    … …  
    372369        lastStep = true;
    373370#ifdef G4DEBUG_FIELD
    374         if (dbg>2)
     371        if (dbg)
    375372        {
    376           int prec= G4cout.precision(12);
    377373          G4cout << "Warning: G4MagIntegratorDriver::AccurateAdvance"
    378374                 << G4endl
    379375                 << "  Integration step 'h' became "
    380                  << h << " due to roundoff. " << G4endl
    381                  << " Calculated as difference of x2= "<< x2 << " and x=" << x
     376                 << h << " due to roundoff " << G4endl
    382377                 << "  Forcing termination of advance." << G4endl;
    383           G4cout.precision(prec);
    384378        }         
    385379#endif
    … …  
    587581                 /  ( sqr(y[9]) + sqr(y[10]) + sqr(y[11]) );
    588582      errspin_sq *= inv_eps_vel_sq;
    589       errmax_sq = std::max( errmax_sq, errspin_sq );
    590    }
     583    }
    591584
    592585    if ( errmax_sq <= 1.0 )  { break; } // Step succeeded.

  • trunk/source/geometry/management/History

    r1228 r1315  
    1 $Id: History,v 1.138 2009/11/27 16:34:53 gcosmo Exp $
     1$Id: History,v 1.139 2010/02/09 16:50:44 gcosmo Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1717     * Reverse chronological order (last date on top), please *
    1818     ----------------------------------------------------------
     19
     20February 9, 2010  G.Cosmo                  geommng-V09-03-00
     21- Fixed initialisation of min/max extent for mother and target volumes in
     22  method BuildNodes() of G4SmartVoxelHeader.
    1923
    2024November 27, 2009  G.Cosmo                 geommng-V09-02-05

  • trunk/source/geometry/management/src/G4SmartVoxelHeader.cc

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4SmartVoxelHeader.cc,v 1.34 2009/10/30 14:05:47 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4SmartVoxelHeader.cc,v 1.35 2010/02/09 16:50:25 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030//
    … …  
    761761                                              EAxis pAxis)
    762762{
    763   G4double motherMinExtent, motherMaxExtent, targetMinExtent, targetMaxExtent;
     763  G4double motherMinExtent= kInfinity, motherMaxExtent= -kInfinity,
     764           targetMinExtent= kInfinity, targetMaxExtent= -kInfinity;
    764765  G4VPhysicalVolume *pDaughter=0;
    765766  G4VPVParameterisation *pParam=0;

  • trunk/source/geometry/navigation/History

    r1228 r1315  
    1 $Id: History,v 1.138 2009/12/11 05:48:53 japost Exp $
     1$Id: History,v 1.139 2010/03/08 13:57:48 gcosmo Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1717     * Reverse chronological order (last date on top), please *
    1818     ----------------------------------------------------------
    19 Dec 11th, 2009 - J.Apostolakis (geomnav-V09-02-11)
    20 ------------------------------
    21 - G4PropagatorInField
    22   * Implementation field by Region
    23       - with extra checking of logical volume
    24 
    25 Nov 30th, 2009 - J.Apostolakis (geomnav-V09-02-10)
    26 ------------------------------
    27 - G4VIntersectionLocator
    28   * Fix problem in ReEstimateEndPoint() behaviour for very small steps.
    29 
    30   * Corrected constructor to ensure that it initialises all data members.
    31   * Labelled the methods (in header) to make noticable those that
    32       must changed at every step.
    33       ( attributes: ChordFinder,  EpsilonStep,  Navigator )
     19
     20March 8th, 2010 - G.Cosmo (geomnav-V09-03-00)
     21-------------------------
     22- Avoid unnecessary creation of string for debug purposes in G4PropagatorInField.
     23  Courtesy of P.Elmer (CMS).
     24- Some printout formatting...
     25
     26December 11th, 2009 - J.Apostolakis (geomnav-V09-02-11)
     27-----------------------------------
     28- Implemented field by region in G4PropagatorInField, with extra checking of
     29  logical volume.
     30
     31November 30th, 2009 - J.Apostolakis (geomnav-V09-02-10)
     32-----------------------------------
     33- Fixes in G4VIntersectionLocator:
     34   * Fixed problem in ReEstimateEndPoint() behaviour for very small steps.
     35   * Corrected constructor to ensure that it initialises all data members.
     36   * Labelled the methods (in header) to make noticable those that must change
     37     at every step (attributes: 'ChordFinder', 'EpsilonStep', 'Navigator').
    3438
    3539Nov 23rd, 2009 - J.Apostolakis
    3640------------------------------
    37 - Fixed G4PropagatorInField.cc:
    38    * Revise condition for flagging ZeroStep to avoid fake triggering.
    39      Small proposed steps (from Physics) could trigger the previous
    40      simple condition for Tiny/Zero Problem steps.
    41    [ This was first commited on branch geomnav-V09-02-08_branch
    42        and tagged geomnav-V09-02-80  ]
    43 
    44 Nov 12th, 2009 - J.Apostolakis (geomnav-V09-02-09)
    45 ------------------------------
    46 - G4PropagatorInField: cleanup of minor methods:
    47    Added new method RefreshIntersectionLocator to update the state of
    48        Use this method to synchronise with IntersectionLocator.
    49    Deleted long obsolete methods: SetAccuraciesWithDeltaOneStep
    50            SetDeltaIntersection,  SetDeltaOneStep
    51    Revised implementation to avoid using the older methods
    52            GetDeltaIntersection   GetDeltaOneStep
    53 
    54 Nov 12th, 2009 - J.Apostolakis (geomnav-V09-02-08)
    55 ------------------------------
     41- Fixes in G4PropagatorInField.cc:
     42   * Revised condition for flagging ZeroStep to avoid fake triggering. Small
     43     proposed steps (from Physics) could trigger the previous simple condition
     44     for Tiny/Zero Problem steps.
     45
     46November 12th, 2009 - J.Apostolakis (geomnav-V09-02-09)
     47-----------------------------------
     48- Cleanup of minor methods in G4PropagatorInField:
     49   * Added new method RefreshIntersectionLocator() to update the state of locator.
     50     It is used to synchronise with G4IntersectionLocator.
     51   * Deleted long obsolete methods: SetAccuraciesWithDeltaOneStep(),
     52     SetDeltaIntersection() and SetDeltaOneStep().
     53   * Revised implementation to avoid using the older methods GetDeltaIntersection()
     54     and GetDeltaOneStep().
     55
     56November 12th, 2009 - J.Apostolakis (geomnav-V09-02-08)
     57-----------------------------------
    5658- Refined G4PropagatorInField.cc:
    5759   * Changed parameters for treating consecutive tiny/zero steps:
    5860     - the value of default decrease factor to 0.25 (from 0.1),
    59      - the value is not used when the step size falls below a certain threshold.s
    60        This threshold is changed. It is expressed as a multiple of fZeroStepThreshold 
    61        instead of kCarTolerance.  The reasons are:
    62          * the threshold were meant to refer to the  ZeroStepThreshold (it was kCarTolerance)
    63          * the values were now much larger than this ZeroStepThreshold.
    64      - The value given to the decrease factor after the first such thresholds is changed to 0.35
    65        instead of 0.25, 0.5. 
     61     - the value is not used when the step size falls below a certain thresholds
     62       This threshold is changed. It is expressed as a multiple of fZeroStepThreshold
     63       instead of kCarTolerance. This, since the thresholds were meant to refer to
     64       the  ZeroStepThreshold (it was kCarTolerance), and the values were now much
     65       larger than 'ZeroStepThreshold'.
     66     - The value given to the decrease factor after the first such thresholds is
     67       changed to 0.35 instead of 0.25, 0.5. 
    6668   * Improved printing of Diagnostic message.
    6769
    68 Nov  3rd, 2009 - J.Apostolakis (geomnav-V09-02-07)
    69 ------------------------------
     70November 3rd, 2009 - J.Apostolakis (geomnav-V09-02-07)
     71----------------------------------
    7072- Refined G4PropagatorInField.cc:
    71    * Added new member fZeroStepThreshold, to enable tuning of threshold for Tiny/Zero steps.
     73   * Added new member fZeroStepThreshold, to enable tuning of threshold for
     74     Tiny/Zero steps.
    7275   * Changed default value from 0.5 * kCarTolerance
    73                           to max (100,000 * kCarTolerance, 0.1 * micron )
    74    * This is first revision to address problem seen at boundary of volumes (Reports of Atlas) .
     76     to max (100,000 * kCarTolerance, 0.1 * micron )
     77   * This is first revision to address problem seen at boundary of volumes
     78     (Reports of ATLAS).
    7579
    7680June 2nd, 2009 - J.Apostolakis (geomnav-V09-02-06)

  • trunk/source/geometry/navigation/src/G4PropagatorInField.cc

    r1228 r1315  
    3636// 17.03.97 John Apostolakis,   renaming new set functions being added
    3737//
    38 // $Id: G4PropagatorInField.cc,v 1.50 2009/12/10 08:41:54 japost Exp $
     38// $Id: G4PropagatorInField.cc,v 1.51 2010/03/08 13:57:34 gcosmo Exp $
    3939// GEANT4 tag $ Name:  $
    4040// ---------------------------------------------------------------------------
    … …  
    129129                G4double&          currentSafety,                // IN/OUT
    130130                G4VPhysicalVolume* pPhysVol)
    131 {
    132   const G4String MethodName("G4PropagatorInField::ComputeStep");
    133    
     131{ 
    134132  // If CurrentProposedStepLength is too small for finding Chords
    135133  // then return with no action (for now - TODO: some action)
    … …  
    238236
    239237#ifdef G4DEBUG_FIELD
    240      G4cout << MethodName << ": "
    241             << " Decreasing step -  "
     238     G4cout << " G4PropagatorInField::ComputeStep(): " << G4endl
     239            << "  Decreasing step -  "
    242240            << " decreaseFactor= " << std::setw(8) << decreaseFactor
    243241            << " stepTrial = "     << std::setw(18) << stepTrial << " "
    … …  
    248246     if( stepTrial == 0.0 )  //  Change to make it < 0.1 * kCarTolerance ??
    249247     {
    250        G4cout << " G4PropagatorInField::ComputeStep "
    251               << " Particle abandoned due to lack of progress in field."
     248       G4cout << " G4PropagatorInField::ComputeStep(): " << G4endl
     249              << "  Particle abandoned due to lack of progress in field."
    252250              << G4endl
    253               << " Properties : " << pFieldTrack << " "
     251              << "  Properties : " << pFieldTrack << " "
    254252              << G4endl;
    255        G4cerr << " G4PropagatorInField::ComputeStep "
    256               << "  ERROR : attempting a zero step= " << stepTrial << G4endl
    257               << " while attempting to progress after " << fNoZeroStep
    258               << " trial steps.  Will abandon step." << G4endl;
     253       G4cerr << " G4PropagatorInField::ComputeStep() - ERROR " << G4endl
     254              << "  Attempting a zero step = " << stepTrial << G4endl
     255              << "  while attempting to progress after " << fNoZeroStep
     256              << " trial steps. Will abandon step." << G4endl;
    259257         fParticleIsLooping= true;
    260258         return 0;  // = stepTrial;
    … …  
    351349    if( (fVerboseLevel > 1) && (do_loop_count > fMax_loop_count-10 )) {
    352350      if( do_loop_count == fMax_loop_count-9 ){
    353         G4cout << "G4PropagatorInField::ComputeStep "
    354                << " Difficult track - taking many sub steps." << G4endl;
     351        G4cout << " G4PropagatorInField::ComputeStep(): " << G4endl
     352               << "  Difficult track - taking many sub steps." << G4endl;
    355353      }
    356354      printStatus( SubStepStartState, CurrentState, CurrentProposedStepLength,
    … …  
    369367    fParticleIsLooping = true;
    370368
    371     if ( fVerboseLevel > 0 ){
    372        G4cout << "G4PropagateInField: Killing looping particle "
     369    if ( fVerboseLevel > 0 )
     370    {
     371       G4cout << " G4PropagateInField::ComputeStep(): " << G4endl
     372              << "  Killing looping particle "
    373373              // << " of " << energy  << " energy "
    374374              << " after " << do_loop_count << " field substeps "
    375375              << " totaling " << StepTaken / mm << " mm " ;
    376376       if( pPhysVol )
    377           G4cout << " in the volume " << pPhysVol->GetName() ;
     377          G4cout << " in volume " << pPhysVol->GetName() ;
    378378       else
    379379         G4cout << " in unknown or null volume. " ;
    … …  
    401401      - End_PointAndTangent.GetCurveLength()) > 3.e-4 * TruePathLength )
    402402  {
    403     G4cerr << " ERROR - G4PropagatorInField::ComputeStep():" << G4endl
    404            << " Curve length mis-match, is advancement wrong ? " << G4endl;
    405     G4cerr << " The curve length of the endpoint should be: "
     403    G4cerr << " G4PropagatorInField::ComputeStep() - ERROR" << G4endl
     404           << "  Curve length mis-match, is advancement wrong ? " << G4endl;
     405    G4cerr << "  The curve length of the endpoint should be: "
    406406           << OriginalState.GetCurveLength() + TruePathLength << G4endl
    407            << " and it is instead: "
     407           << "  and it is instead: "
    408408           << End_PointAndTangent.GetCurveLength() << "." << G4endl
    409            << " A difference of: "
     409           << "  A difference of: "
    410410           << OriginalState.GetCurveLength() + TruePathLength
    411411              - End_PointAndTangent.GetCurveLength() << G4endl;
    412     G4cerr << " Original state= " << OriginalState   << G4endl
    413            << " Proposed state= " << End_PointAndTangent << G4endl;
    414     G4Exception("G4PropagatorInField::ComputeStep()", "IncorrectProposedEndPoint",
    415                 FatalException,
     412    G4cerr << "  Original state = " << OriginalState   << G4endl
     413           << "  Proposed state = " << End_PointAndTangent << G4endl;
     414    G4Exception("G4PropagatorInField::ComputeStep()",
     415                "IncorrectProposedEndPoint", FatalException,
    416416                "Curve length mis-match between original state and proposed endpoint of propagation.");
    417417  }
    … …  
    434434  }
    435435
    436   if( fNoZeroStep > fAbandonThreshold_NoZeroSteps ) {
     436  if( fNoZeroStep > fAbandonThreshold_NoZeroSteps )
     437  {
    437438     fParticleIsLooping = true;
    438      G4cout << " WARNING - G4PropagatorInField::ComputeStep():" << G4endl
    439             << " Zero progress for "  << fNoZeroStep << " attempted steps."
     439     G4cout << " G4PropagatorInField::ComputeStep() - WARNING" << G4endl
     440            << "  Zero progress for "  << fNoZeroStep << " attempted steps."
    440441            << G4endl;
    441      G4cout << "Proposed Step is "<<CurrentProposedStepLength <<" but Step Taken is "<< fFull_CurveLen_of_LastAttempt <<G4endl;
    442      G4cout << "For Particle with Charge ="<<fCharge
    443             << " Momentum="<< fInitialMomentumModulus<<" Mass="<< fMass<<G4endl;
     442     G4cout << "  Proposed Step is " << CurrentProposedStepLength
     443            << " but Step Taken is "<< fFull_CurveLen_of_LastAttempt << G4endl;
     444     G4cout << "  For Particle with Charge = " << fCharge
     445            << " Momentum = "<< fInitialMomentumModulus
     446            << " Mass = " << fMass << G4endl;
    444447       if( pPhysVol )
    445           G4cout << " in the volume " << pPhysVol->GetName() ;
     448         G4cout << " in volume " << pPhysVol->GetName() ;
    446449       else
    447450         G4cout << " in unknown or null volume. " ;
    448451       G4cout << G4endl;
    449452     if ( fVerboseLevel > 2 )
    450        G4cout << " Particle that is stuck will be killed." << G4endl;
     453       G4cout << " Particle is stuck; it will be killed." << G4endl;
    451454     fNoZeroStep = 0;
    452455  }
    … …  
    558561       
    559562     G4cout << "Step taken was " << step_len 
    560             << " out of PhysicalStep= " <<  requestStep << G4endl;
     563            << " out of PhysicalStep = " <<  requestStep << G4endl;
    561564     G4cout << "Final safety is: " << safety << G4endl;
    562565
    … …  
    579582{
    580583#if 0
    581   G4cout << " PiF: NoZeroStep= " << fNoZeroStep
    582          << " CurrentProposedStepLength= " << CurrentProposedStepLength
    583          << " Full_curvelen_last=" << fFull_CurveLen_of_LastAttempt
    584          << " last proposed step-length= " << fLast_ProposedStepLength
     584  G4cout << " PiF: NoZeroStep = " << fNoZeroStep
     585         << " CurrentProposedStepLength = " << CurrentProposedStepLength
     586         << " Full_curvelen_last =" << fFull_CurveLen_of_LastAttempt
     587         << " last proposed step-length = " << fLast_ProposedStepLength
    585588         << " decrease factor = " << decreaseFactor
    586589         << " step trial = " << stepTrial
    … …  
    698701  //
    699702  G4MagInt_Driver* integrDriver= GetChordFinder()->GetIntegrationDriver();
    700   integrDriver->SetVerboseLevel( fVerboseLevel - 2 ); 
     703  integrDriver->SetVerboseLevel( fVerboseLevel - 2 );
    701704  G4cout << "Set Driver verbosity to " << fVerboseLevel - 2 << G4endl;
    702705

  • trunk/source/geometry/solids/Boolean/History

    r831 r1315  
    11
    2 $Id: History,v 1.62 2007/10/24 08:19:11 gcosmo Exp $
     2$Id: History,v 1.64 2010/05/11 15:04:01 gcosmo Exp $
    33-------------------------------------------------------------------
    44
    … …  
    2020     * Reverse chronological order (last date on top), please *
    2121     ----------------------------------------------------------
     22
     23 May 11, 2010  J.Allison                  geom-bool-V09-03-00
     24  - Introduced recursive algorithm in CreatePolyhedron(): it uses
     25    HepPolyhedronProcessor from 'graphics_reps', which tries small shifts
     26    in an attempt to avoid numerical problems in the calculation of the
     27    polyhedron in BooleanProcessor. Recursion allows HepPolyhedronProcessor
     28    to try all permutations, also for Booleans of Booleans.
     29    Helps in reducing the number of cases of "Error in Boolean processor" for
     30    visualization, but still some stubborn cases are left.
    2231
    2332 Oct 24, 2007  V.Grichine                 geom-bool-V09-00-00

  • trunk/source/geometry/solids/Boolean/include/G4BooleanSolid.hh

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4BooleanSolid.hh,v 1.15 2006/10/19 15:34:49 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4BooleanSolid.hh,v 1.17 2010/05/11 15:03:45 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030//
    … …  
    5050#include "G4Transform3D.hh"
    5151
    52 class G4VSolid;
     52class HepPolyhedronProcessor;
    5353
    5454class G4BooleanSolid : public G4VSolid
    … …  
    111111    G4VSolid* fPtrSolidB;
    112112
     113    G4Polyhedron* StackPolyhedron(HepPolyhedronProcessor&,
     114                                  const G4VSolid*) const;
     115      // Stack polyhedra for processing. Return top polyhedron.
     116
    113117  private:
    114118

  • trunk/source/geometry/solids/Boolean/src/G4BooleanSolid.cc

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4BooleanSolid.cc,v 1.21 2006/10/19 15:34:49 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4BooleanSolid.cc,v 1.23 2010/05/11 15:03:45 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// Implementation for the abstract base class for solids created by boolean
    … …  
    4040#include "G4VSolid.hh"
    4141#include "G4Polyhedron.hh"
     42#include "HepPolyhedronProcessor.h"
    4243#include "Randomize.hh"
    4344
    … …  
    231232  return fpPolyhedron;
    232233}
     234
     235//////////////////////////////////////////////////////////////////////////
     236//
     237// Stacks polyhedra for processing. Returns top polyhedron.
     238
     239G4Polyhedron*
     240G4BooleanSolid::StackPolyhedron(HepPolyhedronProcessor& processor,
     241                                const G4VSolid* solid) const
     242{
     243  HepPolyhedronProcessor::Operation operation;
     244  const G4String& type = solid->GetEntityType();
     245  if (type == "G4UnionSolid")
     246    { operation = HepPolyhedronProcessor::UNION; }
     247  else if (type == "G4IntersectionSolid")
     248    { operation = HepPolyhedronProcessor::INTERSECTION; }
     249  else if (type == "G4SubtractionSolid")
     250    { operation = HepPolyhedronProcessor::SUBTRACTION; }
     251  else
     252  {
     253    std::ostringstream oss;
     254    oss << "Solid - " << solid->GetName()
     255        << " - Unrecognised composite solid"
     256        << "\n  Returning NULL !";
     257    G4Exception("StackPolyhedron()", "InvalidSetup",
     258                JustWarning, oss.str().c_str());
     259    return 0;
     260  }
     261
     262  G4Polyhedron* top = 0;
     263  const G4VSolid* solidA = solid->GetConstituentSolid(0);
     264  const G4VSolid* solidB = solid->GetConstituentSolid(1);
     265
     266  if (solidA->GetConstituentSolid(0))
     267  {
     268    top = StackPolyhedron(processor, solidA);
     269  }
     270  else
     271  {
     272    top = solidA->GetPolyhedron();
     273  }
     274  G4Polyhedron* operand = solidB->GetPolyhedron();
     275  processor.push_back (operation, *operand);
     276
     277  return top;
     278}

  • trunk/source/geometry/solids/Boolean/src/G4IntersectionSolid.cc

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4IntersectionSolid.cc,v 1.30 2006/11/08 09:37:41 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4IntersectionSolid.cc,v 1.32 2010/05/11 15:03:45 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// Implementation of methods for the class G4IntersectionSolid
    … …  
    5050#include "G4VGraphicsScene.hh"
    5151#include "G4Polyhedron.hh"
     52#include "HepPolyhedronProcessor.h"
    5253#include "G4NURBS.hh"
    5354// #include "G4NURBSbox.hh"
    … …  
    502503G4IntersectionSolid::CreatePolyhedron () const
    503504{
    504   G4Polyhedron* pA = fPtrSolidA->GetPolyhedron();
    505   G4Polyhedron* pB = fPtrSolidB->GetPolyhedron();
    506   if (pA && pB)
    507   {
    508     G4Polyhedron* resultant = new G4Polyhedron (pA->intersect(*pB));
    509     return resultant;
    510   }
    511   else
    512   {
    513     std::ostringstream oss;
    514     oss << "Solid - " << GetName()
    515         << " - one of the Boolean components has no" << G4endl
    516         << " corresponding polyhedron. Returning NULL !";
    517     G4Exception("G4IntersectionSolid::CreatePolyhedron()", "InvalidSetup",
    518                 JustWarning, oss.str().c_str());
    519     return 0;
    520   }
     505  HepPolyhedronProcessor processor;
     506  // Stack components and components of components recursively
     507  // See G4BooleanSolid::StackPolyhedron
     508  G4Polyhedron* top = StackPolyhedron(processor, this);
     509  G4Polyhedron* result = new G4Polyhedron(*top);
     510  if (processor.execute(*result)) { return result; }
     511  else { return 0; }
    521512}
    522513

  • trunk/source/geometry/solids/Boolean/src/G4SubtractionSolid.cc

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4SubtractionSolid.cc,v 1.31 2007/10/23 14:42:31 grichine Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4SubtractionSolid.cc,v 1.33 2010/05/11 15:03:45 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// Implementation of methods for the class G4IntersectionSolid
    … …  
    4848#include "G4VGraphicsScene.hh"
    4949#include "G4Polyhedron.hh"
     50#include "HepPolyhedronProcessor.h"
    5051#include "G4NURBS.hh"
    5152// #include "G4NURBSbox.hh"
    … …  
    463464G4SubtractionSolid::CreatePolyhedron () const
    464465{
    465   G4Polyhedron* pA = fPtrSolidA->GetPolyhedron();
    466   G4Polyhedron* pB = fPtrSolidB->GetPolyhedron();
    467   if (pA && pB)
    468   {
    469     G4Polyhedron* resultant = new G4Polyhedron (pA->subtract(*pB));
    470     return resultant;
    471   }
    472   else
    473   {
    474     std::ostringstream oss;
    475     oss << "Solid - " << GetName()
    476         << " - one of the Boolean components has no" << G4endl
    477         << " corresponding polyhedron. Returning NULL !";
    478     G4Exception("G4SubtractionSolid::CreatePolyhedron()", "InvalidSetup",
    479                 JustWarning, oss.str().c_str());
    480     return 0;
    481   }
     466  HepPolyhedronProcessor processor;
     467  // Stack components and components of components recursively
     468  // See G4BooleanSolid::StackPolyhedron
     469  G4Polyhedron* top = StackPolyhedron(processor, this);
     470  G4Polyhedron* result = new G4Polyhedron(*top);
     471  if (processor.execute(*result)) { return result; }
     472  else { return 0; }
    482473}
    483474

  • trunk/source/geometry/solids/Boolean/src/G4UnionSolid.cc

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4UnionSolid.cc,v 1.35 2007/10/23 14:42:31 grichine Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4UnionSolid.cc,v 1.37 2010/05/11 15:03:45 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// Implementation of methods for the class G4IntersectionSolid
    … …  
    4949#include "G4VGraphicsScene.hh"
    5050#include "G4Polyhedron.hh"
     51#include "HepPolyhedronProcessor.h"
    5152#include "G4NURBS.hh"
    5253// #include "G4NURBSbox.hh"
    … …  
    454455G4UnionSolid::CreatePolyhedron () const
    455456{
    456   G4Polyhedron* pA = fPtrSolidA->GetPolyhedron();
    457   G4Polyhedron* pB = fPtrSolidB->GetPolyhedron();
    458   if (pA && pB) {
    459     G4Polyhedron* resultant = new G4Polyhedron (pA->add(*pB));
    460     return resultant;
    461   } else {
    462     std::ostringstream oss;
    463     oss << GetName() <<
    464       ": one of the Boolean components has no corresponding polyhedron.";
    465     G4Exception("G4UnionSolid::CreatePolyhedron",
    466                 "", JustWarning, oss.str().c_str());
    467     return 0;
    468   }
     457  HepPolyhedronProcessor processor;
     458  // Stack components and components of components recursively
     459  // See G4BooleanSolid::StackPolyhedron
     460  G4Polyhedron* top = StackPolyhedron(processor, this);
     461  G4Polyhedron* result = new G4Polyhedron(*top);
     462  if (processor.execute(*result)) { return result; }
     463  else { return 0; }
    469464}
    470465

  • trunk/source/geometry/solids/CSG/History

    r1228 r1315  
    1 $Id: History,v 1.119 2009/11/30 10:20:53 gcosmo Exp $
     1$Id: History,v 1.121 2010/05/25 09:19:02 gcosmo Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1717     * Reverse chronological order (last date on top), please *
    1818     ----------------------------------------------------------
     19
     20May 25, 2010  G.Cosmo                    geom-csg-V09-03-01
     21- G4Box: more minor fixes in coding style, following code review discussion.
     22
     23May 18, 2010  G.Cosmo                    geom-csg-V09-03-00
     24- G4Box: first implementation of speed improvements and corrections from joint
     25  code review of G4Box class: introduced half-tolerance constants and better
     26  logic in some conditional statements.
     27- Some code cleanup.
     28
     29Dec 04, 2009  V.Grichine
     30- G4Orb: modified DistanceToIn(p,v) to be more consistent with Inside(p).
    1931
    2032Nov 30, 2009  V.Grichine                 geom-csg-V09-02-08

  • trunk/source/geometry/solids/CSG/include/G4Box.hh

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4Box.hh,v 1.17 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4Box.hh,v 1.18 2010/05/18 10:07:52 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// --------------------------------------------------------------------
    … …  
    101101    std::ostream& StreamInfo(std::ostream& os) const;
    102102
    103   // Functions for visualization
     103  // Utilities for visualization
    104104
    105105    void          DescribeYourselfTo (G4VGraphicsScene& scene) const;
    … …  
    125125
    126126    enum ESide {kUndefined,kPX,kMX,kPY,kMY,kPZ,kMZ};
    127       // Codes for faces (kPX=plus x face,kMY= minus y face etc)
     127      // Codes for faces (kPX= +x face, kMY= -y face, etc...)
    128128
    129129  private:

  • trunk/source/geometry/solids/CSG/src/G4Box.cc

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4Box.cc,v 1.44 2006/10/19 15:33:37 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4Box.cc,v 1.49 2010/05/25 10:14:41 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030//
    … …  
    3232// Implementation for G4Box class
    3333//
    34 //  24.06.98 - V. Grichine: insideEdge in DistanceToIn(p,v)
     34//  24.06.98 - V.Grichine: insideEdge in DistanceToIn(p,v)
    3535//  20.09.98 - V.Grichine: new algorithm of DistanceToIn(p,v)
    3636//  07.05.00 - V.Grichine: d= DistanceToIn(p,v), if d<e/2, d=0
    … …  
    6767  if ( (pX > 2*kCarTolerance)
    6868    && (pY > 2*kCarTolerance)
    69     && (pZ > 2*kCarTolerance) )
     69    && (pZ > 2*kCarTolerance) )  // limit to thickness of surfaces
    7070  {
    7171    fDx = pX ;
    … …  
    105105void G4Box::SetXHalfLength(G4double dx)
    106106{
    107   if(dx > 2*kCarTolerance)
     107  if(dx > 2*kCarTolerance)  // limit to thickness of surfaces
     108  {
    108109    fDx = dx;
     110  }
    109111  else
    110112  {
    … …  
    122124void G4Box::SetYHalfLength(G4double dy)
    123125{
    124   if(dy > 2*kCarTolerance)
     126  if(dy > 2*kCarTolerance)  // limit to thickness of surfaces
     127  {
    125128    fDy = dy;
     129  }
    126130  else
    127131  {
    … …  
    139143void G4Box::SetZHalfLength(G4double dz)
    140144{
    141   if(dz > 2*kCarTolerance)
     145  if(dz > 2*kCarTolerance)  // limit to thickness of surfaces
     146  {
    142147    fDz = dz;
     148  }
    143149  else
    144150  {
    … …  
    153159  fpPolyhedron = 0;
    154160}
    155    
    156 
    157161
    158162////////////////////////////////////////////////////////////////////////
    … …  
    193197    if (pVoxelLimit.IsXLimited())
    194198    {
    195       if ( xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance ||
    196            xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance    ) return false ;
     199      if ((xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance) ||
     200          (xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance)) { return false ; }
    197201      else
    198202      {
    199         if (xMin < pVoxelLimit.GetMinXExtent())
    200         {
    201           xMin = pVoxelLimit.GetMinXExtent() ;
    202         }
    203         if (xMax > pVoxelLimit.GetMaxXExtent())
    204         {
    205           xMax = pVoxelLimit.GetMaxXExtent() ;
    206         }
     203        xMin = std::max(xMin, pVoxelLimit.GetMinXExtent());
     204        xMax = std::min(xMax, pVoxelLimit.GetMaxXExtent());
    207205      }
    208206    }
    … …  
    213211    if (pVoxelLimit.IsYLimited())
    214212    {
    215       if ( yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance ||
    216            yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance   ) return false ;
     213      if ((yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance) ||
     214          (yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance)) { return false ; }
    217215      else
    218216      {
    219         if (yMin < pVoxelLimit.GetMinYExtent())
    220         {
    221           yMin = pVoxelLimit.GetMinYExtent() ;
    222         }
    223         if (yMax > pVoxelLimit.GetMaxYExtent())
    224         {
    225           yMax = pVoxelLimit.GetMaxYExtent() ;
    226         }
     217        yMin = std::max(yMin, pVoxelLimit.GetMinYExtent());
     218        yMax = std::min(yMax, pVoxelLimit.GetMaxYExtent());
    227219      }
    228220    }
    … …  
    233225    if (pVoxelLimit.IsZLimited())
    234226    {
    235       if ( zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance ||
    236            zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance   ) return false ;
     227      if ((zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance) ||
     228          (zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance)) { return false ; }
    237229      else
    238230      {
    239         if (zMin < pVoxelLimit.GetMinZExtent())
    240         {
    241           zMin = pVoxelLimit.GetMinZExtent() ;
    242         }
    243         if (zMax > pVoxelLimit.GetMaxZExtent())
    244         {
    245           zMax = pVoxelLimit.GetMaxZExtent() ;
    246         }
     231        zMin = std::max(zMin, pVoxelLimit.GetMinZExtent());
     232        zMax = std::min(zMax, pVoxelLimit.GetMaxZExtent());
    247233      }
    248234    }
    … …  
    286272    if (pVoxelLimit.IsLimited(pAxis) == false)
    287273    { 
    288       if ( pMin != kInfinity || pMax != -kInfinity )
     274      if ( (pMin != kInfinity) || (pMax != -kInfinity) )
    289275      {
    290276        existsAfterClip = true ;
    … …  
    292278        // Add 2*tolerance to avoid precision troubles
    293279
    294         pMin           -= kCarTolerance;
    295         pMax           += kCarTolerance;
     280        pMin -= kCarTolerance;
     281        pMax += kCarTolerance;
    296282      }
    297283    }     
    … …  
    303289       ( pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
    304290
    305       if ( pMin != kInfinity || pMax != -kInfinity )
     291      if ( (pMin != kInfinity) || (pMax != -kInfinity) )
    306292      {
    307293        existsAfterClip = true ;
    … …  
    356342EInside G4Box::Inside(const G4ThreeVector& p) const
    357343{
     344  static const G4double delta=0.5*kCarTolerance;
    358345  EInside in = kOutside ;
    359 
    360   if ( std::fabs(p.x()) <= fDx - kCarTolerance*0.5 )
    361   {
    362     if (std::fabs(p.y()) <= fDy - kCarTolerance*0.5 )
    363     {
    364       if      (std::fabs(p.z()) <= fDz - kCarTolerance*0.5 ) in = kInside ;
    365       else if (std::fabs(p.z()) <= fDz + kCarTolerance*0.5 ) in = kSurface ;
    366     }
    367     else if (std::fabs(p.y()) <= fDy + kCarTolerance*0.5 )
    368     {
    369       if (std::fabs(p.z()) <= fDz + kCarTolerance*0.5 ) in = kSurface ;
    370     }
    371   }
    372   else if (std::fabs(p.x()) <= fDx + kCarTolerance*0.5 )
    373   {
    374     if (std::fabs(p.y()) <= fDy + kCarTolerance*0.5 )
    375     {
    376       if (std::fabs(p.z()) <= fDz + kCarTolerance*0.5) in = kSurface ;
     346  G4ThreeVector q(std::fabs(p.x()), std::fabs(p.y()), std::fabs(p.z()));
     347
     348  if ( q.x() <= (fDx - delta) )
     349  {
     350    if (q.y() <= (fDy - delta) )
     351    {
     352      if      ( q.z() <= (fDz - delta) ) { in = kInside ;  }
     353      else if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
     354    }
     355    else if ( q.y() <= (fDy + delta) )
     356    {
     357      if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
     358    }
     359  }
     360  else if ( q.x() <= (fDx + delta) )
     361  {
     362    if ( q.y() <= (fDy + delta) )
     363    {
     364      if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
    377365    }
    378366  }
    … …  
    389377{
    390378  G4double distx, disty, distz ;
    391   G4ThreeVector norm ;
     379  G4ThreeVector norm(0.,0.,0.);
    392380
    393381  // Calculate distances as if in 1st octant
    … …  
    400388  // normals
    401389
    402   const G4double delta    = 0.5*kCarTolerance;
     390  static const G4double delta    = 0.5*kCarTolerance;
    403391  const G4ThreeVector nX  = G4ThreeVector( 1.0, 0,0  );
    404392  const G4ThreeVector nmX = G4ThreeVector(-1.0, 0,0  );
    … …  
    415403  {
    416404    noSurfaces ++;
    417     if ( p.x() >= 0.){        // on +X surface
    418       normX= nX ;    // G4ThreeVector( 1.0, 0., 0. );
    419     }else{
    420       normX= nmX;    // G4ThreeVector(-1.0, 0., 0. );
    421     }
     405    if ( p.x() >= 0. )  { normX= nX ; }       // on +X surface : (1,0,0)
     406    else                { normX= nmX; }       //                 (-1,0,0)
    422407    sumnorm= normX;
    423408  }
    … …  
    426411  {
    427412    noSurfaces ++;
    428     if ( p.y() >= 0.){        // on +Y surface
    429       normY= nY;
    430     }else{
    431       normY = nmY;
    432     }
     413    if ( p.y() >= 0. )  { normY= nY;  }       // on +Y surface
     414    else                { normY= nmY; }
    433415    sumnorm += normY;
    434416  }
    … …  
    437419  {
    438420    noSurfaces ++;
    439     if ( p.z() >= 0.){        // on +Z surface
    440       normZ= nZ;
    441     }else{
    442       normZ = nmZ;
    443     }
    444     sumnorm += normZ;
    445   }
    446 
    447   // sumnorm= normX + normY + normZ;
    448   const G4double invSqrt2 = 1.0 / std::sqrt( 2.0);
    449   const G4double invSqrt3 = 1.0 / std::sqrt( 3.0);
    450 
    451   norm= G4ThreeVector( 0., 0., 0.);
     421    if ( p.z() >= 0. )  { normZ= nZ;  }       // on +Z surface
     422    else                { normZ= nmZ; }
     423    sumnorm += normZ;
     424  }
     425
     426  static const G4double invSqrt2 = 1.0 / std::sqrt(2.0);
     427  static const G4double invSqrt3 = 1.0 / std::sqrt(3.0);
     428
    452429  if( noSurfaces > 0 )
    453430  {
    454     if( noSurfaces == 1 ){
     431    if( noSurfaces == 1 )
     432    {
    455433      norm= sumnorm;
    456     }else{
     434    }
     435    else
     436    {
    457437      // norm = sumnorm . unit();
    458       if( noSurfaces == 2 ) {
     438      if( noSurfaces == 2 )
     439      {
    459440        // 2 surfaces -> on edge
    460441        norm = invSqrt2 * sumnorm;
    461       } else {
     442      }
     443      else
     444      {
    462445        // 3 surfaces (on corner)
    463446        norm = invSqrt3 * sumnorm;
    464447      }
    465448    }
    466   }else{
     449  }
     450  else
     451  {
    467452#ifdef G4CSGDEBUG
    468453     G4Exception("G4Box::SurfaceNormal(p)", "Notification", JustWarning,
    … …  
    483468{
    484469  G4double distx, disty, distz ;
    485   G4ThreeVector norm ;
     470  G4ThreeVector norm(0.,0.,0.);
    486471
    487472  // Calculate distances as if in 1st octant
    … …  
    495480    if ( distx <= distz )     // Closest to X
    496481    {
    497       if ( p.x() < 0 ) norm = G4ThreeVector(-1.0,0,0) ;
    498       else             norm = G4ThreeVector( 1.0,0,0) ;
     482      if ( p.x() < 0 ) { norm = G4ThreeVector(-1.0,0,0) ; }
     483      else             { norm = G4ThreeVector( 1.0,0,0) ; }
    499484    }
    500485    else                      // Closest to Z
    501486    {
    502       if ( p.z() < 0 ) norm = G4ThreeVector(0,0,-1.0) ;
    503       else             norm = G4ThreeVector(0,0, 1.0) ;
     487      if ( p.z() < 0 ) { norm = G4ThreeVector(0,0,-1.0) ; }
     488      else             { norm = G4ThreeVector(0,0, 1.0) ; }
    504489    }
    505490  }
    … …  
    508493    if ( disty <= distz )      // Closest to Y
    509494    {
    510       if ( p.y() < 0 ) norm = G4ThreeVector(0,-1.0,0) ;
    511       else             norm = G4ThreeVector(0, 1.0,0) ;
     495      if ( p.y() < 0 ) { norm = G4ThreeVector(0,-1.0,0) ; }
     496      else             { norm = G4ThreeVector(0, 1.0,0) ; }
    512497    }
    513498    else                       // Closest to Z
    514499    {
    515       if ( p.z() < 0 ) norm = G4ThreeVector(0,0,-1.0) ;
    516       else             norm = G4ThreeVector(0,0, 1.0) ;
     500      if ( p.z() < 0 ) { norm = G4ThreeVector(0,0,-1.0) ; }
     501      else             { norm = G4ThreeVector(0,0, 1.0) ; }
    517502    }
    518503  }
    … …  
    541526// shape.
    542527
    543 G4double G4Box::DistanceToIn(const G4ThreeVector& p,const G4ThreeVector& v) const
     528G4double G4Box::DistanceToIn(const G4ThreeVector& p,
     529                             const G4ThreeVector& v) const
    544530{
    545531  G4double safx, safy, safz ;
    … …  
    549535  G4double sOut=kInfinity, sOuty=kInfinity, sOutz=kInfinity ;
    550536
     537  static const G4double delta = 0.5*kCarTolerance;
     538
    551539  safx = std::fabs(p.x()) - fDx ;     // minimum distance to x surface of shape
    552540  safy = std::fabs(p.y()) - fDy ;
    … …  
    558546  // travel is in a direction away from the shape.
    559547
    560   if (    ((p.x()*v.x() >= 0.0) && safx > -kCarTolerance*0.5)
    561        || ((p.y()*v.y() >= 0.0) && safy > -kCarTolerance*0.5)
    562        || ((p.z()*v.z() >= 0.0) && safz > -kCarTolerance*0.5)   )
     548  if (    ((p.x()*v.x() >= 0.0) && (safx > -delta))
     549       || ((p.y()*v.y() >= 0.0) && (safy > -delta))
     550       || ((p.z()*v.z() >= 0.0) && (safz > -delta))   )
    563551  {
    564552    return kInfinity ;  // travel away or parallel within tolerance
    … …  
    568556  // X Planes
    569557
    570   if ( v.x())
     558  if ( v.x() )  // != 0
    571559  {
    572560    stmp = 1.0/std::fabs(v.x()) ;
    … …  
    579567    else
    580568    {
    581       if (v.x() > 0)  sOut = (fDx - p.x())*stmp ;
    582       if (v.x() < 0)  sOut = (fDx + p.x())*stmp ;
     569      if (v.x() < 0)  { sOut = (fDx + p.x())*stmp ; }
     570      else            { sOut = (fDx - p.x())*stmp ; }
    583571    }
    584572  }
    … …  
    586574  // Y Planes
    587575
    588   if ( v.y())
     576  if ( v.y() )  // != 0
    589577  {
    590578    stmp = 1.0/std::fabs(v.y()) ;
    … …  
    595583      smaxy = (fDy+std::fabs(p.y()))*stmp ;
    596584
    597       if (sminy > smin) smin=sminy ;
    598       if (smaxy < smax) smax=smaxy ;
    599 
    600       if (smin >= smax-kCarTolerance*0.5)
     585      if (sminy > smin) { smin=sminy ; }
     586      if (smaxy < smax) { smax=smaxy ; }
     587
     588      if (smin >= (smax-delta))
    601589      {
    602590        return kInfinity ;  // touch XY corner
    … …  
    605593    else
    606594    {
    607       if (v.y() > 0)  sOuty = (fDy - p.y())*stmp ;
    608       if (v.y() < 0)  sOuty = (fDy + p.y())*stmp ;
    609       if( sOuty < sOut ) sOut = sOuty ;
     595      if (v.y() < 0)  { sOuty = (fDy + p.y())*stmp ; }
     596      else            { sOuty = (fDy - p.y())*stmp ; }
     597      if( sOuty < sOut ) { sOut = sOuty ; }
    610598    }     
    611599  }
    … …  
    613601  // Z planes
    614602
    615   if ( v.z() )
     603  if ( v.z() )  // != 0
    616604  {
    617605    stmp = 1.0/std::fabs(v.z()) ;
    618606
    619     if ( safz >= 0.0)
     607    if ( safz >= 0.0 )
    620608    {
    621609      sminz = safz*stmp ;
    622610      smaxz = (fDz+std::fabs(p.z()))*stmp ;
    623611
    624       if (sminz > smin) smin = sminz ;
    625       if (smaxz < smax) smax = smaxz ;
    626 
    627       if (smin >= smax-kCarTolerance*0.5)
     612      if (sminz > smin) { smin = sminz ; }
     613      if (smaxz < smax) { smax = smaxz ; }
     614
     615      if (smin >= (smax-delta))
    628616      {
    629617        return kInfinity ;    // touch ZX or ZY corners
    … …  
    632620    else
    633621    {
    634       if (v.z() > 0)  sOutz = (fDz - p.z())*stmp ;
    635       if (v.z() < 0)  sOutz = (fDz + p.z())*stmp ;
    636       if( sOutz < sOut ) sOut = sOutz ;
    637     }
    638   }
    639 
    640   if ( sOut <= smin + 0.5*kCarTolerance) // travel over edge
     622      if (v.z() < 0)  { sOutz = (fDz + p.z())*stmp ; }
     623      else            { sOutz = (fDz - p.z())*stmp ; }
     624      if( sOutz < sOut ) { sOut = sOutz ; }
     625    }
     626  }
     627
     628  if (sOut <= (smin + delta)) // travel over edge
    641629  {
    642630    return kInfinity ;
    643631  }
    644   if (smin < 0.5*kCarTolerance)  smin = 0.0 ;
     632  if (smin < delta)  { smin = 0.0 ; }
    645633
    646634  return smin ;
    … …  
    662650  safez = std::fabs(p.z()) - fDz ;
    663651
    664   if (safex > safe) safe = safex ;
    665   if (safey > safe) safe = safey ;
    666   if (safez > safe) safe = safez ;
     652  if (safex > safe) { safe = safex ; }
     653  if (safey > safe) { safe = safey ; }
     654  if (safez > safe) { safe = safez ; }
    667655
    668656  return safe ;
    … …  
    671659/////////////////////////////////////////////////////////////////////////
    672660//
    673 // Calcluate distance to surface of box from inside
     661// Calculate distance to surface of box from inside
    674662// by calculating distances to box's x/y/z planes.
    675663// Smallest distance is exact distance to exiting.
    … …  
    682670{
    683671  ESide side = kUndefined ;
    684   G4double pdist,stmp,snxt;
    685 
    686   if (calcNorm) *validNorm = true ; // All normals are valid
    687 
    688   if (v.x() > 0)   // X planes 
     672  G4double pdist,stmp,snxt=kInfinity;
     673
     674  static const G4double delta = 0.5*kCarTolerance;
     675
     676  if (calcNorm) { *validNorm = true ; }  // All normals are valid
     677
     678  if (v.x() > 0)   // X planes
    689679  {
    690680    pdist = fDx - p.x() ;
    691681
    692     if (pdist > kCarTolerance*0.5)
     682    if (pdist > delta)
    693683    {
    694684      snxt = pdist/v.x() ;
    … …  
    697687    else
    698688    {
    699       if (calcNorm) *n    = G4ThreeVector(1,0,0) ;
    700       return         snxt = 0 ;
    701     }
    702   }
    703   else if (v.x() < 0) 
     689      if (calcNorm) { *n   = G4ThreeVector(1,0,0) ; }
     690      return        snxt = 0 ;
     691    }
     692  }
     693  else if (v.x() < 0)
    704694  {
    705695    pdist = fDx + p.x() ;
    706696
    707     if (pdist > kCarTolerance*0.5)
     697    if (pdist > delta)
    708698    {
    709699      snxt = -pdist/v.x() ;
    … …  
    712702    else
    713703    {
    714       if (calcNorm) *n   = G4ThreeVector(-1,0,0) ;
     704      if (calcNorm) { *n   = G4ThreeVector(-1,0,0) ; }
    715705      return        snxt = 0 ;
    716706    }
    717707  }
    718   else snxt = kInfinity ;
    719 
    720   if ( v.y() > 0 )   // Y planes 
    721   {
    722     pdist=fDy-p.y();
    723 
    724     if (pdist>kCarTolerance*0.5)
    725     {
    726       stmp=pdist/v.y();
    727 
    728       if (stmp<snxt)
    729       {
    730         snxt=stmp;
    731         side=kPY;
    732       }
    733     }
    734     else
    735     {
    736       if (calcNorm) *n    = G4ThreeVector(0,1,0) ;
    737       return         snxt = 0 ;
    738     }
    739   }
    740   else if ( v.y() < 0 )
    741   {
    742     pdist = fDy + p.y() ;
    743 
    744     if (pdist > kCarTolerance*0.5)
    745     {
    746       stmp=-pdist/v.y();
    747 
    748       if (stmp<snxt)
    749       {
    750         snxt=stmp;
    751         side=kMY;
    752       }
    753     }
    754     else
    755     {
    756       if (calcNorm) *n    = G4ThreeVector(0,-1,0) ;
    757       return         snxt = 0 ;
    758     }
    759   }
    760   if (v.z()>0)        // Z planes
    761   {
    762     pdist=fDz-p.z();
    763 
    764     if (pdist > kCarTolerance*0.5)
    765     {
    766       stmp=pdist/v.z();
     708
     709  if (v.y() > 0)   // Y planes
     710  {
     711    pdist = fDy-p.y();
     712
     713    if (pdist > delta)
     714    {
     715      stmp = pdist/v.y();
    767716
    768717      if (stmp < snxt)
    769718      {
    770         snxt=stmp;
    771         side=kPZ;
     719        snxt = stmp;
     720        side = kPY;
    772721      }
    773722    }
    774723    else
    775724    {
    776       if (calcNorm) *n    = G4ThreeVector(0,0,1) ;
    777       return         snxt = 0 ;
    778     }
    779   }
    780   else if (v.z()<0)
    781   {
    782     pdist = fDz + p.z() ;
    783 
    784     if (pdist > kCarTolerance*0.5)
    785     {
    786       stmp=-pdist/v.z();
    787 
    788       if (stmp < snxt)
     725      if (calcNorm) { *n   = G4ThreeVector(0,1,0) ; }
     726      return        snxt = 0 ;
     727    }
     728  }
     729  else if (v.y() < 0)
     730  {
     731    pdist = fDy + p.y() ;
     732
     733    if (pdist > delta)
     734    {
     735      stmp = -pdist/v.y();
     736
     737      if ( stmp < snxt )
    789738      {
    790         snxt=stmp;
    791         side=kMZ;
     739        snxt = stmp;
     740        side = kMY;
    792741      }
    793742    }
    794743    else
    795744    {
    796       if (calcNorm) *n    = G4ThreeVector(0,0,-1) ;
    797       return         snxt = 0 ;
    798     }
    799   }
     745      if (calcNorm) { *n   = G4ThreeVector(0,-1,0) ; }
     746      return        snxt = 0 ;
     747    }
     748  }
     749
     750  if (v.z() > 0)        // Z planes
     751  {
     752    pdist = fDz-p.z();
     753
     754    if ( pdist > delta )
     755    {
     756      stmp = pdist/v.z();
     757
     758      if ( stmp < snxt )
     759      {
     760        snxt = stmp;
     761        side = kPZ;
     762      }
     763    }
     764    else
     765    {
     766      if (calcNorm) { *n   = G4ThreeVector(0,0,1) ; }
     767      return        snxt = 0 ;
     768    }
     769  }
     770  else if (v.z() < 0)
     771  {
     772    pdist = fDz + p.z();
     773
     774    if ( pdist > delta )
     775    {
     776      stmp = -pdist/v.z();
     777
     778      if ( stmp < snxt )
     779      {
     780        snxt = stmp;
     781        side = kMZ;
     782      }
     783    }
     784    else
     785    {
     786      if (calcNorm) { *n   = G4ThreeVector(0,0,-1) ; }
     787      return        snxt = 0 ;
     788    }
     789  }
     790
    800791  if (calcNorm)
    801792  {     
    … …  
    875866  // shortest Dist to any boundary now MIN(safx1,safx2,safy1..)
    876867
    877   if (safx2 < safx1) safe = safx2 ;
    878   else               safe = safx1 ;
    879   if (safy1 < safe)  safe = safy1 ;
    880   if (safy2 < safe)  safe = safy2 ;
    881   if (safz1 < safe)  safe = safz1 ;
    882   if (safz2 < safe)  safe = safz2 ;
    883 
    884   if (safe < 0) safe = 0 ;
     868  if (safx2 < safx1) { safe = safx2; }
     869  else               { safe = safx1; }
     870  if (safy1 < safe)  { safe = safy1; }
     871  if (safy2 < safe)  { safe = safy2; }
     872  if (safz1 < safe)  { safe = safz1; }
     873  if (safz2 < safe)  { safe = safz2; }
     874
     875  if (safe < 0) { safe = 0 ; }
    885876  return safe ; 
    886877}
    … …  
    979970    py = -fDy +2*fDy*G4UniformRand();
    980971
    981     if(G4UniformRand() > 0.5) pz =  fDz;
    982     else                      pz = -fDz;
     972    if(G4UniformRand() > 0.5) { pz =  fDz; }
     973    else                      { pz = -fDz; }
    983974  }
    984975  else if ( ( select - Sxy ) < Sxz )
    … …  
    987978    pz = -fDz +2*fDz*G4UniformRand();
    988979
    989     if(G4UniformRand() > 0.5) py =  fDy;
    990     else                      py = -fDy;
     980    if(G4UniformRand() > 0.5) { py =  fDy; }
     981    else                      { py = -fDy; }
    991982  }
    992983  else 
    … …  
    995986    pz = -fDz +2*fDz*G4UniformRand();
    996987
    997     if(G4UniformRand() > 0.5) px =  fDx;
    998     else                      px = -fDx;
     988    if(G4UniformRand() > 0.5) { px =  fDx; }
     989    else                      { px = -fDx; }
    999990  }
    1000991  return G4ThreeVector(px,py,pz);

  • trunk/source/geometry/solids/CSG/src/G4Orb.cc

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4Orb.cc,v 1.30 2009/11/30 10:20:38 gcosmo Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4Orb.cc,v 1.31 2009/12/04 15:39:56 grichine Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// class G4Orb
    … …  
    296296
    297297
    298   rad2 = p.x()*p.x()+p.y()*p.y()+p.z()*p.z() ;
     298  rad2 = p.x()*p.x()+p.y()*p.y()+p.z()*p.z();
    299299
    300300  G4double rad = std::sqrt(rad2);
    … …  
    304304  // sets `in'
    305305 
    306   tolRMax = fRmax - fRmaxTolerance*0.5 ;
     306  tolRMax = fRmax - fRmaxTolerance*0.5;
    307307   
    308   if ( rad <= tolRMax )  { in = kInside ; }
     308  if ( rad <= tolRMax )  { in = kInside; }
    309309  else
    310310  {
    311     tolRMax = fRmax + fRmaxTolerance*0.5 ;       
    312     if ( rad <= tolRMax )  { in = kSurface ; }
    313     else                   { in = kOutside ; }
     311    tolRMax = fRmax + fRmaxTolerance*0.5;       
     312    if ( rad <= tolRMax )  { in = kSurface; }
     313    else                   { in = kOutside; }
    314314  }
    315315  return in;
    … …  
    357357                              const G4ThreeVector& v  ) const
    358358{
    359   G4double snxt = kInfinity ;      // snxt = default return value
    360 
    361   G4double rad2, pDotV3d, tolORMax2, tolIRMax2 ;
    362   G4double c, d2, s = kInfinity ;
     359  G4double snxt = kInfinity;      // snxt = default return value
     360
     361  G4double rad, pDotV3d; // , tolORMax2, tolIRMax2;
     362  G4double c, d2, s = kInfinity;
    363363
    364364  const G4double dRmax = 100.*fRmax;
    … …  
    366366  // General Precalcs
    367367
    368   rad2    = p.x()*p.x() + p.y()*p.y() + p.z()*p.z() ;
    369   pDotV3d = p.x()*v.x() + p.y()*v.y() + p.z()*v.z() ;
     368  rad    = std::sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z());
     369  pDotV3d = p.x()*v.x() + p.y()*v.y() + p.z()*v.z();
    370370
    371371  // Radial Precalcs
    372372
    373   tolORMax2 = (fRmax+fRmaxTolerance*0.5)*(fRmax+fRmaxTolerance*0.5) ;
    374   tolIRMax2 = (fRmax-fRmaxTolerance*0.5)*(fRmax-fRmaxTolerance*0.5) ;
     373  // tolORMax2 = (fRmax+fRmaxTolerance*0.5)*(fRmax+fRmaxTolerance*0.5);
     374  // tolIRMax2 = (fRmax-fRmaxTolerance*0.5)*(fRmax-fRmaxTolerance*0.5);
    375375
    376376  // Outer spherical shell intersection
    … …  
    388388  // => s=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))
    389389
    390 
    391   G4double rad = std::sqrt(rad2);
    392390  c = (rad - fRmax)*(rad + fRmax);
    393391
    394   if ( c > fRmaxTolerance*fRmax )
    395   {
    396     // If outside tolerant boundary of outer G4Orb
    397     // [ should be std::sqrt(rad2) - fRmax > fRmaxTolerance*0.5 ]
    398 
    399     d2 = pDotV3d*pDotV3d - c ;
    400 
    401     if ( d2 >= 0 )
    402     {
    403       s = -pDotV3d - std::sqrt(d2) ;
    404       if ( s >= 0 )
    405       {
    406         if ( s>dRmax ) // Avoid rounding errors due to precision issues seen on
    407         {              // 64 bits systems. Split long distances and recompute
    408           G4double fTerm = s-std::fmod(s,dRmax);
    409           s = fTerm + DistanceToIn(p+fTerm*v,v);
    410         }
    411         return snxt = s;
    412       }
    413     }
    414     else    // No intersection with G4Orb
    415     {
    416       return snxt = kInfinity;
    417     }
    418   }
    419   else
    420   {
    421     if ( c > -fRmaxTolerance*fRmax )  // on surface 
    422     {
    423       d2 = pDotV3d*pDotV3d - c ;             
    424       if ( (d2 < fRmaxTolerance*fRmax) || (pDotV3d >= 0) )
     392  if( rad > fRmax-fRmaxTolerance*0.5 ) // not inside in terms of Inside(p)
     393  {
     394    if ( c > fRmaxTolerance*fRmax )
     395    {
     396      // If outside tolerant boundary of outer G4Orb in terms of c
     397      // [ should be std::sqrt(rad2) - fRmax > fRmaxTolerance*0.5 ]
     398
     399      d2 = pDotV3d*pDotV3d - c;
     400
     401      if ( d2 >= 0 )
     402      {
     403        s = -pDotV3d - std::sqrt(d2);
     404        if ( s >= 0 )
     405        {
     406          if ( s > dRmax ) // Avoid rounding errors due to precision issues seen on
     407          {                // 64 bits systems. Split long distances and recompute
     408            G4double fTerm = s - std::fmod(s,dRmax);
     409            s = fTerm + DistanceToIn(p+fTerm*v,v);
     410          }
     411          return snxt = s;
     412        }
     413      }
     414      else    // No intersection with G4Orb
    425415      {
    426416        return snxt = kInfinity;
    427417      }
    428       else
    429       {
    430         return snxt = 0.;
    431       }
    432     }
     418    }
     419    else // not outside in terms of c
     420    {
     421      if ( c > -fRmaxTolerance*fRmax )  // on surface 
     422      {
     423        d2 = pDotV3d*pDotV3d - c;             
     424        if ( (d2 < fRmaxTolerance*fRmax) || (pDotV3d >= 0) )
     425        {
     426          return snxt = kInfinity;
     427        }
     428        else
     429        {
     430          return snxt = 0.;
     431        }
     432      }
     433    }
     434  }
    433435#ifdef G4CSGDEBUG
    434     else // inside ???
    435     {
     436  else // inside ???
     437  {
    436438      G4Exception("G4Orb::DistanceToIn(p,v)", "Notification",
    437439                  JustWarning, "Point p is inside !?");
    438     }
     440  }
    439441#endif
    440   }
     442
    441443  return snxt;
    442444}
    … …  
    520522        if(calcNorm)
    521523        {
    522           *validNorm = true ;
    523           *n         = G4ThreeVector(p.x()/fRmax,p.y()/fRmax,p.z()/fRmax) ;
     524          *validNorm = true;
     525          *n         = G4ThreeVector(p.x()/fRmax,p.y()/fRmax,p.z()/fRmax);
    524526        }
    525527        return snxt = 0;
    … …  
    528530      {
    529531        snxt = -pDotV3d + std::sqrt(d2);    // second root since inside Rmax
    530         side = kRMax ;
     532        side = kRMax;
    531533      }
    532534    }
    … …  
    596598  if( Inside(p) == kOutside )
    597599  {
    598      G4cout.precision(16) ;
    599      G4cout << G4endl ;
     600     G4cout.precision(16);
     601     G4cout << G4endl;
    600602     DumpInfo();
    601      G4cout << "Position:"  << G4endl << G4endl ;
    602      G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
    603      G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
    604      G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
     603     G4cout << "Position:"  << G4endl << G4endl;
     604     G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl;
     605     G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl;
     606     G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl;
    605607     G4Exception("G4Orb::DistanceToOut(p)", "Notification", JustWarning,
    606608                 "Point p is outside !?" );

  • trunk/source/geometry/solids/specific/History

    r1228 r1315  
    1 $Id: History,v 1.160 2009/11/11 12:25:46 gcosmo Exp $
     1$Id: History,v 1.170 2010/06/11 09:42:59 gcosmo Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1717     * Reverse chronological order (last date on top), please *
    1818     ----------------------------------------------------------
     19
     2011-Jun-2010  T.Nikitina (geom-specific-V09-03-07)
     21- G4GenericTrap: fixed cases of zero dToIn/sToOut when vertices are collapsed
     22  to triangle, line or point. Added new methods handling those specific cases.
     23- Added unit test for G4GenericTrap.
     24
     2510-Jun-2010  I.Hrivnacova
     26- G4GenericTrap: fixed parameter names in CalculateExtent() for the test
     27  case construction through tessellated facets.
     28
     2909-Jun-2010  G.Cosmo
     30- G4GenericTrap: moved internal methods to private section and reordered
     31  in source file. Added missing implementation for IsTwisted() method.
     32
     3303-Jun-2010  T.Nikitina, G.Cosmo (geom-specific-V09-03-06)
     34- G4GenericTrap:
     35  o Fixed initialization of fSurfaceArea and fCubicVolume and
     36    calculation of surface area.
     37  o Fixed error in Inside(p) function, and corrected use std::fabs() instead
     38    of std::abs() for floating point values.
     39  o Added missing initialisation of fpPolyhedron pointer.
     40  o More corrected signatures for use of non-const references for vectors
     41    passed as arguments to functions.
     42
     4302-Jun-2010  G.Cosmo (geom-specific-V09-03-05)
     44- G4GenericTrap: use const reference for vector of vertices passed as argument
     45  in constructor and accessor.
     46
     4727-May-2010  T.Nikitina (geom-specific-V09-03-04)
     48- First implementation of G4GenericTrap shape, a new solid representing an
     49  arbitrary trapezoid with up to 8 vertices standing on two parallel planes
     50  perpendicular to the Z axis.
     51
     5228-Apr-2010  P.R.Truscott (geom-specific-V09-03-03)
     53- Fix in G4TriangularFacet and G4TessellatedSolid to correct treatment of
     54  optical photon transport related to internal reflection at surface.
     55  Addresses problem report #1103.
     56
     5715-Apr-2010  I.Hrivnacova (geom-specific-V09-03-02)
     58- G4ExtrudedSolid: eliminated requirement for clockwise ordering of polygon
     59  vertices. Added a check for vertices ordering; if vertices are defined
     60  anti-clockwise their ordering is reverted.
     61  Fix in polygon facet triangularization for consequent concave vertices.
     62
     6324-Feb-2010  G.Cosmo (geom-specific-V09-03-01)
     64- Adopt caching of Phi in G4PolyconeSide and G4PolyhedraSide to avoid
     65  unnecessary consecutive computations on the same point.
     66
     6710-Feb-2010  G.Cosmo (geom-specific-V09-03-00)
     68- Use kInfinity for initialising minimum and maximum allowed extent for
     69  G4SolidExtentList of faceted solids.
    1970
    207111-Nov-2009  G.Cosmo (geom-specific-V09-02-08)

  • trunk/source/geometry/solids/specific/include/G4ExtrudedSolid.hh

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4ExtrudedSolid.hh,v 1.7 2008/02/27 12:32:48 ivana Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4ExtrudedSolid.hh,v 1.8 2010/04/15 10:23:34 ivana Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030//
    … …  
    4545// const G4String& pName             - solid name
    4646// std::vector<G4TwoVector> polygon  - the vertices of the outlined polygon
    47 //                                     defined in clock-wise order     
     47//                                     defined in clockwise or anti-clockwise order     
    4848// std::vector<ZSection>             - the z-sections defined by
    4949//                                     z position, offset and scale

  • trunk/source/geometry/solids/specific/include/G4PolyconeSide.hh

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4PolyconeSide.hh,v 1.12 2008/05/15 11:41:58 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4PolyconeSide.hh,v 1.13 2010/02/24 11:18:25 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030//
    … …  
    114114  protected:
    115115
     116    G4double DistanceAway( const G4ThreeVector &p, G4bool opposite,
     117                                 G4double &distOutside2, G4double *rzNorm=0 );
     118     
     119    G4bool PointOnCone( const G4ThreeVector &hit, G4double normSign,
     120                        const G4ThreeVector &p,
     121                        const G4ThreeVector &v, G4ThreeVector &normal );
     122
     123    void CopyStuff( const G4PolyconeSide &source );
     124 
     125    static void FindLineIntersect( G4double x1, G4double y1,
     126                                   G4double tx1, G4double ty1,
     127                                   G4double x2, G4double y2,
     128                                 G4double tx2, G4double ty2,
     129                                 G4double &x, G4double &y );
     130
     131    G4double GetPhi( const G4ThreeVector& p );
     132
     133  protected:
     134
    116135    G4double r[2], z[2]; // r, z parameters, in specified order
    117136    G4double startPhi,   // Start phi (0 to 2pi), if phiIsOpen
    … …  
    136155    G4ThreeVector *corners; // The coordinates of the corners (if phiIsOpen)
    137156
    138     G4double DistanceAway( const G4ThreeVector &p, G4bool opposite,
    139                                  G4double &distOutside2, G4double *rzNorm=0 );
    140      
    141     G4bool PointOnCone( const G4ThreeVector &hit, G4double normSign,
    142                         const G4ThreeVector &p,
    143                         const G4ThreeVector &v, G4ThreeVector &normal );
    144 
    145     void CopyStuff( const G4PolyconeSide &source );
    146  
    147     static void FindLineIntersect( G4double x1, G4double y1,
    148                                    G4double tx1, G4double ty1,
    149                                    G4double x2, G4double y2,
    150                                  G4double tx2, G4double ty2,
    151                                  G4double &x, G4double &y );
    152157  private:
    153158
     159    std::pair<G4ThreeVector, G4double> fPhi;  // Cached value for phi
    154160    G4double kCarTolerance; // Geometrical surface thickness
    155161    G4double fSurfaceArea;  // Used for surface calculation

  • trunk/source/geometry/solids/specific/include/G4PolyhedraSide.hh

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4PolyhedraSide.hh,v 1.11 2008/05/15 11:41:59 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4PolyhedraSide.hh,v 1.12 2010/02/24 11:18:25 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030//
    … …  
    166166 
    167167    G4int PhiSegment( G4double phi );
    168  
     168
     169    G4double GetPhi( const G4ThreeVector& p );
     170
    169171    G4double DistanceToOneSide( const G4ThreeVector &p,
    170172                                const G4PolyhedraSideVec &vec,
    … …  
    197199  private:
    198200
     201    std::pair<G4ThreeVector, G4double> fPhi;  // Cached value for phi
    199202    G4double kCarTolerance;  // Geometrical surface thickness
    200203    G4double fSurfaceArea;   // Surface Area

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

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4ExtrudedSolid.cc,v 1.18 2008/10/30 11:47:45 ivana Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4ExtrudedSolid.cc,v 1.19 2010/04/15 10:23:34 ivana Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030//
    … …  
    4040#include <algorithm>
    4141#include <cmath>
     42#include <iomanip>
    4243
    4344#include "G4ExtrudedSolid.hh"
    … …  
    9495        "Z-sections with the same z position are not supported.");
    9596    }
    96   }     
    97 
     97  } 
     98 
     99  // Check if polygon vertices are defined clockwise
     100  // (the area is positive if polygon vertices are defined anti-clockwise)
     101  //
     102  G4double area = 0.;
     103  for ( G4int i=0; i<fNv; ++i ) {
     104    G4int j = i+1;
     105    if ( j == fNv ) j = 0;
     106    area += 0.5 * ( polygon[i].x()*polygon[j].y() - polygon[j].x()*polygon[i].y());
     107  }
     108 
    98109  // Copy polygon
    99110  //
    100   for ( G4int i=0; i<fNv; ++i ) { fPolygon.push_back(polygon[i]); }
     111  if  ( area < 0. ) {   
     112    // Polygon vertices are defined clockwise, we just copy the polygon       
     113    for ( G4int i=0; i<fNv; ++i ) { fPolygon.push_back(polygon[i]); }
     114  }
     115  else {
     116    // Polygon vertices are defined anti-clockwise, we revert them
     117    //G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", errorDescription,
     118    //            JustWarning,
     119    //            "Polygon vertices defined anti-clockwise, reverting polygon");     
     120    for ( G4int i=0; i<fNv; ++i ) { fPolygon.push_back(polygon[fNv-i-1]); }
     121  }
     122   
    101123 
    102124  // Copy z-sections
    … …  
    148170  }
    149171     
     172  // Check if polygon vertices are defined clockwise
     173  // (the area is positive if polygon vertices are defined anti-clockwise)
     174 
     175  G4double area = 0.;
     176  for ( G4int i=0; i<fNv; ++i ) {
     177    G4int j = i+1;
     178    if ( j == fNv ) j = 0;
     179    area += 0.5 * ( polygon[i].x()*polygon[j].y() - polygon[j].x()*polygon[i].y());
     180  }
     181 
    150182  // Copy polygon
    151183  //
    152   for ( G4int i=0; i<fNv; ++i ) { fPolygon.push_back(polygon[i]); }
     184  if  ( area < 0. ) {   
     185    // Polygon vertices are defined clockwise, we just copy the polygon       
     186    for ( G4int i=0; i<fNv; ++i ) { fPolygon.push_back(polygon[i]); }
     187  }
     188  else {
     189    // Polygon vertices are defined anti-clockwise, we revert them
     190    //G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", errorDescription,
     191    //            JustWarning,
     192    //            "Polygon vertices defined anti-clockwise, reverting polygon");     
     193    for ( G4int i=0; i<fNv; ++i ) { fPolygon.push_back(polygon[fNv-i-1]); }
     194  }
    153195 
    154196  // Copy z-sections
    … …  
    453495    //
    454496    G4double angle = GetAngle(c2->first, c3->first, c1->first);
    455 
    456     if ( angle > pi )
     497    //G4cout << "angle " << angle  << G4endl;
     498
     499    G4int counter = 0;
     500    while ( angle > pi )
    457501    {
    458502      // G4cout << "Skipping concave vertex " << c2->second << G4endl;
    … …  
    467511      // G4cout << "Looking at triangle : "
    468512      //        << c1->second << "  " << c2->second
    469       //        << "  " << c3->second << G4endl; 
    470 
     513      //        << "  " << c3->second << G4endl;
     514     
     515      angle = GetAngle(c2->first, c3->first, c1->first);
     516      //G4cout << "angle " << angle  << G4endl;
     517     
     518      counter++;
     519     
     520      if ( counter > fNv) {
     521        G4Exception("G4ExtrudedSolid::AddGeneralPolygonFacets", "InvalidSetup" ,
     522                    FatalException, "Triangularisation has failed.");
     523        break;
     524      } 
    471525    }
    472526
    … …  
    734788}
    735789
    736 
    737790//_____________________________________________________________________________
    738791
    … …  
    751804  for ( G4int i=0; i<fNv; ++i )
    752805  {
    753     os << "   vx = " << fPolygon[i].x()/mm << " mm"
     806    os << std::setw(5) << "#" << i
     807       << "   vx = " << fPolygon[i].x()/mm << " mm"
    754808       << "   vy = " << fPolygon[i].y()/mm << " mm" << G4endl;
    755809  }
    … …  
    764818  }     
    765819
     820/*
     821  // Triangles (for debogging)
     822  os << G4endl;
     823  os << " Triangles:" << G4endl;
     824  os << " Triangle #   vertex1   vertex2   vertex3" << G4endl;
     825
     826  G4int counter = 0;
     827  std::vector< std::vector<G4int> >::const_iterator it;
     828  for ( it = fTriangles.begin(); it != fTriangles.end(); it++ ) {
     829     std::vector<G4int> triangle = *it;
     830     os << std::setw(10) << counter++
     831        << std::setw(10) << triangle[0] << std::setw(10)  << triangle[1]  << std::setw(10)  << triangle[2]
     832        << G4endl;
     833  }         
     834*/
    766835  return os;
    767836} 

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

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4PolyconeSide.cc,v 1.22 2009/11/11 12:23:37 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4PolyconeSide.cc,v 1.24 2010/02/24 11:31:49 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030//
    … …  
    6767  kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
    6868  fSurfaceArea = 0.0;
     69  fPhi.first = G4ThreeVector(0,0,0);
     70  fPhi.second= 0.0;
    6971
    7072  //
    … …  
    492494  if (phiIsOpen)
    493495  {
    494     G4double phi = axis.phi();
     496    G4double phi = GetPhi(axis);
    495497    while( phi < startPhi ) phi += twopi;
    496498   
    … …  
    853855}
    854856
     857//
     858// GetPhi
     859//
     860// Calculate Phi for a given 3-vector (point), if not already cached for the
     861// same point, in the attempt to avoid consecutive computation of the same
     862// quantity
     863//
     864G4double G4PolyconeSide::GetPhi( const G4ThreeVector& p )
     865{
     866  G4double val=0.;
     867
     868  if (fPhi.first != p)
     869  {
     870    val = p.phi();
     871    fPhi.first = p;
     872    fPhi.second = val;
     873  }
     874  else
     875  {
     876    val = fPhi.second;
     877  }
     878  return val;
     879}
    855880
    856881//
    … …  
    922947    // Finally, check phi
    923948    //
    924     G4double phi = p.phi();
     949    G4double phi = GetPhi(p);
    925950    while( phi < startPhi ) phi += twopi;
    926951   
    … …  
    9781003    // PolyPhiFace. See PolyPhiFace::InsideEdgesExact
    9791004    //
    980     G4double phi = hit.phi();
     1005    G4double phi = GetPhi(hit);
    9811006    while( phi < startPhi-phiTolerant ) phi += twopi;
    9821007   

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

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4PolyhedraSide.cc,v 1.15 2008/05/15 11:41:59 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4PolyhedraSide.cc,v 1.17 2010/02/24 11:31:49 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030//
    … …  
    6464                                        G4bool isAllBehind )
    6565{
    66 
    6766  kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
    6867  fSurfaceArea=0.;
     68  fPhi.first = G4ThreeVector(0,0,0);
     69  fPhi.second= 0.0;
     70
    6971  //
    7072  // Record values
    … …  
    561563  // Try the closest phi segment first
    562564  //
    563   G4int iPhi = ClosestPhiSegment( p.phi() );
     565  G4int iPhi = ClosestPhiSegment( GetPhi(p) );
    564566 
    565567  G4ThreeVector pdotc = p - vecs[iPhi].center;
    … …  
    594596  // Which phi segment is closest to this point?
    595597  //
    596   G4int iPhi = ClosestPhiSegment( p.phi() );
     598  G4int iPhi = ClosestPhiSegment( GetPhi(p) );
    597599 
    598600  G4double norm;
    … …  
    624626  // Which phi segment is closest to this point?
    625627  //
    626   G4int iPhi = ClosestPhiSegment( p.phi() );
     628  G4int iPhi = ClosestPhiSegment( GetPhi(p) );
    627629
    628630  //
    … …  
    656658  // Which phi segment, if any, does the axis belong to
    657659  //
    658   iPhi = PhiSegment( axis.phi() );
     660  iPhi = PhiSegment( GetPhi(axis) );
    659661 
    660662  if (iPhi < 0)
    … …  
    971973
    972974//
     975// GetPhi
     976//
     977// Calculate Phi for a given 3-vector (point), if not already cached for the
     978// same point, in the attempt to avoid consecutive computation of the same
     979// quantity
     980//
     981G4double G4PolyhedraSide::GetPhi( const G4ThreeVector& p )
     982{
     983  G4double val=0.;
     984
     985  if (fPhi.first != p)
     986  {
     987    val = p.phi();
     988    fPhi.first = p;
     989    fPhi.second = val;
     990  }
     991  else
     992  {
     993    val = fPhi.second;
     994  }
     995  return val;
     996}
     997
     998
     999//
    9731000// DistanceToOneSide
    9741001//

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

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4SolidExtentList.cc,v 1.5 2007/05/11 13:54:29 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4SolidExtentList.cc,v 1.6 2010/02/10 16:38:37 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030//
    … …  
    5151  axis = kZAxis;
    5252  limited = false;
    53   minLimit = -DBL_MAX;
    54   maxLimit = +DBL_MAX;
     53  minLimit = -kInfinity;
     54  maxLimit = +kInfinity;
    5555}
    5656
    … …  
    7272  else
    7373  {
    74     minLimit = -DBL_MAX;
    75     maxLimit = +DBL_MAX;
     74    minLimit = -kInfinity;
     75    maxLimit = +kInfinity;
    7676  }
    7777}

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

    r1228 r1315  
    2525// ********************************************************************
    2626//
    27 // $Id: G4TessellatedSolid.cc,v 1.19 2009/04/27 08:06:27 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4TessellatedSolid.cc,v 1.20 2010/04/28 16:21:21 flei Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    … …  
    4242// CHANGE HISTORY
    4343// --------------
     44//
     45// 12 April 2010      P R Truscott, QinetiQ, bug fixes to treat optical
     46//                    photon transport, in particular internal reflection
     47//                    at surface.
    4448//
    4549// 14 November 2007   P R Truscott, QinetiQ & Stan Seibert, U Texas
    … …  
    670674    if ((*f)->Intersect(p,v,false,dist,distFromSurface,normal))
    671675    {
     676//
     677//
     678// Set minDist to the new distance to current facet if distFromSurface is in
     679// positive direction and point is not at surface.  If the point is within
     680// 0.5*kCarTolerance of the surface, then force distance to be zero and
     681// leave member function immediately (for efficiency), as proposed by & credit
     682// to Akira Okumura.
     683//
    672684      if (distFromSurface > 0.5*kCarTolerance && dist >= 0.0 && dist < minDist)
    673685      {
    674686        minDist  = dist;
     687      }
     688      else if (-0.5*kCarTolerance <= dist && dist <= 0.5*kCarTolerance)
     689      {
     690        return 0.0;
    675691      }
    676692    }

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

    r1228 r1315  
    2525// ********************************************************************
    2626//
    27 // $Id: G4TriangularFacet.cc,v 1.12 2008/11/13 08:25:07 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4TriangularFacet.cc,v 1.13 2010/04/28 16:21:21 flei Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    … …  
    5656//                  plane of the triangle.
    5757//
     58// 12 April 2010    P R Truscott, QinetiQ, bug fixes to treat optical
     59//                  photon transport, in particular internal reflection
     60//                  at surface.
    5861// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    5962
    … …  
    350353//
    351354//
    352 // Do a heck for rounding errors in the distance-squared.
     355// Do a check for rounding errors in the distance-squared.  It appears that
     356// the conventional methods for calculating sqrDist breaks down when very near
     357// to or at the surface (as required by transport).  We'll therefore also use
     358// the magnitude-squared of the vector displacement.  (Note that I've also
     359// tried to get around this problem by using the existing equations for
     360//
     361//    sqrDist = function(a,b,c,d,s,t)
     362//
     363// and use a more accurate addition process which minimises errors and
     364// breakdown of cummutitivity [where (A+B)+C != A+(B+C)] but this still
     365// doesn't work.
     366// Calculation from u = D + s*E[0] + t*E[1] is less efficient, but appears
     367// more robust.
    353368//
    354369  if (sqrDist < 0.0) { sqrDist = 0.0; }
    355 
    356   return D + s*E[0] + t*E[1];
     370  G4ThreeVector u = D + s*E[0] + t*E[1];
     371  G4double u2     = u.mag2();
     372//
     373//
     374// The following (part of the roundoff error check) is from Oliver Merle's
     375// updates.
     376//
     377  if ( sqrDist > u2 ) sqrDist = u2;
     378
     379  return u;
    357380}
    358381
    … …  
    542565// we intersect.
    543566//
    544       distance = -0.5*kCarTolerance;
     567//      distance = -0.5*kCarTolerance;
     568      distance = 0.0;
    545569      normal   = surfaceNormal;
    546570      return true;
    … …  
    728752  return area;
    729753}
     754////////////////////////////////////////////////////////////////////////
     755//

  • trunk/source/geometry/volumes/History

    r1228 r1315  
    1 $Id: History,v 1.163 2009/11/06 11:35:03 gcosmo Exp $
     1$Id: History,v 1.170 2010/04/23 10:27:49 gcosmo Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1717     * Reverse chronological order (last date on top), please *
    1818     ----------------------------------------------------------
     19
     20Apr 23rd, 2010      G.Cosmo - geomvol-V09-03-06
     21- Corrected initialisation of dummy copy-ctor in G4EnhancedVecAllocator,
     22  fixing compilation problem on WIN32-VC.
     23
     24Apr 22nd, 2010      G.Cosmo - geomvol-V09-03-05
     25- Make use of specialized allocator for handling internal vector in
     26  G4NavigatorHistory, globally controlling the memory pool.
     27  Measured ~2% average run-time speed-up.
     28
     29Apr 15th, 2010      G.Cosmo - geomvol-V09-03-04
     30- Combine use of G4Allocator for vectors in G4NavigationHistory with use
     31  of assign().
     32
     33Apr 13th, 2010      G.Cosmo - geomvol-V09-03-03
     34- Added Reset() method to G4ReflectionFactory for clearing maps of
     35  constituent and reflected volumes.
     36
     37Apr 9th, 2010       G.Cosmo - geomvol-V09-03-02
     38- Use more elegant solution in G4NavigationHistory copy-ctor. Adopt assign().
     39
     40Mar 31th, 2010      G.Cosmo - geomvol-V09-03-01
     41- Restore original vector allocation, but adopt simple loop copy within
     42  copy-constructor. Should provide slight performance improvement and keep
     43  locality.
     44
     45Dec 11th, 2009      G.Cosmo - geomvol-V09-03-00
     46- Use G4Allocator for vectors in G4NavigationHistory, to optimise memory
     47  management and reduce fragmentation.
    1948
    2049Nov 6th, 2009       G.Cosmo - geomvol-V09-02-04

  • trunk/source/geometry/volumes/include/G4NavigationHistory.hh

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4NavigationHistory.hh,v 1.14 2009/11/03 09:15:51 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4NavigationHistory.hh,v 1.18 2010/04/22 09:06:50 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// class G4NavigationHistory
    … …  
    4747#include "geomdefs.hh"
    4848
    49 //#include "G4Allocator.hh"
    5049#include "G4AffineTransform.hh"
    5150#include "G4VPhysicalVolume.hh"
    5251#include "G4NavigationLevel.hh"
     52#include "G4EnhancedVecAllocator.hh"
    5353
    5454#include <vector>
    … …  
    142142 private:
    143143
    144   std::vector<G4NavigationLevel> fNavHistory;
    145   //std::vector<G4NavigationLevel, G4Allocator<G4NavigationLevel> > fNavHistory;
     144  std::vector<G4NavigationLevel,
     145              G4EnhancedVecAllocator<G4NavigationLevel> > fNavHistory;
     146    // The geometrical tree; uses specialized allocator to optimize
     147    // memory handling and reduce possible fragmentation
    146148
    147149  G4int fStackDepth;

  • trunk/source/geometry/volumes/include/G4ReflectionFactory.hh

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4ReflectionFactory.hh,v 1.4 2008/11/13 09:33:20 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4ReflectionFactory.hh,v 1.5 2010/04/13 07:19:01 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030//
    … …  
    181181      // been reflected, after that placement or replication is performed.
    182182
     183    void Reset(); 
     184      // Resets maps of constituent and reflected volumes.
     185      // To be used exclusively when volumes are removed from the stores.
     186
    183187  protected: 
    184188

  • trunk/source/geometry/volumes/src/G4NavigationHistory.cc

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4NavigationHistory.cc,v 1.11 2009/08/03 16:27:37 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4NavigationHistory.cc,v 1.15 2010/04/22 09:06:50 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030//
    … …  
    3838#include "G4ios.hh"
    3939
     40// Initialise static data for the specialized memory pool
     41// for the internal STL vector of histories  ...
     42//
     43G4ChunkIndexType* G4AllocStats::allocStat = 0;
     44G4int             G4AllocStats::totSpace = 0;
     45G4int             G4AllocStats::numCat = 0;
     46
    4047G4NavigationHistory::G4NavigationHistory()
    4148  : fNavHistory(kHistoryMax), fStackDepth(0)
    … …  
    4552
    4653G4NavigationHistory::G4NavigationHistory(const G4NavigationHistory &h)
    47   : fNavHistory(h.fNavHistory), fStackDepth(h.fStackDepth)
     54  : fStackDepth(h.fStackDepth)
    4855{
     56  fNavHistory.assign(h.fNavHistory.begin(),h.fNavHistory.end());
    4957}
    5058

  • trunk/source/geometry/volumes/src/G4ReflectionFactory.cc

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4ReflectionFactory.cc,v 1.9 2008/11/13 09:33:20 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4ReflectionFactory.cc,v 1.10 2010/04/13 07:19:01 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030//
    … …  
    761761//_____________________________________________________________________________
    762762
     763void
     764G4ReflectionFactory::Reset()
     765{
     766  fConstituentLVMap.~map();
     767  fReflectedLVMap.~map();
     768}
     769
     770//_____________________________________________________________________________
     771
    763772void G4ReflectionFactory::PrintConstituentLVMap()
    764773{
Note: See TracChangeset for help on using the changeset viewer.