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 G4MultiLevelLocator implementation 26 // Class G4MultiLevelLocator implementation 27 // 27 // 28 // 27.10.08 - Tatiana Nikitina. 28 // 27.10.08 - Tatiana Nikitina. 29 // 04.10.11 - John Apostolakis, revised conver 29 // 04.10.11 - John Apostolakis, revised convergence to use Surface Normal 30 // ------------------------------------------- 30 // --------------------------------------------------------------------------- 31 31 32 #include <iomanip> 32 #include <iomanip> 33 33 34 #include "G4ios.hh" 34 #include "G4ios.hh" 35 #include "G4MultiLevelLocator.hh" 35 #include "G4MultiLevelLocator.hh" 36 #include "G4LocatorChangeRecord.hh" 36 #include "G4LocatorChangeRecord.hh" 37 #include "G4LocatorChangeLogger.hh" 37 #include "G4LocatorChangeLogger.hh" 38 38 39 G4MultiLevelLocator::G4MultiLevelLocator(G4Nav 39 G4MultiLevelLocator::G4MultiLevelLocator(G4Navigator *theNavigator) 40 : G4VIntersectionLocator(theNavigator) 40 : G4VIntersectionLocator(theNavigator) 41 { 41 { 42 // In case of too slow progress in finding I 42 // In case of too slow progress in finding Intersection Point 43 // intermediates Points on the Track must be 43 // intermediates Points on the Track must be stored. 44 // Initialise the array of Pointers [max_dep 44 // Initialise the array of Pointers [max_depth+1] to do this 45 45 46 G4ThreeVector zeroV(0.0,0.0,0.0); 46 G4ThreeVector zeroV(0.0,0.0,0.0); 47 for (auto & idepth : ptrInterMedFT) << 47 for ( auto idepth=0; idepth<max_depth+1; ++idepth ) 48 { 48 { 49 idepth = new G4FieldTrack( zeroV, zeroV, 0 << 49 ptrInterMedFT[ idepth ] = new G4FieldTrack( zeroV, zeroV, 0., 0., 0., 0.); 50 } 50 } 51 51 52 if (fCheckMode) 52 if (fCheckMode) 53 { 53 { 54 // Trial values Loose Mediu 54 // Trial values Loose Medium Tight 55 // To happen: Infrequent Occas 55 // To happen: Infrequent Occasional Often 56 SetMaxSteps(150); // 300 85 56 SetMaxSteps(150); // 300 85 25 57 SetWarnSteps(80); // 250 60 57 SetWarnSteps(80); // 250 60 15 58 } 58 } 59 } 59 } 60 60 61 G4MultiLevelLocator::~G4MultiLevelLocator() 61 G4MultiLevelLocator::~G4MultiLevelLocator() 62 { 62 { 63 for (auto & idepth : ptrInterMedFT) << 63 for ( auto idepth=0; idepth<max_depth+1; ++idepth ) 64 { 64 { 65 delete idepth; << 65 delete ptrInterMedFT[idepth]; 66 } 66 } 67 #ifdef G4DEBUG_FIELD 67 #ifdef G4DEBUG_FIELD 68 ReportStatistics(); 68 ReportStatistics(); 69 #endif 69 #endif 70 } 70 } 71 71 72 72 73 // ------------------------------------------- 73 // -------------------------------------------------------------------------- 74 // G4bool G4PropagatorInField::LocateIntersect 74 // G4bool G4PropagatorInField::LocateIntersectionPoint( 75 // const G4FieldTrack& CurveStart 75 // const G4FieldTrack& CurveStartPointVelocity, // A 76 // const G4FieldTrack& CurveEndPo 76 // const G4FieldTrack& CurveEndPointVelocity, // B 77 // const G4ThreeVector& TrialPoint 77 // const G4ThreeVector& TrialPoint, // E 78 // G4FieldTrack& Intersecte 78 // G4FieldTrack& IntersectedOrRecalculated // Output 79 // G4bool& recalculat 79 // G4bool& recalculated ) // Out 80 // ------------------------------------------- 80 // -------------------------------------------------------------------------- 81 // 81 // 82 // Function that returns the intersection of t 82 // Function that returns the intersection of the true path with the surface 83 // of the current volume (either the external 83 // of the current volume (either the external one or the inner one with one 84 // of the daughters: 84 // of the daughters: 85 // 85 // 86 // A = Initial point 86 // A = Initial point 87 // B = another point 87 // B = another point 88 // 88 // 89 // Both A and B are assumed to be on the true 89 // Both A and B are assumed to be on the true path: 90 // 90 // 91 // E is the first point of intersection of 91 // E is the first point of intersection of the chord AB with 92 // a volume other than A (on the surface 92 // a volume other than A (on the surface of A or of a daughter) 93 // 93 // 94 // Convention of Use : 94 // Convention of Use : 95 // i) If it returns "true", then Intersect 95 // i) If it returns "true", then IntersectionPointVelocity is set 96 // to the approximate intersection point 96 // to the approximate intersection point. 97 // ii) If it returns "false", no intersecti 97 // ii) If it returns "false", no intersection was found. 98 // Potential reasons: 98 // Potential reasons: 99 // a) no segment found an intersection 99 // a) no segment found an intersection 100 // b) too many steps were required - af 100 // b) too many steps were required - after that it abandoned the effort 101 // and is returning how far it could 101 // and is returning how far it could go. (New - 29 Oct 2015) 102 // (If so, it must set 'recalculated 102 // (If so, it must set 'recalculated' to true.) 103 // TODO/idea: add a new flag: 'unfinish 103 // TODO/idea: add a new flag: 'unfinished' to identify these cases. 104 // 104 // 105 // IntersectedOrRecalculated means diff 105 // IntersectedOrRecalculated means different things: 106 // a) if it is the same curve lenght al 106 // a) if it is the same curve lenght along, it is a revision of the 107 // original enpdoint due to the nee 107 // original enpdoint due to the need for re-integration. 108 // b) if it is at a shorter curve lengt 108 // b) if it is at a shorter curve length, it is 'end of what it could do' 109 // i.e. as far as it could go, becau 109 // i.e. as far as it could go, because it took too many steps! 110 // Note: IntersectedOrRecalculated is v 110 // Note: IntersectedOrRecalculated is valid only if 'recalculated' is 111 // 'true'. 111 // 'true'. 112 // ------------------------------------------- 112 // -------------------------------------------------------------------------- 113 // NOTE: implementation taken from G4Propagato 113 // NOTE: implementation taken from G4PropagatorInField 114 // 114 // 115 G4bool G4MultiLevelLocator::EstimateIntersecti 115 G4bool G4MultiLevelLocator::EstimateIntersectionPoint( 116 const G4FieldTrack& CurveStart 116 const G4FieldTrack& CurveStartPointVelocity, // A 117 const G4FieldTrack& CurveEndPo 117 const G4FieldTrack& CurveEndPointVelocity, // B 118 const G4ThreeVector& TrialPoint 118 const G4ThreeVector& TrialPoint, // E 119 G4FieldTrack& Intersecte 119 G4FieldTrack& IntersectedOrRecalculatedFT, // Output 120 G4bool& recalculat 120 G4bool& recalculatedEndPoint, // Out 121 G4double& previousSa 121 G4double& previousSafety, // In/Out 122 G4ThreeVector& previousSf 122 G4ThreeVector& previousSftOrigin) // In/Out 123 { 123 { 124 // Find Intersection Point ( A, B, E ) of t 124 // Find Intersection Point ( A, B, E ) of true path AB - start at E. 125 const char* MethodName= "G4MultiLevelLocator 125 const char* MethodName= "G4MultiLevelLocator::EstimateIntersectionPoint()"; 126 126 127 G4bool found_approximate_intersection = fals 127 G4bool found_approximate_intersection = false; 128 G4bool there_is_no_intersection = fals 128 G4bool there_is_no_intersection = false; 129 129 130 G4FieldTrack CurrentA_PointVelocity = Curve 130 G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity; 131 G4FieldTrack CurrentB_PointVelocity = Curve 131 G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity; 132 G4ThreeVector CurrentE_Point = TrialPoint; 132 G4ThreeVector CurrentE_Point = TrialPoint; 133 G4bool validNormalAtE = false; 133 G4bool validNormalAtE = false; 134 G4ThreeVector NormalAtEntry; 134 G4ThreeVector NormalAtEntry; 135 135 136 G4FieldTrack ApproxIntersecPointV(CurveEndP 136 G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity); // FT-Def-Construct 137 G4bool validIntersectP= true; // Is 137 G4bool validIntersectP= true; // Is it current ? 138 G4double NewSafety = 0.0; 138 G4double NewSafety = 0.0; 139 G4bool last_AF_intersection = false; 139 G4bool last_AF_intersection = false; 140 140 141 auto integrDriver = GetChordFinderFor()-> 141 auto integrDriver = GetChordFinderFor()->GetIntegrationDriver(); 142 G4bool driverReIntegrates = integrDriver->D 142 G4bool driverReIntegrates = integrDriver->DoesReIntegrate(); 143 143 144 G4bool first_section = true; 144 G4bool first_section = true; 145 recalculatedEndPoint = false; 145 recalculatedEndPoint = false; 146 G4bool restoredFullEndpoint = false; 146 G4bool restoredFullEndpoint = false; 147 147 148 unsigned int substep_no = 0; 148 unsigned int substep_no = 0; 149 149 150 // Statistics for substeps 150 // Statistics for substeps 151 static G4ThreadLocal unsigned int max_no_see 151 static G4ThreadLocal unsigned int max_no_seen= 0; 152 152 153 #ifdef G4DEBUG_FIELD 153 #ifdef G4DEBUG_FIELD 154 unsigned int trigger_substepno_print = 0; 154 unsigned int trigger_substepno_print = 0; 155 const G4double tolerance = 1.0e-8 * CLHEP::m 155 const G4double tolerance = 1.0e-8 * CLHEP::mm; 156 unsigned int biggest_depth = 0; 156 unsigned int biggest_depth = 0; 157 // using kInitialisingCL = G4LocatorChange 157 // using kInitialisingCL = G4LocatorChangeRecord::kInitialisingCL; 158 #endif 158 #endif 159 159 160 // Log the location, iteration of changes in 160 // Log the location, iteration of changes in A,B 161 //------------------------------------------ 161 //---------------------------------------------- 162 static G4ThreadLocal G4LocatorChangeLogger e 162 static G4ThreadLocal G4LocatorChangeLogger endChangeA("StartPointA"), 163 endChangeB("EndPointB"), recApproxPoint("A 163 endChangeB("EndPointB"), recApproxPoint("ApproxPoint"), 164 pointH_logger("Trial points 'E': position, 164 pointH_logger("Trial points 'E': position, normal"); 165 unsigned int eventCount = 0; 165 unsigned int eventCount = 0; 166 166 167 if (fCheckMode) 167 if (fCheckMode) 168 { 168 { 169 // Clear previous call's data 169 // Clear previous call's data 170 endChangeA.clear(); 170 endChangeA.clear(); 171 endChangeB.clear(); 171 endChangeB.clear(); 172 recApproxPoint.clear(); 172 recApproxPoint.clear(); 173 pointH_logger.clear(); 173 pointH_logger.clear(); 174 174 175 // Record the initialisation 175 // Record the initialisation 176 ++eventCount; 176 ++eventCount; 177 endChangeA.AddRecord( G4LocatorChangeRecor 177 endChangeA.AddRecord( G4LocatorChangeRecord::kInitialisingCL, substep_no, 178 eventCount, CurrentA 178 eventCount, CurrentA_PointVelocity ); 179 endChangeB.AddRecord( G4LocatorChangeRecor 179 endChangeB.AddRecord( G4LocatorChangeRecord::kInitialisingCL, substep_no, 180 eventCount, CurrentB 180 eventCount, CurrentB_PointVelocity ); 181 } 181 } 182 182 183 //------------------------------------------ 183 //-------------------------------------------------------------------------- 184 // Algorithm for the case if progress in fo 184 // Algorithm for the case if progress in founding intersection is too slow. 185 // Process is defined too slow if after N=p 185 // Process is defined too slow if after N=param_substeps advances on the 186 // path, it will be only 'fraction_done' of 186 // path, it will be only 'fraction_done' of the total length. 187 // In this case the remaining length is div 187 // In this case the remaining length is divided in two half and 188 // the loop is restarted for each half. 188 // the loop is restarted for each half. 189 // If progress is still too slow, the divis 189 // If progress is still too slow, the division in two halfs continue 190 // until 'max_depth'. 190 // until 'max_depth'. 191 //------------------------------------------ 191 //-------------------------------------------------------------------------- 192 192 193 const G4int param_substeps = 5; // Test val 193 const G4int param_substeps = 5; // Test value for the maximum number 194 // of subst 194 // of substeps 195 const G4double fraction_done = 0.3; 195 const G4double fraction_done = 0.3; 196 196 197 G4bool Second_half = false; // First half 197 G4bool Second_half = false; // First half or second half of divided step 198 198 199 // We need to know this for the 'final_secti 199 // We need to know this for the 'final_section': 200 // real 'final_section' or first half 'final 200 // real 'final_section' or first half 'final_section' 201 // In algorithm it is considered that the 'S 201 // In algorithm it is considered that the 'Second_half' is true 202 // and it becomes false only if we are in th 202 // and it becomes false only if we are in the first-half of level 203 // depthness or if we are in the first secti 203 // depthness or if we are in the first section 204 204 205 unsigned int depth = 0; // Depth counts subd 205 unsigned int depth = 0; // Depth counts subdivisions of initial step made 206 ++fNumCalls; 206 ++fNumCalls; 207 207 208 NormalAtEntry = GetSurfaceNormal(CurrentE_Po 208 NormalAtEntry = GetSurfaceNormal(CurrentE_Point, validNormalAtE); 209 209 210 if (fCheckMode) 210 if (fCheckMode) 211 { 211 { 212 pointH_logger.AddRecord( G4LocatorChangeRe 212 pointH_logger.AddRecord( G4LocatorChangeRecord::kInitialisingCL, 213 substep_no, event 213 substep_no, eventCount, 214 G4FieldTrack( Cur 214 G4FieldTrack( CurrentE_Point,0.,NormalAtEntry,0., 215 0., 215 0., 1.,G4ThreeVector(0.),0.,0.) ); 216 #if (G4DEBUG_FIELD>1) 216 #if (G4DEBUG_FIELD>1) 217 G4ThreeVector StartPosition = CurveStartP 217 G4ThreeVector StartPosition = CurveStartPointVelocity.GetPosition(); 218 if( (TrialPoint - StartPosition).mag2() < 218 if( (TrialPoint - StartPosition).mag2() < sqr(tolerance)) 219 { 219 { 220 ReportImmediateHit( MethodName, StartPo 220 ReportImmediateHit( MethodName, StartPosition, TrialPoint, 221 tolerance, fNumCall 221 tolerance, fNumCalls); 222 } 222 } 223 #endif 223 #endif 224 } 224 } 225 225 226 // Intermediates Points on the Track = Subdi 226 // Intermediates Points on the Track = Subdivided Points must be stored. 227 // Give the initial values to 'InterMedFt' 227 // Give the initial values to 'InterMedFt' 228 // Important is 'ptrInterMedFT[0]', it saves 228 // Important is 'ptrInterMedFT[0]', it saves the 'EndCurvePoint' 229 // 229 // 230 *ptrInterMedFT[0] = CurveEndPointVelocity; 230 *ptrInterMedFT[0] = CurveEndPointVelocity; 231 for ( auto idepth=1; idepth<max_depth+1; ++i 231 for ( auto idepth=1; idepth<max_depth+1; ++idepth ) 232 { 232 { 233 *ptrInterMedFT[idepth] = CurveStartPointVe 233 *ptrInterMedFT[idepth] = CurveStartPointVelocity; 234 } 234 } 235 235 236 // Final_section boolean store 236 // Final_section boolean store 237 // 237 // 238 G4bool fin_section_depth[max_depth]; 238 G4bool fin_section_depth[max_depth]; 239 for (bool & idepth : fin_section_depth) << 239 for ( auto idepth=0; idepth<max_depth; ++idepth ) 240 { 240 { 241 idepth = true; << 241 fin_section_depth[idepth] = true; 242 } 242 } 243 // 'SubStartPoint' is needed to calculate th 243 // 'SubStartPoint' is needed to calculate the length of the divided step 244 // 244 // 245 G4FieldTrack SubStart_PointVelocity = CurveS 245 G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity; 246 246 247 do // Loop checking, 07.10.2016, J.Apostola 247 do // Loop checking, 07.10.2016, J.Apostolakis 248 { 248 { 249 unsigned int substep_no_p = 0; 249 unsigned int substep_no_p = 0; 250 G4bool sub_final_section = false; // the s 250 G4bool sub_final_section = false; // the same as final_section, 251 // but f 251 // but for 'sub_section' 252 SubStart_PointVelocity = CurrentA_PointVel 252 SubStart_PointVelocity = CurrentA_PointVelocity; 253 253 254 do // Loop checking, 07.10.2016, J.Apostol 254 do // Loop checking, 07.10.2016, J.Apostolakis 255 { // REPEAT param 255 { // REPEAT param 256 G4ThreeVector Point_A = CurrentA_PointVe << 256 257 G4ThreeVector Point_B = CurrentB_PointVe << 257 #ifdef G4DEBUG_FIELD 258 << 258 if( CurrentA_PointVelocity.GetCurveLength() >= 259 #ifdef G4DEBUG_FIELD << 259 CurrentB_PointVelocity.GetCurveLength() ) 260 const G4double lenA = CurrentA_PointVelo << 261 const G4double lenB = CurrentB_PointVelo << 262 G4double curv_lenAB = lenB - lenA; << 263 G4double distAB = (Point_B - Point_A << 264 if( curv_lenAB < distAB * ( 1. - 10.*fiE << 265 { 260 { 266 G4cerr << "ERROR> (Start) Point A coin 261 G4cerr << "ERROR> (Start) Point A coincides with or has gone past (end) point B" 267 << "MLL: iters = " << substep_n 262 << "MLL: iters = " << substep_no << G4endl; 268 G4long op=G4cerr.precision(6); << 263 // G4LocatorChangeRecord::ReportVector(G4cerr, "endPointB", endChangeB ); 269 G4cerr << " Difference = " << di << 264 // G4cerr<<"EndPoints A(start) and B(end): combined changes " << G4endl; 270 << " exceeds limit of relative << 265 G4LocatorChangeLogger::ReportEndChanges(G4cerr, endChangeA, endChangeB); 271 << " i.e. limit = " << 10 * fi << 272 G4cerr.precision(9); << 273 G4cerr << " Len A, B = " << len << 274 << " Position A: " << Po << 275 << " Position B: " << Po << 276 G4cerr.precision(op); << 277 // G4LocatorChangeRecord::ReportVector << 278 // G4cerr<<"EndPoints A(start) and B(e << 279 if (fCheckMode) << 280 { << 281 G4LocatorChangeLogger::ReportEndChan << 282 } << 283 } << 284 << 285 if( !validIntersectP ) << 286 { << 287 G4ExceptionDescription errmsg; << 288 errmsg << "Assertion FAILURE - invalid << 289 << substep_no << " call: " << f << 290 if (fCheckMode) << 291 { << 292 G4LocatorChangeRecord::ReportEndChan << 293 } << 294 G4Exception("G4MultiLevelLocator::Esti << 295 JustWarning, errmsg); << 296 } 266 } 297 #endif 267 #endif >> 268 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition(); >> 269 G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition(); 298 270 299 // F = a point on true AB path close to 271 // F = a point on true AB path close to point E 300 // (the closest if possible) 272 // (the closest if possible) 301 // 273 // 302 ApproxIntersecPointV = GetChordFinderFor 274 ApproxIntersecPointV = GetChordFinderFor() 303 ->ApproxCurvePoin 275 ->ApproxCurvePointV( CurrentA_PointVelocity, 304 276 CurrentB_PointVelocity, 305 277 CurrentE_Point, 306 278 GetEpsilonStepFor() ); 307 // The above method is the key & most 279 // The above method is the key & most intuitive part ... 308 280 309 #ifdef G4DEBUG_FIELD 281 #ifdef G4DEBUG_FIELD 310 recApproxPoint.push_back(G4LocatorChange 282 recApproxPoint.push_back(G4LocatorChangeRecord(G4LocatorChangeRecord::kInvalidCL, 311 283 substep_no, eventCount, ApproxIntersecPointV ) ); 312 G4double lenIntsc= ApproxIntersecPointV. 284 G4double lenIntsc= ApproxIntersecPointV.GetCurveLength(); >> 285 G4double lenB = CurrentB_PointVelocity.GetCurveLength(); 313 G4double checkVsEnd= lenB - lenIntsc; 286 G4double checkVsEnd= lenB - lenIntsc; 314 287 315 if( lenIntsc > lenB ) 288 if( lenIntsc > lenB ) 316 { 289 { 317 std::ostringstream errmsg; << 290 std::ostringstream errmsg; 318 errmsg.precision(17); << 291 errmsg.precision(17); 319 G4double ratio = checkVsEnd / lenB; << 292 G4double ratio = checkVsEnd / lenB; 320 G4double ratioTol = std::fabs(ratio) / << 293 G4double ratioTol = std::fabs(ratio) / tolerance; 321 errmsg << "Intermediate F point is pas << 294 errmsg << "Intermediate F point is past end B point" << G4endl 322 << " l( intersection ) = " << lenInt << 295 << " l( intersection ) = " << lenIntsc << G4endl 323 << " l( endpoint ) = " << lenB << 296 << " l( endpoint ) = " << lenB << G4endl; 324 errmsg.precision(8); << 297 errmsg.precision(8); 325 errmsg << " l_end - l_inters = " << << 298 errmsg << " l_end - l_inters = " << checkVsEnd << G4endl 326 << " / l_end = " << rati << 299 << " / l_end = " << ratio << G4endl 327 << " ratio / tolerance = " << rati << 300 << " ratio / tolerance = " << ratioTol << G4endl; 328 if( ratioTol < 1.0 ) << 301 if( ratioTol < 1.0 ) 329 G4Exception(MethodName, "GeomNav0003 << 302 G4Exception(MethodName, "GeomNav0003", JustWarning, errmsg ); 330 else << 303 else 331 G4Exception(MethodName, "GeomNav0003 << 304 G4Exception(MethodName, "GeomNav0003", FatalException, errmsg ); 332 } 305 } 333 #endif 306 #endif 334 307 335 G4ThreeVector CurrentF_Point= ApproxInte 308 G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition(); 336 309 337 // First check whether EF is small - the 310 // First check whether EF is small - then F is a good approx. point 338 // Calculate the length and direction of 311 // Calculate the length and direction of the chord AF 339 // 312 // 340 G4ThreeVector ChordEF_Vector = CurrentF 313 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point; 341 314 342 G4ThreeVector NewMomentumDir = ApproxIn 315 G4ThreeVector NewMomentumDir = ApproxIntersecPointV.GetMomentumDir(); 343 G4double MomDir_dot_Norm = NewMome 316 G4double MomDir_dot_Norm = NewMomentumDir.dot( NormalAtEntry ); 344 317 345 #ifdef G4DEBUG_FIELD 318 #ifdef G4DEBUG_FIELD 346 if( fVerboseLevel > 3 ) 319 if( fVerboseLevel > 3 ) 347 { 320 { 348 G4ThreeVector ChordAB = Po 321 G4ThreeVector ChordAB = Point_B - Point_A; 349 G4double ABchord_length = Ch 322 G4double ABchord_length = ChordAB.mag(); 350 G4double MomDir_dot_ABchord; 323 G4double MomDir_dot_ABchord; 351 MomDir_dot_ABchord = (1.0 / ABchord_l 324 MomDir_dot_ABchord = (1.0 / ABchord_length) 352 * NewMomentumDir.d 325 * NewMomentumDir.dot( ChordAB ); 353 G4VIntersectionLocator::ReportTrialSt 326 G4VIntersectionLocator::ReportTrialStep( substep_no, ChordAB, 354 ChordEF_Vector, NewMomentumDir, 327 ChordEF_Vector, NewMomentumDir, NormalAtEntry, validNormalAtE ); 355 G4cout << " dot( MomentumDir, ABchord 328 G4cout << " dot( MomentumDir, ABchord_unit ) = " 356 << MomDir_dot_ABchord << G4end 329 << MomDir_dot_ABchord << G4endl; 357 } 330 } 358 #endif 331 #endif 359 G4bool adequate_angle = 332 G4bool adequate_angle = 360 ( MomDir_dot_Norm >= 0.0 ) // Can 333 ( MomDir_dot_Norm >= 0.0 ) // Can use ( > -epsilon) instead 361 || (! validNormalAtE ); // Inv 334 || (! validNormalAtE ); // Invalid, cannot use this criterion 362 G4double EF_dist2 = ChordEF_Vector.mag2( 335 G4double EF_dist2 = ChordEF_Vector.mag2(); 363 if ( ( EF_dist2 <= sqr(fiDeltaIntersecti 336 if ( ( EF_dist2 <= sqr(fiDeltaIntersection) && ( adequate_angle ) ) 364 || ( EF_dist2 <= kCarTolerance*kCarTol 337 || ( EF_dist2 <= kCarTolerance*kCarTolerance ) ) 365 { 338 { 366 found_approximate_intersection = true; 339 found_approximate_intersection = true; 367 340 368 // Create the "point" return value 341 // Create the "point" return value 369 // 342 // 370 IntersectedOrRecalculatedFT = ApproxIn 343 IntersectedOrRecalculatedFT = ApproxIntersecPointV; 371 IntersectedOrRecalculatedFT.SetPositio 344 IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point ); 372 345 373 if ( GetAdjustementOfFoundIntersection 346 if ( GetAdjustementOfFoundIntersection() ) 374 { 347 { 375 // Try to Get Correction of Intersec 348 // Try to Get Correction of IntersectionPoint using SurfaceNormal() 376 // 349 // 377 G4ThreeVector IP; 350 G4ThreeVector IP; 378 G4ThreeVector MomentumDir=ApproxInte 351 G4ThreeVector MomentumDir=ApproxIntersecPointV.GetMomentumDirection(); 379 G4bool goodCorrection = AdjustmentOf 352 G4bool goodCorrection = AdjustmentOfFoundIntersection(Point_A, 380 CurrentE_P 353 CurrentE_Point, CurrentF_Point, MomentumDir, 381 last_AF_in 354 last_AF_intersection, IP, NewSafety, 382 previousSa 355 previousSafety, previousSftOrigin ); 383 if ( goodCorrection ) 356 if ( goodCorrection ) 384 { 357 { 385 IntersectedOrRecalculatedFT = Appr 358 IntersectedOrRecalculatedFT = ApproxIntersecPointV; 386 IntersectedOrRecalculatedFT.SetPos 359 IntersectedOrRecalculatedFT.SetPosition(IP); 387 } 360 } 388 } 361 } 389 // Note: in order to return a point on 362 // Note: in order to return a point on the boundary, 390 // we must return E. But it is F 363 // we must return E. But it is F on the curve. 391 // So we must "cheat": we are us 364 // So we must "cheat": we are using the position at point E 392 // and the velocity at point F ! 365 // and the velocity at point F !!! 393 // 366 // 394 // This must limit the length we can a 367 // This must limit the length we can allow for displacement! 395 } 368 } 396 else // E is NOT close enough to the cu 369 else // E is NOT close enough to the curve (ie point F) 397 { 370 { 398 // Check whether any volumes are encou 371 // Check whether any volumes are encountered by the chord AF 399 // ----------------------------------- 372 // --------------------------------------------------------- 400 // First relocate to restore any Voxel 373 // First relocate to restore any Voxel etc information 401 // in the Navigator before calling Com 374 // in the Navigator before calling ComputeStep() 402 // 375 // 403 GetNavigatorFor()->LocateGlobalPointWi 376 GetNavigatorFor()->LocateGlobalPointWithinVolume( Point_A ); 404 377 405 G4ThreeVector PointG; // Candidate i 378 G4ThreeVector PointG; // Candidate intersection point 406 G4double stepLengthAF; 379 G4double stepLengthAF; 407 G4bool Intersects_FB = false; 380 G4bool Intersects_FB = false; 408 G4bool Intersects_AF = IntersectChord( 381 G4bool Intersects_AF = IntersectChord( Point_A, CurrentF_Point, 409 382 NewSafety, previousSafety, 410 383 previousSftOrigin, 411 384 stepLengthAF, 412 385 PointG ); 413 last_AF_intersection = Intersects_AF; 386 last_AF_intersection = Intersects_AF; 414 if( Intersects_AF ) 387 if( Intersects_AF ) 415 { 388 { 416 // G is our new Candidate for the in 389 // G is our new Candidate for the intersection point. 417 // It replaces "E" and we will see 390 // It replaces "E" and we will see if it's good enough. 418 CurrentB_PointVelocity = ApproxInter 391 CurrentB_PointVelocity = ApproxIntersecPointV; // B <- F 419 CurrentE_Point = PointG; 392 CurrentE_Point = PointG; // E <- G 420 393 421 validIntersectP = true; // 'E' ha 394 validIntersectP = true; // 'E' has been updated. 422 395 423 G4bool validNormalLast; 396 G4bool validNormalLast; 424 NormalAtEntry = GetSurfaceNormal( P 397 NormalAtEntry = GetSurfaceNormal( PointG, validNormalLast ); 425 validNormalAtE = validNormalLast; 398 validNormalAtE = validNormalLast; 426 399 427 // As we move point B, must take car 400 // As we move point B, must take care in case the current 428 // AF has no intersection to try cur 401 // AF has no intersection to try current FB!! 429 fin_section_depth[depth] = false; 402 fin_section_depth[depth] = false; 430 403 431 if (fCheckMode) 404 if (fCheckMode) 432 { 405 { 433 ++eventCount; 406 ++eventCount; 434 endChangeB.push_back( 407 endChangeB.push_back( 435 G4LocatorChangeRecord(G4LocatorC 408 G4LocatorChangeRecord(G4LocatorChangeRecord::kIntersectsAF, 436 substep_no, event 409 substep_no, eventCount, CurrentB_PointVelocity) ); 437 } 410 } 438 #ifdef G4VERBOSE 411 #ifdef G4VERBOSE 439 if( fVerboseLevel > 3 ) 412 if( fVerboseLevel > 3 ) 440 { 413 { 441 G4cout << "G4PiF::LI> Investigatin 414 G4cout << "G4PiF::LI> Investigating intermediate point" 442 << " at s=" << ApproxInters 415 << " at s=" << ApproxIntersecPointV.GetCurveLength() 443 << " on way to full s=" 416 << " on way to full s=" 444 << CurveEndPointVelocity.Ge 417 << CurveEndPointVelocity.GetCurveLength() << G4endl; 445 } 418 } 446 #endif 419 #endif 447 } 420 } 448 else // not Intersects_AF 421 else // not Intersects_AF 449 { 422 { 450 // In this case: 423 // In this case: 451 // There is NO intersection of AF wi 424 // There is NO intersection of AF with a volume boundary. 452 // We must continue the search in th 425 // We must continue the search in the segment FB! 453 // 426 // 454 GetNavigatorFor()->LocateGlobalPoint 427 GetNavigatorFor()->LocateGlobalPointWithinVolume( CurrentF_Point ); 455 428 456 G4double stepLengthFB; 429 G4double stepLengthFB; 457 G4ThreeVector PointH; 430 G4ThreeVector PointH; 458 431 459 // Check whether any volumes are enc 432 // Check whether any volumes are encountered by the chord FB 460 // --------------------------------- 433 // --------------------------------------------------------- 461 434 462 Intersects_FB = IntersectChord( Curr 435 Intersects_FB = IntersectChord( CurrentF_Point, Point_B, 463 NewS 436 NewSafety, previousSafety, 464 prev 437 previousSftOrigin, 465 step 438 stepLengthFB, 466 Poin 439 PointH ); 467 if( Intersects_FB ) 440 if( Intersects_FB ) 468 { 441 { 469 // There is an intersection of FB 442 // There is an intersection of FB with a volume boundary 470 // H <- First Intersection of Chor 443 // H <- First Intersection of Chord FB 471 444 472 // H is our new Candidate for the 445 // H is our new Candidate for the intersection point. 473 // It replaces "E" and we will re 446 // It replaces "E" and we will repeat the test to see if 474 // it is a good enough approximate 447 // it is a good enough approximate point for us. 475 448 476 // Note that F must be in volume v 449 // Note that F must be in volume volA (the same as A) 477 // (otherwise AF would meet a volu 450 // (otherwise AF would meet a volume boundary!) 478 // A <- F 451 // A <- F 479 // E <- H 452 // E <- H 480 // 453 // 481 CurrentA_PointVelocity = ApproxInt 454 CurrentA_PointVelocity = ApproxIntersecPointV; 482 CurrentE_Point = PointH; 455 CurrentE_Point = PointH; 483 456 484 validIntersectP = true; // 'E' 457 validIntersectP = true; // 'E' has been updated. 485 458 486 G4bool validNormalH; 459 G4bool validNormalH; 487 NormalAtEntry = GetSurfaceNormal( 460 NormalAtEntry = GetSurfaceNormal( PointH, validNormalH ); 488 validNormalAtE = validNormalH; 461 validNormalAtE = validNormalH; 489 462 490 if (fCheckMode) 463 if (fCheckMode) 491 { 464 { 492 ++eventCount; 465 ++eventCount; 493 endChangeA.push_back( 466 endChangeA.push_back( 494 G4LocatorChangeRecord(G4Locat 467 G4LocatorChangeRecord(G4LocatorChangeRecord::kIntersectsFB, 495 substep_no, event 468 substep_no, eventCount, CurrentA_PointVelocity) ); 496 G4FieldTrack intersectH_pn('0'); 469 G4FieldTrack intersectH_pn('0'); // Point and normal 497 470 // nothing else will be valid 498 intersectH_pn.SetPosition( Point 471 intersectH_pn.SetPosition( PointH ); 499 intersectH_pn.SetMomentum( Norma 472 intersectH_pn.SetMomentum( NormalAtEntry ); 500 pointH_logger.AddRecord(G4Locato 473 pointH_logger.AddRecord(G4LocatorChangeRecord::kIntersectsFB, 501 substep_no 474 substep_no, eventCount, intersectH_pn ); 502 } 475 } 503 } 476 } 504 else // not Intersects_FB 477 else // not Intersects_FB 505 { 478 { 506 validIntersectP = false; // Int << 507 if( fin_section_depth[depth] ) 479 if( fin_section_depth[depth] ) 508 { 480 { 509 // If B is the original endpoint 481 // If B is the original endpoint, this means that whatever 510 // volume(s) intersected the ori 482 // volume(s) intersected the original chord, none touch the 511 // smaller chords we have used. 483 // smaller chords we have used. 512 // The value of 'IntersectedOrRe 484 // The value of 'IntersectedOrRecalculatedFT' returned is 513 // likely not valid 485 // likely not valid 514 486 515 // Check on real final_section o 487 // Check on real final_section or SubEndSection 516 // 488 // 517 if( ((Second_half)&&(depth==0)) 489 if( ((Second_half)&&(depth==0)) || (first_section) ) 518 { 490 { 519 there_is_no_intersection = tru 491 there_is_no_intersection = true; // real final_section 520 } 492 } 521 else 493 else 522 { 494 { 523 // end of subsection, not real 495 // end of subsection, not real final section 524 // exit from the and go to the 496 // exit from the and go to the depth-1 level 525 substep_no_p = param_substeps+ 497 substep_no_p = param_substeps+2; // exit from the loop 526 498 527 // but 'Second_half' is still 499 // but 'Second_half' is still true because we need to find 528 // the 'CurrentE_point' for th 500 // the 'CurrentE_point' for the next loop 529 Second_half = true; 501 Second_half = true; 530 sub_final_section = true; 502 sub_final_section = true; 531 } 503 } 532 } 504 } 533 else 505 else 534 { 506 { 535 CurrentA_PointVelocity = Current 507 CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B 536 CurrentB_PointVelocity = (depth= 508 CurrentB_PointVelocity = (depth==0) ? CurveEndPointVelocity 537 509 : *ptrInterMedFT[depth] ; 538 SubStart_PointVelocity = Current 510 SubStart_PointVelocity = CurrentA_PointVelocity; 539 restoredFullEndpoint = true; 511 restoredFullEndpoint = true; 540 512 541 validIntersectP = false; // 'E' 513 validIntersectP = false; // 'E' has NOT been updated. 542 514 543 if (fCheckMode) 515 if (fCheckMode) 544 { 516 { 545 ++eventCount; 517 ++eventCount; 546 endChangeA.push_back( 518 endChangeA.push_back( 547 G4LocatorChangeRecord( 519 G4LocatorChangeRecord( 548 G4LocatorChangeRecord::kNo 520 G4LocatorChangeRecord::kNoIntersectAForFB, 549 substep_no, event 521 substep_no, eventCount, CurrentA_PointVelocity) ); 550 endChangeB.push_back( 522 endChangeB.push_back( 551 G4LocatorChangeRecord ( 523 G4LocatorChangeRecord ( 552 G4LocatorChangeRecord::kNo 524 G4LocatorChangeRecord::kNoIntersectAForFB, 553 substep_no, event 525 substep_no, eventCount, CurrentB_PointVelocity) ); 554 } 526 } 555 } 527 } 556 } // Endif (Intersects_FB) 528 } // Endif (Intersects_FB) 557 } // Endif (Intersects_AF) 529 } // Endif (Intersects_AF) 558 530 559 G4int errorEndPt = 0; // Default: no e 531 G4int errorEndPt = 0; // Default: no error (if not calling CheckAnd... 560 532 561 G4bool recalculatedB= false; 533 G4bool recalculatedB= false; 562 if( driverReIntegrates ) 534 if( driverReIntegrates ) 563 { 535 { 564 G4FieldTrack RevisedB_FT = CurrentB_ 536 G4FieldTrack RevisedB_FT = CurrentB_PointVelocity; 565 recalculatedB= CheckAndReEstimateEnd 537 recalculatedB= CheckAndReEstimateEndpoint(CurrentA_PointVelocity, 566 538 CurrentB_PointVelocity, 567 539 RevisedB_FT, 568 540 errorEndPt ); 569 if( recalculatedB ) 541 if( recalculatedB ) 570 { 542 { 571 CurrentB_PointVelocity = RevisedB_ 543 CurrentB_PointVelocity = RevisedB_FT; // Use it ! 572 // Do not invalidate intersection 544 // Do not invalidate intersection F -- it is still roughly OK. 573 // 545 // 574 // The best course would be to inv 546 // The best course would be to invalidate (reset validIntersectP) 575 // BUT if we invalidate it, we mus 547 // BUT if we invalidate it, we must re-estimate it somewhere! E.g. 576 // validIntersectP = false; // 548 // validIntersectP = false; // 'E' has NOT been updated. 577 549 578 if ( (fin_section_depth[depth]) 550 if ( (fin_section_depth[depth]) // real final section 579 &&( first_section || ((Second_h 551 &&( first_section || ((Second_half)&&(depth==0)) ) ) 580 { 552 { 581 recalculatedEndPoint = true; 553 recalculatedEndPoint = true; 582 IntersectedOrRecalculatedFT = Re 554 IntersectedOrRecalculatedFT = RevisedB_FT; 583 // So that we can return it, if 555 // So that we can return it, if it is the endpoint! 584 } 556 } 585 // else 557 // else 586 // Move forward the other points 558 // Move forward the other points 587 // - or better flag it, so that 559 // - or better flag it, so that they are re-computed when next used 588 // [ Implementation: a counter 560 // [ Implementation: a counter for # of recomputations 589 // => avoids extra work] 561 // => avoids extra work] 590 } 562 } >> 563 // else >> 564 // Move forward the other points >> 565 // - or better flag it, so that they are re-computed when next used >> 566 // [ Implementation: a counter for # of recomputations >> 567 // => avoids extra work] 591 if (fCheckMode) 568 if (fCheckMode) 592 { 569 { 593 ++eventCount; 570 ++eventCount; 594 endChangeB.push_back( 571 endChangeB.push_back( 595 G4LocatorChangeRecord( G4Locator 572 G4LocatorChangeRecord( G4LocatorChangeRecord::kRecalculatedB, 596 substep_n 573 substep_no, eventCount, RevisedB_FT ) ); 597 } 574 } 598 } 575 } 599 else 576 else 600 { 577 { 601 if( CurrentB_PointVelocity.GetCurveL << 578 if( CurrentB_PointVelocity.GetCurveLength() < CurrentA_PointVelocity.GetCurveLength() ) 602 < CurrentA_PointVelocity.GetCurveL << 603 { << 604 errorEndPt = 2; 579 errorEndPt = 2; 605 } << 606 } 580 } 607 << 581 608 if( errorEndPt > 1 ) // errorEndPt = 582 if( errorEndPt > 1 ) // errorEndPt = 1 is milder, just: len(B)=len(A) 609 { 583 { 610 std::ostringstream errmsg; << 584 std::ostringstream errmsg; 611 << 585 ReportReversedPoints(errmsg, 612 ReportReversedPoints(errmsg, << 586 CurveStartPointVelocity, CurveEndPointVelocity, 613 CurveStartPoint << 587 NewSafety, fiEpsilonStep, 614 NewSafety, fiEp << 588 CurrentA_PointVelocity, CurrentB_PointVelocity, 615 CurrentA_PointV << 589 SubStart_PointVelocity, CurrentE_Point, 616 SubStart_PointV << 590 ApproxIntersecPointV, substep_no, substep_no_p, depth); 617 ApproxIntersecP << 591 errmsg << G4endl << " * Location: " << MethodName 618 << 592 << "- After EndIf(Intersects_AF)" << G4endl; 619 if (fCheckMode) << 593 errmsg << " * Bool flags: Recalculated = " << recalculatedB 620 { << 594 << " Intersects_AF = " << Intersects_AF 621 G4LocatorChangeRecord::ReportEndCh << 595 << " Intersects_FB = " << Intersects_FB << G4endl; 622 } << 596 errmsg << " * Number of calls to MLL:EIP= " << fNumCalls << G4endl; 623 << 597 G4Exception(MethodName, "GeomNav0003", FatalException, errmsg); 624 errmsg << G4endl << " * Location: " << 625 << "- After EndIf(Intersects_ << 626 errmsg << " * Bool flags: Recalcula << 627 << " Intersects_AF = " << I << 628 << " Intersects_FB = " << I << 629 errmsg << " * Number of calls to MLL << 630 G4Exception(MethodName, "GeomNav0003 << 631 } 598 } 632 if( restoredFullEndpoint ) 599 if( restoredFullEndpoint ) 633 { 600 { 634 fin_section_depth[depth] = restoredF 601 fin_section_depth[depth] = restoredFullEndpoint; 635 restoredFullEndpoint = false; 602 restoredFullEndpoint = false; 636 } 603 } 637 } // EndIf ( E is close enough to the cu 604 } // EndIf ( E is close enough to the curve, ie point F. ) 638 // tests ChordAF_Vector.mag() <= maxim 605 // tests ChordAF_Vector.mag() <= maximum_lateral_displacement 639 606 640 #ifdef G4DEBUG_FIELD 607 #ifdef G4DEBUG_FIELD 641 if( trigger_substepno_print == 0) 608 if( trigger_substepno_print == 0) 642 { 609 { 643 trigger_substepno_print= fWarnSteps - 610 trigger_substepno_print= fWarnSteps - 20; 644 } 611 } 645 612 646 if( substep_no >= trigger_substepno_prin 613 if( substep_no >= trigger_substepno_print ) 647 { 614 { 648 G4cout << "Difficulty in converging in 615 G4cout << "Difficulty in converging in " << MethodName 649 << " Substep no = " << subst 616 << " Substep no = " << substep_no << G4endl; 650 if( substep_no == trigger_substepno_pr 617 if( substep_no == trigger_substepno_print ) 651 { 618 { 652 G4cout << " Start: "; 619 G4cout << " Start: "; 653 printStatus( CurveStartPointVelocity 620 printStatus( CurveStartPointVelocity, CurveEndPointVelocity, 654 -1.0, NewSafety, 0 ); 621 -1.0, NewSafety, 0 ); 655 if( fCheckMode ) { << 622 656 G4LocatorChangeRecord::ReportEndCh << 623 G4cout << " ** Change records: " << G4endl; 657 } else { << 624 G4cout << "endPoints A (start) and B (end): combined changes of AB intervals" << G4endl; 658 G4cout << " ** For more informatio << 625 G4LocatorChangeRecord::ReportEndChanges(G4cout, endChangeA, endChangeB ); 659 << "-- (it saves and can ou << 660 } << 661 } 626 } 662 G4cout << " Point A: "; 627 G4cout << " Point A: "; 663 printStatus( CurrentA_PointVelocity, C 628 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity, 664 -1.0, NewSafety, substep_ 629 -1.0, NewSafety, substep_no-1 ); 665 G4cout << " Point B: "; 630 G4cout << " Point B: "; 666 printStatus( CurrentA_PointVelocity, C 631 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity, 667 -1.0, NewSafety, substep_ 632 -1.0, NewSafety, substep_no ); 668 } 633 } 669 #endif 634 #endif 670 ++substep_no; 635 ++substep_no; 671 ++substep_no_p; 636 ++substep_no_p; 672 637 673 } while ( ( ! found_approximate_intersect 638 } while ( ( ! found_approximate_intersection ) 674 && ( ! there_is_no_intersection ) 639 && ( ! there_is_no_intersection ) 675 && validIntersectP // New c << 676 && ( substep_no_p <= param_substep 640 && ( substep_no_p <= param_substeps) ); // UNTIL found or 677 641 // failed param substep 678 642 679 if( (!found_approximate_intersection) && ( 643 if( (!found_approximate_intersection) && (!there_is_no_intersection) ) 680 { 644 { 681 G4double did_len = std::abs( CurrentA_Po 645 G4double did_len = std::abs( CurrentA_PointVelocity.GetCurveLength() 682 - SubStart_PointVelocit 646 - SubStart_PointVelocity.GetCurveLength()); 683 G4double all_len = std::abs( CurrentB_Po 647 G4double all_len = std::abs( CurrentB_PointVelocity.GetCurveLength() 684 - SubStart_PointVelocit 648 - SubStart_PointVelocity.GetCurveLength()); 685 649 686 G4double distAB = -1; 650 G4double distAB = -1; 687 // 651 // 688 // Is progress is too slow, and is it po 652 // Is progress is too slow, and is it possible to go deeper? 689 // If so, then *halve the step* 653 // If so, then *halve the step* 690 // ============== 654 // ============== 691 if( (did_len < fraction_done*all_len) 655 if( (did_len < fraction_done*all_len) 692 && (depth<max_depth) && (!sub_final_s 656 && (depth<max_depth) && (!sub_final_section) ) 693 { 657 { 694 #ifdef G4DEBUG_FIELD 658 #ifdef G4DEBUG_FIELD 695 static G4ThreadLocal unsigned int numS 659 static G4ThreadLocal unsigned int numSplits=0; // For debugging only 696 biggest_depth = std::max(depth, bigges 660 biggest_depth = std::max(depth, biggest_depth); 697 ++numSplits; 661 ++numSplits; 698 #endif 662 #endif 699 Second_half = false; 663 Second_half = false; 700 ++depth; 664 ++depth; 701 first_section = false; 665 first_section = false; 702 666 703 G4double Sub_len = (all_len-did_len)/( 667 G4double Sub_len = (all_len-did_len)/(2.); 704 G4FieldTrack midPoint = CurrentA_Point 668 G4FieldTrack midPoint = CurrentA_PointVelocity; 705 G4bool fullAdvance= 669 G4bool fullAdvance= 706 integrDriver->AccurateAdvance(midPo 670 integrDriver->AccurateAdvance(midPoint, Sub_len, fiEpsilonStep); 707 671 708 ++fNumAdvanceTrials; 672 ++fNumAdvanceTrials; 709 if( fullAdvance ) { ++fNumAdvanceFull 673 if( fullAdvance ) { ++fNumAdvanceFull; } 710 674 711 G4double lenAchieved= 675 G4double lenAchieved= 712 midPoint.GetCurveLength()-CurrentA_ 676 midPoint.GetCurveLength()-CurrentA_PointVelocity.GetCurveLength(); 713 677 714 const G4double adequateFraction = (1.0 678 const G4double adequateFraction = (1.0-CLHEP::perThousand); 715 G4bool goodAdvance = (lenAchieved >= a 679 G4bool goodAdvance = (lenAchieved >= adequateFraction * Sub_len); 716 if ( goodAdvance ) { ++fNumAdvanceGoo 680 if ( goodAdvance ) { ++fNumAdvanceGood; } 717 681 718 #ifdef G4DEBUG_FIELD 682 #ifdef G4DEBUG_FIELD 719 else // !goodAdvance 683 else // !goodAdvance 720 { 684 { 721 G4cout << "MLL> AdvanceChordLimited 685 G4cout << "MLL> AdvanceChordLimited not full at depth=" << depth 722 << " total length achieved 686 << " total length achieved = " << lenAchieved << " of " 723 << Sub_len << " fraction= " 687 << Sub_len << " fraction= "; 724 if (Sub_len != 0.0 ) { G4cout << le 688 if (Sub_len != 0.0 ) { G4cout << lenAchieved / Sub_len; } 725 else { G4cout << "D 689 else { G4cout << "DivByZero"; } 726 G4cout << " Good-enough fraction = 690 G4cout << " Good-enough fraction = " << adequateFraction; 727 G4cout << " Number of call to mll 691 G4cout << " Number of call to mll = " << fNumCalls 728 << " iteration " << subste 692 << " iteration " << substep_no 729 << " inner = " << substep_ 693 << " inner = " << substep_no_p << G4endl; 730 G4cout << " Epsilon of integration 694 G4cout << " Epsilon of integration = " << fiEpsilonStep << G4endl; 731 G4cout << " State at start is 695 G4cout << " State at start is = " << CurrentA_PointVelocity 732 << G4endl 696 << G4endl 733 << " at end (midpoint 697 << " at end (midpoint)= " << midPoint << G4endl; 734 G4cout << " Particle mass = " << m 698 G4cout << " Particle mass = " << midPoint.GetRestMass() << G4endl; 735 699 736 G4EquationOfMotion *equation = inte 700 G4EquationOfMotion *equation = integrDriver->GetEquationOfMotion(); 737 ReportFieldValue( CurrentA_PointVel 701 ReportFieldValue( CurrentA_PointVelocity, "start", equation ); 738 ReportFieldValue( midPoint, "midPoi 702 ReportFieldValue( midPoint, "midPoint", equation ); 739 G4cout << " Original Start = " 703 G4cout << " Original Start = " 740 << CurveStartPointVelocity < 704 << CurveStartPointVelocity << G4endl; 741 G4cout << " Original End = " 705 G4cout << " Original End = " 742 << CurveEndPointVelocity << 706 << CurveEndPointVelocity << G4endl; 743 G4cout << " Original TrialPoint = 707 G4cout << " Original TrialPoint = " 744 << TrialPoint << G4endl; 708 << TrialPoint << G4endl; 745 G4cout << " (this is STRICT mode) 709 G4cout << " (this is STRICT mode) " 746 << " num Splits= " << numSp 710 << " num Splits= " << numSplits; 747 G4cout << G4endl; 711 G4cout << G4endl; 748 } 712 } 749 #endif 713 #endif 750 714 751 *ptrInterMedFT[depth] = midPoint; 715 *ptrInterMedFT[depth] = midPoint; 752 CurrentB_PointVelocity = midPoint; 716 CurrentB_PointVelocity = midPoint; 753 717 754 if (fCheckMode) 718 if (fCheckMode) 755 { 719 { 756 ++eventCount; 720 ++eventCount; 757 endChangeB.push_back( 721 endChangeB.push_back( 758 G4LocatorChangeRecord( G4LocatorCh 722 G4LocatorChangeRecord( G4LocatorChangeRecord::kInsertingMidPoint, 759 substep_no, 723 substep_no, eventCount, midPoint) ); 760 } 724 } 761 725 762 // Adjust 'SubStartPoint' to calculate 726 // Adjust 'SubStartPoint' to calculate the 'did_length' in next loop 763 // 727 // 764 SubStart_PointVelocity = CurrentA_Poin 728 SubStart_PointVelocity = CurrentA_PointVelocity; 765 729 766 // Find new trial intersection point n 730 // Find new trial intersection point needed at start of the loop 767 // 731 // 768 G4ThreeVector Point_A = CurrentA_Point 732 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition(); 769 G4ThreeVector Point_B = CurrentB_Point 733 G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition(); 770 734 771 G4ThreeVector PointGe; 735 G4ThreeVector PointGe; 772 GetNavigatorFor()->LocateGlobalPointWi 736 GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_A); 773 G4bool Intersects_AB = IntersectChord( 737 G4bool Intersects_AB = IntersectChord(Point_A, Point_B, 774 738 NewSafety, previousSafety, 775 739 previousSftOrigin, distAB, 776 740 PointGe); 777 if( Intersects_AB ) 741 if( Intersects_AB ) 778 { 742 { 779 last_AF_intersection = Intersects_AB 743 last_AF_intersection = Intersects_AB; 780 CurrentE_Point = PointGe; 744 CurrentE_Point = PointGe; 781 fin_section_depth[depth] = true; 745 fin_section_depth[depth] = true; 782 746 783 validIntersectP = true; // 'E' has 747 validIntersectP = true; // 'E' has been updated. 784 748 785 G4bool validNormalAB; 749 G4bool validNormalAB; 786 NormalAtEntry = GetSurfaceNormal( P 750 NormalAtEntry = GetSurfaceNormal( PointGe, validNormalAB ); 787 validNormalAtE = validNormalAB; 751 validNormalAtE = validNormalAB; 788 } 752 } 789 else 753 else 790 { 754 { 791 // No intersection found for first p 755 // No intersection found for first part of curve 792 // (CurrentA,InterMedPoint[depth]). 756 // (CurrentA,InterMedPoint[depth]). Go to the second part 793 // 757 // 794 Second_half = true; 758 Second_half = true; 795 759 796 validIntersectP= false; // No new 760 validIntersectP= false; // No new 'E' chord intersection found 797 } 761 } 798 } // if did_len 762 } // if did_len 799 763 >> 764 unsigned int levelPops = 0; >> 765 800 G4bool unfinished = Second_half; 766 G4bool unfinished = Second_half; 801 while ( unfinished && (depth>0) ) // Lo 767 while ( unfinished && (depth>0) ) // Loop checking, 07.10.2016, JA 802 { 768 { 803 // Second part of curve (InterMed[dept 769 // Second part of curve (InterMed[depth],Intermed[depth-1])) 804 // On the depth-1 level normally we ar 770 // On the depth-1 level normally we are on the 'second_half' 805 771 >> 772 ++levelPops; >> 773 806 // Find new trial intersection point 774 // Find new trial intersection point needed at start of the loop 807 // 775 // 808 SubStart_PointVelocity = *ptrInterMedF 776 SubStart_PointVelocity = *ptrInterMedFT[depth]; 809 CurrentA_PointVelocity = *ptrInterMedF 777 CurrentA_PointVelocity = *ptrInterMedFT[depth]; 810 CurrentB_PointVelocity = *ptrInterMedF 778 CurrentB_PointVelocity = *ptrInterMedFT[depth-1]; 811 779 812 if (fCheckMode) 780 if (fCheckMode) 813 { 781 { 814 ++eventCount; 782 ++eventCount; 815 G4LocatorChangeRecord chngPop_a( G4L 783 G4LocatorChangeRecord chngPop_a( G4LocatorChangeRecord::kLevelPop, 816 substep_no, even 784 substep_no, eventCount, CurrentA_PointVelocity); 817 endChangeA.push_back( chngPop_a ); 785 endChangeA.push_back( chngPop_a ); 818 G4LocatorChangeRecord chngPop_b( G4L 786 G4LocatorChangeRecord chngPop_b( G4LocatorChangeRecord::kLevelPop, 819 substep_no, even 787 substep_no, eventCount, CurrentB_PointVelocity); 820 endChangeB.push_back( chngPop_b ); 788 endChangeB.push_back( chngPop_b ); 821 } 789 } 822 790 823 // Ensure that the new endpoints are n 791 // Ensure that the new endpoints are not further apart in space 824 // than on the curve due to different 792 // than on the curve due to different errors in the integration 825 // 793 // 826 G4int errorEndPt = -1; 794 G4int errorEndPt = -1; 827 G4bool recalculatedB= false; 795 G4bool recalculatedB= false; 828 if( driverReIntegrates ) 796 if( driverReIntegrates ) 829 { 797 { 830 // Ensure that the new endpoints are 798 // Ensure that the new endpoints are not further apart in space 831 // than on the curve due to differen 799 // than on the curve due to different errors in the integration 832 // 800 // 833 G4FieldTrack RevisedEndPointFT = Cur 801 G4FieldTrack RevisedEndPointFT = CurrentB_PointVelocity; 834 recalculatedB = 802 recalculatedB = 835 CheckAndReEstimateEndpoint( Curre 803 CheckAndReEstimateEndpoint( CurrentA_PointVelocity, 836 Curre 804 CurrentB_PointVelocity, 837 Revis 805 RevisedEndPointFT, 838 error 806 errorEndPt ); 839 if( recalculatedB ) 807 if( recalculatedB ) 840 { 808 { 841 CurrentB_PointVelocity = RevisedEn 809 CurrentB_PointVelocity = RevisedEndPointFT; // Use it ! 842 810 843 if ( depth == 1 ) 811 if ( depth == 1 ) 844 { 812 { 845 recalculatedEndPoint = true; 813 recalculatedEndPoint = true; 846 IntersectedOrRecalculatedFT = R 814 IntersectedOrRecalculatedFT = RevisedEndPointFT; 847 // So that we can return it, if 815 // So that we can return it, if it is the endpoint! 848 } 816 } 849 } 817 } 850 else 818 else 851 { 819 { 852 if( CurrentB_PointVelocity.GetCurv << 820 if( CurrentB_PointVelocity.GetCurveLength() < CurrentA_PointVelocity.GetCurveLength() ) 853 < CurrentA_PointVelocity.GetCurv << 854 { << 855 errorEndPt = 2; 821 errorEndPt = 2; 856 } << 857 } 822 } 858 823 859 if (fCheckMode) 824 if (fCheckMode) 860 { 825 { 861 ++eventCount; 826 ++eventCount; 862 endChangeB.push_back( 827 endChangeB.push_back( 863 G4LocatorChangeRecord(G4LocatorC 828 G4LocatorChangeRecord(G4LocatorChangeRecord::kRecalculatedBagn, 864 substep_no 829 substep_no, eventCount, RevisedEndPointFT)); 865 } 830 } 866 } 831 } 867 if( errorEndPt > 1 ) // errorEndPt = 832 if( errorEndPt > 1 ) // errorEndPt = 1 is milder, just: len(B)=len(A) 868 { 833 { 869 std::ostringstream errmsg; 834 std::ostringstream errmsg; 870 ReportReversedPoints(errmsg, 835 ReportReversedPoints(errmsg, 871 CurveStartPointVelocity, C 836 CurveStartPointVelocity, CurveEndPointVelocity, 872 NewSafety, fiEpsilonStep, 837 NewSafety, fiEpsilonStep, 873 CurrentA_PointVelocity, Cu 838 CurrentA_PointVelocity, CurrentB_PointVelocity, 874 SubStart_PointVelocity, Cu 839 SubStart_PointVelocity, CurrentE_Point, 875 ApproxIntersecPointV, subs 840 ApproxIntersecPointV, substep_no, substep_no_p, depth); 876 errmsg << " * Location: " << Method 841 errmsg << " * Location: " << MethodName << "- Second-Half" << G4endl; 877 errmsg << " * Recalculated = " << re 842 errmsg << " * Recalculated = " << recalculatedB << G4endl; // false 878 G4Exception(MethodName, "GeomNav0003 843 G4Exception(MethodName, "GeomNav0003", FatalException, errmsg); 879 } 844 } 880 845 881 G4ThreeVector Point_A = CurrentA_Point 846 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition(); 882 G4ThreeVector Point_B = CurrentB_Point 847 G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition(); 883 G4ThreeVector PointGi; 848 G4ThreeVector PointGi; 884 GetNavigatorFor()->LocateGlobalPointWi 849 GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_A); 885 G4bool Intersects_AB = IntersectChord( 850 G4bool Intersects_AB = IntersectChord(Point_A, Point_B, NewSafety, 886 851 previousSafety, 887 852 previousSftOrigin, distAB, 888 853 PointGi); 889 if( Intersects_AB ) 854 if( Intersects_AB ) 890 { 855 { 891 last_AF_intersection = Intersects_AB 856 last_AF_intersection = Intersects_AB; 892 CurrentE_Point = PointGi; 857 CurrentE_Point = PointGi; 893 858 894 validIntersectP = true; // 'E' has 859 validIntersectP = true; // 'E' has been updated. 895 NormalAtEntry = GetSurfaceNormal( P 860 NormalAtEntry = GetSurfaceNormal( PointGi, validNormalAtE ); 896 } 861 } 897 else 862 else 898 { 863 { 899 validIntersectP = false; // No new 864 validIntersectP = false; // No new 'E' chord intersection found 900 if( depth == 1) 865 if( depth == 1) 901 { 866 { 902 there_is_no_intersection = true; 867 there_is_no_intersection = true; 903 } 868 } 904 } 869 } 905 depth--; 870 depth--; 906 fin_section_depth[depth] = true; 871 fin_section_depth[depth] = true; 907 unfinished = !validIntersectP; 872 unfinished = !validIntersectP; 908 } 873 } 909 #ifdef G4DEBUG_FIELD 874 #ifdef G4DEBUG_FIELD 910 if( ! ( validIntersectP || there_is_no_i 875 if( ! ( validIntersectP || there_is_no_intersection ) ) 911 { 876 { 912 // What happened ?? 877 // What happened ?? 913 G4cout << "MLL - WARNING Potential FA 878 G4cout << "MLL - WARNING Potential FAILURE: Conditions not met!" 914 << G4endl 879 << G4endl 915 << " Depth = " << depth << G4e 880 << " Depth = " << depth << G4endl >> 881 << " Levels popped = " << levelPops 916 << " Num Substeps= " << subste 882 << " Num Substeps= " << substep_no << G4endl; 917 G4cout << " Found intersection= " << 883 G4cout << " Found intersection= " << found_approximate_intersection 918 << G4endl; 884 << G4endl; 919 G4cout << " Progress report: -- " << 885 G4cout << " Progress report: -- " << G4endl; 920 ReportProgress(G4cout, 886 ReportProgress(G4cout, 921 CurveStartPointVelocit 887 CurveStartPointVelocity, CurveEndPointVelocity, 922 substep_no, CurrentA_P 888 substep_no, CurrentA_PointVelocity, 923 CurrentB_PointVelocity 889 CurrentB_PointVelocity, 924 NewSafety, depth ); 890 NewSafety, depth ); 925 G4cout << G4endl; 891 G4cout << G4endl; 926 } 892 } 927 #endif 893 #endif 928 } // if(!found_aproximate_intersection) 894 } // if(!found_aproximate_intersection) 929 895 930 assert( validIntersectP || there_is_no_int 896 assert( validIntersectP || there_is_no_intersection 931 || found_approxima 897 || found_approximate_intersection); 932 898 933 } while ( ( ! found_approximate_intersection 899 } while ( ( ! found_approximate_intersection ) 934 && ( ! there_is_no_intersection ) 900 && ( ! there_is_no_intersection ) 935 && ( substep_no <= fMaxSteps) ); / 901 && ( substep_no <= fMaxSteps) ); // UNTIL found or failed 936 902 937 if( substep_no > max_no_seen ) 903 if( substep_no > max_no_seen ) 938 { 904 { 939 max_no_seen = substep_no; 905 max_no_seen = substep_no; 940 #ifdef G4DEBUG_FIELD 906 #ifdef G4DEBUG_FIELD 941 if( max_no_seen > fWarnSteps ) 907 if( max_no_seen > fWarnSteps ) 942 { 908 { 943 trigger_substepno_print = max_no_seen-20 909 trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps 944 } 910 } 945 #endif 911 #endif 946 } 912 } 947 913 948 if( !there_is_no_intersection && !found_appr 914 if( !there_is_no_intersection && !found_approximate_intersection ) 949 { 915 { 950 if( substep_no >= fMaxSteps) 916 if( substep_no >= fMaxSteps) 951 { 917 { 952 // Since we cannot go further (yet), w 918 // Since we cannot go further (yet), we return as far as we have gone 953 919 954 recalculatedEndPoint = true; 920 recalculatedEndPoint = true; 955 IntersectedOrRecalculatedFT = CurrentA 921 IntersectedOrRecalculatedFT = CurrentA_PointVelocity; 956 found_approximate_intersection = false 922 found_approximate_intersection = false; 957 923 958 std::ostringstream message; 924 std::ostringstream message; 959 message << G4endl; 925 message << G4endl; 960 message << "Convergence is requiring t 926 message << "Convergence is requiring too many substeps: " 961 << substep_no << " ( limit = 927 << substep_no << " ( limit = "<< fMaxSteps << ")" 962 << G4endl 928 << G4endl 963 << " Abandoning effort to int 929 << " Abandoning effort to intersect. " << G4endl << G4endl; 964 message << " Number of calls to MLL 930 message << " Number of calls to MLL: " << fNumCalls; 965 message << " iteration = " << substep 931 message << " iteration = " << substep_no <<G4endl << G4endl; 966 932 967 message.precision( 12 ); 933 message.precision( 12 ); 968 G4double done_len = CurrentA_PointVelo 934 G4double done_len = CurrentA_PointVelocity.GetCurveLength(); 969 G4double full_len = CurveEndPointVeloc 935 G4double full_len = CurveEndPointVelocity.GetCurveLength(); 970 message << " Undertaken only le 936 message << " Undertaken only length: " << done_len 971 << " out of " << full_len << " 937 << " out of " << full_len << " required." << G4endl 972 << " Remaining length = 938 << " Remaining length = " << full_len - done_len; 973 939 974 message << " Start and end-point o 940 message << " Start and end-point of requested Step:" << G4endl; 975 printStatus( CurveStartPointVelocity, 941 printStatus( CurveStartPointVelocity, CurveEndPointVelocity, 976 -1.0, NewSafety, 0, m 942 -1.0, NewSafety, 0, message, -1 ); 977 message << " Start and end-poin 943 message << " Start and end-point of current Sub-Step:" << G4endl; 978 printStatus( CurrentA_PointVelocity, C 944 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity, 979 -1.0, NewSafety, substep_ 945 -1.0, NewSafety, substep_no-1, message, -1 ); 980 printStatus( CurrentA_PointVelocity, C 946 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity, 981 -1.0, NewSafety, substep_ 947 -1.0, NewSafety, substep_no, message, -1 ); 982 948 983 G4Exception(MethodName, "GeomNav0003", 949 G4Exception(MethodName, "GeomNav0003", JustWarning, message); 984 } 950 } 985 else if( substep_no >= fWarnSteps) 951 else if( substep_no >= fWarnSteps) 986 { 952 { 987 std::ostringstream message; 953 std::ostringstream message; 988 message << "Many substeps while trying 954 message << "Many substeps while trying to locate intersection." 989 << G4endl 955 << G4endl 990 << " Undertaken lengt 956 << " Undertaken length: " 991 << CurrentB_PointVelocity.GetC 957 << CurrentB_PointVelocity.GetCurveLength() 992 << " - Needed: " << substep_n 958 << " - Needed: " << substep_no << " substeps." << G4endl 993 << " Warning number = 959 << " Warning number = " << fWarnSteps 994 << " and maximum substeps = " 960 << " and maximum substeps = " << fMaxSteps; 995 G4Exception(MethodName, "GeomNav1002", 961 G4Exception(MethodName, "GeomNav1002", JustWarning, message); 996 } 962 } 997 } 963 } 998 964 999 return (!there_is_no_intersection) && found 965 return (!there_is_no_intersection) && found_approximate_intersection; 1000 // Success or failure 966 // Success or failure 1001 } 967 } 1002 968 1003 void G4MultiLevelLocator::ReportStatistics() 969 void G4MultiLevelLocator::ReportStatistics() 1004 { 970 { 1005 G4cout << " Number of calls = " << fNumCal 971 G4cout << " Number of calls = " << fNumCalls << G4endl; 1006 G4cout << " Number of split level ('advanc 972 G4cout << " Number of split level ('advances'): " 1007 << fNumAdvanceTrials << G4endl; 973 << fNumAdvanceTrials << G4endl; 1008 G4cout << " Number of full advances: 974 G4cout << " Number of full advances: " 1009 << fNumAdvanceGood << G4endl; 975 << fNumAdvanceGood << G4endl; 1010 G4cout << " Number of good advances: 976 G4cout << " Number of good advances: " 1011 << fNumAdvanceFull << G4endl; 977 << fNumAdvanceFull << G4endl; 1012 } 978 } 1013 979 1014 void G4MultiLevelLocator::ReportFieldValue( c 980 void G4MultiLevelLocator::ReportFieldValue( const G4FieldTrack& locationPV, 1015 c 981 const char* nameLoc, 1016 c 982 const G4EquationOfMotion* equation ) 1017 { 983 { 1018 enum { maxNumFieldComp = 24 }; 984 enum { maxNumFieldComp = 24 }; 1019 985 1020 G4ThreeVector position = locationPV.GetPos 986 G4ThreeVector position = locationPV.GetPosition(); 1021 G4double startPoint[4] = { position.x(), p 987 G4double startPoint[4] = { position.x(), position.y(), position.z(), 1022 locationPV.GetL 988 locationPV.GetLabTimeOfFlight() }; 1023 G4double FieldVec[maxNumFieldComp]; // 24 989 G4double FieldVec[maxNumFieldComp]; // 24 ; 1024 for (double & i : FieldVec) << 990 for (auto i=0; i<maxNumFieldComp; ++i ) 1025 { 991 { 1026 i = 0.0; << 992 FieldVec[i] = 0.0; 1027 } 993 } 1028 equation->GetFieldValue( startPoint, Field 994 equation->GetFieldValue( startPoint, FieldVec); 1029 G4cout << " B-field value (" << nameLoc < 995 G4cout << " B-field value (" << nameLoc << ")= " 1030 << FieldVec[0] << " " << FieldVec[1 996 << FieldVec[0] << " " << FieldVec[1] << " " << FieldVec[2]; 1031 G4double Emag2= G4ThreeVector( FieldVec[3] 997 G4double Emag2= G4ThreeVector( FieldVec[3], 1032 FieldVec[4] 998 FieldVec[4], 1033 FieldVec[5] 999 FieldVec[5] ).mag2(); 1034 if( Emag2 > 0.0 ) 1000 if( Emag2 > 0.0 ) 1035 { 1001 { 1036 G4cout << " Electric = " << FieldVec[3] 1002 G4cout << " Electric = " << FieldVec[3] << " " 1037 << FieldVec[4] 1003 << FieldVec[4] << " " 1038 << FieldVec[5] 1004 << FieldVec[5]<< G4endl; 1039 } 1005 } 1040 return; 1006 return; 1041 } 1007 } 1042 1008