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