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