| 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 G4PropagatorInField Implementation 26 // class G4PropagatorInField Implementation
27 // 27 //
28 // This class implements an algorithm to trac 28 // This class implements an algorithm to track a particle in a
29 // non-uniform magnetic field. It utilises an 29 // non-uniform magnetic field. It utilises an ODE solver (with
30 // the Runge - Kutta method) to evolve the pa 30 // the Runge - Kutta method) to evolve the particle, and drives it
31 // until the particle has traveled a set dist 31 // until the particle has traveled a set distance or it enters a new
32 // volume. 32 // volume.
33 // 33 //
34 // 14.10.96 John Apostolakis, design and imple 34 // 14.10.96 John Apostolakis, design and implementation
35 // 17.03.97 John Apostolakis, renaming new set 35 // 17.03.97 John Apostolakis, renaming new set functions being added
36 // ------------------------------------------- 36 // ---------------------------------------------------------------------------
37 37
38 #include <iomanip> 38 #include <iomanip>
39 39
40 #include "G4PropagatorInField.hh" 40 #include "G4PropagatorInField.hh"
41 #include "G4ios.hh" 41 #include "G4ios.hh"
42 #include "G4SystemOfUnits.hh" 42 #include "G4SystemOfUnits.hh"
43 #include "G4ThreeVector.hh" 43 #include "G4ThreeVector.hh"
44 #include "G4Material.hh" 44 #include "G4Material.hh"
45 #include "G4VPhysicalVolume.hh" 45 #include "G4VPhysicalVolume.hh"
46 #include "G4Navigator.hh" 46 #include "G4Navigator.hh"
47 #include "G4GeometryTolerance.hh" 47 #include "G4GeometryTolerance.hh"
48 #include "G4VCurvedTrajectoryFilter.hh" 48 #include "G4VCurvedTrajectoryFilter.hh"
49 #include "G4ChordFinder.hh" 49 #include "G4ChordFinder.hh"
50 #include "G4MultiLevelLocator.hh" 50 #include "G4MultiLevelLocator.hh"
51 51
52 <<
53 // ------------------------------------------- 52 // ---------------------------------------------------------------------------
54 // Constructors and destructor 53 // Constructors and destructor
55 // 54 //
56 G4PropagatorInField::G4PropagatorInField( G4Na 55 G4PropagatorInField::G4PropagatorInField( G4Navigator* theNavigator,
57 G4Fi 56 G4FieldManager* detectorFieldMgr,
58 G4VI 57 G4VIntersectionLocator* vLocator )
59 : fDetectorFieldMgr(detectorFieldMgr), 58 : fDetectorFieldMgr(detectorFieldMgr),
60 fNavigator(theNavigator), 59 fNavigator(theNavigator),
61 fCurrentFieldMgr(detectorFieldMgr), 60 fCurrentFieldMgr(detectorFieldMgr),
62 End_PointAndTangent(G4ThreeVector(0.,0.,0. 61 End_PointAndTangent(G4ThreeVector(0.,0.,0.),
63 G4ThreeVector(0.,0.,0. 62 G4ThreeVector(0.,0.,0.),0.0,0.0,0.0,0.0,0.0)
64 { 63 {
65 fEpsilonStep = (fDetectorFieldMgr != nullptr 64 fEpsilonStep = (fDetectorFieldMgr != nullptr)
66 ? fDetectorFieldMgr->GetMaximum 65 ? fDetectorFieldMgr->GetMaximumEpsilonStep() : 1.0e-5;
67 << 66 fLargestAcceptableStep = 1000.0 * meter;
68 67
69 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) 68 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.);
70 kCarTolerance = G4GeometryTolerance::GetInst 69 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
71 fZeroStepThreshold = std::max( 1.0e5 * kCarT 70 fZeroStepThreshold = std::max( 1.0e5 * kCarTolerance, 1.0e-1 * micrometer );
72 71
73 fLargestAcceptableStep = 100.0 * meter; // <<
74 fMaxStepSizeMultiplier= 0.1 ; // 0.1 in <<
75 fMinBigDistance= 100. * CLHEP::mm; <<
76 #ifdef G4DEBUG_FIELD 72 #ifdef G4DEBUG_FIELD
77 G4cout << " PiF: Zero Step Threshold set to 73 G4cout << " PiF: Zero Step Threshold set to "
78 << fZeroStepThreshold / millimeter 74 << fZeroStepThreshold / millimeter
79 << " mm." << G4endl; 75 << " mm." << G4endl;
80 G4cout << " PiF: Value of kCarTolerance = 76 G4cout << " PiF: Value of kCarTolerance = "
81 << kCarTolerance / millimeter 77 << kCarTolerance / millimeter
82 << " mm. " << G4endl; 78 << " mm. " << G4endl;
83 fVerboseLevel = 2; 79 fVerboseLevel = 2;
84 fVerbTracePiF = true; 80 fVerbTracePiF = true;
85 #endif 81 #endif
86 82
87 // Defining Intersection Locator and his par 83 // Defining Intersection Locator and his parameters
88 if ( vLocator == nullptr ) 84 if ( vLocator == nullptr )
89 { 85 {
90 fIntersectionLocator = new G4MultiLevelLoc 86 fIntersectionLocator = new G4MultiLevelLocator(theNavigator);
91 fAllocatedLocator = true; 87 fAllocatedLocator = true;
92 } 88 }
93 else 89 else
94 { 90 {
95 fIntersectionLocator = vLocator; 91 fIntersectionLocator = vLocator;
96 fAllocatedLocator = false; 92 fAllocatedLocator = false;
97 } 93 }
98 RefreshIntersectionLocator(); // Copy all 94 RefreshIntersectionLocator(); // Copy all relevant parameters
99 } 95 }
100 96
101 // ------------------------------------------- 97 // ---------------------------------------------------------------------------
102 // 98 //
103 G4PropagatorInField::~G4PropagatorInField() 99 G4PropagatorInField::~G4PropagatorInField()
104 { 100 {
105 if(fAllocatedLocator) { delete fIntersecti 101 if(fAllocatedLocator) { delete fIntersectionLocator; }
106 } 102 }
107 103
108 // ------------------------------------------- 104 // ---------------------------------------------------------------------------
109 // Update the IntersectionLocator with current 105 // Update the IntersectionLocator with current parameters
110 // 106 //
111 void G4PropagatorInField::RefreshIntersectionL 107 void G4PropagatorInField::RefreshIntersectionLocator()
112 { 108 {
113 fIntersectionLocator->SetEpsilonStepFor(fEps 109 fIntersectionLocator->SetEpsilonStepFor(fEpsilonStep);
114 fIntersectionLocator->SetDeltaIntersectionFo 110 fIntersectionLocator->SetDeltaIntersectionFor(fCurrentFieldMgr->GetDeltaIntersection());
115 fIntersectionLocator->SetChordFinderFor(GetC 111 fIntersectionLocator->SetChordFinderFor(GetChordFinder());
116 fIntersectionLocator->SetSafetyParametersFor 112 fIntersectionLocator->SetSafetyParametersFor( fUseSafetyForOptimisation);
117 } 113 }
118 114
119 // ------------------------------------------- 115 // ---------------------------------------------------------------------------
120 // Compute the next geometric Step 116 // Compute the next geometric Step
121 // 117 //
122 G4double G4PropagatorInField::ComputeStep( 118 G4double G4PropagatorInField::ComputeStep(
123 G4FieldTrack& pFieldTrack 119 G4FieldTrack& pFieldTrack,
124 G4double CurrentProp 120 G4double CurrentProposedStepLength,
125 G4double& currentSafe 121 G4double& currentSafety, // IN/OUT
126 G4VPhysicalVolume* pPhysVol, 122 G4VPhysicalVolume* pPhysVol,
127 G4bool canRelaxDel 123 G4bool canRelaxDeltaChord)
128 { 124 {
129 GetChordFinder()->OnComputeStep(&pFieldTrack << 125 GetChordFinder()->OnComputeStep();
130 const G4double deltaChord = GetChordFinder() 126 const G4double deltaChord = GetChordFinder()->GetDeltaChord();
131 127
132 // If CurrentProposedStepLength is too small 128 // If CurrentProposedStepLength is too small for finding Chords
133 // then return with no action (for now - TOD 129 // then return with no action (for now - TODO: some action)
134 // 130 //
135 const char* methodName = "G4PropagatorInFiel 131 const char* methodName = "G4PropagatorInField::ComputeStep";
136 if (CurrentProposedStepLength<kCarTolerance) 132 if (CurrentProposedStepLength<kCarTolerance)
137 { 133 {
138 return kInfinity; 134 return kInfinity;
139 } 135 }
140 136
141 // Introducing smooth trajectory display (ja 137 // Introducing smooth trajectory display (jacek 01/11/2002)
142 // 138 //
143 if (fpTrajectoryFilter != nullptr) << 139 if (fpTrajectoryFilter)
144 { 140 {
145 fpTrajectoryFilter->CreateNewTrajectorySeg 141 fpTrajectoryFilter->CreateNewTrajectorySegment();
146 } 142 }
147 143
148 fFirstStepInVolume = fNewTrack ? true : fLas 144 fFirstStepInVolume = fNewTrack ? true : fLastStepInVolume;
149 fLastStepInVolume = false; 145 fLastStepInVolume = false;
150 fNewTrack = false; 146 fNewTrack = false;
151 147
152 if( fVerboseLevel > 2 ) 148 if( fVerboseLevel > 2 )
153 { 149 {
154 G4cout << methodName << " called" << G4end 150 G4cout << methodName << " called" << G4endl;
155 G4cout << " Starting FT: " << pFieldTrac 151 G4cout << " Starting FT: " << pFieldTrack;
156 G4cout << " Requested length = " << Curr 152 G4cout << " Requested length = " << CurrentProposedStepLength << G4endl;
157 G4cout << " PhysVol = "; 153 G4cout << " PhysVol = ";
158 if( pPhysVol != nullptr ) << 154 if( pPhysVol )
159 { <<
160 G4cout << pPhysVol->GetName() << G4endl 155 G4cout << pPhysVol->GetName() << G4endl;
161 } <<
162 else 156 else
163 { <<
164 G4cout << " N/A "; 157 G4cout << " N/A ";
165 } <<
166 G4cout << G4endl; 158 G4cout << G4endl;
167 } 159 }
168 160
169 // Parameters for adaptive Runge-Kutta integ 161 // Parameters for adaptive Runge-Kutta integration
170 162
171 G4double h_TrialStepSize; // 1st Step 163 G4double h_TrialStepSize; // 1st Step Size
172 G4double TruePathLength = CurrentProposedSte 164 G4double TruePathLength = CurrentProposedStepLength;
173 G4double StepTaken = 0.0; 165 G4double StepTaken = 0.0;
174 G4double s_length_taken, epsilon; 166 G4double s_length_taken, epsilon;
175 G4bool intersects; 167 G4bool intersects;
176 G4bool first_substep = true; 168 G4bool first_substep = true;
177 169
178 G4double NewSafety; 170 G4double NewSafety;
179 fParticleIsLooping = false; 171 fParticleIsLooping = false;
180 172
181 // If not yet done, 173 // If not yet done,
182 // Set the field manager to the local one 174 // Set the field manager to the local one if the volume has one,
183 // or to the global one 175 // or to the global one if not
184 // 176 //
185 if( !fSetFieldMgr ) 177 if( !fSetFieldMgr )
186 { 178 {
187 fCurrentFieldMgr = FindAndSetFieldManager( 179 fCurrentFieldMgr = FindAndSetFieldManager( pPhysVol );
188 } 180 }
189 fSetFieldMgr = false; // For next call, the 181 fSetFieldMgr = false; // For next call, the field manager must be set again
190 182
191 G4FieldTrack CurrentState(pFieldTrack); 183 G4FieldTrack CurrentState(pFieldTrack);
192 G4FieldTrack OriginalState = CurrentState; 184 G4FieldTrack OriginalState = CurrentState;
193 185
194 // If the Step length is "infinite", then an 186 // If the Step length is "infinite", then an approximate-maximum Step
195 // length (used to calculate the relative ac 187 // length (used to calculate the relative accuracy) must be guessed
196 // 188 //
197 if( CurrentProposedStepLength >= fLargestAcc 189 if( CurrentProposedStepLength >= fLargestAcceptableStep )
198 { 190 {
199 G4ThreeVector StartPointA, VelocityUnit; 191 G4ThreeVector StartPointA, VelocityUnit;
200 StartPointA = pFieldTrack.GetPosition(); 192 StartPointA = pFieldTrack.GetPosition();
201 VelocityUnit = pFieldTrack.GetMomentumDir( 193 VelocityUnit = pFieldTrack.GetMomentumDir();
202 194
203 G4double trialProposedStep = fMaxStepSizeM << 195 G4double trialProposedStep = 1.e2 * ( 10.0 * cm +
204 fNavigator->GetWorldVolume()->GetLogical 196 fNavigator->GetWorldVolume()->GetLogicalVolume()->
205 GetSolid()->DistanceToOut(St 197 GetSolid()->DistanceToOut(StartPointA, VelocityUnit) );
206 CurrentProposedStepLength = std::min( tria 198 CurrentProposedStepLength = std::min( trialProposedStep,
207 fLar 199 fLargestAcceptableStep );
208 } 200 }
209 epsilon = fCurrentFieldMgr->GetDeltaOneStep( 201 epsilon = fCurrentFieldMgr->GetDeltaOneStep() / CurrentProposedStepLength;
210 G4double epsilonMin= fCurrentFieldMgr->GetMi 202 G4double epsilonMin= fCurrentFieldMgr->GetMinimumEpsilonStep();
211 G4double epsilonMax= fCurrentFieldMgr->GetMa 203 G4double epsilonMax= fCurrentFieldMgr->GetMaximumEpsilonStep();
212 if( epsilon < epsilonMin ) { epsilon = epsi 204 if( epsilon < epsilonMin ) { epsilon = epsilonMin; }
213 if( epsilon > epsilonMax ) { epsilon = epsi 205 if( epsilon > epsilonMax ) { epsilon = epsilonMax; }
214 SetEpsilonStep( epsilon ); 206 SetEpsilonStep( epsilon );
215 207
216 // Values for Intersection Locator has to be 208 // Values for Intersection Locator has to be updated on each call for the
217 // case that CurrentFieldManager has changed 209 // case that CurrentFieldManager has changed from the one of previous step
218 // 210 //
219 RefreshIntersectionLocator(); 211 RefreshIntersectionLocator();
220 212
221 // Shorten the proposed step in case of earl 213 // Shorten the proposed step in case of earlier problems (zero steps)
222 // 214 //
223 if( fNoZeroStep > fActionThreshold_NoZeroSte 215 if( fNoZeroStep > fActionThreshold_NoZeroSteps )
224 { 216 {
225 G4double stepTrial; 217 G4double stepTrial;
226 218
227 stepTrial = fFull_CurveLen_of_LastAttempt; 219 stepTrial = fFull_CurveLen_of_LastAttempt;
228 if( (stepTrial <= 0.0) && (fLast_ProposedS 220 if( (stepTrial <= 0.0) && (fLast_ProposedStepLength > 0.0) )
229 { 221 {
230 stepTrial = fLast_ProposedStepLength; 222 stepTrial = fLast_ProposedStepLength;
231 } 223 }
232 224
233 G4double decreaseFactor = 0.9; // Unused d 225 G4double decreaseFactor = 0.9; // Unused default
234 if( (fNoZeroStep < fSevereActionThreshol 226 if( (fNoZeroStep < fSevereActionThreshold_NoZeroSteps)
235 && (stepTrial > 100.0*fZeroStepThreshol 227 && (stepTrial > 100.0*fZeroStepThreshold) )
236 { 228 {
237 // Attempt quick convergence 229 // Attempt quick convergence
238 // 230 //
239 decreaseFactor= 0.25; 231 decreaseFactor= 0.25;
240 } 232 }
241 else 233 else
242 { 234 {
243 // We are in significant difficulties, p 235 // We are in significant difficulties, probably at a boundary that
244 // is either geometrically sharp or betw 236 // is either geometrically sharp or between very different materials.
245 // Careful decreases to cope with tolera 237 // Careful decreases to cope with tolerance are required
246 // 238 //
247 if( stepTrial > 100.0*fZeroStepThreshold << 239 if( stepTrial > 100.0*fZeroStepThreshold )
248 decreaseFactor = 0.35; // Try decr 240 decreaseFactor = 0.35; // Try decreasing slower
249 } else if( stepTrial > 30.0*fZeroStepThr << 241 else if( stepTrial > 30.0*fZeroStepThreshold )
250 decreaseFactor= 0.5; // Try yet 242 decreaseFactor= 0.5; // Try yet slower decrease
251 } else if( stepTrial > 10.0*fZeroStepThr << 243 else if( stepTrial > 10.0*fZeroStepThreshold )
252 decreaseFactor= 0.75; // Try even 244 decreaseFactor= 0.75; // Try even slower decreases
253 } else { << 245 else
254 decreaseFactor= 0.9; // Try very 246 decreaseFactor= 0.9; // Try very slow decreases
255 } <<
256 } 247 }
257 stepTrial *= decreaseFactor; 248 stepTrial *= decreaseFactor;
258 249
259 #ifdef G4DEBUG_FIELD 250 #ifdef G4DEBUG_FIELD
260 if( fVerboseLevel > 2 251 if( fVerboseLevel > 2
261 || (fNoZeroStep >= fSevereActionThreshol 252 || (fNoZeroStep >= fSevereActionThreshold_NoZeroSteps) )
262 { 253 {
263 G4cerr << " " << methodName 254 G4cerr << " " << methodName
264 << " Decreasing step after " < 255 << " Decreasing step after " << fNoZeroStep << " zero steps "
265 << " - in volume " << pPhysVol; 256 << " - in volume " << pPhysVol;
266 if( pPhysVol ) 257 if( pPhysVol )
267 G4cerr << " with name " << pPhysVol 258 G4cerr << " with name " << pPhysVol->GetName();
268 else 259 else
269 G4cerr << " i.e. *unknown* volume." 260 G4cerr << " i.e. *unknown* volume.";
270 G4cerr << G4endl; 261 G4cerr << G4endl;
271 PrintStepLengthDiagnostic(CurrentPropo 262 PrintStepLengthDiagnostic(CurrentProposedStepLength, decreaseFactor,
272 stepTrial, p 263 stepTrial, pFieldTrack);
273 } 264 }
274 #endif 265 #endif
275 if( stepTrial == 0.0 ) // Change to mak 266 if( stepTrial == 0.0 ) // Change to make it < 0.1 * kCarTolerance ??
276 { 267 {
277 std::ostringstream message; 268 std::ostringstream message;
278 message << "Particle abandoned due to l 269 message << "Particle abandoned due to lack of progress in field."
279 << G4endl 270 << G4endl
280 << " Properties : " << pFieldT 271 << " Properties : " << pFieldTrack << G4endl
281 << " Attempting a zero step = 272 << " Attempting a zero step = " << stepTrial << G4endl
282 << " while attempting to progr 273 << " while attempting to progress after " << fNoZeroStep
283 << " trial steps. Will abandon 274 << " trial steps. Will abandon step.";
284 G4Exception(methodName, "GeomNav1002", 275 G4Exception(methodName, "GeomNav1002", JustWarning, message);
285 fParticleIsLooping = true; 276 fParticleIsLooping = true;
286 return 0; // = stepTrial; 277 return 0; // = stepTrial;
287 } 278 }
288 if( stepTrial < CurrentProposedStepLength 279 if( stepTrial < CurrentProposedStepLength )
289 { 280 {
290 CurrentProposedStepLength = stepTrial; 281 CurrentProposedStepLength = stepTrial;
291 } 282 }
292 } 283 }
293 fLast_ProposedStepLength = CurrentProposedSt 284 fLast_ProposedStepLength = CurrentProposedStepLength;
294 285
295 G4int do_loop_count = 0; 286 G4int do_loop_count = 0;
296 do // Loop checking, 07.10.2016, JA 287 do // Loop checking, 07.10.2016, JA
297 { 288 {
298 G4FieldTrack SubStepStartState = CurrentSt 289 G4FieldTrack SubStepStartState = CurrentState;
299 G4ThreeVector SubStartPoint = CurrentState 290 G4ThreeVector SubStartPoint = CurrentState.GetPosition();
300 291
301 if(!first_substep) 292 if(!first_substep)
302 { 293 {
303 if( fVerboseLevel > 4 ) 294 if( fVerboseLevel > 4 )
304 { 295 {
305 G4cout << " PiF: Calling Nav/Locate Gl 296 G4cout << " PiF: Calling Nav/Locate Global Point within-Volume "
306 << G4endl; 297 << G4endl;
307 } 298 }
308 fNavigator->LocateGlobalPointWithinVolum 299 fNavigator->LocateGlobalPointWithinVolume( SubStartPoint );
309 } 300 }
310 301
311 // How far to attempt to move the particle 302 // How far to attempt to move the particle !
312 // 303 //
313 h_TrialStepSize = CurrentProposedStepLengt 304 h_TrialStepSize = CurrentProposedStepLength - StepTaken;
314 305
315 if (canRelaxDeltaChord && 306 if (canRelaxDeltaChord &&
316 fIncreaseChordDistanceThreshold > 0 & 307 fIncreaseChordDistanceThreshold > 0 &&
317 do_loop_count > fIncreaseChordDistance 308 do_loop_count > fIncreaseChordDistanceThreshold &&
318 do_loop_count % fIncreaseChordDistance 309 do_loop_count % fIncreaseChordDistanceThreshold == 0)
319 { 310 {
320 GetChordFinder()->SetDeltaChord( 311 GetChordFinder()->SetDeltaChord(
321 GetChordFinder()->GetDeltaChord() * 312 GetChordFinder()->GetDeltaChord() * 2.0
322 ); 313 );
323 } 314 }
324 315
325 // Integrate as far as "chord miss" rule a 316 // Integrate as far as "chord miss" rule allows.
326 // 317 //
327 s_length_taken = GetChordFinder()->Advance 318 s_length_taken = GetChordFinder()->AdvanceChordLimited(
328 CurrentState, 319 CurrentState, // Position & velocity
329 h_TrialStepSize, 320 h_TrialStepSize,
330 fEpsilonStep, 321 fEpsilonStep,
331 fPreviousSftOrigi 322 fPreviousSftOrigin,
332 fPreviousSafety ) 323 fPreviousSafety );
333 // CurrentState is now updated with the 324 // CurrentState is now updated with the final position and velocity
334 325
335 fFull_CurveLen_of_LastAttempt = s_length_t 326 fFull_CurveLen_of_LastAttempt = s_length_taken;
336 327
337 G4ThreeVector EndPointB = CurrentState.Get 328 G4ThreeVector EndPointB = CurrentState.GetPosition();
338 G4ThreeVector InterSectionPointE; 329 G4ThreeVector InterSectionPointE;
339 G4double LinearStepLength; 330 G4double LinearStepLength;
340 331
341 // Intersect chord AB with geometry 332 // Intersect chord AB with geometry
342 // 333 //
343 intersects= IntersectChord( SubStartPoint, 334 intersects= IntersectChord( SubStartPoint, EndPointB,
344 NewSafety, Lin 335 NewSafety, LinearStepLength,
345 InterSectionPo 336 InterSectionPointE );
346 // E <- Intersection Point of chord AB a 337 // E <- Intersection Point of chord AB and either volume A's surface
347 // or a 338 // or a daughter volume's surface ..
348 339
349 if( first_substep ) 340 if( first_substep )
350 { 341 {
351 currentSafety = NewSafety; 342 currentSafety = NewSafety;
352 } // Updating safety in other steps is pot 343 } // Updating safety in other steps is potential future extention
353 344
354 if( intersects ) 345 if( intersects )
355 { 346 {
356 G4FieldTrack IntersectPointVelct_G(Curr 347 G4FieldTrack IntersectPointVelct_G(CurrentState); // FT-Def-Construct
357 348
358 // Find the intersection point of AB tr 349 // Find the intersection point of AB true path with the surface
359 // of vol(A), if it exists. Start wit 350 // of vol(A), if it exists. Start with point E as first "estimate".
360 G4bool recalculatedEndPt = false; 351 G4bool recalculatedEndPt = false;
361 352
362 G4bool found_intersection = fIntersecti 353 G4bool found_intersection = fIntersectionLocator->
363 EstimateIntersectionPoint( SubStepSta 354 EstimateIntersectionPoint( SubStepStartState, CurrentState,
364 InterSecti 355 InterSectionPointE, IntersectPointVelct_G,
365 recalculat 356 recalculatedEndPt, fPreviousSafety,
366 fPreviousS 357 fPreviousSftOrigin);
367 intersects = found_intersection; 358 intersects = found_intersection;
368 if( found_intersection ) 359 if( found_intersection )
369 { 360 {
370 End_PointAndTangent= IntersectPointV 361 End_PointAndTangent= IntersectPointVelct_G; // G is our EndPoint ...
371 StepTaken = TruePathLength = Interse 362 StepTaken = TruePathLength = IntersectPointVelct_G.GetCurveLength()
372 - Origina 363 - OriginalState.GetCurveLength();
373 } 364 }
374 else 365 else
375 { 366 {
376 // Either "minor" chords do not inte 367 // Either "minor" chords do not intersect
377 // or else stopped (due to too many 368 // or else stopped (due to too many steps)
378 // 369 //
379 if( recalculatedEndPt ) 370 if( recalculatedEndPt )
380 { 371 {
381 G4double endAchieved = IntersectP 372 G4double endAchieved = IntersectPointVelct_G.GetCurveLength();
382 G4double endExpected = CurrentSta 373 G4double endExpected = CurrentState.GetCurveLength();
383 374
384 // Detect failure - due to too ma 375 // Detect failure - due to too many steps
385 G4bool shortEnd = endAchieved 376 G4bool shortEnd = endAchieved
386 < (endExpected*(1 377 < (endExpected*(1.0-CLHEP::perMillion));
387 378
388 G4double stepAchieved = endAchiev 379 G4double stepAchieved = endAchieved
389 - SubStepSt 380 - SubStepStartState.GetCurveLength();
390 381
391 // Update remaining state - must 382 // Update remaining state - must work for 'full' step or
392 // abandonned intersection 383 // abandonned intersection
393 // 384 //
394 CurrentState = IntersectPointVelc 385 CurrentState = IntersectPointVelct_G;
395 s_length_taken = stepAchieved; 386 s_length_taken = stepAchieved;
396 if( shortEnd ) 387 if( shortEnd )
397 { 388 {
398 fParticleIsLooping = true; 389 fParticleIsLooping = true;
399 } 390 }
400 } 391 }
401 } 392 }
402 } 393 }
403 if( !intersects ) 394 if( !intersects )
404 { 395 {
405 StepTaken += s_length_taken; 396 StepTaken += s_length_taken;
406 397
407 if (fpTrajectoryFilter != nullptr) // Fo << 398 if (fpTrajectoryFilter) // For smooth trajectory display (jacek 1/11/2002)
408 { 399 {
409 fpTrajectoryFilter->TakeIntermediatePo 400 fpTrajectoryFilter->TakeIntermediatePoint(CurrentState.GetPosition());
410 } 401 }
411 } 402 }
412 first_substep = false; 403 first_substep = false;
413 404
414 #ifdef G4DEBUG_FIELD 405 #ifdef G4DEBUG_FIELD
415 if( fNoZeroStep > fActionThreshold_NoZeroS 406 if( fNoZeroStep > fActionThreshold_NoZeroSteps )
416 { 407 {
417 if( fNoZeroStep > fSevereActionThreshold 408 if( fNoZeroStep > fSevereActionThreshold_NoZeroSteps )
418 G4cout << " Above 'Severe Action' thre 409 G4cout << " Above 'Severe Action' threshold -- for Zero steps. ";
419 else 410 else
420 G4cout << " Above 'action' threshold - 411 G4cout << " Above 'action' threshold -- for Zero steps. ";
421 G4cout << " Number of zero steps = " << 412 G4cout << " Number of zero steps = " << fNoZeroStep << G4endl;
422 printStatus( SubStepStartState, // or O 413 printStatus( SubStepStartState, // or OriginalState,
423 CurrentState, CurrentPropos 414 CurrentState, CurrentProposedStepLength,
424 NewSafety, do_loop_count, p 415 NewSafety, do_loop_count, pPhysVol );
425 } 416 }
426 if( (fVerboseLevel > 1) && (do_loop_count 417 if( (fVerboseLevel > 1) && (do_loop_count > fMax_loop_count-10 ))
427 { 418 {
428 if( do_loop_count == fMax_loop_count-9 ) 419 if( do_loop_count == fMax_loop_count-9 )
429 { 420 {
430 G4cout << " G4PropagatorInField::Compu 421 G4cout << " G4PropagatorInField::ComputeStep(): " << G4endl
431 << " Difficult track - taking 422 << " Difficult track - taking many sub steps." << G4endl;
432 printStatus( SubStepStartState, SubSte 423 printStatus( SubStepStartState, SubStepStartState, CurrentProposedStepLength,
433 NewSafety, 0, pPhysVol ); 424 NewSafety, 0, pPhysVol );
434 } 425 }
435 printStatus( SubStepStartState, CurrentS 426 printStatus( SubStepStartState, CurrentState, CurrentProposedStepLength,
436 NewSafety, do_loop_count, p 427 NewSafety, do_loop_count, pPhysVol );
437 } 428 }
438 #endif 429 #endif
439 430
440 ++do_loop_count; 431 ++do_loop_count;
441 432
442 } while( (!intersects ) 433 } while( (!intersects )
443 && (!fParticleIsLooping) 434 && (!fParticleIsLooping)
444 && (StepTaken + kCarTolerance < Curren 435 && (StepTaken + kCarTolerance < CurrentProposedStepLength)
445 && ( do_loop_count < fMax_loop_count ) 436 && ( do_loop_count < fMax_loop_count ) );
446 437
447 if( do_loop_count >= fMax_loop_count 438 if( do_loop_count >= fMax_loop_count
448 && (StepTaken + kCarTolerance < CurrentPro 439 && (StepTaken + kCarTolerance < CurrentProposedStepLength) )
449 { 440 {
450 fParticleIsLooping = true; 441 fParticleIsLooping = true;
451 } 442 }
452 if ( ( fParticleIsLooping ) && (fVerboseLeve 443 if ( ( fParticleIsLooping ) && (fVerboseLevel > 0) )
453 { 444 {
454 ReportLoopingParticle( do_loop_count, Step 445 ReportLoopingParticle( do_loop_count, StepTaken,
455 CurrentProposedStep 446 CurrentProposedStepLength, methodName,
456 CurrentState.GetMom 447 CurrentState.GetMomentum(), pPhysVol );
457 } 448 }
458 449
459 if( !intersects ) 450 if( !intersects )
460 { 451 {
461 // Chord AB or "minor chords" do not inter 452 // Chord AB or "minor chords" do not intersect
462 // B is the endpoint Step of the current S 453 // B is the endpoint Step of the current Step.
463 // 454 //
464 End_PointAndTangent = CurrentState; 455 End_PointAndTangent = CurrentState;
465 TruePathLength = StepTaken; // Original 456 TruePathLength = StepTaken; // Original code
466 457
467 // Tried the following to avoid potential 458 // Tried the following to avoid potential issue with round-off error
468 // - but has issues... Suppressing this ch 459 // - but has issues... Suppressing this change JA 2015/05/02
469 // TruePathLength = CurrentProposedStepLen 460 // TruePathLength = CurrentProposedStepLength;
470 } 461 }
471 fLastStepInVolume = intersects; 462 fLastStepInVolume = intersects;
472 463
473 // Set pFieldTrack to the return value 464 // Set pFieldTrack to the return value
474 // 465 //
475 pFieldTrack = End_PointAndTangent; 466 pFieldTrack = End_PointAndTangent;
476 467
477 #ifdef G4VERBOSE 468 #ifdef G4VERBOSE
478 // Check that "s" is correct 469 // Check that "s" is correct
479 // 470 //
480 if( std::fabs(OriginalState.GetCurveLength() 471 if( std::fabs(OriginalState.GetCurveLength() + TruePathLength
481 - End_PointAndTangent.GetCurveLength()) 472 - End_PointAndTangent.GetCurveLength()) > 3.e-4 * TruePathLength )
482 { 473 {
483 std::ostringstream message; 474 std::ostringstream message;
484 message << "Curve length mis-match between 475 message << "Curve length mis-match between original state "
485 << "and proposed endpoint of propa 476 << "and proposed endpoint of propagation." << G4endl
486 << " The curve length of the endp 477 << " The curve length of the endpoint should be: "
487 << OriginalState.GetCurveLength() 478 << OriginalState.GetCurveLength() + TruePathLength << G4endl
488 << " and it is instead: " 479 << " and it is instead: "
489 << End_PointAndTangent.GetCurveLen 480 << End_PointAndTangent.GetCurveLength() << "." << G4endl
490 << " A difference of: " 481 << " A difference of: "
491 << OriginalState.GetCurveLength() 482 << OriginalState.GetCurveLength() + TruePathLength
492 - End_PointAndTangent.GetCurveL 483 - End_PointAndTangent.GetCurveLength() << G4endl
493 << " Original state = " << Origin 484 << " Original state = " << OriginalState << G4endl
494 << " Proposed state = " << End_Po 485 << " Proposed state = " << End_PointAndTangent;
495 G4Exception(methodName, "GeomNav0003", Fat 486 G4Exception(methodName, "GeomNav0003", FatalException, message);
496 } 487 }
497 #endif 488 #endif
498 489
499 if( TruePathLength+kCarTolerance >= CurrentP 490 if( TruePathLength+kCarTolerance >= CurrentProposedStepLength )
500 { 491 {
501 fNoZeroStep = 0; 492 fNoZeroStep = 0;
502 } 493 }
503 else 494 else
504 { 495 {
505 // In particular anomalous cases, we can 496 // In particular anomalous cases, we can get repeated zero steps
506 // We identify these cases and take corre 497 // We identify these cases and take corrective action when they occur.
507 // 498 //
508 if( TruePathLength < std::max( fZeroStepT 499 if( TruePathLength < std::max( fZeroStepThreshold, 0.5*kCarTolerance ) )
509 { 500 {
510 ++fNoZeroStep; 501 ++fNoZeroStep;
511 } 502 }
512 else 503 else
513 { 504 {
514 fNoZeroStep = 0; 505 fNoZeroStep = 0;
515 } 506 }
516 } 507 }
517 if( fNoZeroStep > fAbandonThreshold_NoZeroSt 508 if( fNoZeroStep > fAbandonThreshold_NoZeroSteps )
518 { 509 {
519 fParticleIsLooping = true; 510 fParticleIsLooping = true;
520 ReportStuckParticle( fNoZeroStep, Current 511 ReportStuckParticle( fNoZeroStep, CurrentProposedStepLength,
521 fFull_CurveLen_of_La 512 fFull_CurveLen_of_LastAttempt, pPhysVol );
522 fNoZeroStep = 0; 513 fNoZeroStep = 0;
523 } 514 }
524 515
525 GetChordFinder()->SetDeltaChord(deltaChord); 516 GetChordFinder()->SetDeltaChord(deltaChord);
526 return TruePathLength; 517 return TruePathLength;
527 } 518 }
528 519
529 // ------------------------------------------- 520 // ---------------------------------------------------------------------------
530 // Dumps status of propagator 521 // Dumps status of propagator
531 // 522 //
532 void 523 void
533 G4PropagatorInField::printStatus( const G4Fiel 524 G4PropagatorInField::printStatus( const G4FieldTrack& StartFT,
534 const G4Fiel 525 const G4FieldTrack& CurrentFT,
535 G4doub 526 G4double requestStep,
536 G4doub 527 G4double safety,
537 G4int 528 G4int stepNo,
538 G4VPhy 529 G4VPhysicalVolume* startVolume)
539 { 530 {
540 const G4int verboseLevel = fVerboseLevel; 531 const G4int verboseLevel = fVerboseLevel;
541 const G4ThreeVector StartPosition = St 532 const G4ThreeVector StartPosition = StartFT.GetPosition();
542 const G4ThreeVector StartUnitVelocity = St 533 const G4ThreeVector StartUnitVelocity = StartFT.GetMomentumDir();
543 const G4ThreeVector CurrentPosition = Cu 534 const G4ThreeVector CurrentPosition = CurrentFT.GetPosition();
544 const G4ThreeVector CurrentUnitVelocity = Cu 535 const G4ThreeVector CurrentUnitVelocity = CurrentFT.GetMomentumDir();
545 536
546 G4double step_len = CurrentFT.GetCurveLength 537 G4double step_len = CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
547 538
548 G4long oldprec; // cout/cerr precision set << 539 G4int oldprec; // cout/cerr precision settings
549 540
550 if( ((stepNo == 0) && (verboseLevel <3)) || 541 if( ((stepNo == 0) && (verboseLevel <3)) || (verboseLevel >= 3) )
551 { 542 {
552 oldprec = G4cout.precision(4); 543 oldprec = G4cout.precision(4);
553 G4cout << std::setw( 5) << "Step#" 544 G4cout << std::setw( 5) << "Step#"
554 << std::setw(10) << " s " << " " 545 << std::setw(10) << " s " << " "
555 << std::setw(10) << "X(mm)" << " " 546 << std::setw(10) << "X(mm)" << " "
556 << std::setw(10) << "Y(mm)" << " " 547 << std::setw(10) << "Y(mm)" << " "
557 << std::setw(10) << "Z(mm)" << " " 548 << std::setw(10) << "Z(mm)" << " "
558 << std::setw( 7) << " N_x " << " " 549 << std::setw( 7) << " N_x " << " "
559 << std::setw( 7) << " N_y " << " " 550 << std::setw( 7) << " N_y " << " "
560 << std::setw( 7) << " N_z " << " " 551 << std::setw( 7) << " N_z " << " " ;
561 G4cout << std::setw( 7) << " Delta|N|" << 552 G4cout << std::setw( 7) << " Delta|N|" << " "
562 << std::setw( 9) << "StepLen" << " 553 << std::setw( 9) << "StepLen" << " "
563 << std::setw(12) << "StartSafety" < 554 << std::setw(12) << "StartSafety" << " "
564 << std::setw( 9) << "PhsStep" << " 555 << std::setw( 9) << "PhsStep" << " ";
565 if( startVolume != nullptr ) 556 if( startVolume != nullptr )
566 { G4cout << std::setw(18) << "NextVolume 557 { G4cout << std::setw(18) << "NextVolume" << " "; }
567 G4cout.precision(oldprec); 558 G4cout.precision(oldprec);
568 G4cout << G4endl; 559 G4cout << G4endl;
569 } 560 }
570 if((stepNo == 0) && (verboseLevel <=3)) 561 if((stepNo == 0) && (verboseLevel <=3))
571 { 562 {
572 // Recurse to print the start values 563 // Recurse to print the start values
573 // 564 //
574 printStatus( StartFT, StartFT, -1.0, safet 565 printStatus( StartFT, StartFT, -1.0, safety, -1, startVolume);
575 } 566 }
576 if( verboseLevel <= 3 ) 567 if( verboseLevel <= 3 )
577 { 568 {
578 if( stepNo >= 0) 569 if( stepNo >= 0)
579 { G4cout << std::setw( 4) << stepNo << " 570 { G4cout << std::setw( 4) << stepNo << " "; }
580 else 571 else
581 { G4cout << std::setw( 5) << "Start" ; } 572 { G4cout << std::setw( 5) << "Start" ; }
582 oldprec = G4cout.precision(8); 573 oldprec = G4cout.precision(8);
583 G4cout << std::setw(10) << CurrentFT.GetCu 574 G4cout << std::setw(10) << CurrentFT.GetCurveLength() << " ";
584 G4cout.precision(8); 575 G4cout.precision(8);
585 G4cout << std::setw(10) << CurrentPosition 576 G4cout << std::setw(10) << CurrentPosition.x() << " "
586 << std::setw(10) << CurrentPosition 577 << std::setw(10) << CurrentPosition.y() << " "
587 << std::setw(10) << CurrentPosition 578 << std::setw(10) << CurrentPosition.z() << " ";
588 G4cout.precision(4); 579 G4cout.precision(4);
589 G4cout << std::setw( 7) << CurrentUnitVelo 580 G4cout << std::setw( 7) << CurrentUnitVelocity.x() << " "
590 << std::setw( 7) << CurrentUnitVelo 581 << std::setw( 7) << CurrentUnitVelocity.y() << " "
591 << std::setw( 7) << CurrentUnitVelo 582 << std::setw( 7) << CurrentUnitVelocity.z() << " ";
592 G4cout.precision(3); 583 G4cout.precision(3);
593 G4cout << std::setw( 7) 584 G4cout << std::setw( 7)
594 << CurrentFT.GetMomentum().mag()-St 585 << CurrentFT.GetMomentum().mag()-StartFT.GetMomentum().mag() << " ";
595 G4cout << std::setw( 9) << step_len << " " 586 G4cout << std::setw( 9) << step_len << " ";
596 G4cout << std::setw(12) << safety << " "; 587 G4cout << std::setw(12) << safety << " ";
597 if( requestStep != -1.0 ) 588 if( requestStep != -1.0 )
598 { G4cout << std::setw( 9) << requestStep 589 { G4cout << std::setw( 9) << requestStep << " "; }
599 else 590 else
600 { G4cout << std::setw( 9) << "Init/NotKn 591 { G4cout << std::setw( 9) << "Init/NotKnown" << " "; }
601 if( startVolume != nullptr) << 592 if( startVolume != 0)
602 { G4cout << std::setw(12) << startVolume 593 { G4cout << std::setw(12) << startVolume->GetName() << " "; }
603 G4cout.precision(oldprec); 594 G4cout.precision(oldprec);
604 G4cout << G4endl; 595 G4cout << G4endl;
605 } 596 }
606 else // if( verboseLevel > 3 ) 597 else // if( verboseLevel > 3 )
607 { 598 {
608 // Multi-line output 599 // Multi-line output
609 600
610 G4cout << "Step taken was " << step_len 601 G4cout << "Step taken was " << step_len
611 << " out of PhysicalStep = " << re 602 << " out of PhysicalStep = " << requestStep << G4endl;
612 G4cout << "Final safety is: " << safety << 603 G4cout << "Final safety is: " << safety << G4endl;
613 G4cout << "Chord length = " << (CurrentPos 604 G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
614 << G4endl; 605 << G4endl;
615 G4cout << G4endl; 606 G4cout << G4endl;
616 } 607 }
617 } 608 }
618 609
619 // ------------------------------------------- 610 // ---------------------------------------------------------------------------
620 // Prints Step diagnostics 611 // Prints Step diagnostics
621 // 612 //
622 void 613 void
623 G4PropagatorInField::PrintStepLengthDiagnostic 614 G4PropagatorInField::PrintStepLengthDiagnostic(
624 G4double CurrentProp 615 G4double CurrentProposedStepLength,
625 G4double decreaseFac 616 G4double decreaseFactor,
626 G4double stepTrial, 617 G4double stepTrial,
627 const G4FieldTrack& ) 618 const G4FieldTrack& )
628 { 619 {
629 G4long iprec= G4cout.precision(8); << 620 G4int iprec= G4cout.precision(8);
630 G4cout << " " << std::setw(12) << " PiF: NoZ 621 G4cout << " " << std::setw(12) << " PiF: NoZeroStep "
631 << " " << std::setw(20) << " CurrentP 622 << " " << std::setw(20) << " CurrentProposed len "
632 << " " << std::setw(18) << " Full_cur 623 << " " << std::setw(18) << " Full_curvelen_last"
633 << " " << std::setw(18) << " last pro 624 << " " << std::setw(18) << " last proposed len "
634 << " " << std::setw(18) << " decrease 625 << " " << std::setw(18) << " decrease factor "
635 << " " << std::setw(15) << " step tri 626 << " " << std::setw(15) << " step trial "
636 << G4endl; 627 << G4endl;
637 628
638 G4cout << " " << std::setw(10) << fNoZeroSte 629 G4cout << " " << std::setw(10) << fNoZeroStep << " "
639 << " " << std::setw(20) << CurrentPro 630 << " " << std::setw(20) << CurrentProposedStepLength
640 << " " << std::setw(18) << fFull_Curv 631 << " " << std::setw(18) << fFull_CurveLen_of_LastAttempt
641 << " " << std::setw(18) << fLast_Prop 632 << " " << std::setw(18) << fLast_ProposedStepLength
642 << " " << std::setw(18) << decreaseFa 633 << " " << std::setw(18) << decreaseFactor
643 << " " << std::setw(15) << stepTrial 634 << " " << std::setw(15) << stepTrial
644 << G4endl; 635 << G4endl;
645 G4cout.precision( iprec ); 636 G4cout.precision( iprec );
646 } 637 }
647 638
648 // Access the points which have passed through 639 // Access the points which have passed through the filter. The
649 // points are stored as ThreeVectors for the i 640 // points are stored as ThreeVectors for the initial impelmentation
650 // only (jacek 30/10/2002) 641 // only (jacek 30/10/2002)
651 // Responsibility for deleting the points lies 642 // Responsibility for deleting the points lies with
652 // SmoothTrajectoryPoint, which is the points' 643 // SmoothTrajectoryPoint, which is the points' final
653 // destination. The points pointer is set to N 644 // destination. The points pointer is set to NULL, to ensure that
654 // the points are not re-used in subsequent st 645 // the points are not re-used in subsequent steps, therefore THIS
655 // METHOD MUST BE CALLED EXACTLY ONCE PER STEP 646 // METHOD MUST BE CALLED EXACTLY ONCE PER STEP. (jacek 08/11/2002)
656 647
657 std::vector<G4ThreeVector>* 648 std::vector<G4ThreeVector>*
658 G4PropagatorInField::GimmeTrajectoryVectorAndF 649 G4PropagatorInField::GimmeTrajectoryVectorAndForgetIt() const
659 { 650 {
660 // NB, GimmeThePointsAndForgetThem really fo 651 // NB, GimmeThePointsAndForgetThem really forgets them, so it can
661 // only be called (exactly) once for each st 652 // only be called (exactly) once for each step.
662 653
663 if (fpTrajectoryFilter != nullptr) 654 if (fpTrajectoryFilter != nullptr)
664 { 655 {
665 return fpTrajectoryFilter->GimmeThePointsA 656 return fpTrajectoryFilter->GimmeThePointsAndForgetThem();
666 } 657 }
667 return nullptr; << 658 else
>> 659 {
>> 660 return nullptr;
>> 661 }
668 } 662 }
669 663
670 // ------------------------------------------- 664 // ---------------------------------------------------------------------------
671 // 665 //
672 void 666 void
673 G4PropagatorInField::SetTrajectoryFilter(G4VCu 667 G4PropagatorInField::SetTrajectoryFilter(G4VCurvedTrajectoryFilter* filter)
674 { 668 {
675 fpTrajectoryFilter = filter; 669 fpTrajectoryFilter = filter;
676 } 670 }
677 671
678 // ------------------------------------------- 672 // ---------------------------------------------------------------------------
679 // 673 //
680 void G4PropagatorInField::ClearPropagatorState 674 void G4PropagatorInField::ClearPropagatorState()
681 { 675 {
682 // Goal: Clear all memory of previous steps, 676 // Goal: Clear all memory of previous steps, cached information
683 677
684 fParticleIsLooping = false; 678 fParticleIsLooping = false;
685 fNoZeroStep = 0; 679 fNoZeroStep = 0;
686 680
687 fSetFieldMgr = false; // Has field-manager <<
688 fEpsilonStep= 1.0e-5; // Relative accuracy <<
689 <<
690 End_PointAndTangent= G4FieldTrack( G4ThreeVe 681 End_PointAndTangent= G4FieldTrack( G4ThreeVector(0.,0.,0.),
691 G4ThreeVe 682 G4ThreeVector(0.,0.,0.),
692 0.0,0.0,0 683 0.0,0.0,0.0,0.0,0.0);
693 fFull_CurveLen_of_LastAttempt = -1; 684 fFull_CurveLen_of_LastAttempt = -1;
694 fLast_ProposedStepLength = -1; 685 fLast_ProposedStepLength = -1;
695 686
696 fPreviousSftOrigin= G4ThreeVector(0.,0.,0.); 687 fPreviousSftOrigin= G4ThreeVector(0.,0.,0.);
697 fPreviousSafety= 0.0; 688 fPreviousSafety= 0.0;
698 <<
699 fNewTrack = true; <<
700 } 689 }
701 690
702 // ------------------------------------------- 691 // ---------------------------------------------------------------------------
703 // 692 //
704 G4FieldManager* G4PropagatorInField:: 693 G4FieldManager* G4PropagatorInField::
705 FindAndSetFieldManager( G4VPhysicalVolume* pCu 694 FindAndSetFieldManager( G4VPhysicalVolume* pCurrentPhysicalVolume )
706 { 695 {
707 G4FieldManager* currentFieldMgr; 696 G4FieldManager* currentFieldMgr;
708 697
709 currentFieldMgr = fDetectorFieldMgr; 698 currentFieldMgr = fDetectorFieldMgr;
710 if( pCurrentPhysicalVolume != nullptr ) 699 if( pCurrentPhysicalVolume != nullptr )
711 { 700 {
712 G4FieldManager *pRegionFieldMgr = nullptr 701 G4FieldManager *pRegionFieldMgr = nullptr, *localFieldMgr = nullptr;
713 G4LogicalVolume* pLogicalVol = pCurrentPh 702 G4LogicalVolume* pLogicalVol = pCurrentPhysicalVolume->GetLogicalVolume();
714 703
715 if( pLogicalVol != nullptr ) 704 if( pLogicalVol != nullptr )
716 { 705 {
717 // Value for Region, if any, overrides 706 // Value for Region, if any, overrides
718 // 707 //
719 G4Region* pRegion = pLogicalVol->GetR 708 G4Region* pRegion = pLogicalVol->GetRegion();
720 if( pRegion != nullptr ) 709 if( pRegion != nullptr )
721 { 710 {
722 pRegionFieldMgr = pRegion->GetField 711 pRegionFieldMgr = pRegion->GetFieldManager();
723 if( pRegionFieldMgr != nullptr ) 712 if( pRegionFieldMgr != nullptr )
724 { 713 {
725 currentFieldMgr= pRegionFieldMgr 714 currentFieldMgr= pRegionFieldMgr;
726 } 715 }
727 } 716 }
728 717
729 // 'Local' Value from logical volume, 718 // 'Local' Value from logical volume, if any, overrides
730 // 719 //
731 localFieldMgr = pLogicalVol->GetFieldM 720 localFieldMgr = pLogicalVol->GetFieldManager();
732 if ( localFieldMgr != nullptr ) 721 if ( localFieldMgr != nullptr )
733 { 722 {
734 currentFieldMgr = localFieldMgr; 723 currentFieldMgr = localFieldMgr;
735 } 724 }
736 } 725 }
737 } 726 }
738 fCurrentFieldMgr = currentFieldMgr; 727 fCurrentFieldMgr = currentFieldMgr;
739 728
740 // Flag that field manager has been set 729 // Flag that field manager has been set
741 // 730 //
742 fSetFieldMgr = true; 731 fSetFieldMgr = true;
743 732
744 return currentFieldMgr; 733 return currentFieldMgr;
745 } 734 }
746 735
747 // ------------------------------------------- 736 // ---------------------------------------------------------------------------
748 // 737 //
749 G4int G4PropagatorInField::SetVerboseLevel( G4 738 G4int G4PropagatorInField::SetVerboseLevel( G4int level )
750 { 739 {
751 G4int oldval = fVerboseLevel; 740 G4int oldval = fVerboseLevel;
752 fVerboseLevel = level; 741 fVerboseLevel = level;
753 742
754 // Forward the verbose level 'reduced' to Ch 743 // Forward the verbose level 'reduced' to ChordFinder,
755 // MagIntegratorDriver ... ? 744 // MagIntegratorDriver ... ?
756 // 745 //
757 auto integrDriver = GetChordFinder()->GetInt 746 auto integrDriver = GetChordFinder()->GetIntegrationDriver();
758 integrDriver->SetVerboseLevel( fVerboseLevel 747 integrDriver->SetVerboseLevel( fVerboseLevel - 2 );
759 G4cout << "Set Driver verbosity to " << fVer 748 G4cout << "Set Driver verbosity to " << fVerboseLevel - 2 << G4endl;
760 749
761 return oldval; 750 return oldval;
762 } 751 }
763 752
764 // ------------------------------------------- 753 // ---------------------------------------------------------------------------
765 // 754 //
766 void G4PropagatorInField::ReportLoopingParticl 755 void G4PropagatorInField::ReportLoopingParticle( G4int count,
767 756 G4double StepTaken,
768 757 G4double StepRequested,
769 758 const char* methodName,
770 << 759 G4ThreeVector momentumVec,
771 760 G4VPhysicalVolume* pPhysVol )
772 { 761 {
773 std::ostringstream message; 762 std::ostringstream message;
774 G4double fraction = StepTaken / StepRequest 763 G4double fraction = StepTaken / StepRequested;
775 message << " Unfinished integration of trac 764 message << " Unfinished integration of track (likely looping particle) "
776 << " of momentum " << momentumVec < 765 << " of momentum " << momentumVec << " ( magnitude = "
777 << momentumVec.mag() << " ) " << G4 766 << momentumVec.mag() << " ) " << G4endl
778 << " after " << count << " field su 767 << " after " << count << " field substeps "
779 << " totaling " << std::setprecisio 768 << " totaling " << std::setprecision(12) << StepTaken / mm << " mm "
780 << " out of requested step " << std 769 << " out of requested step " << std::setprecision(12)
781 << StepRequested / mm << " mm "; 770 << StepRequested / mm << " mm ";
782 message << " a fraction of "; 771 message << " a fraction of ";
783 G4int prec = 4; 772 G4int prec = 4;
784 if( fraction > 0.99 ) 773 if( fraction > 0.99 )
785 { 774 {
786 prec = 7; 775 prec = 7;
787 } 776 }
788 else 777 else
789 { 778 {
790 if (fraction > 0.97 ) { prec = 5; } 779 if (fraction > 0.97 ) { prec = 5; }
791 } 780 }
792 message << std::setprecision(prec) 781 message << std::setprecision(prec)
793 << 100. * StepTaken / StepRequested 782 << 100. * StepTaken / StepRequested << " % " << G4endl ;
794 if( pPhysVol != nullptr ) << 783 if( pPhysVol )
795 { 784 {
796 message << " in volume " << pPhysVol->Get << 785 message << " in volume " << pPhysVol->GetName() ;
797 auto material = pPhysVol->GetLogicalVolum << 786 auto material = pPhysVol->GetLogicalVolume()->GetMaterial();
798 if( material != nullptr ) << 787 if( material != nullptr )
799 { << 788 message << " with material " << material->GetName()
800 message << " with material " << materia << 789 << " ( density = "
801 << " ( density = " << 790 << material->GetDensity() / ( g/(cm*cm*cm) ) << " g / cm^3 ) ";
802 << material->GetDensity() / ( g <<
803 } <<
804 } 791 }
805 else 792 else
806 { 793 {
807 message << " in unknown (null) volume. " << 794 message << " in unknown (null) volume. " ;
808 } 795 }
809 G4Exception(methodName, "GeomNav1002", Just 796 G4Exception(methodName, "GeomNav1002", JustWarning, message);
810 } 797 }
811 798
812 // ------------------------------------------- 799 // ---------------------------------------------------------------------------
813 // 800 //
814 void G4PropagatorInField::ReportStuckParticle( 801 void G4PropagatorInField::ReportStuckParticle( G4int noZeroSteps,
815 802 G4double proposedStep,
816 803 G4double lastTriedStep,
817 804 G4VPhysicalVolume* physVol )
818 { 805 {
819 std::ostringstream message; 806 std::ostringstream message;
820 message << "Particle is stuck; it will be k 807 message << "Particle is stuck; it will be killed." << G4endl
821 << " Zero progress for " << noZero 808 << " Zero progress for " << noZeroSteps << " attempted steps."
822 << G4endl 809 << G4endl
823 << " Proposed Step is " << propose 810 << " Proposed Step is " << proposedStep
824 << " but Step Taken is "<< lastTrie 811 << " but Step Taken is "<< lastTriedStep << G4endl;
825 if( physVol != nullptr ) 812 if( physVol != nullptr )
826 { <<
827 message << " in volume " << physVol->Get 813 message << " in volume " << physVol->GetName() ;
828 } <<
829 else 814 else
830 { <<
831 message << " in unknown or null volume. 815 message << " in unknown or null volume. " ;
832 } <<
833 G4Exception("G4PropagatorInField::ComputeSt 816 G4Exception("G4PropagatorInField::ComputeStep()",
834 "GeomNav1002", JustWarning, mes 817 "GeomNav1002", JustWarning, message);
835 } <<
836 <<
837 // ------------------------------------------- <<
838 <<
839 // ------------------------------------------- <<
840 // Methods to alter Parameters <<
841 // ------------------------------------------- <<
842 <<
843 // Was a data member (of an object) -- now mov <<
844 G4double G4PropagatorInField::GetLargestAccep <<
845 { <<
846 return fLargestAcceptableStep; <<
847 } <<
848 <<
849 // ------------------------------------------- <<
850 // <<
851 void G4PropagatorInField::SetLargestAcceptable <<
852 { <<
853 if( fLargestAcceptableStep>0.0 ) <<
854 { <<
855 fLargestAcceptableStep = newBigDist; <<
856 } <<
857 } <<
858 <<
859 // ------------------------------------------- <<
860 <<
861 G4double G4PropagatorInField::GetMaxStepSizeMu <<
862 { <<
863 return fMaxStepSizeMultiplier; <<
864 } <<
865 <<
866 // ------------------------------------------- <<
867 <<
868 void G4PropagatorInField::SetMaxStepSizeMu <<
869 { <<
870 fMaxStepSizeMultiplier=vm; <<
871 } <<
872 <<
873 // ------------------------------------------- <<
874 <<
875 G4double G4PropagatorInField::GetMinBigDistanc <<
876 { <<
877 return fMinBigDistance; <<
878 } <<
879 <<
880 // ------------------------------------------- <<
881 <<
882 void G4PropagatorInField::SetMinBigDistanc <<
883 { <<
884 fMinBigDistance= val; <<
885 } 818 }
886 819