Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // Class G4SimpleLocator implementation 27 // 28 // 27.10.08 - Tatiana Nikitina, extracted from 29 // 04.10.11 - John Apostolakis, revised conver 30 // ------------------------------------------- 31 32 #include <iomanip> 33 34 #include "G4ios.hh" 35 #include "G4SimpleLocator.hh" 36 37 G4SimpleLocator::G4SimpleLocator(G4Navigator * 38 : G4VIntersectionLocator(theNavigator) 39 { 40 } 41 42 G4SimpleLocator::~G4SimpleLocator() = default; 43 44 // ------------------------------------------- 45 // G4bool G4PropagatorInField::LocateIntersect 46 // const G4FieldTrack& CurveStart 47 // const G4FieldTrack& CurveEndPo 48 // const G4ThreeVector& TrialPoint 49 // G4FieldTrack& Intersecte 50 // G4bool& recalculat 51 // ------------------------------------------- 52 // 53 // Function that returns the intersection of t 54 // of the current volume (either the external 55 // of the daughters: 56 // 57 // A = Initial point 58 // B = another point 59 // 60 // Both A and B are assumed to be on the true 61 // 62 // E is the first point of intersection of 63 // a volume other than A (on the surface 64 // 65 // Convention of Use : 66 // i) If it returns "true", then Intersect 67 // to the approximate intersection point 68 // ii) If it returns "false", no intersecti 69 // The validity of IntersectedOrRecal 70 // a) if latter is false, then Intersec 71 // b) if latter is true, then Intersec 72 // the new endpoint, due to a re-i 73 // ------------------------------------------- 74 // NOTE: implementation taken from G4Propagato 75 // 76 G4bool G4SimpleLocator::EstimateIntersectionPo 77 const G4FieldTrack& CurveStart 78 const G4FieldTrack& CurveEndPo 79 const G4ThreeVector& TrialPoint 80 G4FieldTrack& Intersecte 81 G4bool& recalculat 82 G4double &fPrevious 83 G4ThreeVector &fPrevious 84 { 85 // Find Intersection Point ( A, B, E ) of t 86 87 G4bool found_approximate_intersection = fals 88 G4bool there_is_no_intersection = fals 89 90 G4FieldTrack CurrentA_PointVelocity = Curve 91 G4FieldTrack CurrentB_PointVelocity = Curve 92 G4ThreeVector CurrentE_Point = TrialPoint; 93 G4bool validNormalAtE = false; 94 G4ThreeVector NormalAtEntry; 95 96 G4FieldTrack ApproxIntersecPointV(CurveEndP 97 G4double NewSafety = 0.0; 98 G4bool last_AF_intersection = false; 99 G4bool final_section = true; // Shows wheth 100 // (i.e. B=ful 101 recalculatedEndPoint = false; 102 103 G4bool restoredFullEndpoint = false; 104 105 G4int substep_no = 0; 106 107 // Limits for substep number 108 // 109 const G4int max_substeps = 100000000; // T 110 const G4int warn_substeps = 1000; // 111 112 // Statistics for substeps 113 // 114 static G4ThreadLocal G4int max_no_seen= -1; 115 116 NormalAtEntry = GetSurfaceNormal( CurrentE_P 117 118 #ifdef G4DEBUG_FIELD 119 const G4double tolerance = 1.0e-8; 120 G4ThreeVector StartPosition= CurveStartPoin 121 if( (TrialPoint - StartPosition).mag() < tol 122 { 123 G4Exception("G4SimpleLocator::EstimateInt 124 "GeomNav1002", JustWarning, 125 "Intersection point F is exac 126 } 127 #endif 128 129 do 130 { 131 G4ThreeVector Point_A = CurrentA_PointVel 132 G4ThreeVector Point_B = CurrentB_PointVel 133 134 // F = a point on true AB path close to p 135 // (the closest if possible) 136 // 137 ApproxIntersecPointV = GetChordFinderFor( 138 ->ApproxCurvePoint 139 140 141 142 // The above method is the key & most 143 144 #ifdef G4DEBUG_FIELD 145 if( ApproxIntersecPointV.GetCurveLength() 146 CurrentB_PointVelocity.GetCurveLength 147 { 148 G4Exception("G4SimpleLocator::EstimateI 149 "GeomNav0003", FatalExcepti 150 "Intermediate F point is pa 151 } 152 #endif 153 154 G4ThreeVector CurrentF_Point = ApproxInte 155 156 // First check whether EF is small - then 157 // Calculate the length and direction of 158 // 159 G4ThreeVector ChordEF_Vector = CurrentF_P 160 161 G4ThreeVector NewMomentumDir = ApproxInte 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 { 180 found_approximate_intersection = true; 181 182 // Create the "point" return value 183 // 184 IntersectedOrRecalculatedFT = ApproxInt 185 IntersectedOrRecalculatedFT.SetPosition 186 187 if ( GetAdjustementOfFoundIntersection( 188 { 189 // Try to Get Correction of Intersect 190 // 191 G4ThreeVector IP; 192 G4ThreeVector MomentumDir= ApproxInte 193 G4bool goodCorrection = AdjustmentOfF 194 CurrentE_Po 195 last_AF_int 196 fPreviousSa 197 198 if ( goodCorrection ) 199 { 200 IntersectedOrRecalculatedFT = Appro 201 IntersectedOrRecalculatedFT.SetPosi 202 } 203 } 204 205 // Note: in order to return a point on 206 // we must return E. But it is F 207 // So we must "cheat": we are usi 208 // and the velocity at point F !! 209 // 210 // This must limit the length we can al 211 } 212 else // E is NOT close enough to the cur 213 { 214 // Check whether any volumes are encoun 215 // ------------------------------------ 216 // First relocate to restore any Voxel 217 // in the Navigator before calling Comp 218 // 219 GetNavigatorFor()->LocateGlobalPointWit 220 221 G4ThreeVector PointG; // Candidate in 222 G4double stepLengthAF; 223 G4bool usedNavigatorAF = false; 224 G4bool Intersects_AF = IntersectChord( 225 226 227 228 229 230 231 232 last_AF_intersection = Intersects_AF; 233 if( Intersects_AF ) 234 { 235 // G is our new Candidate for the int 236 // It replaces "E" and we will repea 237 // it is a good enough approximate po 238 // B <- F 239 // E <- G 240 241 CurrentB_PointVelocity = ApproxInters 242 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 254 // By moving point B, must take care 255 // AF has no intersection to try curr 256 // 257 final_section = false; 258 259 #ifdef G4VERBOSE 260 if( fVerboseLevel > 3 ) 261 { 262 G4cout << "G4PiF::LI> Investigating 263 << " at s=" << ApproxInterse 264 << " on way to full s=" 265 << CurveEndPointVelocity.Get 266 } 267 #endif 268 } 269 else // not Intersects_AF 270 { 271 // In this case: 272 // There is NO intersection of AF wit 273 // We must continue the search in the 274 // 275 GetNavigatorFor()->LocateGlobalPointW 276 277 G4double stepLengthFB; 278 G4ThreeVector PointH; 279 G4bool usedNavigatorFB = false; 280 281 // Check whether any volumes are enco 282 // ---------------------------------- 283 284 G4bool Intersects_FB = IntersectChord 285 286 287 288 289 if( Intersects_FB ) 290 { 291 // There is an intersection of FB w 292 // H <- First Intersection of Chord 293 294 // H is our new Candidate for the i 295 // It replaces "E" and we will rep 296 // it is a good enough approximate 297 298 // Note that F must be in volume vo 299 // (otherwise AF would meet a volum 300 // A <- F 301 // E <- H 302 // 303 CurrentA_PointVelocity = ApproxInte 304 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 } 316 else // not Intersects_FB 317 { 318 // There is NO intersection of FB w 319 320 if( final_section ) 321 { 322 // If B is the original endpoint, 323 // volume(s) intersected the orig 324 // smaller chords we have used. 325 // The value of 'IntersectedOrRec 326 // likely not valid 327 328 there_is_no_intersection = true; 329 } 330 else 331 { 332 // We must restore the original e 333 334 CurrentA_PointVelocity = CurrentB 335 CurrentB_PointVelocity = CurveEnd 336 restoredFullEndpoint = true; 337 } 338 } // Endif (Intersects_FB) 339 } // Endif (Intersects_AF) 340 341 // Ensure that the new endpoints are no 342 // than on the curve due to different e 343 // 344 G4double linDistSq, curveDist; 345 linDistSq = ( CurrentB_PointVelocity.Ge 346 - CurrentA_PointVelocity.Ge 347 curveDist = CurrentB_PointVelocity.GetC 348 - CurrentA_PointVelocity.Ge 349 350 // Change this condition for very stric 351 // 352 if( curveDist*curveDist*(1+2* GetEpsilo 353 { 354 // Re-integrate to obtain a new B 355 // 356 G4FieldTrack newEndPointFT = 357 ReEstimateEndpoint( CurrentA_ 358 CurrentB_ 359 linDistSq 360 curveDist 361 G4FieldTrack oldPointVelB = CurrentB_ 362 CurrentB_PointVelocity = newEndPointF 363 364 if( (final_section)) // real final se 365 { 366 recalculatedEndPoint = true; 367 IntersectedOrRecalculatedFT = newEn 368 // So that we can return it, if i 369 } 370 } 371 if( curveDist < 0.0 ) 372 { 373 fVerboseLevel = 5; // Print out a max 374 printStatus( CurrentA_PointVelocity, 375 -1.0, NewSafety, subste 376 std::ostringstream message; 377 message << "Error in advancing propag 378 << " Point A (start) i 379 << G4endl 380 << " Point B (end) i 381 << G4endl 382 << " Curve distance is 383 << G4endl 384 << "The final curve point is 385 << " than the original!" << G 386 387 if( recalculatedEndPoint ) 388 { 389 message << "Recalculation of EndPoi 390 << GetEpsilonStepFor() << G 391 } 392 message.precision(20); 393 message << " Point A (Curve start) 394 << G4endl 395 << " Point B (Curve end) 396 << G4endl 397 << " Point A (Current start) 398 << G4endl 399 << " Point B (Current end) 400 << G4endl 401 << " Point E (Trial Point) 402 << G4endl 403 << " Point F (Intersection) 404 << G4endl 405 << " LocateIntersectio 406 << substep_no; 407 408 G4Exception("G4SimpleLocator::Estimat 409 "GeomNav0003", FatalExcep 410 } 411 412 if ( restoredFullEndpoint ) 413 { 414 final_section = restoredFullEndpoint; 415 restoredFullEndpoint = false; 416 } 417 } // EndIf ( E is close enough to the cur 418 // tests ChordAF_Vector.mag() <= maximu 419 420 #ifdef G4DEBUG_LOCATE_INTERSECTION 421 G4int trigger_substepno_print= warn_subst 422 423 if( substep_no >= trigger_substepno_print 424 { 425 G4cout << "Difficulty in converging in 426 << "G4SimpleLocator::EstimateInt 427 << G4endl 428 << " Substep no = " << subste 429 if( substep_no == trigger_substepno_pri 430 { 431 printStatus( CurveStartPointVelocity, 432 -1.0, NewSafety, 0); 433 } 434 G4cout << " State of point A: "; 435 printStatus( CurrentA_PointVelocity, Cu 436 -1.0, NewSafety, substep_n 437 G4cout << " State of point B: "; 438 printStatus( CurrentA_PointVelocity, Cu 439 -1.0, NewSafety, substep_n 440 } 441 #endif 442 ++substep_no; 443 444 } while ( ( ! found_approximate_intersection 445 && ( ! there_is_no_intersection ) 446 && ( substep_no <= max_substeps) ); 447 448 if( substep_no > max_no_seen ) 449 { 450 max_no_seen = substep_no; 451 #ifdef G4DEBUG_LOCATE_INTERSECTION 452 if( max_no_seen > warn_substeps ) 453 { 454 trigger_substepno_print = max_no_seen-20 455 } 456 #endif 457 } 458 459 if( ( substep_no >= max_substeps) 460 && !there_is_no_intersection 461 && !found_approximate_intersection ) 462 { 463 G4cout << "ERROR - G4SimpleLocator::Estima 464 << " Start and Endpoint of R 465 printStatus( CurveStartPointVelocity, Curv 466 -1.0, NewSafety, 0); 467 G4cout << G4endl 468 << " Start and end-point of 469 printStatus( CurrentA_PointVelocity, Curre 470 -1.0, NewSafety, substep_no-1 471 printStatus( CurrentA_PointVelocity, Curre 472 -1.0, NewSafety, substep_no); 473 474 std::ostringstream message; 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 484 G4double full_len = CurveEndPointVelocity. 485 message << " Undertaken only leng 486 << " out of " << full_len << " req 487 << " Remaining length = " 488 489 G4Exception("G4SimpleLocator::EstimateInte 490 "GeomNav0003", FatalException, 491 } 492 else if( substep_no >= warn_substeps ) 493 { 494 std::ostringstream message; 495 message.precision(10); 496 497 message << "Many substeps while trying to 498 << " Undertaken length: " 499 << CurrentB_PointVelocity.GetCurve 500 << " - Needed: " << substep_no << 501 << " Warning level = " << 502 << " and maximum substeps = " << m 503 G4Exception("G4SimpleLocator::EstimateInte 504 "GeomNav1002", JustWarning, me 505 } 506 return !there_is_no_intersection; // Succe 507 } 508