Geant4-11
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// class G4NavigationLogger Implementation
27//
28// Author: G.Cosmo, 2010
29// --------------------------------------------------------------------
30
31#include <iomanip>
33
34#include "G4NavigationLogger.hh"
36
38
40 : fId(id)
41{
42}
43
45{
46}
47
48// ********************************************************************
49// PreComputeStepLog
50// ********************************************************************
51//
52void
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//
133void
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//
257void
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 {
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 constexpr G4double eps= 1.0e-10;
317 G4bool DaughterEntryIsOutside = SuspiciousDaughterDist
318 && ( (sampleStep * (1.0+eps) < distToReEntry) || (insideMother == kOutside ) );
319 G4bool EntryIsMotherExit = std::fabs(sampleStep-motherStep) < kCarTolerance;
320
321 // Check for more subtle error - is exit point of daughter correct ?
322 G4ThreeVector sampleEntryPoint = samplePoint+sampleStep*sampleDirection;
323 G4double sampleCrossingDist = sampleSolid->DistanceToOut( sampleEntryPoint,
324 sampleDirection );
325 G4double sampleExitDist = sampleStep+sampleCrossingDist;
326 G4ThreeVector sampleExitPoint = samplePoint+sampleExitDist*sampleDirection;
327
328 G4bool TransitProblem = ( (sampleStep < motherStep)
329 && (sampleExitDist > motherStep + kCarTolerance) )
330 || ( EntryIsMotherExit && (sampleCrossingDist > kCarTolerance) );
331
332 if( DaughterEntryIsOutside
333 || TransitProblem
334 || (SuspiciousDaughterDist && (fVerbose > 3) ) )
335 {
337 msg.precision(16);
338
339 if( DaughterEntryIsOutside )
340 {
341 msg << "WARNING> Intersection distance to Daughter volume is further"
342 << " than the distance to boundary." << G4endl
343 << " It appears that part of the daughter volume is *outside*"
344 << " this mother. " << G4endl;
345 msg << " One of the following checks signaled a problem:" << G4endl
346 << " -sampleStep (dist to daugh) < mother-exit dist + distance "
347 << "to ReEntry point for mother " << G4endl
348 << " -position of daughter intersection is outside mother volume."
349 << G4endl;
350 }
351 else if( TransitProblem )
352 {
353 G4double protrusion = sampleExitDist - motherStep;
354
355 msg << "WARNING> Daughter volume extends beyond mother boundary. "
356 << G4endl;
357 if ( ( sampleStep < motherStep )
358 && (sampleExitDist > motherStep + kCarTolerance ) )
359 {
360 // 1st Issue with Daughter
361 msg << " Crossing distance in the daughter causes is to extend"
362 << " beyond the mother exit. " << G4endl;
363 msg << " Length protruding = " << protrusion << G4endl;
364 }
365 if( EntryIsMotherExit )
366 {
367 // 1st Issue with Daughter
368 msg << " Intersection distance to Daughter is within "
369 << " tolerance of the distance" << G4endl;
370 msg << " to the mother boundary * and * " << G4endl;
371 msg << " the crossing distance in the daughter is > tolerance."
372 << G4endl;
373 }
374 }
375 else
376 {
377 msg << "NearMiss> Intersection to Daughter volume is in extension past the"
378 << " current exit point of the mother volume." << G4endl;
379 msg << " This is not an error - just an unusual occurrence,"
380 << " possible in the case of concave volume. " << G4endl;
381 }
382 msg << "---- Information about intersection with daughter, mother: "
383 << G4endl;
384 msg << " sampleStep (daughter) = " << sampleStep << G4endl
385 << " motherStep = " << motherStep << G4endl
386 << " distToRentry(mother) = " << distToReEntry << G4endl
387 << " Inside(entry pnt daug): " << solidResponse << G4endl
388 << " dist across daughter = " << sampleCrossingDist << G4endl;
389 msg << " Mother Name (Solid) : " << motherSolid->GetName() << G4endl
390 << " In local (mother) coordinates: " << G4endl
391 << " Starting Point = " << localPoint << G4endl
392 << " Direction = " << localDirection << G4endl
393 << " Exit Point (mother)= " << localExitMotherPos << G4endl
394 << " Entry Point (daughter)= " << localPoint+sampleStep*localDirection
395 << G4endl;
396 if( distToReEntry < kInfinity )
397 {
398 msg << " ReEntry Point (mother)= " << localReEntryPoint << G4endl;
399 }
400 else
401 {
402 msg << " No ReEntry - track does not encounter mother volume again! "
403 << G4endl;
404 }
405 msg << " Daughter Name (Solid): " << sampleSolid->GetName() << G4endl
406 << " In daughter coordinates: " << G4endl
407 << " Starting Point = " << samplePoint << G4endl
408 << " Direction = " << sampleDirection << G4endl
409 << " Entry Point (daughter)= " << sampleEntryPoint
410 << G4endl;
411 msg << " Description of mother solid: " << G4endl
412 << *motherSolid << G4endl
413 << " Description of daughter solid: " << G4endl
414 << *sampleSolid << G4endl;
415 G4String fType = fId + "::ComputeStep()";
416
417 if( DaughterEntryIsOutside || TransitProblem )
418 {
419 G4Exception(fType, "GeomNav0003", JustWarning, msg);
420 }
421 else
422 {
423 G4cout << fType
424 << " -- Checked distance of Entry to daughter vs exit of mother"
425 << G4endl;
426 G4cout << msg.str();
427 G4cout << G4endl;
428 }
429 }
430}
431
432// ********************************************************************
433// PostComputeStepLog
434// ********************************************************************
435//
436void
438 const G4ThreeVector& localPoint,
439 const G4ThreeVector& localDirection,
440 G4double motherStep,
441 G4double motherSafety) const
442{
443 if ( fVerbose == 1 || fVerbose > 4 )
444 {
445 G4cout << " Mother "
446 << std::setw(15) << motherSafety << " "
447 << std::setw(15) << motherStep << " " << localPoint << " - "
448 << motherSolid->GetEntityType() << ": " << motherSolid->GetName()
449 << G4endl;
450 }
451 if( ( motherStep < 0.0 ) || ( motherStep >= kInfinity) )
452 {
453 G4String fType = fId + "::ComputeStep()";
454 G4int oldPrOut = G4cout.precision(16);
455 G4int oldPrErr = G4cerr.precision(16);
456 std::ostringstream message;
457 message << "Current point is outside the current solid !" << G4endl
458 << " Problem in Navigation" << G4endl
459 << " Point (local coordinates): "
460 << localPoint << G4endl
461 << " Local Direction: " << localDirection << G4endl
462 << " Solid: " << motherSolid->GetName();
463 motherSolid->DumpInfo();
464 G4Exception(fType, "GeomNav0003", FatalException, message);
465 G4cout.precision(oldPrOut);
466 G4cerr.precision(oldPrErr);
467 }
468 if ( fVerbose > 1 )
469 {
470 static const G4int precVerf = 20; // Precision
471 G4int oldprec = G4cout.precision(precVerf);
472 G4cout << " Mother " << std::setw(12) << motherSolid->GetName() << " "
473 << std::setw(4+precVerf) << localPoint << " "
474 << std::setw(4+precVerf) << motherSafety << " "
475 << std::setw(4+precVerf) << motherStep << " "
476 << std::setw(16) << "distanceToOut" << " "
477 << std::setw(4+precVerf) << localDirection << " "
478 << G4endl;
479 G4cout.precision(oldprec);
480 }
481}
482
483// ********************************************************************
484// ComputeSafetyLog
485// ********************************************************************
486//
487void
489 const G4ThreeVector& point,
490 G4double safety,
491 G4bool isMotherVolume,
492 G4int banner) const
493{
494 if( banner < 0 )
495 {
496 banner = isMotherVolume;
497 }
498 if( fVerbose >= 1 )
499 {
500 G4String volumeType = isMotherVolume ? " Mother " : "Daughter";
501 if (banner)
502 {
503 G4cout << "************** " << fId << "::ComputeSafety() ****************"
504 << G4endl;
505 G4cout << " VolType "
506 << std::setw(15) << "Safety/mm" << " "
507 << std::setw(52) << "Position (local coordinates)"
508 << " - Solid" << G4endl;
509 }
510 G4cout << volumeType
511 << std::setw(15) << safety << " " << point << " - "
512 << solid->GetEntityType() << ": " << solid->GetName() << G4endl;
513 }
514}
515
516// ********************************************************************
517// PrintDaughterLog
518// ********************************************************************
519//
520void
522 const G4ThreeVector& samplePoint,
523 G4double sampleSafety,
524 G4bool withStep,
525 const G4ThreeVector& sampleDirection, G4double sampleStep ) const
526{
527 if ( fVerbose >= 1 )
528 {
529 G4int oldPrec = G4cout.precision(8);
530 G4cout << "Daughter "
531 << std::setw(15) << sampleSafety << " ";
532 if (withStep) // (sampleStep != -1.0 )
533 {
534 G4cout << std::setw(15) << sampleStep << " ";
535 }
536 else
537 {
538 G4cout << std::setw(15) << "Not-Available" << " ";
539 }
540 G4cout << samplePoint << " - "
541 << sampleSolid->GetEntityType() << ": " << sampleSolid->GetName();
542 if( withStep )
543 {
544 G4cout << " dir= " << sampleDirection;
545 }
546 G4cout << G4endl;
547 G4cout.precision(oldPrec);
548 }
549}
550
551// ********************************************************************
552// CheckAndReportBadNormal
553// ********************************************************************
554//
555G4bool
557CheckAndReportBadNormal(const G4ThreeVector& unitNormal,
558 const G4ThreeVector& localPoint,
559 const G4ThreeVector& localDirection,
560 G4double step,
561 const G4VSolid* solid,
562 const char* msg ) const
563{
564 G4double normMag2 = unitNormal.mag2();
565 G4bool badLength = ( std::fabs ( normMag2 - 1.0 ) > CLHEP::perMillion );
566
567 if( badLength )
568 {
569 G4double normMag = std::sqrt(normMag2);
571 message.precision(10);
572 message << "============================================================"
573 << G4endl;
574 message << " WARNING> Normal is not a unit vector. "
575 << " - but |normal| = " << normMag
576 << " - and |normal|^2 = " << normMag2 << G4endl
577 << " which differ from 1.0 by: " << G4endl
578 << " |normal|-1 = " << normMag - 1.0
579 << " and |normal|^2 - 1 = " << normMag2 - 1.0 << G4endl
580 << " n = " << unitNormal << G4endl;
581 message << " Info string: " << msg << G4endl;
582 message << "============================================================"
583 << G4endl;
584
585 message.precision(16);
586
587 message << " Information on call to DistanceToOut: " << G4endl;
588 message << " Position = " << localPoint << G4endl
589 << " Direction = " << localDirection << G4endl;
590 message << " Obtained> distance = " << step << G4endl;
591 message << " > Exit position = " << localPoint+step*localDirection
592 << G4endl;
593 message << " Parameters of solid: " << G4endl;
594 message << *solid;
595 message << "============================================================";
596
597 G4String fMethod = fId + "::ComputeStep()";
598 G4Exception( fMethod, "GeomNav0003", JustWarning, message);
599 }
600 return badLength;
601}
602
603// ********************************************************************
604// CheckAndReportBadNormal - due to Rotation Matrix
605// ********************************************************************
606//
607G4bool
609CheckAndReportBadNormal(const G4ThreeVector& rotatedNormal,
610 const G4ThreeVector& originalNormal,
611 const G4RotationMatrix& rotationM,
612 const char* msg ) const
613{
614 G4double normMag2 = rotatedNormal.mag2();
615 G4bool badLength = ( std::fabs ( normMag2 - 1.0 ) > CLHEP::perMillion );
616
617 if( badLength )
618 {
619 G4double normMag = std::sqrt(normMag2);
621 message.precision(10);
622 message << "============================================================"
623 << G4endl;
624 message << " WARNING> Rotated n(ormal) is not a unit vector. " << G4endl
625 << " |normal| = " << normMag
626 << " and |normal|^2 = " << normMag2 << G4endl
627 << " Diff from 1.0: " << G4endl
628 << " |normal|-1 = " << normMag - 1.0
629 << " and |normal|^2 - 1 = " << normMag2 - 1.0 << G4endl;
630 message << " Rotated n = (" << rotatedNormal.x() << "," << rotatedNormal.y() << ","
631 << rotatedNormal.z() << ")" << G4endl;
632 message << " Original n = (" << originalNormal.x() << "," << originalNormal.y() << ","
633 << originalNormal.z() << ")" << G4endl;
634 message << " Info string: " << msg << G4endl;
635 message << "============================================================"
636 << G4endl;
637
638 message.precision(16);
639
640 message << " Information on RotationMatrix : " << G4endl;
641 message << " Original: " << G4endl;
642 message << rotationM << G4endl;
643 message << " Inverse (used in transformation): " << G4endl;
644 message << rotationM.inverse() << G4endl;
645 message << "============================================================";
646
647 G4String fMethod = fId + "::ComputeStep()";
648 G4Exception( fMethod, "GeomNav0003", JustWarning, message);
649 }
650 return badLength;
651}
652
653// ********************************************************************
654// ReportOutsideMother
655// ********************************************************************
656//
657// Report Exception if point is outside mother.
658// Fatal exception will be used if either 'checkMode error is > triggerDist
659//
660void
662 const G4ThreeVector& localDirection,
663 const G4VPhysicalVolume* physical,
664 G4double triggerDist) const
665{
666 const G4LogicalVolume* logicalVol = physical
667 ? physical->GetLogicalVolume() : nullptr;
668 const G4VSolid* solid = logicalVol
669 ? logicalVol->GetSolid() : nullptr;
670
671 G4String fMethod = fId + "::ComputeStep()";
672
673 if( solid == nullptr )
674 {
675 G4Exception(fMethod, "GeomNav0003", FatalException,
676 "Erroneous call to ReportOutsideMother: no Solid is available");
677 return;
678 }
679 const G4double kCarTolerance = solid->GetTolerance();
680
681 // Double check reply - it should be kInfinity
682 const G4double distToOut = solid->DistanceToOut(localPoint, localDirection);
683 const EInside inSolid = solid->Inside(localPoint);
684 const G4double safetyToIn = solid->DistanceToIn(localPoint);
685 const G4double safetyToOut = solid->DistanceToOut(localPoint);
686 // const G4double distToInPos =
687 // solid->DistanceToIn(localPoint, localDirection);
688
689 // 1. Check consistency between Safety obtained and report from distance
690 // We must ensure that (mother)Safety <= 0.0
691 // in the case that the point is outside!
692 // [ If this is not the case, this problem can easily go undetected,
693 // except in Check mode ! ]
694 if( safetyToOut > kCarTolerance
695 && ( distToOut < 0.0 || distToOut >= kInfinity ) )
696 {
698 // fNavClerk->ReportBadSafety(localPoint, localDirection,
699 // motherPhysical, motherSafety, motherStep);
700 msg1 << " Dangerous inconsistency in response of solid." << G4endl
701 << " Solid type: " << solid->GetEntityType()
702 << " Name= " << solid->GetName() << G4endl;
703 msg1 << " Mother volume gives safety > 0 despite being called for *Outside* point "
704 << G4endl
705 << " Location = " << localPoint << G4endl
706 << " Direction= " << localDirection << G4endl
707 << " - Safety (Isotropic d) = " << safetyToOut << G4endl
708 << " - Intersection Distance= " << distToOut << G4endl
709 << G4endl;
710 G4Exception( fMethod, "GeomNav0123", JustWarning, msg1);
711 }
712
713 // 2. Inconsistency - Too many distances are zero (or will be rounded to zero)
714
715// if( std::fabs(distToOut) < kCarTolerance
716// && std::fabs(distToInPos) < kCarTolerance )
717// {
718 // If both distanceToIn and distanceToOut (p,v) are zero for
719 // one direction, the particle could get stuck!
720// }
721
723 msg.precision(10);
724
725 if( std::fabs(distToOut) < kCarTolerance )
726 {
727 // 3. Soft error - safety is not rounded to zero
728 // Report nothing - except in 'loud' mode
730 {
731 msg << " Warning> DistanceToOut(p,v): "
732 << "Distance from surface is not rounded to zero" << G4endl;
733 }
734 else
735 {
736 return;
737 }
738 }
739 else
740 {
741 // 4. General message - complain that the point is outside!
742 // and provide all information about the current location,
743 // direction and the answers of the solid
744 msg << "============================================================"
745 << G4endl;
746 msg << " WARNING> Current Point appears to be Outside mother volume !! "
747 << G4endl;
748 msg << " Response of DistanceToOut was negative or kInfinity"
749 << " when called in " << fMethod << G4endl;
750 }
751
752 // Generate and 'print'/stream all the information needed
753 this->ReportVolumeAndIntersection(msg, localPoint, localDirection, physical);
754
755 // Default for distance of 'major' error
756 if( triggerDist <= 0.0 )
757 {
758 triggerDist = std::max ( 1.0e+6 * kCarTolerance, // Well beyond tolerance
760 }
761
762 G4bool majorError = inSolid == kOutside
763 ? ( safetyToIn > triggerDist )
764 : ( safetyToOut > triggerDist );
765
766 G4ExceptionSeverity exceptionType = JustWarning;
767 if ( majorError )
768 {
769 exceptionType = FatalException;
770 }
771
772 G4Exception( fMethod, "GeomNav0003", exceptionType, msg);
773}
774
776{
777 const G4String EInsideNames[3] = { "kOutside", "kSurface", "kInside" };
778}
779
781ReportVolumeAndIntersection( std::ostream& os,
782 const G4ThreeVector& localPoint,
783 const G4ThreeVector& localDirection,
784 const G4VPhysicalVolume* physical ) const
785{
786 G4String fMethod = fId + "::ComputeStep()";
787 const G4LogicalVolume* logicalVol = physical
788 ? physical->GetLogicalVolume() : nullptr;
789 const G4VSolid* solid = logicalVol
790 ? logicalVol->GetSolid() : nullptr;
791 if( solid == nullptr )
792 {
793 os << " ERROR> Solid is not available. Logical Volume = "
794 << logicalVol << std::endl;
795 return;
796 }
797 const G4double kCarTolerance = solid->GetTolerance();
798
799 // Double check reply - it should be kInfinity
800 const G4double distToOut = solid->DistanceToOut(localPoint, localDirection);
801 const G4double distToOutNeg = solid->DistanceToOut(localPoint,
802 -localDirection);
803 const EInside inSolid = solid->Inside(localPoint);
804 const G4double safetyToIn = solid->DistanceToIn(localPoint);
805 const G4double safetyToOut = solid->DistanceToOut(localPoint);
806
807 const G4double distToInPos = solid->DistanceToIn(localPoint, localDirection);
808 const G4double distToInNeg = solid->DistanceToIn(localPoint, -localDirection);
809
810 const G4ThreeVector exitNormal = solid->SurfaceNormal(localPoint);
811
812 // Double check whether points nearby are in/surface/out
813 const G4double epsilonDist = 1000.0 * kCarTolerance;
814 const G4ThreeVector PointPlusDir = localPoint + epsilonDist * localDirection;
815 const G4ThreeVector PointMinusDir = localPoint - epsilonDist * localDirection;
816 const G4ThreeVector PointPlusNorm = localPoint + epsilonDist * exitNormal;
817 const G4ThreeVector PointMinusNorm = localPoint - epsilonDist * exitNormal;
818
819 const EInside inPlusDir = solid->Inside(PointPlusDir);
820 const EInside inMinusDir = solid->Inside(PointMinusDir);
821 const EInside inPlusNorm = solid->Inside(PointPlusNorm);
822 const EInside inMinusNorm = solid->Inside(PointMinusNorm);
823
824 // Basic information
825 os << " Current physical volume = " << physical->GetName() << G4endl;
826 os << " Position (loc) = " << localPoint << G4endl
827 << " Direction (dir) = " << localDirection << G4endl;
828 os << " For confirmation:" << G4endl;
829 os << " Response of DistanceToOut (loc, +dir)= " << distToOut << G4endl;
830 os << " Response of DistanceToOut (loc, -dir)= " << distToOutNeg << G4endl;
831
832 os << " Inside responds = " << inSolid << " , ie: ";
833 if( inSolid == kOutside )
834 {
835 os << " Outside -- a problem, as observed in " << fMethod << G4endl;
836 }
837 else if( inSolid == kSurface )
838 {
839 os << " Surface -- unexpected / inconsistent response ! " << G4endl;
840 }
841 else
842 {
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/direction." << G4endl
855 << " Solid responses: " << G4endl
856 << " +eps in direction : "
858 << " +eps in Normal : "
860 << " -eps in direction : "
862 << " -eps in Normal : "
864
865 os << " Parameters of solid: " << G4endl;
866 os << *solid;
867 os << "============================================================";
868}
const G4double kCarTolerance
static const G4double eps
G4ExceptionSeverity
@ JustWarning
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:40
static constexpr double millimeter
Definition: G4SIunits.hh:66
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
G4GLOB_DLL std::ostream G4cerr
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
double z() const
double x() const
double mag2() const
double y() const
double dot(const Hep3Vector &) const
HepRotation inverse() const
static G4GeometryTolerance * GetInstance()
G4VSolid * GetSolid() const
void PreComputeStepLog(const G4VPhysicalVolume *motherPhysical, G4double motherSafety, const G4ThreeVector &localPoint) const
G4NavigationLogger(const G4String &id)
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
void PrintDaughterLog(const G4VSolid *sampleSolid, const G4ThreeVector &samplePoint, G4double sampleSafety, G4bool onlySafety, const G4ThreeVector &sampleDirection, G4double sampleStep) const
G4bool CheckAndReportBadNormal(const G4ThreeVector &unitNormal, const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, G4double step, const G4VSolid *solid, const char *msg) const
void ComputeSafetyLog(const G4VSolid *solid, const G4ThreeVector &point, G4double safety, G4bool isMotherVolume, G4int banner=-1) const
void ReportVolumeAndIntersection(std::ostream &ostrm, const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, const G4VPhysicalVolume *physical) const
void PostComputeStepLog(const G4VSolid *motherSolid, const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, G4double motherStep, G4double motherSafety) const
void ReportOutsideMother(const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, const G4VPhysicalVolume *motherPV, G4double tDist=30.0 *CLHEP::cm) const
void AlongComputeStepLog(const G4VSolid *sampleSolid, const G4ThreeVector &samplePoint, const G4ThreeVector &sampleDirection, const G4ThreeVector &localDirection, G4double sampleSafety, G4double sampleStep) const
G4LogicalVolume * GetLogicalVolume() const
const G4String & GetName() const
G4String GetName() const
G4double GetTolerance() const
virtual EInside Inside(const G4ThreeVector &p) const =0
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=nullptr, G4ThreeVector *n=nullptr) const =0
void DumpInfo() const
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const =0
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const =0
virtual G4GeometryType GetEntityType() const =0
EInside
Definition: geomdefs.hh:67
@ kInside
Definition: geomdefs.hh:70
@ kOutside
Definition: geomdefs.hh:68
@ kSurface
Definition: geomdefs.hh:69
static const G4double kInfinity
Definition: geomdefs.hh:41
static constexpr double perMillion
static constexpr double millimeter
Definition: SystemOfUnits.h:63
T max(const T t1, const T t2)
brief Return the largest of the two arguments