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$ >> 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 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 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 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