source: trunk/source/geometry/navigation/src/G4NormalNavigation.cc @ 1315

Last change on this file since 1315 was 1228, checked in by garnier, 14 years ago

update geant4.9.3 tag
File size: 17.1 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26//
27// $Id: G4NormalNavigation.cc,v 1.9 2007/05/11 13:43:59 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-03 $
29//
30//
31// class G4NormalNavigation Implementation
32//
33// Author: P.Kent, 1996
34//
35// --------------------------------------------------------------------
36
37#include "G4NormalNavigation.hh"
38#include "G4AffineTransform.hh"
39
40#include <iomanip>
41
42// ********************************************************************
43// Constructor
44// ********************************************************************
45//
46G4NormalNavigation::G4NormalNavigation()
47  : fCheck(false), fVerbose(0)
48{
49}
50
51// ********************************************************************
52// Destructor
53// ********************************************************************
54//
55G4NormalNavigation::~G4NormalNavigation()
56{;}
57
58// ********************************************************************
59// ComputeStep
60// ********************************************************************
61//
62G4double
63G4NormalNavigation::ComputeStep(const G4ThreeVector &localPoint,
64                                const G4ThreeVector &localDirection,
65                                const G4double currentProposedStepLength,
66                                      G4double &newSafety,
67                                      G4NavigationHistory &history,
68                                      G4bool &validExitNormal,
69                                      G4ThreeVector &exitNormal,
70                                      G4bool &exiting,
71                                      G4bool &entering,
72                                      G4VPhysicalVolume *(*pBlockedPhysical),
73                                      G4int &blockedReplicaNo)
74{
75  G4VPhysicalVolume *motherPhysical, *samplePhysical, *blockedExitedVol=0;
76  G4LogicalVolume *motherLogical;
77  G4VSolid *motherSolid;
78  G4ThreeVector sampleDirection;
79  G4double ourStep=currentProposedStepLength, motherSafety, ourSafety;
80  G4int localNoDaughters, sampleNo;
81
82  motherPhysical = history.GetTopVolume();
83  motherLogical  = motherPhysical->GetLogicalVolume();
84  motherSolid    = motherLogical->GetSolid();
85
86  // Compute mother safety
87  //
88  motherSafety = motherSolid->DistanceToOut(localPoint);
89  ourSafety = motherSafety; // Working isotropic safety
90 
91#ifdef G4VERBOSE
92  static G4int precVerf= 20;  // Precision
93  if ( fCheck )
94  {
95    if( fVerbose == 1 )
96    {
97      G4cout << "*** G4NormalNavigation::ComputeStep(): ***" << G4endl
98             << "    Invoked DistanceToOut(p) for mother solid: "
99             << motherSolid->GetName()
100             << ". Solid replied: " << motherSafety << G4endl
101             << "    For local point p: " << localPoint << G4endl
102             << "    To be considered as 'mother safety'." << G4endl;
103    }
104    if ( motherSafety < 0.0 )
105    {
106      G4cerr << "ERROR - G4NormalNavigation::ComputeStep()" << G4endl
107             << "        Current solid " << motherSolid->GetName()
108             << " gave negative safety: " << motherSafety << G4endl
109             << "        for the current (local) point " << localPoint
110             << G4endl;
111      motherSolid->DumpInfo();
112      G4Exception("G4NormalNavigation::ComputeStep()",
113                  "NegativeSafetyMotherVol", FatalException,
114                  "Negative Safety In Voxel Navigation !" );
115    }
116    if( motherSolid->Inside(localPoint)==kOutside )
117    {
118      G4cout << "WARNING - G4NormalNavigation::ComputeStep()" << G4endl
119             << "          Point " << localPoint
120             << " is outside current volume " << motherPhysical->GetName()
121             << G4endl;
122      G4double estDistToSolid= motherSolid->DistanceToIn(localPoint); 
123      G4cout << "          Estimated isotropic distance to solid (distToIn)= " 
124             << estDistToSolid << G4endl;
125      if( estDistToSolid > 100.0 * motherSolid->GetTolerance() )
126      {
127        motherSolid->DumpInfo();
128        G4Exception("G4NormalNavigation::ComputeStep()",
129                    "FarOutsideCurrentVolume", FatalException,
130                    "Point is far outside Current Volume !" ); 
131      }
132      else
133        G4Exception("G4NormalNavigation::ComputeStep()", "OutsideCurrentVolume", 
134                    JustWarning, "Point is a little outside Current Volume."); 
135    }
136
137    // Verification / verbosity
138    //
139    if ( fVerbose > 1 )
140    {
141      G4int oldprec = G4cout.precision(precVerf);
142      G4cout << " G4NormalNavigation::ComputeStep()"
143             << " - Information on mother / key daughters ..." << G4endl;
144      G4cout << " Type   " << std::setw(12) << "Solid-Name"   << " " 
145             << std::setw(3*(6+precVerf))   << " local point" << " "
146             << std::setw(4+precVerf)       << "solid-Safety" << " "
147             << std::setw(4+precVerf)       << "solid-Step"   << " "
148             << std::setw(17)               << "distance Method "
149             << std::setw(3*(6+precVerf))   << " local direction" << " "
150             << G4endl;
151      G4cout << " Mother " << std::setw(12) << motherSolid->GetName() << " "
152             << std::setw(4+precVerf)       << localPoint   << " "
153             << std::setw(4+precVerf)       << motherSafety << " "
154             << G4endl;
155      G4cout.precision(oldprec);
156    }
157
158  }
159#endif
160
161  //
162  // Compute daughter safeties & intersections
163  //
164
165  // Exiting normal optimisation
166  //
167  if ( exiting&&validExitNormal )
168  {
169    if ( localDirection.dot(exitNormal)>=kMinExitingNormalCosine )
170    {
171      // Block exited daughter volume
172      //
173      blockedExitedVol =* pBlockedPhysical;
174      ourSafety = 0;
175    }
176  }
177  exiting  = false;
178  entering = false;
179
180  localNoDaughters = motherLogical->GetNoDaughters();
181  for ( sampleNo=localNoDaughters-1; sampleNo>=0; sampleNo--)
182  {
183    samplePhysical = motherLogical->GetDaughter(sampleNo);
184    if ( samplePhysical!=blockedExitedVol )
185    {
186      G4AffineTransform sampleTf(samplePhysical->GetRotation(),
187                                 samplePhysical->GetTranslation());
188      sampleTf.Invert();
189      const G4ThreeVector samplePoint =
190              sampleTf.TransformPoint(localPoint);
191      const G4VSolid *sampleSolid =
192              samplePhysical->GetLogicalVolume()->GetSolid();
193      const G4double sampleSafety =
194              sampleSolid->DistanceToIn(samplePoint);
195      if ( sampleSafety<ourSafety )
196      {
197        ourSafety=sampleSafety;
198      }
199      if ( sampleSafety<=ourStep )
200      {
201        sampleDirection = sampleTf.TransformAxis(localDirection);
202        const G4double sampleStep =
203                sampleSolid->DistanceToIn(samplePoint,sampleDirection);
204#ifdef G4VERBOSE
205        if(( fCheck ) && ( fVerbose == 1 ))
206        {
207          G4cout << "*** G4NormalNavigation::ComputeStep(): ***" << G4endl
208                 << "    Invoked DistanceToIn(p,v) for daughter solid: "
209                 << sampleSolid->GetName()
210                 << ". Solid replied: " << sampleStep << G4endl
211                 << "    For local point p: " << samplePoint << G4endl
212                 << "    Direction v: " << sampleDirection
213                 << ", to be considered as 'daughter step'." << G4endl;
214        }
215#endif
216        if ( sampleStep<=ourStep )
217        {
218          ourStep  = sampleStep;
219          entering = true;
220          exiting  = false;
221          *pBlockedPhysical = samplePhysical;
222          blockedReplicaNo  = -1;
223
224#ifdef G4VERBOSE
225          // Check to see that the resulting point is indeed in/on volume.
226          // This check could eventually be made only for successful candidate.
227
228          if ( ( fCheck ) && ( sampleStep < kInfinity ) )
229          {
230            G4ThreeVector intersectionPoint;
231            intersectionPoint= samplePoint + sampleStep * sampleDirection;
232            EInside insideIntPt= sampleSolid->Inside(intersectionPoint); 
233            G4String solidResponse = "-kInside-";
234            if (insideIntPt == kOutside)
235              solidResponse = "-kOutside-";
236            else if (insideIntPt == kSurface)
237              solidResponse = "-kSurface-";
238            if( fVerbose == 1 )
239            {
240              G4cout << "*** G4NormalNavigation::ComputeStep(): ***" << G4endl
241                     << "    Invoked Inside() for solid: "
242                     << sampleSolid->GetName()
243                     << ". Solid replied: " << solidResponse << G4endl
244                     << "    For point p: " << intersectionPoint
245                     << ", considered as 'intersection' point." << G4endl;
246            }
247            if ( insideIntPt != kSurface )
248            {
249              G4int oldcoutPrec = G4cout.precision(16); 
250              G4cout << "WARNING - G4NormalNavigation::ComputeStep()" << G4endl
251                     << "          Inaccurate DistanceToIn for solid "
252                     << sampleSolid->GetName() << G4endl;
253              G4cout << "          Solid gave DistanceToIn = " << sampleStep
254                     << " yet returns " << solidResponse
255                     << " for this point !" << G4endl; 
256              G4cout << "          Point = " << intersectionPoint << G4endl;
257              if ( insideIntPt != kInside )
258                G4cout << "        DistanceToIn(p) = " 
259                       << sampleSolid->DistanceToIn(intersectionPoint)
260                       << G4endl;
261              if ( insideIntPt != kOutside ) 
262                G4cout << "        DistanceToOut(p) = " 
263                       << sampleSolid->DistanceToOut(intersectionPoint)
264                       << G4endl;
265              G4Exception("G4NormalNavigation::ComputeStep()", 
266                          "InaccurateDistanceToIn", JustWarning,
267                          "Navigator gets conflicting response from Solid."); 
268              G4cout.precision(oldcoutPrec);
269            }
270          }
271
272          // Verification / verbosity
273          //
274          if ( fVerbose > 1 )
275          {
276            G4int oldprec = G4cout.precision(precVerf);
277            G4cout << " Daught "
278                   << std::setw(12)         << sampleSolid->GetName() << " "
279                   << std::setw(4+precVerf) << samplePoint  << " "
280                   << std::setw(4+precVerf) << sampleSafety << " "
281                   << std::setw(4+precVerf) << sampleStep   << " "
282                   << std::setw(16)         << "distanceToIn" << " "
283                   << std::setw(4+precVerf) << localDirection << " "
284                   << G4endl;
285            G4cout.precision(oldprec);
286          }
287#endif
288        }
289      }
290    }
291  }
292  if ( currentProposedStepLength<ourSafety )
293  {
294    // Guaranteed physics limited
295    //
296    entering = false;
297    exiting  = false;
298    *pBlockedPhysical = 0;
299    ourStep = kInfinity;
300  }
301  else
302  {
303    // Compute mother intersection if required
304    //
305    if ( motherSafety<=ourStep )
306    {
307      G4double motherStep = motherSolid->DistanceToOut(localPoint,
308                                                       localDirection,
309                                                       true,
310                                                       &validExitNormal,
311                                                       &exitNormal);
312#ifdef G4VERBOSE
313      if ( fCheck )
314      {
315        if( fVerbose == 1 )
316        {
317          G4cout << "*** G4NormalNavigation::ComputeStep(): ***" << G4endl
318                 << "    Invoked DistanceToOut(p,v,...) for mother solid: "
319                 << motherSolid->GetName()
320                 << ". Solid replied: " << motherStep << G4endl
321                 << "    For local point p: " << localPoint << G4endl
322                 << "    Direction v: " << localDirection
323                 << ", to be considered as 'mother step'." << G4endl;
324        }
325        if( ( motherStep < 0.0 ) || ( motherStep >= kInfinity) )
326        {
327          G4cerr << "ERROR - G4NormalNavigation::ComputeStep()" << G4endl
328                 << "        Problem in Navigation"  << G4endl
329                 << "        Point (local coordinates): "
330                 << localPoint << G4endl
331                 << "        Local Direction: " << localDirection << G4endl
332                 << "        Solid: " << motherSolid->GetName() << G4endl; 
333          motherSolid->DumpInfo();
334          G4Exception("G4NormalNavigation::ComputeStep()",
335                      "PointDistOutInvalid", FatalException,
336                      "Current point is outside the current solid !");
337        }
338      }
339      if ( fVerbose > 1 )
340      {
341        G4int oldprec = G4cout.precision(precVerf);
342        G4cout << " Mother " << std::setw(12) << motherSolid->GetName() << " "
343               << std::setw(4+precVerf)       << localPoint   << " "
344               << std::setw(4+precVerf)       << motherSafety << " "
345               << std::setw(4+precVerf)       << motherStep   << " "
346               << std::setw(16)               << "distanceToOut" << " "
347               << std::setw(4+precVerf)       << localDirection << " "
348               << G4endl;
349        G4cout.precision(oldprec);     
350      }
351#endif
352
353      if ( motherStep<=ourStep )
354      {
355        ourStep  = motherStep;
356        exiting  = true;
357        entering = false;
358        if ( validExitNormal )
359        {
360          const G4RotationMatrix *rot = motherPhysical->GetRotation();
361          if (rot)
362          {
363            exitNormal *= rot->inverse();
364          }
365        }
366      }
367      else
368      {
369        validExitNormal = false;
370      }
371    }
372  }
373  newSafety = ourSafety;
374  return ourStep;
375}
376
377// ********************************************************************
378// ComputeSafety
379// ********************************************************************
380//
381G4double G4NormalNavigation::ComputeSafety(const G4ThreeVector &localPoint,
382                                           const G4NavigationHistory &history,
383                                           const G4double)
384{
385  G4VPhysicalVolume *motherPhysical, *samplePhysical;
386  G4LogicalVolume *motherLogical;
387  G4VSolid *motherSolid;
388  G4double motherSafety, ourSafety;
389  G4int localNoDaughters, sampleNo;
390
391  motherPhysical = history.GetTopVolume();
392  motherLogical  = motherPhysical->GetLogicalVolume();
393  motherSolid    = motherLogical->GetSolid();
394
395  // Compute mother safety
396  //
397  motherSafety = motherSolid->DistanceToOut(localPoint);
398  ourSafety = motherSafety; // Working isotropic safety
399
400#ifdef G4VERBOSE
401  if(( fCheck ) && ( fVerbose == 1 ))
402  {
403    G4cout << "*** G4NormalNavigation::ComputeSafety(): ***" << G4endl
404           << "    Invoked DistanceToOut(p) for mother solid: "
405           << motherSolid->GetName()
406           << ". Solid replied: " << motherSafety << G4endl
407           << "    For local point p: " << localPoint
408           << ", to be considered as 'mother safety'." << G4endl;
409  }
410#endif
411
412  // Compute daughter safeties
413  //
414  localNoDaughters = motherLogical->GetNoDaughters();
415  for ( sampleNo=localNoDaughters-1; sampleNo>=0; sampleNo-- )
416  {
417    samplePhysical = motherLogical->GetDaughter(sampleNo);
418    G4AffineTransform sampleTf(samplePhysical->GetRotation(),
419                               samplePhysical->GetTranslation());
420    sampleTf.Invert();
421    const G4ThreeVector samplePoint =
422            sampleTf.TransformPoint(localPoint);
423    const G4VSolid *sampleSolid =
424            samplePhysical->GetLogicalVolume()->GetSolid();
425    const G4double sampleSafety =
426            sampleSolid->DistanceToIn(samplePoint);
427    if ( sampleSafety<ourSafety )
428    {
429      ourSafety = sampleSafety;
430    }
431#ifdef G4VERBOSE
432    if(( fCheck ) && ( fVerbose == 1 ))
433    {
434      G4cout << "*** G4NormalNavigation::ComputeSafety(): ***" << G4endl
435             << "    Invoked DistanceToIn(p) for daughter solid: "
436             << sampleSolid->GetName()
437             << ". Solid replied: " << sampleSafety << G4endl
438             << "    For local point p: " << samplePoint
439             << ", to be considered as 'daughter safety'." << G4endl;
440    }
441#endif
442  }
443  return ourSafety;
444}
Note: See TracBrowser for help on using the repository browser.