Geant4  10.02.p02
G4NavigationLogger.cc
Go to the documentation of this file.
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: G4NavigationLogger.cc 97507 2016-06-03 12:48:42Z gcosmo $
28 //
29 //
30 // class G4NavigationLogger Implementation
31 //
32 // Author: G.Cosmo, 2010
33 //
34 // --------------------------------------------------------------------
35 
36 #include <iomanip>
37 #include <CLHEP/Units/SystemOfUnits.h>
38 
39 #include "G4NavigationLogger.hh"
40 #include "G4GeometryTolerance.hh"
41 
42 // const double millimeter = CLHEP::Units::millimeter;
43 using CLHEP::millimeter;
44 
45 G4NavigationLogger::G4NavigationLogger(const G4String& id)
46  : fId(id), fVerbose(0), fMinTriggerDistance( DBL_MAX ), fReportSoftWarnings( false )
47 {
48 }
49 
50 G4NavigationLogger::~G4NavigationLogger()
51 {
52 }
53 
54 // ********************************************************************
55 // PreComputeStepLog
56 // ********************************************************************
57 //
58 void
59 G4NavigationLogger::PreComputeStepLog(const G4VPhysicalVolume* motherPhysical,
60  G4double motherSafety,
61  const G4ThreeVector& localPoint) const
62 {
63  G4VSolid* motherSolid = motherPhysical->GetLogicalVolume()->GetSolid();
64  G4String fType = fId + "::ComputeStep()";
65  // const double millimeter = CLHEP::millimeter;
66 
67  if ( fVerbose == 1 || fVerbose > 4 )
68  {
69  G4cout << "*************** " << fType << " *****************" << G4endl
70  << " VolType "
71  << std::setw(15) << "Safety/mm" << " "
72  << std::setw(15) << "Distance/mm" << " "
73  << std::setw(52) << "Position (local coordinates)"
74  << " - Solid" << G4endl;
75  G4cout << " Mother "
76  << std::setw(15) << motherSafety / millimeter << " "
77  << std::setw(15) << "N/C" << " " << localPoint << " - "
78  << motherSolid->GetEntityType() << ": " << motherSolid->GetName()
79  << G4endl;
80  }
81  if ( motherSafety < 0.0 )
82  {
83  std::ostringstream message;
84  message << "Negative Safety In Voxel Navigation !" << G4endl
85  << " Current solid " << motherSolid->GetName()
86  << " gave negative safety: " << motherSafety / millimeter << G4endl
87  << " for the current (local) point " << localPoint;
88  message << " Solid info: " << *motherSolid << G4endl;
89  G4Exception(fType, "GeomNav0003", FatalException, message);
90  }
91  if( motherSolid->Inside(localPoint)==kOutside )
92  {
93  std::ostringstream message;
94  message << "Point is outside Current Volume - " << G4endl
95  << " Point " << localPoint / millimeter
96  << " is outside current volume '" << motherPhysical->GetName()
97  << "'" << G4endl;
98  G4double estDistToSolid= motherSolid->DistanceToIn(localPoint);
99  message << " Estimated isotropic distance to solid (distToIn)= "
100  << estDistToSolid << G4endl;
101  if( estDistToSolid > 100.0 * motherSolid->GetTolerance() )
102  {
103  message << " Solid info: " << *motherSolid << G4endl;
104  G4Exception(fType, "GeomNav0003", JustWarning, message,
105  "Point is far outside Current Volume !" );
106  }
107  else
108  {
109  G4Exception(fType, "GeomNav1001", JustWarning, message,
110  "Point is a little outside Current Volume.");
111  }
112  }
113 
114  // Verification / verbosity
115  //
116  if ( fVerbose > 1 )
117  {
118  static const G4int precVerf= 16; // Precision
119  G4int oldprec = G4cout.precision(precVerf);
120  G4cout << " - Information on mother / key daughters ..." << G4endl;
121  G4cout << " Type " << std::setw(12) << "Solid-Name" << " "
122  << std::setw(3*(6+precVerf)) << " local point" << " "
123  << std::setw(4+precVerf) << "solid-Safety" << " "
124  << std::setw(4+precVerf) << "solid-Step" << " "
125  << std::setw(17) << "distance Method "
126  << std::setw(3*(6+precVerf)) << " local direction" << " "
127  << G4endl;
128  G4cout << " Mother " << std::setw(12) << motherSolid->GetName() << " "
129  << std::setw(4+precVerf) << localPoint << " "
130  << std::setw(4+precVerf) << motherSafety << " "
131  << G4endl;
132  G4cout.precision(oldprec);
133  }
134 }
135 
136 // ********************************************************************
137 // AlongComputeStepLog
138 // ********************************************************************
139 //
140 void
141 G4NavigationLogger::AlongComputeStepLog(const G4VSolid* sampleSolid,
142  const G4ThreeVector& samplePoint,
143  const G4ThreeVector& sampleDirection,
144  const G4ThreeVector& localDirection,
145  G4double sampleSafety,
146  G4double sampleStep) const
147 {
148  // Check to see that the resulting point is indeed in/on volume.
149  // This check could eventually be made only for successful candidate.
150 
151  if ( sampleStep < kInfinity )
152  {
153  G4ThreeVector intersectionPoint;
154  intersectionPoint= samplePoint + sampleStep * sampleDirection;
155  EInside insideIntPt= sampleSolid->Inside(intersectionPoint);
156  G4String fType = fId + "::ComputeStep()";
157 
158  G4String solidResponse = "-kInside-";
159  if (insideIntPt == kOutside)
160  { solidResponse = "-kOutside-"; }
161  else if (insideIntPt == kSurface)
162  { solidResponse = "-kSurface-"; }
163 
164  if ( fVerbose == 1 || fVerbose > 4 )
165  {
166  G4cout << " Invoked Inside() for solid: "
167  << sampleSolid->GetName()
168  << ". Solid replied: " << solidResponse << G4endl
169  << " For point p: " << intersectionPoint
170  << ", considered as 'intersection' point." << G4endl;
171  }
172 
173  G4double safetyIn= -1, safetyOut= -1; // Set to invalid values
174  G4double newDistIn= -1, newDistOut= -1;
175  if( insideIntPt != kInside )
176  {
177  safetyIn= sampleSolid->DistanceToIn(intersectionPoint);
178  newDistIn= sampleSolid->DistanceToIn(intersectionPoint,
179  sampleDirection);
180  }
181  if( insideIntPt != kOutside )
182  {
183  safetyOut= sampleSolid->DistanceToOut(intersectionPoint);
184  newDistOut= sampleSolid->DistanceToOut(intersectionPoint,
185  sampleDirection);
186  }
187  if( insideIntPt != kSurface )
188  {
189  std::ostringstream message;
190  message.precision(16);
191  message << "Conflicting response from Solid." << G4endl
192  << " Inaccurate solid DistanceToIn"
193  << " for solid " << sampleSolid->GetName() << G4endl
194  << " Solid gave DistanceToIn = "
195  << sampleStep << " yet returns " << solidResponse
196  << " for this point !" << G4endl
197  << " Original Point = " << samplePoint << G4endl
198  << " Original Direction = " << sampleDirection << G4endl
199  << " Intersection Point = " << intersectionPoint << G4endl
200  << " Safety values: " << G4endl;
201  if ( insideIntPt != kInside )
202  {
203  message << " DistanceToIn(p) = " << safetyIn;
204  }
205  if ( insideIntPt != kOutside )
206  {
207  message << " DistanceToOut(p) = " << safetyOut;
208  }
209  message << G4endl;
210  message << " Solid Parameters: " << *sampleSolid;
211  G4Exception(fType, "GeomNav1001", JustWarning, message);
212  }
213  else
214  {
215  // If it is on the surface, *ensure* that either DistanceToIn
216  // or DistanceToOut returns a finite value ( >= Tolerance).
217  //
218  if( std::max( newDistIn, newDistOut ) <=
219  G4GeometryTolerance::GetInstance()->GetSurfaceTolerance() )
220  {
221  std::ostringstream message;
222  message << "Zero from both Solid DistanceIn and Out(p,v)." << G4endl
223  << " Identified point for which the solid "
224  << sampleSolid->GetName() << G4endl
225  << " has MAJOR problem: " << G4endl
226  << " --> Both DistanceToIn(p,v) and DistanceToOut(p,v) "
227  << "return Zero, an equivalent value or negative value."
228  << G4endl
229  << " Solid: " << sampleSolid << G4endl
230  << " Point p= " << intersectionPoint << G4endl
231  << " Direction v= " << sampleDirection << G4endl
232  << " DistanceToIn(p,v) = " << newDistIn << G4endl
233  << " DistanceToOut(p,v,..) = " << newDistOut << G4endl
234  << " Safety values: " << G4endl
235  << " DistanceToIn(p) = " << safetyIn << G4endl
236  << " DistanceToOut(p) = " << safetyOut;
237  G4Exception(fType, "GeomNav0003", FatalException, message);
238  }
239  }
240 
241  // Verification / verbosity
242  //
243  if ( fVerbose > 1 )
244  {
245  static const G4int precVerf= 20; // Precision
246  G4int oldprec = G4cout.precision(precVerf);
247  G4cout << "Daughter "
248  << std::setw(12) << sampleSolid->GetName() << " "
249  << std::setw(4+precVerf) << samplePoint << " "
250  << std::setw(4+precVerf) << sampleSafety << " "
251  << std::setw(4+precVerf) << sampleStep << " "
252  << std::setw(16) << "distanceToIn" << " "
253  << std::setw(4+precVerf) << localDirection << " "
254  << G4endl;
255  G4cout.precision(oldprec);
256  }
257  }
258 }
259 
260 // ********************************************************************
261 // CheckDaughterEntryPoint
262 // ********************************************************************
263 //
264 void
265 G4NavigationLogger::CheckDaughterEntryPoint(const G4VSolid* sampleSolid,
266  const G4ThreeVector& samplePoint,
267  const G4ThreeVector& sampleDirection,
268  const G4VSolid* motherSolid,
269  const G4ThreeVector& localPoint,
270  const G4ThreeVector& localDirection,
271  G4double motherStep,
272  G4double sampleStep) const
273 {
274  const G4double kCarTolerance = motherSolid->GetTolerance();
275 
276  // Double check the expected condition of being called
277  //
278  G4bool SuspiciousDaughterDist = ( sampleStep >= motherStep )
279  && ( sampleStep < kInfinity );
280 
281  if( sampleStep >= kInfinity )
282  {
283  G4ExceptionDescription msg;
284  msg.precision(12);
285  msg << " WARNING - Called with 'infinite' step. " << G4endl;
286  msg << " Checks have no meaning if daughter step is infinite." << G4endl;
287  msg << " kInfinity = " << kInfinity / millimeter << G4endl;
288  msg << " sampleStep = " << sampleStep / millimeter << G4endl;
289  msg << " sampleStep < kInfinity " << (sampleStep<kInfinity) << G4endl;
290  msg << " kInfinity - sampleStep " << (kInfinity-sampleStep) / millimeter << G4endl;
291  msg << " Returning immediately.";
292  G4Exception("G4NavigationLogger::CheckDaughterEntryPoint()",
293  "GeomNav0003", JustWarning, msg);
294  return;
295  }
296 
297  // The intersection point with the daughter is after the exit point
298  // from the mother volume !!
299  // This is legal / allowed to occur only if the mother is concave
300  // ****************************************************************
301  // If mother is convex the daughter volume must be encountered
302  // before the exit from the current volume!
303 
304  // Check #1) whether the track will re-enter the current mother
305  // in the extension past its current exit point
306  G4ThreeVector localExitMotherPos= localPoint+motherStep*localDirection;
307  G4double distExitToReEntry= motherSolid->DistanceToIn(localExitMotherPos,
308  localDirection);
309 
310  // Check #2) whether the 'entry' point in the daughter is inside the mother
311  //
312  G4ThreeVector localEntryInDaughter = localPoint+sampleStep*localDirection;
313  EInside insideMother= motherSolid->Inside( localEntryInDaughter );
314 
315  G4String solidResponse = "-kInside-";
316  if (insideMother == kOutside) { solidResponse = "-kOutside-"; }
317  else if (insideMother == kSurface) { solidResponse = "-kSurface-"; }
318 
319  G4double distToReEntry = distExitToReEntry + motherStep;
320  G4ThreeVector localReEntryPoint = localPoint+distToReEntry*localDirection;
321 
322  // Clear error -- Daughter entry point is bad
323  G4bool DaughterEntryIsOutside = SuspiciousDaughterDist
324  && ( (sampleStep < distToReEntry) || (insideMother == kOutside ) );
325  G4bool EntryIsMotherExit = std::fabs(sampleStep-motherStep) < kCarTolerance;
326 
327  // Check for more subtle error - is exit point of daughter correct ?
328  G4ThreeVector sampleEntryPoint= samplePoint+sampleStep*sampleDirection;
329  G4double sampleCrossingDist= sampleSolid->DistanceToOut( sampleEntryPoint,
330  sampleDirection );
331  G4double sampleExitDist = sampleStep+sampleCrossingDist;
332  G4ThreeVector sampleExitPoint= samplePoint+sampleExitDist*sampleDirection;
333 
334  G4bool TransitProblem = ( (sampleStep < motherStep)
335  && (sampleExitDist > motherStep + kCarTolerance) )
336  || ( EntryIsMotherExit && (sampleCrossingDist > kCarTolerance) );
337 
338  if( DaughterEntryIsOutside
339  || TransitProblem
340  || (SuspiciousDaughterDist && (fVerbose > 3) ) )
341  {
342  G4ExceptionDescription msg;
343  msg.precision(16);
344 
345  if( DaughterEntryIsOutside )
346  {
347  msg << "WARNING> Intersection distance to Daughter volume is further"
348  << " than the distance to boundary." << G4endl
349  << " It appears that part of the daughter volume is *outside*"
350  << " this mother. " << G4endl;
351  msg << " One of the following checks signaled a problem:" << G4endl
352  << " -sampleStep (dist to daugh) < mother-exit dist + distance "
353  << "to ReEntry point for mother " << G4endl
354  << " -position of daughter intersection is outside mother volume."
355  << G4endl;
356  }
357  else if( TransitProblem )
358  {
359  G4double protrusion = sampleExitDist - motherStep;
360 
361  msg << "WARNING> Daughter volume extends beyond mother boundary. "
362  << G4endl;
363  if ( ( sampleStep < motherStep )
364  && (sampleExitDist > motherStep + kCarTolerance ) )
365  {
366  // 1st Issue with Daughter
367  msg << " Crossing distance in the daughter causes is to extend"
368  << " beyond the mother exit. " << G4endl;
369  msg << " Length protruding = " << protrusion << G4endl;
370  }
371  if( EntryIsMotherExit )
372  {
373  // 1st Issue with Daughter
374  msg << " Intersection distance to Daughter is within "
375  << " tolerance of the distance" << G4endl;
376  msg << " to the mother boundary * and * " << G4endl;
377  msg << " the crossing distance in the daughter is > tolerance."
378  << G4endl;
379  }
380  }
381  else
382  {
383  msg << "NearMiss> Intersection to Daughter volume is in extension past the"
384  << " current exit point of the mother volume." << G4endl;
385  msg << " This is not an error - just an unusual occurence,"
386  << " possible in the case of concave volume. " << G4endl;
387  }
388  msg << "---- Information about intersection with daughter, mother: "
389  << G4endl;
390  msg << " sampleStep (daughter) = " << sampleStep << G4endl
391  << " motherStep = " << motherStep << G4endl
392  << " distToRentry(mother) = " << distToReEntry << G4endl
393  << " Inside(entry pnt daug): " << solidResponse << G4endl
394  << " dist across daughter = " << sampleCrossingDist << G4endl;
395  msg << " Mother Name (Solid) : " << motherSolid->GetName() << G4endl
396  << " In local (mother) coordinates: " << G4endl
397  << " Starting Point = " << localPoint << G4endl
398  << " Direction = " << localDirection << G4endl
399  << " Exit Point (mother)= " << localExitMotherPos << G4endl
400  << " Entry Point (daughter)= " << localPoint+sampleStep*localDirection
401  << G4endl;
402  if( distToReEntry < kInfinity )
403  {
404  msg << " ReEntry Point (mother)= " << localReEntryPoint << G4endl;
405  }
406  else
407  {
408  msg << " No ReEntry - track does not encounter mother volume again! "
409  << G4endl;
410  }
411  msg << " Daughter Name (Solid): " << sampleSolid->GetName() << G4endl
412  << " In daughter coordinates: " << G4endl
413  << " Starting Point = " << samplePoint << G4endl
414  << " Direction = " << sampleDirection << G4endl
415  << " Entry Point (daughter)= " << sampleEntryPoint
416  << G4endl;
417  msg << " Description of mother solid: " << G4endl
418  << *motherSolid << G4endl
419  << " Description of daughter solid: " << G4endl
420  << *sampleSolid << G4endl;
421  G4String fType = fId + "::ComputeStep()";
422 
423  if( DaughterEntryIsOutside || TransitProblem )
424  {
425  G4Exception(fType, "GeomNav0003", JustWarning, msg);
426  }
427  else
428  {
429  G4cout << fType
430  << " -- Checked distance of Entry to daughter vs exit of mother"
431  << G4endl;
432  G4cout << msg.str();
433  G4cout << G4endl;
434  }
435  }
436 }
437 
438 // ********************************************************************
439 // PostComputeStepLog
440 // ********************************************************************
441 //
442 void
443 G4NavigationLogger::PostComputeStepLog(const G4VSolid* motherSolid,
444  const G4ThreeVector& localPoint,
445  const G4ThreeVector& localDirection,
446  G4double motherStep,
447  G4double motherSafety) const
448 {
449  if ( fVerbose == 1 || fVerbose > 4 )
450  {
451  G4cout << " Mother "
452  << std::setw(15) << motherSafety << " "
453  << std::setw(15) << motherStep << " " << localPoint << " - "
454  << motherSolid->GetEntityType() << ": " << motherSolid->GetName()
455  << G4endl;
456  }
457  if( ( motherStep < 0.0 ) || ( motherStep >= kInfinity) )
458  {
459  G4String fType = fId + "::ComputeStep()";
460  G4int oldPrOut= G4cout.precision(16);
461  G4int oldPrErr= G4cerr.precision(16);
462  std::ostringstream message;
463  message << "Current point is outside the current solid !" << G4endl
464  << " Problem in Navigation" << G4endl
465  << " Point (local coordinates): "
466  << localPoint << G4endl
467  << " Local Direction: " << localDirection << G4endl
468  << " Solid: " << motherSolid->GetName();
469  motherSolid->DumpInfo();
470  G4Exception(fType, "GeomNav0003", FatalException, message);
471  G4cout.precision(oldPrOut);
472  G4cerr.precision(oldPrErr);
473  }
474  if ( fVerbose > 1 )
475  {
476  static const G4int precVerf= 20; // Precision
477  G4int oldprec = G4cout.precision(precVerf);
478  G4cout << " Mother " << std::setw(12) << motherSolid->GetName() << " "
479  << std::setw(4+precVerf) << localPoint << " "
480  << std::setw(4+precVerf) << motherSafety << " "
481  << std::setw(4+precVerf) << motherStep << " "
482  << std::setw(16) << "distanceToOut" << " "
483  << std::setw(4+precVerf) << localDirection << " "
484  << G4endl;
485  G4cout.precision(oldprec);
486  }
487 }
488 
489 // ********************************************************************
490 // ComputeSafetyLog
491 // ********************************************************************
492 //
493 void
494 G4NavigationLogger::ComputeSafetyLog(const G4VSolid* solid,
495  const G4ThreeVector& point,
496  G4double safety,
497  G4bool isMotherVolume,
498  G4int banner) const
499 {
500  if( banner < 0 )
501  {
502  banner = isMotherVolume;
503  }
504  if( fVerbose >= 1 )
505  {
506  G4String volumeType = isMotherVolume ? " Mother " : "Daughter";
507  if (banner)
508  {
509  G4cout << "************** " << fId << "::ComputeSafety() ****************"
510  << G4endl;
511  G4cout << " VolType "
512  << std::setw(15) << "Safety/mm" << " "
513  << std::setw(52) << "Position (local coordinates)"
514  << " - Solid" << G4endl;
515  }
516  G4cout << volumeType
517  << std::setw(15) << safety << " " << point << " - "
518  << solid->GetEntityType() << ": " << solid->GetName() << G4endl;
519  }
520 }
521 
522 // ********************************************************************
523 // PrintDaughterLog
524 // ********************************************************************
525 //
526 void
527 G4NavigationLogger::PrintDaughterLog (const G4VSolid* sampleSolid,
528  const G4ThreeVector& samplePoint,
529  G4double sampleSafety,
530  G4bool withStep,
531  const G4ThreeVector& sampleDirection, G4double sampleStep ) const
532 {
533  if ( fVerbose >= 1 )
534  {
535  G4int oldPrec= G4cout.precision(8);
536  G4cout << "Daughter "
537  << std::setw(15) << sampleSafety << " ";
538  if (withStep) // (sampleStep != -1.0 )
539  {
540  G4cout << std::setw(15) << sampleStep << " ";
541  }
542  else
543  {
544  G4cout << std::setw(15) << "Not-Available" << " ";
545  }
546  G4cout << samplePoint << " - "
547  << sampleSolid->GetEntityType() << ": " << sampleSolid->GetName();
548  if( withStep )
549  {
550  G4cout << " dir= " << sampleDirection;
551  }
552  G4cout << G4endl;
553  G4cout.precision(oldPrec);
554  }
555 }
556 
557 // ********************************************************************
558 // CheckAndReportBadNormal
559 // ********************************************************************
560 //
561 G4bool
562 G4NavigationLogger::
563 CheckAndReportBadNormal(const G4ThreeVector& unitNormal,
564  const G4ThreeVector& localPoint,
565  const G4ThreeVector& localDirection,
566  G4double step,
567  const G4VSolid* solid,
568  const char* msg ) const
569 {
570  G4double normMag2 = unitNormal.mag2();
571  G4bool badLength = ( std::fabs ( normMag2 - 1.0 ) > CLHEP::perMillion );
572 
573  if( badLength )
574  {
575  G4double normMag= std::sqrt(normMag2);
576  G4ExceptionDescription message;
577  message.precision(10);
578  message << "============================================================"
579  << G4endl;
580  message << " WARNING> Normal is not a unit vector. "
581  << " - but |normal| = " << normMag
582  << " - and |normal|^2 = " << normMag2 << G4endl
583  << " which differ from 1.0 by: " << G4endl
584  << " |normal|-1 = " << normMag - 1.0
585  << " and |normal|^2 - 1 = " << normMag2 - 1.0 << G4endl
586  << " n = " << unitNormal << G4endl;
587  message << " Info string: " << msg << G4endl;
588  message << "============================================================"
589  << G4endl;
590 
591  message.precision(16);
592 
593  message << " Information on call to DistanceToOut: " << G4endl;
594  message << " Position = " << localPoint << G4endl
595  << " Direction = " << localDirection << G4endl;
596  message << " Obtained> distance = " << step << G4endl;
597  message << " > Exit position = " << localPoint+step*localDirection
598  << G4endl;
599  message << " Parameters of solid: " << G4endl;
600  message << *solid;
601  message << "============================================================";
602 
603  G4String fMethod = fId + "::ComputeStep()";
604  G4Exception( fMethod, "GeomNav0003", JustWarning, message);
605  }
606  return badLength;
607 }
608 
609 // ********************************************************************
610 // CheckAndReportBadNormal - due to Rotation Matrix
611 // ********************************************************************
612 //
613 G4bool
614 G4NavigationLogger::
615 CheckAndReportBadNormal(const G4ThreeVector& rotatedNormal,
616  const G4ThreeVector& originalNormal,
617  const G4RotationMatrix& rotationM,
618  const char* msg ) const
619 {
620  G4double normMag2 = rotatedNormal.mag2();
621  G4bool badLength = ( std::fabs ( normMag2 - 1.0 ) > CLHEP::perMillion );
622 
623  if( badLength )
624  {
625  G4double normMag= std::sqrt(normMag2);
626  G4ExceptionDescription message;
627  message.precision(10);
628  message << "============================================================"
629  << G4endl;
630  message << " WARNING> Rotated n(ormal) is not a unit vector. " << G4endl
631  << " |normal| = " << normMag
632  << " and |normal|^2 = " << normMag2 << G4endl
633  << " Diff from 1.0: " << G4endl
634  << " |normal|-1 = " << normMag - 1.0
635  << " and |normal|^2 - 1 = " << normMag2 - 1.0 << G4endl;
636  message << " Rotated n = (" << rotatedNormal.x() << "," << rotatedNormal.y() << ","
637  << rotatedNormal.z() << ")" << G4endl;
638  message << " Original n = (" << originalNormal.x() << "," << originalNormal.y() << ","
639  << originalNormal.z() << ")" << G4endl;
640  message << " Info string: " << msg << G4endl;
641  message << "============================================================"
642  << G4endl;
643 
644  message.precision(16);
645 
646  message << " Information on RotationMatrix : " << G4endl;
647  message << " Original: " << G4endl;
648  message << rotationM << G4endl;
649  message << " Inverse (used in transformation): " << G4endl;
650  message << rotationM.inverse() << G4endl;
651  message << "============================================================";
652 
653  G4String fMethod = fId + "::ComputeStep()";
654  G4Exception( fMethod, "GeomNav0003", JustWarning, message);
655  }
656  return badLength;
657 }
658 
659 // ********************************************************************
660 // ReportOutsideMother
661 // ********************************************************************
662 //
663 // Report Exception if point is outside mother.
664 // Fatal exception will be used if either 'checkMode error is > triggerDist
665 //
666 void
667 G4NavigationLogger::ReportOutsideMother(const G4ThreeVector& localPoint,
668  const G4ThreeVector& localDirection,
669  const G4VPhysicalVolume* physical,
670  G4double triggerDist) const
671 {
672  const G4LogicalVolume* logicalVol = physical
673  ? physical->GetLogicalVolume() : 0;
674  const G4VSolid* solid = logicalVol
675  ? logicalVol->GetSolid() : 0;
676 
677  G4String fMethod = fId + "::ComputeStep()";
678 
679  if( solid == 0 )
680  {
681  G4Exception(fMethod, "GeomNav0003", FatalException,
682  "Erroneous call to ReportOutsideMother: no Solid is available");
683  return;
684  }
685  const G4double kCarTolerance = solid->GetTolerance();
686 
687  // Double check reply - it should be kInfinity
688  const G4double distToOut = solid->DistanceToOut(localPoint, localDirection);
689  // const G4double distToOutNeg = solid->DistanceToOut(localPoint, -localDirection);
690 
691  const EInside inSolid = solid->Inside(localPoint);
692  const G4double safetyToIn = solid->DistanceToIn(localPoint);
693  const G4double safetyToOut = solid->DistanceToOut(localPoint);
694 
695  const G4double distToInPos = solid->DistanceToIn(localPoint, localDirection);
696  // const G4double distToInNeg = solid->DistanceToIn(localPoint, -localDirection);
697 
698  // 1. Check consistency between Safety obtained and report from distance
699  // We must ensure that (mother)Safety <= 0.0
700  // in the case that the point is outside!
701  // [ If this is not the case, this problem can easily go undetected,
702  // except in Check mode ! ]
703  if( safetyToOut > kCarTolerance
704  && ( distToOut < 0.0 || distToOut >= kInfinity ) )
705  {
706  G4ExceptionDescription msg1;
707  // fNavClerk->ReportBadSafety(localPoint, localDirection,
708  // motherPhysical, motherSafety, motherStep);
709  msg1 << " Dangerous inconsistency in response of solid." << G4endl
710  << " Solid type: " << solid->GetEntityType()
711  << " Name= " << solid->GetName() << G4endl;
712  msg1 << " Mother volume gives safety > 0 despite being called for *Outside* point "
713  << G4endl
714  << " Location = " << localPoint << G4endl
715  << " Direction= " << localDirection << G4endl
716  << " - Safety (Isotropic d) = " << safetyToOut << G4endl
717  << " - Intersection Distance= " << distToOut << G4endl
718  << G4endl;
719  G4Exception( fMethod, "GeomNav0123", JustWarning, msg1);
720  }
721 
722  // 2. Inconsistency - Too many distances are zero (or will be rounded to zero)
723 
724  G4ExceptionDescription msg2;
725  if( std::fabs(distToOut) < kCarTolerance && std::fabs(distToInPos) < kCarTolerance )
726  {
727  // If both distanceToIn and distanceToOut (p,v) are zero for one direction,
728  // the particle could get stuck!
729  }
730 
731  G4ExceptionDescription msg;
732  msg.precision(10);
733  // G4bool reportIssue= true;
734 
735  if( std::fabs(distToOut) < kCarTolerance )
736  {
737  // 3. Soft error - safety is not rounded to zero
738  // Report nothing - except in 'loud' mode
739  if( fReportSoftWarnings )
740  {
741  // reportIssue= true;
742  msg << " Warning> DistanceToOut(p,v): Distance from surface is not rounded to zero" << G4endl;
743  } else {
744  // reportIssue= false;
745  return;
746  }
747  }
748  else
749  {
750  // 4. General message - complain that the point is outside!
751  // and provide all information about the current location,
752  // direction and the answers of the solid
753  msg << "============================================================" << G4endl;
754  msg << " WARNING> Current Point appears to be Outside mother volume !! " << G4endl;
755  msg << " Response of DistanceToOut was negative or kInfinity when called in "
756  << fMethod << G4endl;
757  }
758 
759  // Generate and 'print'/stream all the information needed
760  this->ReportVolumeAndIntersection( msg, localPoint, localDirection, physical );
761 
762  // Default for distance of 'major' error
763  if( triggerDist <= 0.0 ) {
764  // triggerDist = 1.e+6 * kCarTolerance; // Well beyond tolerance
765  triggerDist = std::max ( 1.0e+6 * kCarTolerance, // Well beyond tolerance
766  fMinTriggerDistance );
767  }
768 
769  G4bool majorError = inSolid==kOutside
770  ? ( safetyToIn > triggerDist )
771  : ( safetyToOut > triggerDist );
772 
773  G4ExceptionSeverity exceptionType = JustWarning;
774  if ( majorError )
775  {
776  exceptionType = FatalException;
777  }
778 
779  G4Exception( fMethod, "GeomNav0003", exceptionType, msg);
780 }
781 
782 namespace G4NavigationLogger_Namespace {
783  const G4String EInsideNames[3] = { "kOutside", "kSurface", "kInside" };
784 }
785 
786 void
787 G4NavigationLogger::ReportVolumeAndIntersection( std::ostream& os,
788  const G4ThreeVector& localPoint,
789  const G4ThreeVector& localDirection,
790  const G4VPhysicalVolume* physical ) const
791 {
792  G4String fMethod = fId + "::ComputeStep()";
793  const G4LogicalVolume* logicalVol = physical
794  ? physical->GetLogicalVolume() : 0;
795  const G4VSolid* solid = logicalVol
796  ? logicalVol->GetSolid() : 0;
797  if( solid == 0 )
798  {
799  os << " ERROR> Solid is not available. Logical Volume = " << logicalVol << std::endl;
800  return;
801  }
802  const G4double kCarTolerance = solid->GetTolerance();
803 
804  // Double check reply - it should be kInfinity
805  const G4double distToOut = solid->DistanceToOut(localPoint, localDirection);
806  const G4double distToOutNeg = solid->DistanceToOut(localPoint, -localDirection);
807 
808  const EInside inSolid = solid->Inside(localPoint);
809  const G4double safetyToIn = solid->DistanceToIn(localPoint);
810  const G4double safetyToOut = solid->DistanceToOut(localPoint);
811 
812  const G4double distToInPos = solid->DistanceToIn(localPoint, localDirection);
813  const G4double distToInNeg = solid->DistanceToIn(localPoint, -localDirection);
814 
815  const G4ThreeVector exitNormal = solid->SurfaceNormal(localPoint);
816 
817  // Double check whether points nearby are in/surface/out
818  const G4double epsilonDist = 1000.0 * kCarTolerance;
819  const G4ThreeVector PointPlusDir = localPoint + epsilonDist * localDirection;
820  const G4ThreeVector PointMinusDir = localPoint - epsilonDist * localDirection;
821  const G4ThreeVector PointPlusNorm = localPoint + epsilonDist * exitNormal;
822  const G4ThreeVector PointMinusNorm = localPoint - epsilonDist * exitNormal;
823 
824  const EInside inPlusDir= solid->Inside(PointPlusDir);
825  const EInside inMinusDir= solid->Inside(PointMinusDir);
826  const EInside inPlusNorm= solid->Inside(PointPlusNorm);
827  const EInside inMinusNorm= solid->Inside(PointMinusNorm);
828 
829  // Basic information
830  os << " Current physical volume = " << physical->GetName() << G4endl;
831  os << " Position (loc) = " << localPoint << G4endl
832  << " Direction (dir) = " << localDirection << G4endl;
833  os << " For confirmation:" << G4endl;
834  os << " Response of DistanceToOut (loc, +dir)= " << distToOut << G4endl;
835  os << " Response of DistanceToOut (loc, -dir)= " << distToOutNeg << G4endl;
836 
837  os << " Inside responds = " << inSolid << " , ie: ";
838  if( inSolid == kOutside ) {
839  os << " Outside -- a problem, as observed in " << fMethod << G4endl;
840  } else if( inSolid == kSurface ) {
841  os << " Surface -- unexpected / inconsistent response ! " << G4endl;
842  } else {
843  os << " Inside -- unexpected / inconsistent response ! " << G4endl;
844  }
845  os << " Obtain safety(ToIn) = " << safetyToIn << G4endl;
846  os << " Obtain safety(ToOut) = " << safetyToOut << G4endl;
847  os << " Response of DistanceToIn (loc, +dir)= " << distToInPos << G4endl;
848  os << " Response of DistanceToIn (loc, -dir)= " << distToInNeg << G4endl;
849 
850  os << " Exit Normal at loc = " << exitNormal << G4endl;
851  os << " Dir . Normal = " << exitNormal.dot( localDirection );
852  os << G4endl;
853 
854  os << " Checking points moved from position by distance/dir. Solid responses: " << G4endl
855  << " +eps in direction : " << G4NavigationLogger_Namespace::EInsideNames[inPlusDir]
856  << " +eps in Normal : " << G4NavigationLogger_Namespace::EInsideNames[inPlusNorm] << G4endl
857  << " -eps in direction : " << G4NavigationLogger_Namespace::EInsideNames[inMinusDir]
858  << " -eps in Normal : " << G4NavigationLogger_Namespace::EInsideNames[inMinusNorm] << G4endl;
859 
860  os << " Parameters of solid: " << G4endl;
861  os << *solid;
862  os << "============================================================";
863 }
G4NavigationLogger::AlongComputeStepLog
void AlongComputeStepLog(const G4VSolid *sampleSolid, const G4ThreeVector &samplePoint, const G4ThreeVector &sampleDirection, const G4ThreeVector &localDirection, G4double sampleSafety, G4double sampleStep) const
Definition: G4NavigationLogger.cc:141
G4VSolid::GetName
G4String GetName() const
G4NavigationLogger::PrintDaughterLog
void PrintDaughterLog(const G4VSolid *sampleSolid, const G4ThreeVector &samplePoint, G4double sampleSafety, G4bool onlySafety, const G4ThreeVector &sampleDirection, G4double sampleStep) const
Definition: G4NavigationLogger.cc:527
kInfinity
static const G4double kInfinity
Definition: geomdefs.hh:42
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4RotationMatrix
CLHEP::HepRotation G4RotationMatrix
Definition: G4RotationMatrix.hh:43
G4VSolid::GetTolerance
G4double GetTolerance() const
G4VSolid::GetEntityType
virtual G4GeometryType GetEntityType() const =0
JustWarning
Definition: G4ExceptionSeverity.hh:64
G4int
int G4int
Definition: G4Types.hh:78
G4VSolid::DumpInfo
void DumpInfo() const
G4NavigationLogger::G4NavigationLogger
G4NavigationLogger(const G4String &id)
Definition: G4NavigationLogger.cc:45
G4NavigationLogger::ReportOutsideMother
void ReportOutsideMother(const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, const G4VPhysicalVolume *motherPV, G4double tDist=30.0 *CLHEP::cm) const
Definition: G4NavigationLogger.cc:667
G4cout
G4GLOB_DLL std::ostream G4cout
G4VPhysicalVolume::GetName
const G4String & GetName() const
G4VSolid::Inside
virtual EInside Inside(const G4ThreeVector &p) const =0
G4NavigationLogger.hh
G4NavigationLogger::ReportVolumeAndIntersection
void ReportVolumeAndIntersection(std::ostream &ostrm, const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, const G4VPhysicalVolume *physical) const
Definition: G4NavigationLogger.cc:787
G4bool
bool G4bool
Definition: G4Types.hh:79
G4VSolid::SurfaceNormal
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const =0
G4NavigationLogger::CheckAndReportBadNormal
G4bool CheckAndReportBadNormal(const G4ThreeVector &unitNormal, const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, G4double step, const G4VSolid *solid, const char *msg) const
Definition: G4NavigationLogger.cc:563
G4VSolid::DistanceToIn
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const =0
G4NavigationLogger_Namespace::EInsideNames
const G4String EInsideNames[3]
Definition: G4NavigationLogger.cc:783
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
perMillion
static const double perMillion
Definition: G4SIunits.hh:331
G4ExceptionSeverity
G4ExceptionSeverity
Definition: G4ExceptionSeverity.hh:59
G4VPhysicalVolume::GetLogicalVolume
G4LogicalVolume * GetLogicalVolume() const
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112
EInside
EInside
Definition: geomdefs.hh:58
G4NavigationLogger::ComputeSafetyLog
void ComputeSafetyLog(const G4VSolid *solid, const G4ThreeVector &point, G4double safety, G4bool isMotherVolume, G4int banner=-1) const
Definition: G4NavigationLogger.cc:494
G4NavigationLogger::PreComputeStepLog
void PreComputeStepLog(const G4VPhysicalVolume *motherPhysical, G4double motherSafety, const G4ThreeVector &localPoint) const
Definition: G4NavigationLogger.cc:59
G4NavigationLogger_Namespace
Definition: G4NavigationLogger.cc:782
FatalException
Definition: G4ExceptionSeverity.hh:60
millimeter
static const double millimeter
Definition: G4SIunits.hh:85
G4endl
#define G4endl
Definition: G4ios.hh:61
G4NavigationLogger::PostComputeStepLog
void PostComputeStepLog(const G4VSolid *motherSolid, const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, G4double motherStep, G4double motherSafety) const
Definition: G4NavigationLogger.cc:443
G4NavigationLogger::CheckDaughterEntryPoint
void CheckDaughterEntryPoint(const G4VSolid *sampleSolid, const G4ThreeVector &samplePoint, const G4ThreeVector &sampleDirection, const G4VSolid *motherSolid, const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, G4double motherStep, G4double sampleStep) const
Definition: G4NavigationLogger.cc:265
G4double
double G4double
Definition: G4Types.hh:76
G4VSolid::DistanceToOut
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const =0
G4GeometryTolerance.hh
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83
G4GeometryTolerance::GetInstance
static G4GeometryTolerance * GetInstance()
Definition: G4GeometryTolerance.cc:73
G4LogicalVolume::GetSolid
G4VSolid * GetSolid() const
G4cerr
G4GLOB_DLL std::ostream G4cerr