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 // $Id: G4SimpleLocator.cc,v 1.5 2008/12/11 10:27:58 tnikitin Exp $ >> 27 // GEANT4 tag $Name: geant4-09-02 $ >> 28 // 26 // Class G4SimpleLocator implementation 29 // Class G4SimpleLocator implementation 27 // 30 // 28 // 27.10.08 - Tatiana Nikitina, extracted from << 31 // 27.10.08 - Tatiana Nikitina. 29 // 04.10.11 - John Apostolakis, revised conver << 30 // ------------------------------------------- 32 // --------------------------------------------------------------------------- 31 33 32 #include <iomanip> 34 #include <iomanip> 33 35 34 #include "G4ios.hh" 36 #include "G4ios.hh" 35 #include "G4SimpleLocator.hh" 37 #include "G4SimpleLocator.hh" 36 38 37 G4SimpleLocator::G4SimpleLocator(G4Navigator * 39 G4SimpleLocator::G4SimpleLocator(G4Navigator *theNavigator) 38 : G4VIntersectionLocator(theNavigator) 40 : G4VIntersectionLocator(theNavigator) 39 { 41 { 40 } 42 } 41 43 42 G4SimpleLocator::~G4SimpleLocator() = default; << 44 G4SimpleLocator::~G4SimpleLocator() >> 45 { >> 46 } 43 47 44 // ------------------------------------------- 48 // -------------------------------------------------------------------------- 45 // G4bool G4PropagatorInField::LocateIntersect 49 // G4bool G4PropagatorInField::LocateIntersectionPoint( 46 // const G4FieldTrack& CurveStart 50 // const G4FieldTrack& CurveStartPointVelocity, // A 47 // const G4FieldTrack& CurveEndPo 51 // const G4FieldTrack& CurveEndPointVelocity, // B 48 // const G4ThreeVector& TrialPoint 52 // const G4ThreeVector& TrialPoint, // E 49 // G4FieldTrack& Intersecte 53 // G4FieldTrack& IntersectedOrRecalculated // Output 50 // G4bool& recalculat 54 // G4bool& recalculated ) // Out 51 // ------------------------------------------- 55 // -------------------------------------------------------------------------- 52 // 56 // 53 // Function that returns the intersection of t 57 // Function that returns the intersection of the true path with the surface 54 // of the current volume (either the external 58 // of the current volume (either the external one or the inner one with one 55 // of the daughters: 59 // of the daughters: 56 // 60 // 57 // A = Initial point 61 // A = Initial point 58 // B = another point 62 // B = another point 59 // 63 // 60 // Both A and B are assumed to be on the true 64 // Both A and B are assumed to be on the true path: 61 // 65 // 62 // E is the first point of intersection of 66 // E is the first point of intersection of the chord AB with 63 // a volume other than A (on the surface 67 // a volume other than A (on the surface of A or of a daughter) 64 // 68 // 65 // Convention of Use : 69 // Convention of Use : 66 // i) If it returns "true", then Intersect 70 // i) If it returns "true", then IntersectionPointVelocity is set 67 // to the approximate intersection point 71 // to the approximate intersection point. 68 // ii) If it returns "false", no intersecti 72 // ii) If it returns "false", no intersection was found. 69 // The validity of IntersectedOrRecal 73 // The validity of IntersectedOrRecalculated depends on 'recalculated' 70 // a) if latter is false, then Intersec 74 // a) if latter is false, then IntersectedOrRecalculated is invalid. 71 // b) if latter is true, then Intersec 75 // b) if latter is true, then IntersectedOrRecalculated is 72 // the new endpoint, due to a re-i 76 // the new endpoint, due to a re-integration. 73 // ------------------------------------------- 77 // -------------------------------------------------------------------------- 74 // NOTE: implementation taken from G4Propagato 78 // NOTE: implementation taken from G4PropagatorInField 75 // 79 // 76 G4bool G4SimpleLocator::EstimateIntersectionPo 80 G4bool G4SimpleLocator::EstimateIntersectionPoint( 77 const G4FieldTrack& CurveStart 81 const G4FieldTrack& CurveStartPointVelocity, // A 78 const G4FieldTrack& CurveEndPo 82 const G4FieldTrack& CurveEndPointVelocity, // B 79 const G4ThreeVector& TrialPoint 83 const G4ThreeVector& TrialPoint, // E 80 G4FieldTrack& Intersecte 84 G4FieldTrack& IntersectedOrRecalculatedFT, // Output 81 G4bool& recalculat 85 G4bool& recalculatedEndPoint, // Out 82 G4double &fPrevious 86 G4double &fPreviousSafety, //In/Out 83 G4ThreeVector &fPrevious 87 G4ThreeVector &fPreviousSftOrigin) //In/Out 84 { 88 { 85 // Find Intersection Point ( A, B, E ) of t 89 // Find Intersection Point ( A, B, E ) of true path AB - start at E. 86 90 87 G4bool found_approximate_intersection = fals 91 G4bool found_approximate_intersection = false; 88 G4bool there_is_no_intersection = fals 92 G4bool there_is_no_intersection = false; 89 93 90 G4FieldTrack CurrentA_PointVelocity = Curve 94 G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity; 91 G4FieldTrack CurrentB_PointVelocity = Curve 95 G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity; 92 G4ThreeVector CurrentE_Point = TrialPoint; 96 G4ThreeVector CurrentE_Point = TrialPoint; 93 G4bool validNormalAtE = false; << 94 G4ThreeVector NormalAtEntry; << 95 << 96 G4FieldTrack ApproxIntersecPointV(CurveEndP 97 G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity); // FT-Def-Construct 97 G4double NewSafety = 0.0; 98 G4double NewSafety = 0.0; 98 G4bool last_AF_intersection = false; 99 G4bool last_AF_intersection = false; 99 G4bool final_section = true; // Shows wheth 100 G4bool final_section = true; // Shows whether current section is last 100 // (i.e. B=ful 101 // (i.e. B=full end) 101 recalculatedEndPoint = false; 102 recalculatedEndPoint = false; 102 103 103 G4bool restoredFullEndpoint = false; 104 G4bool restoredFullEndpoint = false; 104 105 105 G4int substep_no = 0; 106 G4int substep_no = 0; 106 107 107 // Limits for substep number 108 // Limits for substep number 108 // 109 // 109 const G4int max_substeps = 100000000; // T 110 const G4int max_substeps = 100000000; // Test 120 (old value 100 ) 110 const G4int warn_substeps = 1000; // 111 const G4int warn_substeps = 1000; // 100 111 112 112 // Statistics for substeps 113 // Statistics for substeps 113 // 114 // 114 static G4ThreadLocal G4int max_no_seen= -1; << 115 static G4int max_no_seen= -1; 115 << 116 static G4int trigger_substepno_print= warn_substeps - 20; 116 NormalAtEntry = GetSurfaceNormal( CurrentE_P << 117 117 << 118 #ifdef G4DEBUG_FIELD 118 #ifdef G4DEBUG_FIELD 119 const G4double tolerance = 1.0e-8; << 119 static G4double tolerance= 1.0e-8; 120 G4ThreeVector StartPosition= CurveStartPoin 120 G4ThreeVector StartPosition= CurveStartPointVelocity.GetPosition(); 121 if( (TrialPoint - StartPosition).mag() < tol << 121 if( (TrialPoint - StartPosition).mag() < tolerance * mm ) 122 { 122 { >> 123 G4cerr << "WARNING - G4SimpleLocator::EstimateIntersectionPoint()" >> 124 << G4endl >> 125 << " Intermediate F point is on top of starting point A." >> 126 << G4endl; 123 G4Exception("G4SimpleLocator::EstimateInt 127 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()", 124 "GeomNav1002", JustWarning, << 128 "IntersectionPointIsAtStart", JustWarning, 125 "Intersection point F is exac 129 "Intersection point F is exactly at start point A." ); 126 } 130 } 127 #endif 131 #endif 128 132 129 do 133 do 130 { 134 { 131 G4ThreeVector Point_A = CurrentA_PointVel 135 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition(); 132 G4ThreeVector Point_B = CurrentB_PointVel 136 G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition(); 133 137 134 // F = a point on true AB path close to p 138 // F = a point on true AB path close to point E 135 // (the closest if possible) 139 // (the closest if possible) 136 // 140 // 137 ApproxIntersecPointV = GetChordFinderFor( 141 ApproxIntersecPointV = GetChordFinderFor() 138 ->ApproxCurvePoint 142 ->ApproxCurvePointV( CurrentA_PointVelocity, 139 143 CurrentB_PointVelocity, 140 144 CurrentE_Point, 141 145 GetEpsilonStepFor()); 142 // The above method is the key & most 146 // The above method is the key & most intuitive part ... 143 147 144 #ifdef G4DEBUG_FIELD 148 #ifdef G4DEBUG_FIELD 145 if( ApproxIntersecPointV.GetCurveLength() 149 if( ApproxIntersecPointV.GetCurveLength() > 146 CurrentB_PointVelocity.GetCurveLength 150 CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) ) 147 { 151 { >> 152 G4cerr << "ERROR - G4SimpleLocator::EstimateIntersectionPoint()" >> 153 << G4endl >> 154 << " Intermediate F point is more advanced than" >> 155 << " endpoint B." << G4endl; 148 G4Exception("G4SimpleLocator::EstimateI 156 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()", 149 "GeomNav0003", FatalExcepti << 157 "IntermediatePointConfusion", FatalException, 150 "Intermediate F point is pa << 158 "Intermediate F point is past end B point" ); 151 } 159 } 152 #endif 160 #endif 153 161 154 G4ThreeVector CurrentF_Point = ApproxInte << 162 G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition(); 155 163 156 // First check whether EF is small - then 164 // First check whether EF is small - then F is a good approx. point 157 // Calculate the length and direction of 165 // Calculate the length and direction of the chord AF 158 // 166 // 159 G4ThreeVector ChordEF_Vector = CurrentF_P << 167 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point; 160 168 161 G4ThreeVector NewMomentumDir = ApproxInte << 169 if ( ChordEF_Vector.mag2() <= sqr(GetDeltaIntersectionFor()) ) 162 G4double MomDir_dot_Norm = NewMomentumDir << 163 << 164 G4ThreeVector ChordAB = Point_B - Point_A << 165 << 166 #ifdef G4DEBUG_FIELD << 167 G4VIntersectionLocator:: << 168 ReportTrialStep( substep_no, ChordAB, C << 169 NewMomentumDir, Normal << 170 #endif << 171 // Check Sign is always exiting !! TODO << 172 // Could ( > -epsilon) be used instead? << 173 // << 174 G4bool adequate_angle = ( MomDir_dot_Norm << 175 || (! validNormalAtE << 176 G4double EF_dist2= ChordEF_Vector.mag2(); << 177 if ( ( EF_dist2 <= sqr(fiDeltaIntersecti << 178 || ( EF_dist2 <= kCarTolerance*kCarTole << 179 { 170 { 180 found_approximate_intersection = true; 171 found_approximate_intersection = true; 181 172 182 // Create the "point" return value 173 // Create the "point" return value 183 // 174 // >> 175 184 IntersectedOrRecalculatedFT = ApproxInt 176 IntersectedOrRecalculatedFT = ApproxIntersecPointV; 185 IntersectedOrRecalculatedFT.SetPosition 177 IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point ); 186 178 187 if ( GetAdjustementOfFoundIntersection( 179 if ( GetAdjustementOfFoundIntersection() ) 188 { 180 { 189 // Try to Get Correction of Intersect 181 // Try to Get Correction of IntersectionPoint using SurfaceNormal() 190 // 182 // 191 G4ThreeVector IP; 183 G4ThreeVector IP; 192 G4ThreeVector MomentumDir= ApproxInte 184 G4ThreeVector MomentumDir= ApproxIntersecPointV.GetMomentumDirection(); 193 G4bool goodCorrection = AdjustmentOfF 185 G4bool goodCorrection = AdjustmentOfFoundIntersection( Point_A, 194 CurrentE_Po 186 CurrentE_Point, CurrentF_Point, MomentumDir, 195 last_AF_int 187 last_AF_intersection, IP, NewSafety, 196 fPreviousSa 188 fPreviousSafety, fPreviousSftOrigin ); 197 189 198 if ( goodCorrection ) << 190 if(goodCorrection) 199 { 191 { 200 IntersectedOrRecalculatedFT = Appro 192 IntersectedOrRecalculatedFT = ApproxIntersecPointV; 201 IntersectedOrRecalculatedFT.SetPosi 193 IntersectedOrRecalculatedFT.SetPosition(IP); 202 } 194 } 203 } 195 } 204 196 205 // Note: in order to return a point on 197 // Note: in order to return a point on the boundary, 206 // we must return E. But it is F 198 // we must return E. But it is F on the curve. 207 // So we must "cheat": we are usi 199 // So we must "cheat": we are using the position at point E 208 // and the velocity at point F !! 200 // and the velocity at point F !!! 209 // 201 // 210 // This must limit the length we can al 202 // This must limit the length we can allow for displacement! 211 } 203 } 212 else // E is NOT close enough to the cur 204 else // E is NOT close enough to the curve (ie point F) 213 { 205 { 214 // Check whether any volumes are encoun 206 // Check whether any volumes are encountered by the chord AF 215 // ------------------------------------ 207 // --------------------------------------------------------- 216 // First relocate to restore any Voxel 208 // First relocate to restore any Voxel etc information 217 // in the Navigator before calling Comp 209 // in the Navigator before calling ComputeStep() 218 // 210 // 219 GetNavigatorFor()->LocateGlobalPointWit 211 GetNavigatorFor()->LocateGlobalPointWithinVolume( Point_A ); 220 212 221 G4ThreeVector PointG; // Candidate in 213 G4ThreeVector PointG; // Candidate intersection point 222 G4double stepLengthAF; 214 G4double stepLengthAF; 223 G4bool usedNavigatorAF = false; << 215 G4bool Intersects_AF = IntersectChord( Point_A, CurrentF_Point, 224 G4bool Intersects_AF = IntersectChord( << 216 NewSafety,fPreviousSafety, 225 << 226 << 227 << 228 217 fPreviousSftOrigin, 229 218 stepLengthAF, 230 << 219 PointG ); 231 << 232 last_AF_intersection = Intersects_AF; 220 last_AF_intersection = Intersects_AF; 233 if( Intersects_AF ) 221 if( Intersects_AF ) 234 { 222 { 235 // G is our new Candidate for the int 223 // G is our new Candidate for the intersection point. 236 // It replaces "E" and we will repea 224 // It replaces "E" and we will repeat the test to see if 237 // it is a good enough approximate po 225 // it is a good enough approximate point for us. 238 // B <- F 226 // B <- F 239 // E <- G 227 // E <- G 240 228 241 CurrentB_PointVelocity = ApproxInters 229 CurrentB_PointVelocity = ApproxIntersecPointV; 242 CurrentE_Point = PointG; << 230 CurrentE_Point = PointG; 243 << 244 // Need to recalculate the Exit Norma << 245 // Relies on a call to Navigator::Com << 246 // If the safety was adequate (for th << 247 // But this must not occur, no inters << 248 // so this branch, ie if( Intersects_ << 249 // << 250 G4bool validNormalLast; << 251 NormalAtEntry = GetSurfaceNormal( Po << 252 validNormalAtE = validNormalLast; << 253 231 254 // By moving point B, must take care 232 // By moving point B, must take care if current 255 // AF has no intersection to try curr 233 // AF has no intersection to try current FB!! 256 // 234 // 257 final_section = false; << 235 final_section= false; 258 236 259 #ifdef G4VERBOSE 237 #ifdef G4VERBOSE 260 if( fVerboseLevel > 3 ) 238 if( fVerboseLevel > 3 ) 261 { 239 { 262 G4cout << "G4PiF::LI> Investigating 240 G4cout << "G4PiF::LI> Investigating intermediate point" 263 << " at s=" << ApproxInterse 241 << " at s=" << ApproxIntersecPointV.GetCurveLength() 264 << " on way to full s=" 242 << " on way to full s=" 265 << CurveEndPointVelocity.Get 243 << CurveEndPointVelocity.GetCurveLength() << G4endl; 266 } 244 } 267 #endif 245 #endif 268 } 246 } 269 else // not Intersects_AF 247 else // not Intersects_AF 270 { 248 { 271 // In this case: 249 // In this case: 272 // There is NO intersection of AF wit 250 // There is NO intersection of AF with a volume boundary. 273 // We must continue the search in the 251 // We must continue the search in the segment FB! 274 // 252 // 275 GetNavigatorFor()->LocateGlobalPointW 253 GetNavigatorFor()->LocateGlobalPointWithinVolume( CurrentF_Point ); 276 254 277 G4double stepLengthFB; 255 G4double stepLengthFB; 278 G4ThreeVector PointH; 256 G4ThreeVector PointH; 279 G4bool usedNavigatorFB = false; << 280 257 281 // Check whether any volumes are enco 258 // Check whether any volumes are encountered by the chord FB 282 // ---------------------------------- 259 // --------------------------------------------------------- 283 260 284 G4bool Intersects_FB = IntersectChord 261 G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B, 285 262 NewSafety,fPreviousSafety, 286 263 fPreviousSftOrigin, 287 264 stepLengthFB, 288 << 265 PointH ); 289 if( Intersects_FB ) 266 if( Intersects_FB ) 290 { 267 { 291 // There is an intersection of FB w 268 // There is an intersection of FB with a volume boundary 292 // H <- First Intersection of Chord 269 // H <- First Intersection of Chord FB 293 270 294 // H is our new Candidate for the i 271 // H is our new Candidate for the intersection point. 295 // It replaces "E" and we will rep 272 // It replaces "E" and we will repeat the test to see if 296 // it is a good enough approximate 273 // it is a good enough approximate point for us. 297 274 298 // Note that F must be in volume vo 275 // Note that F must be in volume volA (the same as A) 299 // (otherwise AF would meet a volum 276 // (otherwise AF would meet a volume boundary!) 300 // A <- F 277 // A <- F 301 // E <- H 278 // E <- H 302 // 279 // 303 CurrentA_PointVelocity = ApproxInte 280 CurrentA_PointVelocity = ApproxIntersecPointV; 304 CurrentE_Point = PointH; 281 CurrentE_Point = PointH; 305 << 306 // Need to recalculate the Exit Nor << 307 // Relies on call to Navigator::Com << 308 // If safety was adequate (for the << 309 // But this must not occur, no inte << 310 // so this branch, i.e. if( Interse << 311 // << 312 G4bool validNormalLast; << 313 NormalAtEntry = GetSurfaceNormal( << 314 validNormalAtE = validNormalLast; << 315 } 282 } 316 else // not Intersects_FB 283 else // not Intersects_FB 317 { 284 { 318 // There is NO intersection of FB w 285 // There is NO intersection of FB with a volume boundary 319 286 320 if( final_section ) 287 if( final_section ) 321 { 288 { 322 // If B is the original endpoint, 289 // If B is the original endpoint, this means that whatever 323 // volume(s) intersected the orig 290 // volume(s) intersected the original chord, none touch the 324 // smaller chords we have used. 291 // smaller chords we have used. 325 // The value of 'IntersectedOrRec 292 // The value of 'IntersectedOrRecalculatedFT' returned is 326 // likely not valid 293 // likely not valid 327 294 328 there_is_no_intersection = true; 295 there_is_no_intersection = true; // real final_section 329 } 296 } 330 else 297 else 331 { 298 { 332 // We must restore the original e 299 // We must restore the original endpoint 333 300 334 CurrentA_PointVelocity = CurrentB 301 CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B 335 CurrentB_PointVelocity = CurveEnd 302 CurrentB_PointVelocity = CurveEndPointVelocity; 336 restoredFullEndpoint = true; 303 restoredFullEndpoint = true; 337 } 304 } 338 } // Endif (Intersects_FB) 305 } // Endif (Intersects_FB) 339 } // Endif (Intersects_AF) 306 } // Endif (Intersects_AF) 340 307 341 // Ensure that the new endpoints are no 308 // Ensure that the new endpoints are not further apart in space 342 // than on the curve due to different e 309 // than on the curve due to different errors in the integration 343 // 310 // 344 G4double linDistSq, curveDist; 311 G4double linDistSq, curveDist; 345 linDistSq = ( CurrentB_PointVelocity.Ge 312 linDistSq = ( CurrentB_PointVelocity.GetPosition() 346 - CurrentA_PointVelocity.Ge 313 - CurrentA_PointVelocity.GetPosition() ).mag2(); 347 curveDist = CurrentB_PointVelocity.GetC 314 curveDist = CurrentB_PointVelocity.GetCurveLength() 348 - CurrentA_PointVelocity.Ge 315 - CurrentA_PointVelocity.GetCurveLength(); 349 316 350 // Change this condition for very stric 317 // Change this condition for very strict parameters of propagation 351 // 318 // 352 if( curveDist*curveDist*(1+2* GetEpsilo 319 if( curveDist*curveDist*(1+2* GetEpsilonStepFor()) < linDistSq ) 353 { 320 { 354 // Re-integrate to obtain a new B 321 // Re-integrate to obtain a new B 355 // 322 // 356 G4FieldTrack newEndPointFT = 323 G4FieldTrack newEndPointFT = 357 ReEstimateEndpoint( CurrentA_ 324 ReEstimateEndpoint( CurrentA_PointVelocity, 358 CurrentB_ 325 CurrentB_PointVelocity, 359 linDistSq 326 linDistSq, // to avoid recalculation 360 curveDist 327 curveDist ); 361 G4FieldTrack oldPointVelB = CurrentB_ 328 G4FieldTrack oldPointVelB = CurrentB_PointVelocity; 362 CurrentB_PointVelocity = newEndPointF 329 CurrentB_PointVelocity = newEndPointFT; 363 330 364 if( (final_section)) // real final se 331 if( (final_section)) // real final section 365 { 332 { 366 recalculatedEndPoint = true; 333 recalculatedEndPoint = true; 367 IntersectedOrRecalculatedFT = newEn 334 IntersectedOrRecalculatedFT = newEndPointFT; 368 // So that we can return it, if i 335 // So that we can return it, if it is the endpoint! 369 } 336 } 370 } 337 } 371 if( curveDist < 0.0 ) 338 if( curveDist < 0.0 ) 372 { 339 { >> 340 G4cerr << "ERROR - G4SimpleLocator::EstimateIntersectionPoint()" >> 341 << G4endl >> 342 << " Error in advancing propagation." << G4endl; 373 fVerboseLevel = 5; // Print out a max 343 fVerboseLevel = 5; // Print out a maximum of information 374 printStatus( CurrentA_PointVelocity, 344 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity, 375 -1.0, NewSafety, subste 345 -1.0, NewSafety, substep_no ); 376 std::ostringstream message; << 346 G4cerr << " Point A (start) is " << CurrentA_PointVelocity 377 message << "Error in advancing propag << 347 << G4endl; 378 << " Point A (start) i << 348 G4cerr << " Point B (end) is " << CurrentB_PointVelocity 379 << G4endl << 349 << G4endl; 380 << " Point B (end) i << 350 G4cerr << " Curve distance is " << curveDist << G4endl; 381 << G4endl << 351 G4cerr << G4endl 382 << " Curve distance is << 352 << "The final curve point is not further along" 383 << G4endl << 353 << " than the original!" << G4endl; 384 << "The final curve point is << 385 << " than the original!" << G << 386 354 387 if( recalculatedEndPoint ) 355 if( recalculatedEndPoint ) 388 { 356 { 389 message << "Recalculation of EndPoi << 357 G4cerr << "Recalculation of EndPoint was called with fEpsStep= " 390 << GetEpsilonStepFor() << G << 358 << GetEpsilonStepFor() << G4endl; 391 } 359 } 392 message.precision(20); << 360 G4cerr.precision(20); 393 message << " Point A (Curve start) << 361 G4cerr << " Point A (Curve start) is " << CurveStartPointVelocity 394 << G4endl << 362 << G4endl; 395 << " Point B (Curve end) << 363 G4cerr << " Point B (Curve end) is " << CurveEndPointVelocity 396 << G4endl << 364 << G4endl; 397 << " Point A (Current start) << 365 G4cerr << " Point A (Current start) is " << CurrentA_PointVelocity 398 << G4endl << 366 << G4endl; 399 << " Point B (Current end) << 367 G4cerr << " Point B (Current end) is " << CurrentB_PointVelocity 400 << G4endl << 368 << G4endl; 401 << " Point E (Trial Point) << 369 G4cerr << " Point E (Trial Point) is " << CurrentE_Point 402 << G4endl << 370 << G4endl; 403 << " Point F (Intersection) << 371 G4cerr << " Point F (Intersection) is " << ApproxIntersecPointV 404 << G4endl << 372 << G4endl; 405 << " LocateIntersectio << 373 G4cerr << " LocateIntersection parameters are : Substep no= " 406 << substep_no; << 374 << substep_no << G4endl; 407 375 408 G4Exception("G4SimpleLocator::Estimat 376 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()", 409 "GeomNav0003", FatalExcep << 377 "FatalError", FatalException, >> 378 "Error in advancing propagation."); 410 } 379 } 411 380 412 if ( restoredFullEndpoint ) << 381 if(restoredFullEndpoint) 413 { 382 { 414 final_section = restoredFullEndpoint; 383 final_section = restoredFullEndpoint; 415 restoredFullEndpoint = false; 384 restoredFullEndpoint = false; 416 } 385 } 417 } // EndIf ( E is close enough to the cur 386 } // EndIf ( E is close enough to the curve, ie point F. ) 418 // tests ChordAF_Vector.mag() <= maximu 387 // tests ChordAF_Vector.mag() <= maximum_lateral_displacement 419 388 420 #ifdef G4DEBUG_LOCATE_INTERSECTION 389 #ifdef G4DEBUG_LOCATE_INTERSECTION 421 G4int trigger_substepno_print= warn_subst << 422 << 423 if( substep_no >= trigger_substepno_print 390 if( substep_no >= trigger_substepno_print ) 424 { 391 { 425 G4cout << "Difficulty in converging in 392 G4cout << "Difficulty in converging in " 426 << "G4SimpleLocator::EstimateInt 393 << "G4SimpleLocator::EstimateIntersectionPoint():" 427 << G4endl 394 << G4endl 428 << " Substep no = " << subste 395 << " Substep no = " << substep_no << G4endl; 429 if( substep_no == trigger_substepno_pri 396 if( substep_no == trigger_substepno_print ) 430 { 397 { 431 printStatus( CurveStartPointVelocity, 398 printStatus( CurveStartPointVelocity, CurveEndPointVelocity, 432 -1.0, NewSafety, 0); 399 -1.0, NewSafety, 0); 433 } 400 } 434 G4cout << " State of point A: "; 401 G4cout << " State of point A: "; 435 printStatus( CurrentA_PointVelocity, Cu 402 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity, 436 -1.0, NewSafety, substep_n 403 -1.0, NewSafety, substep_no-1, 0); 437 G4cout << " State of point B: "; 404 G4cout << " State of point B: "; 438 printStatus( CurrentA_PointVelocity, Cu 405 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity, 439 -1.0, NewSafety, substep_n 406 -1.0, NewSafety, substep_no); 440 } 407 } 441 #endif 408 #endif 442 ++substep_no; << 409 substep_no++; 443 410 444 } while ( ( ! found_approximate_intersection 411 } while ( ( ! found_approximate_intersection ) 445 && ( ! there_is_no_intersection ) 412 && ( ! there_is_no_intersection ) 446 && ( substep_no <= max_substeps) ); 413 && ( substep_no <= max_substeps) ); // UNTIL found or failed 447 414 448 if( substep_no > max_no_seen ) 415 if( substep_no > max_no_seen ) 449 { 416 { 450 max_no_seen = substep_no; 417 max_no_seen = substep_no; 451 #ifdef G4DEBUG_LOCATE_INTERSECTION << 452 if( max_no_seen > warn_substeps ) 418 if( max_no_seen > warn_substeps ) 453 { 419 { 454 trigger_substepno_print = max_no_seen-20 420 trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps 455 } 421 } 456 #endif << 457 } 422 } 458 423 459 if( ( substep_no >= max_substeps) 424 if( ( substep_no >= max_substeps) 460 && !there_is_no_intersection 425 && !there_is_no_intersection 461 && !found_approximate_intersection ) 426 && !found_approximate_intersection ) 462 { 427 { 463 G4cout << "ERROR - G4SimpleLocator::Estima << 428 G4cerr << "WARNING - G4SimpleLocator::EstimateIntersectionPoint()" 464 << " Start and Endpoint of R << 429 << G4endl >> 430 << " Convergence is requiring too many substeps: " >> 431 << substep_no << G4endl; >> 432 G4cerr << " Abandoning effort to intersect. " << G4endl; >> 433 G4cerr << " Information on start & current step follows in cout." >> 434 << G4endl; >> 435 G4cout << "WARNING - G4SimpleLocator::EstimateIntersectionPoint()" >> 436 << G4endl >> 437 << " Convergence is requiring too many substeps: " >> 438 << substep_no << G4endl; >> 439 G4cout << " Found intersection = " >> 440 << found_approximate_intersection << G4endl >> 441 << " Intersection exists = " >> 442 << !there_is_no_intersection << G4endl; >> 443 G4cout << " Start and Endpoint of Requested Step:" << G4endl; 465 printStatus( CurveStartPointVelocity, Curv 444 printStatus( CurveStartPointVelocity, CurveEndPointVelocity, 466 -1.0, NewSafety, 0); 445 -1.0, NewSafety, 0); 467 G4cout << G4endl << 446 G4cout << G4endl; 468 << " Start and end-point of << 447 G4cout << " 'Bracketing' starting and endpoint of current Sub-Step" >> 448 << G4endl; 469 printStatus( CurrentA_PointVelocity, Curre 449 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity, 470 -1.0, NewSafety, substep_no-1 450 -1.0, NewSafety, substep_no-1); 471 printStatus( CurrentA_PointVelocity, Curre 451 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity, 472 -1.0, NewSafety, substep_no); 452 -1.0, NewSafety, substep_no); 473 << 453 G4cout << G4endl; 474 std::ostringstream message; << 454 G4cout.precision( 10 ); 475 message << "Convergence is requiring too m << 476 << substep_no << G4endl << 477 << " Abandoning effort to << 478 << " Found intersection = << 479 << found_approximate_intersection << 480 << " Intersection exists << 481 << !there_is_no_intersection << G4 << 482 message.precision(10); << 483 G4double done_len = CurrentA_PointVelocity 455 G4double done_len = CurrentA_PointVelocity.GetCurveLength(); 484 G4double full_len = CurveEndPointVelocity. 456 G4double full_len = CurveEndPointVelocity.GetCurveLength(); 485 message << " Undertaken only leng << 457 G4cout << "ERROR - G4SimpleLocator::EstimateIntersectionPoint()" 486 << " out of " << full_len << " req << 458 << G4endl 487 << " Remaining length = " << 459 << " Undertaken only length: " << done_len >> 460 << " out of " << full_len << " required." << G4endl; >> 461 G4cout << " Remaining length = " << full_len - done_len << G4endl; 488 462 489 G4Exception("G4SimpleLocator::EstimateInte 463 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()", 490 "GeomNav0003", FatalException, << 464 "UnableToLocateIntersection", FatalException, >> 465 "Too many substeps while trying to locate intersection."); 491 } 466 } 492 else if( substep_no >= warn_substeps ) 467 else if( substep_no >= warn_substeps ) 493 { 468 { 494 std::ostringstream message; << 469 G4int oldprc= G4cout.precision( 10 ); 495 message.precision(10); << 470 G4cout << "WARNING - G4SimpleLocator::EstimateIntersectionPoint()" 496 << 471 << G4endl 497 message << "Many substeps while trying to << 472 << " Undertaken length: " 498 << " Undertaken length: " << 473 << CurrentB_PointVelocity.GetCurveLength(); 499 << CurrentB_PointVelocity.GetCurve << 474 G4cout << " - Needed: " << substep_no << " substeps." << G4endl 500 << " - Needed: " << substep_no << << 475 << " Warning level = " << warn_substeps 501 << " Warning level = " << << 476 << " and maximum substeps = " << max_substeps << G4endl; 502 << " and maximum substeps = " << m << 503 G4Exception("G4SimpleLocator::EstimateInte 477 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()", 504 "GeomNav1002", JustWarning, me << 478 "DifficultyToLocateIntersection", JustWarning, >> 479 "Many substeps while trying to locate intersection."); >> 480 G4cout.precision( oldprc ); 505 } 481 } 506 return !there_is_no_intersection; // Succe 482 return !there_is_no_intersection; // Success or failure 507 } 483 } 508 484