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