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