Geant4 10.7.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: G.Cosmo, 2010
29// --------------------------------------------------------------------
30
31#include <iomanip>
33
34#include "G4NavigationLogger.hh"
36
37using CLHEP::millimeter;
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 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 {
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//
435void
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//
486void
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//
519void
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//
554G4bool
556CheckAndReportBadNormal(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);
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//
606G4bool
608CheckAndReportBadNormal(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);
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//
659void
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 {
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
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
775{
776 const G4String EInsideNames[3] = { "kOutside", "kSurface", "kInside" };
777}
778
780ReportVolumeAndIntersection( 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 : "
857 << " +eps in Normal : "
859 << " -eps in direction : "
861 << " -eps in Normal : "
863
864 os << " Parameters of solid: " << G4endl;
865 os << *solid;
866 os << "============================================================";
867}
const G4double kCarTolerance
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
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