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