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