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